典型光学窗口流场的气动光学效应数值模拟

谭小童，许和勇，*，田仁治

(1.西北工业大学 翼型、叶栅空气动力学国家级重点实验室，西安 710072；2.代尔夫特理工大学 航空航天工程学院，荷兰 代尔夫特 999025)

0 引言

近些年，随着激光武器和超声速飞行器的发展，气动光学受到越来越多的关注。根据光线的传输路径，气动光学系统可以分为主动光学系统（激光发射器）和被动光学系统（红外寻的器）。在主动光学系统中，气动光学效应导致高能激光发射器发出的激光束能量减弱，失去攻击力或聚焦偏离攻击目标；在被动光学系统中，气动光学效应导致光学接收器上的图像变形或者光学信号减弱，使跟踪方向发生偏差。无论是主动光学系统还是被动光学系统都分为近场“气动光学”区域和远场“大气光学”区域。

迄今为止，人们已经对湍流边界层[1-2]、自由剪切层[3-4]、分离剪切层和湍流尾迹[5-6]等几种典型流动进行了广泛研究。截至目前，广大研究人员已经设计出许多气动光学系统装置，其中圆柱半球转塔是目前激光武器的首选发射装置[7]。针对圆柱半球转塔的气动光学，科研工作者对具有平坦或共形的光学窗口的转塔进行了许多风洞实验和数值研究[8-10]。美国开展了ALL、ABL、ATL、SHIELD 等一系列项目来验证各种转塔系统的气动光学性能，并做了一些流动控制研究[11]。Wang 等[12]、Jumper 和Gordeyev[13]也对近年来的计算和实验研究工作进行了总结。流动控制也被应用于减轻三维转塔的流动分离和相关气动光学畸变[14-16]。

进入21 世纪后，我国的气动光学研究开始起步，2003 年殷兴良[17]出版的《气动光学原理》一书奠定了我国气动光学的基础。李桂春[18]在其著作《气动光学》一书中对气动光学原理和实验测量方法进行了详细论述。史可天等[19]对新近的一些计算气动光学方法进行了归纳总结。孙喜万等[20]基于WCNS-E-5 对斜激波、超声速混合剪切层、亚声速混合流进行了广泛的数值模拟验证研究，并系统回顾了近年来气动光学研究的进展[21]。丁浩林等[2]对超声速湍流边界层的气动光学效应开展了风洞实验研究；此外，他们还开展了带超声速气膜的高超声速光学头罩的气动光学效应抑制实验[22-23]。董航等[24]对圆柱半球转塔模型进行了气动光学效应时空特性分析。

为了方便研究，通常将绕流流场分为平均流场分量和脉动流场分量，因此气动光学效应也分为平均流场效应和脉动流场效应。平均流场效应决定平移、倾斜、离焦以及像散等低阶项，容易通过自适应光学系统矫正，但是其绝对值较大；脉动流场效应影响气动光学效应中的高阶项，使用自适应光学系统矫正很困难[25-26]。

目前，采用Zernike 多项式拟合波前畸变并结合自适应光学思想矫正低阶项的研究方法还较少，而且将近场气动光学和远场衍射结合的研究方法也不多。本文主要基于DDES 方法和光线追迹算法计算3D 共形窗口转塔和3D 平面窗口转塔的气动光学效应；采用Zernike 多项式拟合波前畸变，并结合自适应光学方法初步分析了两种光学窗口在不同转角下的气动光学效应；基于光线方程和点扩散函数，模拟光线经过近场气动光学畸变后的远场自由衍射光强。

1 物理模型和数值计算方法

1.1 流场计算方法

为了提高湍流模型在转塔后部流动分离区域的预测能力，本文采用延迟脱体涡模拟（DDES），该方法是将雷诺平均Navier-Stokes（RANS）模拟和大涡模拟（LES）相结合。Strelets[27]在1997 年提出的SSTDES 模型，即在壁面附近采用SST 模型求解RANS方程模拟小尺度湍流结构；在分离区域转化为LES求解大尺度湍流结构，并引入亚格子模型模化该区域的小尺度湍流。而RANS/LES 的转换则是通过当地的湍流尺度和网格尺度的相互对应关系来确定。单一的LES 方法求解边界层区域需要耗费大量的计算资源，而DES 方法只需要在边界层区域的网格尺度达到RANS 的精度即可，因此可以大大地减少计算量。

