计算气动光学研究进展

史可天, 马汉东

(中国航天空气动力技术研究院, 北京 100074)

0 引 言

带有光学成像探测制导系统的飞行器在大气层内高速飞行时，其光学窗口与来流之间发生剧烈的相互作用，形成激波、边界层、混合层等复杂流场，对光学成像探测系统造成图像传输干扰，引起目标图像偏移、抖动、模糊以及能量散布，这种效应称为气动光学效应[1-3]。气动光学效应的存在会严重降低光学成像探测系统的性能，已成为红外成像制导等高速流动光学传输系统研制的关键技术之一。

计算气动光学是采用数值计算的方法研究空气动力流场对光波传输和光学成像影响及其校正的一门交叉学科[4-5]，包括流场数值计算和光学传输数值计算两大方面。光学折射率n和流场密度ρ通过Gladstone-Dale关系式相联系：

n=1+KGDρ(1)

其中KGD为Gladstone-Dale系数，其值与入射光波波长和流体特性相关。

计算气动光学的发展主要经历了三个阶段，一是早期的基于CFD简化方法的统计计算，二是基于RANS流场的光学传输计算，三是基于LES/DNS瞬态流场的光学传输计算。三种方法相辅相成，前两种方法可以快速地预测气动光学效应，第三种方法可以对前两种方法进行校核和验证，并可以细致分析气动光学效应，研究其形成机理。

基于CFD简化方法的统计估算方法采用流场简化模型计算方法获得光学传输模拟所需的流场信息，然后用统计光学原理构建脉动流场的光学畸变统计模型，以波面均方值、光学传递函数平均值作为气动光学效应参数，对光学畸变进行平均或统计描述。

基于RANS流场的光学传输计算方法针对流场平均密度分布，采用几何光线追迹法和Fourier光学相结合的方法，计算平均流场引起的光学畸变。通过波面光程差分布、像面光强(点扩散函数)、瞄视误差、Strehl比、环围能量曲线等多种气动光学效应参数对流场光学传输特性进行定量描述。由于RANS计算采用的Reynolds平均抹平了脉动运动的细节，所以难以直接反映拟序结构运动对光学传输的影响。

基于LES/DNS流场的光学传输计算方法采用几何光学/波动光学和Fourier光学相结合的方法，对流场瞬态拟序结构的光学传输特性进行计算，能够了解光学畸变和流场结构特征变化的时间历程，更有利于建立二者间的内在关联。

以上三种方法都是目前气动光学研究正在应用的计算方法，但后两种方法更为多见。这三种方法涉及的计算流体力学和光学传输方法的关系，如表1所示。

表1 计算气动光学涉及的流场和光传输方法Table 1 CFD and optical transmission method in the computational aero-optics

本文对这三种方法的内涵、研究发展历程以及新近的一些研究进展进行综述，最后指出未来的研究重点和方向。

1 基于CFD简化方法的光学统计估算

早在1952年Liepmann[6]就开展了气动光学研究的开创性工作，他为了确定纹影系统应用于高速流场时的极限灵敏度，应用几何光学原理推导了小孔径光束通过湍流流场后，偏折角均方值〈θ2〉的计算公式：

〈(∂v/∂y)2〉Rv(|y-ζ|)dydζ(2)

其中：折射率脉动值v=KGDρ′/n0，Rv(y)表示折射率变化相关函数。这一计算公式首次建立了光学畸变参数与流场参数之间的统计关联，并且还在计算中考虑了沿光线传播方向变化的相关函数。Liepmann提出的光束偏折角均方值的计算公式获得了实验证实[7]，实验测量结果与Liepmann的计算十分吻合。

Sutton[8]根据统计光学理论提出了近场光波相位均方值的计算模型，首次建立了大孔径条件下波面畸变参数与流场参数之间的定量联系：

在这一计算公式中同样考虑了密度脉动相关尺度对光学畸变的影响。

在随后的10多年中，大部分的气动光学研究工作都围绕Sutton统计模型开展，主要是基于这一统计模型，应用直接或间接的方法计算近场光波相位均方值[9]。

这一时期的计算气动光学研究，其所需的流场密度脉动信息主要是通过计算流体力学简化模型方法获得的，包括点涡法、无粘位势流方法等[10-11]。

虽然Sutton统计模型在工程应用中取得了一定的成功，但仍有学者对统计方法在高超声速飞行器气动光学效应研究中的适应性提出质疑。Havaner[12]通过理论分析认为Sutton统计模型仅适用于完全气体状态下的二维亚声速均匀各向同性湍流流动，对于具有多尺度涡结构的高速复杂流场，如超声速湍流混合层、两相流等非均匀各向同性流动，并不适用。