Strelets 提出的SST-DES 模型的原理是，当RANS模型预测的湍流长度Lt大于局部网格间距的区域时，从SST 模式切换到LES 模式。在这种情况下，湍动能方程中计算耗散项所用的长度尺度由局部网格间距替代，具体的转换原理如下：

在DES 的公式中选择当地网格最大边长 δ是因为模型在边界层中应使用RANS 模式，最大长度是符合要求的最安全估计。F1是用于区分k-ω和k-ε方程的混合函数。但是SST-DES 方法产生的实际问题便是当地网格的最大边长 δ小于边界层厚度时，RANS 模式提前过渡到LES 模式，导致模化应力不足（MSD），容易引起网格诱导分离问题（GIS）。所以，本文在SSTDES 的基础上对其耗散项进行修正，得到SSTDDES 模型[28]，将耗散项修正为如下的形式：

该修正推迟了RANS 模式转换到LES 模式，有效地避免了GIS 问题。

1.2 物理模型及边界条件

为了研究转塔附近不同流态对气动光学效应的影响，本文选取两种不同的光学窗口（平面光学窗口和共形光学窗口），并对比研究了两者在不同转角下的气动光学效应。转塔模型的几何模型如图1 所示，底部圆柱半径和顶部半球的半径R=0.152 4 m，圆柱高H=0.114 3 m，光学窗口的平面半径r=0.137 2 m，光线的发射仰角EL=30°，而光线的转角AL=0°、90°、180°，其中共形光学窗口模型在三个转角下的流场计算模型相同。

图1 转塔几何模型示意图及光学窗口视角的定义Fig.1 Schematic diagram of the turret geometry model and visual angle definition of the optical window

图2 展示了转角为 0◦时，平面光学窗口转塔的流场计算域和边界条件设置，其余工况的流场计算域和边界条件与该工况相同。为了稳健地计算流场，将上游距离圆心O7.5R处设置为计算域入口，采用速度入口边界条件；而将下游距离圆心O20R处设置为计算域出口，采用静压出口边界条件；上表面及两侧面采用绝热自由滑移边界条件，而下表面和转塔表面则采用绝热无滑移边界条件[29]。

图2 转塔的流场计算域及其边界条件示意图Fig.2 Computational domain and boundary conditions for the flow simulation of a turret

数值计算选取的马赫数为Ma∞=0.4，雷诺数为Re∞=1.43×106，其余计算条件为标准海平面条件。图3 展示了两个模型的表面结构化面网格分布，为了保证壁面第一层网格y+≤1，在无滑移边界条件上的第一层网格高度为1 .67×10−5m。平面光学窗口转塔的总网格量为1 .4×107，而共形光学窗口转塔的网格量为1 .35×107。

图3 转塔表面面网格分布Fig.3 Turret surface grid distribution

1.3 光线追迹方法

通过GladStone-Dale 关系式可知，折射率和空气中的密度直接相关，其表达式如下：

式中：r=xi+y j+zk，表示光线流场中的位置，而n(r) 和 ρ(r)分别为该位置处的折射率和密度；KGD为GladStone-Dale系数，在可见光条件下，KGD近似取值2.27×10−4m3/kg。在光线光学中，光线沿路径对折射率的积分定义为光程（optical path length,OPL），其表达式如下：

式中，（x′,y′）为光学窗口平面中的光学坐标，z′为光线的发射方向。而光线在不均匀流场中的传播导致了各条光线光程的不同，被定义为光程差（optical path difference,OPD），其表达式如下：

其中〈OPL(x′,y′)〉表示 OPL(x′,y′)的空间平均。

根据波动方程的射线近似，可以得到光线在不均匀介质中的传播规律即光线方程，该方程的表达式如下：

由于该方程是典型的二阶常微分方程，所以运用四级四阶Runge-Kutta 方法[30]求解该方程。

为了验证本文追迹算法的准确性，本文对螺旋线进行光线追迹[30]，该光线的表达式如下：

其中n0和α为常系数，表征了径向变折射率介质中折射率分布的变化情况。上述径向变折射率分布的介质是有解析解的，其光线轨迹的解析解表达式为：

其中，(x0,y0,z0)为光线初射点的坐标，l0、p0、q0分别为初始光点的折射率在对应的x、y、z轴的光学方向余弦的乘积。图4 展示了光线追迹结果和解析解的对比情况，追迹得到的光线和理想光线的轨迹基本重合，所以本文将使用该追迹算法模拟光线在流场中的传播。

图4 光线追迹算法验证Fig.4 Verification of the ray tracing algorithm

对于实际的光在流场中的传播，初始平面光波穿过流场后不再是平面，这种实际波面与理想波面之间的偏差称为波前像差。波前像差可以由一系列多项式的线性组合来表示。通常描述波前像差的多项式为Zernike 多项式[25]。光学系统像差、大气湍流像差等静态和动态像差都可以用Zernike多项式来表示。Zernike 多项式每一项都有明确的物理意义，并且在单位圆内正交。关于Zernike 多项式的描述都在极坐标下进行。

为了验证Zernike 多项式拟合波前的合理性，图5 给出了波前拟合Zernike 多项式前36 项的图像和原始波前图像的对比，其中Nx和Ny代表网格点数。两者的误差在 1×10−10量级，相比波前OPD 在1×10−7量级，误差在千分之一，该误差可以忽略不计。

图5 原始波前和Zernike 拟合后的OPD 云图对比Fig.5 Comparison of OPD contours between the original wavefront and the Zernike polynomial fitted wave-front

在后续的研究中将使用OPD的均方根 OPDrms评估气动光学效应。在自适应光学中，通常可以消除平移、倾斜、离焦和像散等项来减弱气动光学效应。所以本文采用自适应光学的方法，利用波前 OPD拟合36 项Zernike 多项式，并去掉平移、倾斜、离焦和像散项来消除低阶项对气动光学效应的影响，并将去掉低阶项后的OPD 的均方根表示为。

1.4 光波空间自由传输理论

波阵面穿过湍流尾迹后，波前相位畸变通常表现为一阶或更高，而量化远场畸变的方法是直接计算远场光强图。在气动光学区域，不均匀的传播介质使得光线传输产生了光程差（OPD(r)），考虑光的波动性，光程差将使得波阵面上产生相位差，导致相位畸变，其表达式如下：

式中，λ表示光线波长，∆φ表示光线在z平面的相位差。流场引起的光波相位畸变会严重影响光波的复振幅分布，光线的光波复振幅为E，则有：

假设光线经过的远场大气光学区域介质均匀，则光线的传播遵循标量波动方程[18]。由波动光学可知，在强光或受到强扰动的条件下，光波的独立传播原理将不再满足，光的衍射等波动现象将严重影响光波的光强辐照度（optical intensity）分布。其中光强为复振幅模的平方，即：

2 计算结果及其分析

2.1 流场计算结果

为了验证网格拓扑方法和CFD 计算方法的准确性，本文对1 000 个时间步长内的流场进行了统计平均，其在物理时间上对应0.01 s。同时，每隔10 个时间步长统计一次瞬时流场，总共统计100 个瞬时流场用于后续气动光学效应的计算。本文选取了文献[31]中共形窗口转塔的流场实验结果与本文数值模拟结果进行对比验证。数值模拟条件与实验条件相同，马赫数为Ma∞=0.4，雷诺数为Re∞=1.43×106，其余计算条件为标准海平面条件。图6 展示了本文计算的共形窗口转塔顶部半球中心截面压力系数和实验值的对比。本文的数值模拟结果在附着流区域和实验值吻合较好；而在分离流区域，压力系数与实验值还存在微小差距。这是因为本文采用的是DDES 算法，在壁面位置处为了减少计算量采用RANS 求解，在分离区对流场壁面的解析能力不足，导致了分离区位置处的压力系数与实验值有所差距；但是DDES 相较于LES 可以大大减少计算量、计算时间，同时本文数值模拟结果与实验值的误差在可接受范围内，所以仍然采用该种方法。