殷兴良[13]从光学传递函数基本定义出发，基于密度脉动的空间平稳性假设，提出了三维超声速流场均匀脉动的光学畸变统计计算模型：

其中密度脉动相关函数C(x,y,z)常采用von Karman模型、指数模型和Gauss模型等形式。Gauss模型：

史可天等[14]放弃了密度脉动的空间平稳性假设，仅认为密度脉动随机过程是Gauss随机过程，发展了一种适用于三维超声速流场非均匀脉动的光学畸变统计模型:

[Rρ(ξ+x,η+y,ζ+z,0,0-z)/2+

Rρ(ξ,η,ζ;0,0,z)/2-

Rρ(ξ,η,ζ;x,y,z)]dzdζ}dξdη(2)

这一计算模型可用于复杂外形光学窗口脉动流场的光学传递函数平均值计算。如果密度脉动分布具有空间平稳性，式(6)中最外层积分号内的表达式将与ξ和η无关，公式则可以退化为殷兴良提出的统计计算模型公式(4)。

2 基于RANS流场的光学传输计算

RANS方法认为湍流的脉动运动是随机和无规则的，通过湍流模型描述流场脉动运动对平均运动的作用[15-16]。RANS方法可以针对复杂外形飞行器的绕流流场进行计算，提供流场平均密度分布，并在一定精度上对脉动密度进行重构。

基于RANS方法获得流场平均密度后，可以采用几何光线追迹法计算光瞳范围内每一条光线的传输路径，由此得到入射平面光波通过流场后产生的光程差，获得光瞳函数；然后采用Fourier光学方法计算畸变波面通过成像系统到达像平面的光强分布，即点扩散函数，再由点扩散函数计算瞄视误差、Strehl比、环围能量曲线等气动光学效应参数对平均流场引起的成像偏移、模糊等效应进行定量描述[17]。

光线追迹法根据几何光学原理、光线微分方程对光线在气体介质内的传输路径进行求解[18]：

首先给定光线传输路径上一点处的位置向量和光线单位切向量，然后根据该点处的密度以及密度梯度，计算得到传输路径下一点的位置向量和光线单位切向量。为了提高光线方程的求解精度，可将光线方程化为1阶微分方程组，然后采用4步4阶Runge-Kutta方法显式推进[18]：

获得每条光线传播路径后，即可求出光程和光瞳函数[19]:

Fourier光学方法对光瞳函数进行Fourier变换，得到像平面光波复振幅分布，进而计算出点扩散函数(即像平面光强分布)和光学传递函数[20-21]：

PSF(x′,y′)=|U(x′,y′)|2(11)

图1给出光学传输效应计算采用的坐标系。

图1 光学传输效应计算坐标系示意图[17]Fig.1 The coordinate for the optical effects simulation[17]

RANS方法除了提供光学传输计算所需的平均密度，还可以根据平均密度分布应用工程经验公式对脉动密度进行估算，以此作为计算输入条件，应用统计方法对脉动密度引起的光学畸变进行数值计算。

传统的脉动密度预测模型基于混合长度理论[22]:

其中:C为加权常数，对亚声速流动C=0.32，对超声速流动C=0.57。

潘宏禄等[23-24]对预测模型中的加权常数C进行了修正。其采用高阶精度的大涡模拟方法，对超声速剪切湍流场进行计算，由流场瞬态密度分布统计得到脉动密度，给出了混合长度模型中系数C的最适取值，见图2。在转捩区内，C在0.1～0.55间剧烈变化，而在充分发展湍流区内C在0.45～0.42间变化，波动范围较小，建议取值0.43。

图2 平板边界层模型系数分布[23]Fig.2 Coefficients for the plane boundary layer[23]

闫溟等[25]针对气动光学效应的RANS方法进行了分析，对常用的几种湍流模型进行了评估，将RANS方程和折射率输运方程联合求解，对折射率脉动均方值g=〈n′2〉进行了预测：

图3为闫溟等计算得到的湍流密度脉动重构结果。

(a) Ma=1.72

(b) Ma=3.56

应用几何光学、Fourier光学方法，针对RANS平均流场数据可以计算得到平均流场的光学传输特性；应用统计光学方法，针对密度脉动重构数据，可以获得脉动流场的光学传输特性。如果将平均流场和脉动流场分别看作空间不变的线性光学系统，就可以通过卷积方法将这二者的点扩散函数结合，获得湍流流场气动光学传输特性[26]：

PSF(x′,y′)=PSFM(x′,y′)*PSFT(x′,y′)(14)