图6 共形窗口转塔顶部半球中心截面的时均Cp 分布Fig.6 Time-averaged pressure coefficient distribution along the central plane of the conformal window turret dome

2.2 近场气动光学效应

为了研究共形光学窗口和平面光学窗口对气动光学效应的影响，本文分析了两种光学窗口分别在0°、90°和180°三个转角下的气动光学效应。图7（a）展示了光学网格示意图，光学传输的坐标轴分别为（x′,y′,z′），光学网格尺寸为 0.9R×0.9R×4R，光学网格数为 95×95×151，光学网格的解析依赖于流场的解析结果。图7（b）展示了光线的传输过程，光线从转塔内部发射，经过近场气动光学区域，得到畸变的波前分布，然后远场衍射至接收平面。

图7 光线在流场中的传输过程Fig.7 Transmission of light in the flow field

基于Gladstone-Dale 关系式（式（6）），图8 展示了两种光学窗口转塔在不同转角下的折射率分布云图。在本文的计算中，光波波长选取 λ=0.75 μm，则Gladstone-Dale 系数KGD=2.27×10−4m3/kg。为了更加真实地模拟光线发射过程，本文将光线追迹的起始位置设置在转塔内部，如图8 所示，其中红色矩形框表示光学网格。在流场区域内的光学网格使用当地的折射率插值，而在流场区域外的光学网格使用自由来流的折射率插值。为了考虑湍流脉动效应的影响，本文将计算100 个瞬时流场的OPD，再将其平均，并计算平均后的OPD 的均方根 OPDrms。

图8 不同转角下的中心截面瞬时折射率分布Fig.8 Instantaneous refraction distributions in the central plane for different angles (the optical grids are indicated by red rectangles)

图9 展示了共形光学窗口和平面光学窗口在转角为0°时的 OPDrms和随着追迹距离的变化趋势。图9（a）表明随着追迹距离的增加，两种光学窗口的 OPDrms都是先增加到1.5R处，然后再缓慢下降。但是，共形光学窗口的 OPDrms始终大于平面窗口的OPDrms；且随着追迹距离的增加，两者的差距还在不断增大。这是因为共形窗口本身的曲面形状会对气动光学效应有一定的影响，而平面窗口的形状则对气动光学效应没有影响。图9（b）则表明去掉平移、倾斜、离焦和像散等低阶项后，两种光学窗口的随着追迹距离的增加而增加至1.5R位置，后保持不变，且相差不大。结合图8（a）和图8（d），光线发射路径都经过了转塔后部的分离区，由于在该仰角下，共形转塔和平面转塔的分离位置相近，所以后部分离区也较为相似，从而使得两者的折射率分布差异较小。因此，在去掉低阶项后两者的气动光学效应相近。

图9 两种光学窗口在0°转角下的O PDrms和 Fig.9 O PDrms and of the two optical windows at the angle of 0°

图10 展示了共形光学窗口和平面光学窗口在转角为90°时的 OPDrms和随着追迹距离的变化趋势。图10（a）表明随着追迹距离的增加，平面光学窗口的 OPDrms随着追迹距离的增加先增加至1.5R处，之后缓慢减小；而共形光学窗口随着追迹距离的增加先增加后保持不变，但是整体上平面光学窗口的OPDrms大于共形窗口的 OPDrms。图10（b）则表明去掉平移、倾斜、离焦和像散后，仍旧是平面窗口气动光学效应更强。对比图8（b）和图8（c），平面窗口转塔附近的折射率梯度明显大于共形窗口转塔的折射率变化率，所以在此转角下，共形窗口的气动光学效应更小。

图10 两种光学窗口在90°转角下的O PDrms和 Fig.10 O PDrms and of the two optical windows at the angle of 90°