其中:PSFM(x′,y′)为平均流场产生的点扩散函数，PSFT(x′,y′)为脉动流场产生的点扩散函数，*代表卷积。

在式(14)对光学畸变的计算中，根据红外成像探测系统的成像积分时间，确定截止频率，可以将密度脉动区分为低频脉动和高频脉动分别进行处理，低频脉动引起成像抖动，高频脉动引起成像模糊。截止频率kFZ根据流场特征速度VT和内时间尺度τin确定：

kFZ=VT·τin(15)

基于RANS方法的气动光学效应计算是气动光学效应工程预测中常用的方法，可以针对复杂外形光学窗口绕流流动进行光学传输效应预测，定量描述平均流场和脉动流场引起的光学畸变，为光学系统设计提供依据。

陈澄等[27]、史可天等 [28-29]、陈勇等 [30]、张黎等 [31]、万士正等 [32]针对不同形式的光学窗口(凹窗、凸窗及头罩等)、不同飞行弹道进行了高速飞行器气动光学效应计算分析，对光学传输效应的影响因素进行了分析，包括飞行高度、飞行Mach数、光线入射角、成像系统工作波长及成像积分时间等参数，在工程设计中取得了成功的应用。

由于RANS计算采用的Reynolds平均方法抹平了脉动运动的细节，丢失了很多包含在脉动运动中的重要信息[33]，所以难以直接反映拟序结构运动等湍流特性对光学传输的影响，在气动光学效应机理研究中有一定的局限性。

3 基于LES/DNS瞬态流场的光学传输计算

随着计算流体力学的发展和计算能力的提高，大涡数值模拟(LES)和直接数值模拟(DNS)以其具备的气动流场瞬态描述能力，为分析波前畸变和流场结构变化的时间历程，研究气动光学效应机理提供了重要的手段。

DNS虽然能直接求解流场各种尺度的涡结构瞬时信息，但对计算机的运算速度和容量要求很高，在现有的计算条件下，大多是对低Reynolds数下的简单外形湍流问题或局部湍流流动的模拟[34]。例如时间发展剪切层的观察点随流体运动，模化的是观察点附近流场的时间演化过程(图4)，与真实的物理过程并不完全等价[35]。

图4 时间发展剪切层涡结构DNS结果[35]Fig.4 Turbulence vortex structure from temporal developing mixing layer DNS[35]

LES是介于RANS和DNS之间的一种数值模拟方法，它通过空间滤波将流场结构分为大小两种尺度，直接模拟尺度大于计算网格的流场结构，以亚格子模型近似小尺度涡的耗散行为[36-37]。大涡模拟方法可以获得流场拟序结构分布(图5)，为光学传输计算提供密度信息[38-39]。

图5 可压缩混合层LES获得的湍流涡结构[37]Fig.5 Turbulence vortex structure from compressible mixing layer LES[37]

获得瞬态流场后，应用几何光学、Fourier光学方法即可获得每一时刻流场的点扩散函数等光学畸变定量特性。由于DNS和LES方法获得丰富的流场结构信息，如果光学计算采用异于流场计算的独立网格，需要确保流场与光场间数据转换的密度信息精确度。

由于DNS和LES方法对流动细节的分辨能力远强于RANS方法，可以反映一些较小尺度的流动结构，光波在这些小尺度结构的传输过程中，有可能表现出一定的波动效应，需要采用波动光学的方法来进行高精度光学传输计算。在气动流场的光学传输计算中，只要求解傍轴近似条件下的光波方程即可获得足够精确的结果[40]：

White等[41-42]针对三维可压缩剪切层流场进行LES计算，然后比较了几何光学和波动光学的计算精度，计算结果表明由两种方法获得的光程分布差别甚微，仅有0.03λ。White认为目前大多数光学传输计算中，高阶精度的几何光学方法，如4步4阶Runge-Kutta推进的光线追迹法，可以满足计算精度需求。

LES为气动光学机理研究提供了有力的研究手段。基于LES的气动光学机理研究目前主要集中在可压缩混合层流动[43-44]、超声速边界层流动[45]、后台阶流动[46]等高速流动的光学传输效应计算分析。这些流动虽然几何外形简单，却具有可压缩湍流的重要特征。对这些流动的光学传输计算分析，可以认识和了解大尺度拟序结构的空间发展变化、激波/激波干扰、激波/边界层干扰等现象对光学传输的影响。