图11 展示了共形光学窗口和平面光学窗口在转角为180°时的 OPDrms和。图11（a）表明两种光学窗口的 OPDrms随着追迹距离的增加一直增大。结合图8（a～f），此时光线经过非均匀折射率的流场区域远远大于0°和90°转角下的非均匀流场区域；而且该转角下折射率梯度也更大，导致了该转角下OPDrms比前两个转角下的 OPDrms都更大。此时共形光学窗口的 OPDrms大于平面窗口的 OPDrms。对比图8（c）和图8（f），在光学窗口附近共形窗口转塔的折射率梯度大于平面窗口转塔；加之共形窗口形状本身对光线传输的影响，两者共同作用导致了共形窗口的 OPDrms更大。在图11（b）中，去掉平移、倾斜、离焦和像散等低阶项后，两种光学窗口的都比0°和90°转角下的小。这是因为在180°转角下都是附着流，脉动效应较小，对气动光学效应的低阶项影响较大，而对高阶项影响较小。同时，去掉低阶项后在该角度下平面窗口的则会更大。

图11 两种光学窗口在180°转角下的O PDrms和 Fig.11 O PDrms and of the two optical windows at the angle of 180°

表1 不同转角下追迹末端的高阶项OPD'rmsTable 1 OPD'rms at the end of the tracing for different angles

2.3 远场衍射结果

图12（c～h）展示了两种光学窗口在不同转角下无自适应光学矫正时，光线远场衍射30 000R距离（真空）下的光强分布云图；并将其与光学窗口出射处光线（图12（a））以及未畸变光线传播相同距离的理论光强分布（图12（b））进行了比较。在0°转角下，平面光学窗口和共形光学窗口发射的光线都出现了光斑的位置偏移，但是平面光学窗口还有明显的光斑分散；在90°转角下，平面光学窗口和共形光学窗口发射的光线也都出现了光斑的位置偏移，共形光学窗口还出现了光斑的抖动现象；而在180°转角下，两种光学窗口都出现了显著的光斑偏移。

图12 两种光学窗口不同转角下（0°、90°、180°）的远场衍射光强分布云图Fig.12 Far-field diffraction intensity contours at the angles of 0°,90° and 180° for the two optical windows

图13 展示了在与图12 相同的衍射距离下，不同光学窗口的远场光强分布。在0°转角下，两种光学窗口的光强峰值都出现了一定的位置偏移，但共形光学窗口的光强峰值大于平面光学窗口；在90°转角下，平面光学窗口的光强峰值发生了较严重的位置偏移且光强峰值有所增大，而共形光学窗口则是出现了光强峰值的减小；在180°转角下，两种光学窗口的光强峰值都发生了大幅的位置偏移和光强数值的增加。值得注意的是，在转角为90°和180°时，畸变波面远场衍射后的光强大于未畸变波面远场衍射的光强。根据波动光学中的惠更斯原理，畸变波面可能使得子波向一方偏折和聚集，导致光强峰值增加和位置偏移。

图13 不同转角下（0°、90°、180°）的远场衍射光强曲线图Fig.13 Far-field diffraction intensity profiles at the angle of 0°,90° and 180°

3 结论

在本文的研究中，无自适应光学矫正下，共形光学窗口和平面光学窗口在各个角度下的气动光学效应各有优劣。然而在利用自适应光学方法去掉平移、倾斜、离焦和像散等低阶项后，流场的脉动效应决定了其高阶项的大小，共形窗口和平面光学窗口在0°转角下发射的光线经过了相似的分离区域，导致了此时两种光学窗口的气动光学效应相近。而在90°和180°转角下，平面光学窗口的平面形状加剧了附近流场的脉动效应，导致了平面光学窗口的气动光学效应远大于共形光学窗口。在3 种不同转角下，光线传输区域位于分离区时（0°）的气动光学效应大于位于附着流区域时（90°、180°）的气动光学效应。综上，共形光学窗口发射的光线的传输性能优于平面光学窗口。

对两种光学窗口发射的光线进行远场衍射计算，结果表明经过近场畸变的光线远场自由衍射后可能比未畸变的光线远场自由衍射的光强峰值更大，且此时峰值位置偏移也更为严重。