通过对可压缩湍流流场的气动光学效应计算分析，得到了一些光学畸变的影响规律。例如，混合层流动从层流到转捩直至完全湍流全过程的光学传输计算表明，随着转捩的发生，流场引起严重的光学畸变，而转捩完成以后，光学畸变程度有所减弱，但仍然严重于转捩以前的情况(图6)；亚声速边界层气动光学效应的计算发现，黏性底层和缓冲层对光学传输的影响很小，光学畸变主要是由对数区流动引起的。这些光学传输的影响规律为气动光学效应校正提供了理论依据。

图6 剪切层波面畸变空间分布[44]Fig.6 Wavefront distortion in the mixing layer[44]

在开展典型湍流流动气动光学效应机理研究的同时，有学者采用LES方法对真实光学窗口外形进行流场数值模拟，对其瞬态光学传输效应进行分析，获得了更精确的高速飞行器光学系统气动光学效应数据。其中，最值得关注的是美国空军气动光学飞行实验室(Airborne Aero-Optics Laboratory)凸台激光系统的气动光学效应计算研究[47]。

在这一研究中，Mathews等[48]采用带壁面模型的LES方法对凸台绕流场进行模拟，计算了流场引起的激光束分散特性。他们首先通过与试验数据的对比验证，确认了波面畸变均方值计算的准确性(图7)，然后对不同高度角和周向角开展了共计390个光线入射角的光学传输计算，分析了流场涡结构对不同口径激光束的影响规律，给出了激光束失效的入射角度范围，图8为其计算获得的绕流流线分布。为了获得精准的气动光学效应数据，他们甚至考虑了凸台系统安装缝隙对流场结构和光学传输效应的影响[49]。

图7 光程差均方根值计算与试验结果比较[48]Fig.7 Comparison of OPDrms between computation and wind tunnel measurement[48]

根据数值计算和试验测量获得的基础数据，美国空军研究实验室采用流体控制和光学补偿方法，研制了航空自适应光束控制台系统，使激光武器可以在高速移动平台上从任意角度精确打击敌方目标。如果没有这一系统，由于机体周围强烈的湍流会造成激光束的严重分散，激光武器只能对正前方的目标实施精确打击。

图8 绕流流线分布[48]Fig.8 Streamlines of the velocity field[48]

4 结 论

在高速流场气动光学效应数值模拟中，三种计算方法相互补充，基于CFD简化方法的光学统计估算和基于RANS流场的光学传输计算具有较高的效率，可以快速提供流场的光学畸变数据，适用于复杂外形光学系统设计的气动光学传输效应分析，而基于LES/DNS瞬态流场的光学传输计算虽然对计算资源要求较高，却可以提供更精细的光学畸变信息。通过对三种方法的总结，得到以下结论：

(1)基于CFD简化方法的光学统计估算方法通过近几年的研究，从细光束偏折估算发展到大孔径平面光波相位均方值、传递函数均值计算，适用范围和计算精度都获得了大幅提升。但是，计算过程中采用的密度相关函数在各向异性复杂非均匀流动条件下还有待进一步发展和验证。

(2)基于RANS流场的气动光学效应计算方法可以针对复杂外形飞行器的绕流流场进行气动光学计算，提供平均流场的成像偏移和模糊特性；通过脉动流场重构方法，还可以对脉动密度引起的成像模糊进行评估，是气动光学效应应用研究中常用的一种方法。由于光学传输计算需要考虑光波通过流场全路径上的密度及其脉动的影响，所以如何精准描述流场的空间分布是RANS方法应用于气动光学效应计算时有待解决的问题。此外，密度脉动重构方法的精度还需要通过风洞试验进行确认。

(3)基于LES/DNS流场的光学传输计算方法能够提供光学畸变的时间变化历程，更有利于分析其与流场特性间的关联。目前DNS方法受到计算条件限制，适用于局部湍流流动的气动光学效应机理分析。而LES方法的计算量相对DNS较小，同时其获得的大尺度结构又正是造成光学相位畸变的主要原因，因此LES方法在当前及未来相当长时期内，都将是计算气动光学机理和应用分析的有力工具。为了能提供更精确的流场信息，LES方法在高精度数值格式、亚格子模型适用性等方面还需进一步研究。

(4)气动光学效应瞬态计算虽然获得了光学畸变的实时数据，但其与流场特性之间的关联多数研究还是定性的，例如发现转捩区的大尺度结构引起严重的光学畸变效应等，目前尚未建立二者之间的定量关联式。下一步计算气动光学的研究中，采用一些流场结构定量分析方法，如湍流结构系综统计动力学方法等[50-51]，才有可能对流场结构及其光学传输特性进行定量分析，揭示湍流流动结构的光学传输机理，发现流场结构对光学畸变的定量影响规律，为最终实现气动光学自适应校正提供理论依据。