摘要:对于含激波等大密度脉动结构气动光学流场, 沿着直线路径对折射率积分会带来较大的光程误差. 因此, 数值求解光线方程进行光线追迹是必要的. 与数值求解光线方程不同, 元胞自动机(cellular automata, CA)通过给定光线位置和方向变换规则来模拟光线在介质中传输路径. 本文基于已有的实验测量、数值仿真所获得的高超声速流场密度场数据, 分别采用数值求解光线方程法和CA光线追迹算法进行光线追迹, 进而得到光线出射流场后的光程差. 结果表明, CA算法对于二维气动光学流场中光线追迹的适用性, 且较数值求解光线方程方法具有更高的效率.
关键词: 气动光学/
元胞自动机/
光线追迹/
光程差

English Abstract

全文HTML

--> --> -->

气动光学是研究稳态或非稳态气体流场与在其中传输的光束相互作用及传输规律的学科^[1,2]. 流场中的激波、边界层/剪切层和尾流等大梯度、非均匀折射率场分布会使在其中传输的光束的波前产生畸变, 这就是气动光学效应^[3,4].
对于折射率脉动是微小量(10^–6)的流场, 流场特征长度远大于传输光波长, 因此光束在流场中的传输可认为满足傍轴近似条件, 进而可认为光线在流场中的传输路径近似为直线并忽略流场对振幅的衰减, 只需考虑非均匀流场对相位的畸变效应. 因此在这种情况下, 可以对直线光路径上各点的折射率直接积分获取光程(optical path length, OPL). 但是对于大脉动量(由非定常流动和湍流大尺度结构产生)流场^[5,6], 例如超声速光学头罩绕流流场, 由于流场激波的存在, 简单的对折射率场沿直线路径进行积分计算OPL会存在较大误差. 因此就需要光线追迹得到光线在流场中的完整传输路径, 综合多条光线的信息得到光波穿过流场后的OPL值, 进而得到光程差(optical path difference, OPD)、斯特列尔比(Strehl ratio, SR)和光学传递函数(optical transmission function, OTF)等光学量并据此分析波前畸变.
Montagnino^[7]在1968年通过泰勒级数展开法数值求解了光线方程, 为光在非均匀介质中传输的数值计算提供了基础. 近年来, 一些光线追迹的新方法不断涌现. Chang等^[8]基于三维Snell定律提出了一种在三维气动光学流场中进行光线追迹的方法, 该方法能够有效地计算气动光学的波前参数. Xu等^[9]开发了一种反向光线追迹的方法, 该方法以探测器为追迹初始位置, 以远距离的目标为追迹末位置, 可明显简化传统的气动光学光线追迹计算.
本文提出了一种基于元胞自动机(cellular automata, CA)的气动光学流场光线追迹方法. CA是结构简单但具有复杂自组织特性的离散动力学系统, 是包含大量相同组分在局部相互作用的复杂自然系统的数学模型. 用CA来模拟一个物理过程的优点在于避免了用微分方程作为控制方程, 而直接通过制定变换规则来模拟非线性物理现象^[10]. CA模型中, 空间被离散成许多元胞, 这些元胞只在有限的状态集中取值, 元胞的状态更新遵循一定的变换法则. CA目前已被广泛应用于如流体流动^[11]、模式识别^[12,13]、相位变化^[14]和人口动力学^[15,16]等多学科的研究中. 对于气动光学光线追迹, 光线在流场中传输可以看成是光子在满足变换规则网格间的传输, 因此如果能通过制定合理的光子移动规则, 利用CA可实现气动光学流场中光线追迹计算.
本文基于纳米粒子示踪平面激光散射技术(nano-tracer-based planar laser scattering, NPLS)^[17,18]实验得到的超声速二维剪切层流场和脱体涡(detached eddy simulation, DES)^[19]模拟得到的含激波的超声速光学头罩绕流流场的二维流场, 计算、对比了CA光线追迹算法与数值求解光线方程光线追迹法(numerical solving the ray equation, NSRE)得到的光线路径以及对折射率沿路径积分得到的OPL数据, 通过OPL得到了OPD的曲线. 结果表明, CA算法对二维随机流场光线追迹非常有效, 并具有相对于光线方程数值求解方法更高的计算效率.

任意折射率分布介质中的光线传输路径可用光线方程进行描述

$\frac{{\rm{d}}}{{{\rm{d}}s}}\left( {n\frac{{{\rm{d}}{r}}}{{{\rm{d}}s}}} \right) = \nabla n,$

其中, s为光线传输路径上的弧长, r为光线传输路径的位置矢量, n为折射率, $\nabla n$为折射率梯度. 该方程仅当流场折射率满足特定分布时才有解析解. 对于折射率随机变化的气动光学流场, 光线传输路径的计算通常采用NSRE算法, 该方法的基本思想是: 对过定点的一条光线矢量, 基于光线方程, 依次计算出下一个光线矢量^[20,21], 进而得到流场区域内光线的路径. 该方法能有效地处理激波处的折射率突变, 具有较高的精度^[22,23]. 然而由于该方法在处理非均匀网格时首先需要进行插值运算, 且每一次迭代都要重置初始值并求解非线性微分方程, 因此存在计算量大的不足. 本文所采用的数值求解光线方程的差分方法为三阶龙格库塔算法^[24], 该算法主要通过求解光线位置和方向同时改变的复杂差分方程来获得光线完整路径.
CA算法中, 当光线传输至某元胞(网格)时, 该元胞的状态为1, 反之则为0. 邻居类型采用的是摩尔型邻居^[25], 当前时刻的元胞状态取决于上一时刻及其邻居元胞的状态. 元胞的坐标和折射率可以分别用该元胞的节点坐标(j, k)和折射率n_{j, k}来表示. 光子在元胞间移动遵循一定位置变换规则和方向变换规则. 位置变换规则用来获得光线矢量矢端点的位置信息和判断本次迭代后是否跨越网格, 如果跨越网格, 则需要根据方向变换规则计算光线偏转角, 反之则保持方向进入下一次位置变换. 可见, CA追迹算法与差分数值求解光线方程追迹法最大的不同在于: 将光线位置的改变与光线方向的改变拆分为两个相互独立的过程. 避免了光线位置改变和方向改变的同时计算, 因此编程更容易实现且具有较高的效率. CA光线追迹算法的具体步骤如下:
1)给定初始光线矢量${{r}_0} = (0, 0)$. ${{V}_1}$表示初始元胞(j₀, k₀)处归一化光速, $\left| {{{V}_1}} \right| = 1$. θ_x,1 和 θ_y,1分别表示${{V}_1}$和x, y轴之间的夹角.
2)位置变换规则: 对于第i次迭代, 光线矢量${{r}_i} = {{r}_{i - 1}} + \Delta {{r}_i} \left( {i = 1, 2, 3 \ldots } \right)$, 其中$\Delta {{r}_i} = {{V}_i}\Delta t$($\Delta t = 1$为迭代时间步长)表示每次迭代光线矢量${{r}_i}$的改变量; ${{V}_i} = (\cos {\theta _{x, i}}{\hat j} + \cos {\theta _{y, i}}{\hat k})/{n_{j, }}_k\left(\left| {{{V}_i}} \right| = 1/{n_{j, }}_k \right),$ ${\theta _{x, i}}$和${\theta _{y, i}}$分别表示${{V}_i}$和x, y轴之间的夹角. ${r_{x, i}}$和${r_{y, i}}$为${{r}_i}$在x, y轴方向的分量. 定义$\left[ {{r_{x, i}}} \right]$和$\left[ {{r_{y, i}}} \right]$分别表示对${r_{x, i}}$和${r_{y, i}}$向整数取整, 并且定义${r_{{\rm{c}}x}} = \left[ {{r_{x, i}}} \right] - \left[ {{r_{x, i - 1}}} \right]$和${r_{{\rm{c}}y}} = \left[ {{r_{y, i}}} \right] - \left[ {{r_{y, i - 1}}} \right]$来判断每一次迭代后是否跨越网格, 从而判断${{V}_{i + 1}}$的大小和方向. 第i次迭代后:

$\left| {{{V}_{i + 1}}} \right| = \begin{cases}{1/{n_{j,k}},}&{\left( {{r_{{\rm{c}}x}} = 0\;{{\rm{and}}}\;{r_{{\rm{c}}y}} = 0} \right)}\\{1/{n_{j + 1,k}},}&{\left( {{r_{{\rm{c}}x}} = 1\;{\rm{and}}\;{r_{{\rm{c}}y}} = 0} \right)}\\{1/{n_{j,k + 1}},}&{\left( {{r_{{\rm{c}}x}} = 0\;{\rm{and}}\;{r_{{\rm{c}}y}} = 1} \right)}\\{1/{n_{j + 1,k + 1}},}&{\left( {{r_{{\rm{c}}x}} = 1\;{\rm{and}}\;{r_{{\rm{c}}y}} = 1} \right)}\end{cases}.$

如(2)式所示, 当${r_{{\rm{c}}x}}$和${r_{{\rm{c}}y}}$中有一个为1时, 表明元胞路径向该方向推进一个网格; 当${r_{{\rm{c}}x}}$和${r_{{\rm{c}}y}}$同时为1时, 表明元胞路径向倾斜方向推进一个网格. 上述两种情况需要根据(5)式来计算偏角${\theta _{x, i + 1}}$和${\theta _{y, i + 1}}$, 否则, ${{V}_{i + 1}}$的大小和方向保持不变.
3)方向变换规则: 认为s傍近元胞路径L, 因此(1)式可以写为

$\frac{{\rm{d}}}{{{\rm{d}}s}}\left( {n\frac{{{\rm{d}}{r}}}{{{\rm{d}}L}}} \right) = \nabla n,$

进而(3)式可写成x和y方向的分量的形式:

$\left\{ {\begin{aligned}&{\frac{{{\rm{d}}{n_{j,k}}}}{{{\rm{d}}{L_{x,i}}}}\frac{{{\rm{d}}x}}{{{\rm{d}}{L_{x,i}}}} + {n_{j,k}}\frac{{{{\rm{d}}^2}x}}{{{\rm{d}}s{\rm{d}}{L_{x,i}}}} = \frac{{\partial {n_{j,k}}}}{{\partial x}},}\\&{\frac{{{\rm{d}}{n_{j,k}}}}{{{\rm{d}}{L_{y,i}}}}\frac{{{\rm{d}}y}}{{{\rm{d}}{L_{y,i}}}} + {n_{j,k}}\frac{{{{\rm{d}}^2}y}}{{{\rm{d}}s{\rm{d}}{L_{y,i}}}} = \frac{{\partial {n_{j,k}}}}{{\partial y}},}\end{aligned}} \right.$

其中, ${L_{x, i}}$和${L_{y, i}}$分别表示L在x和y方向的分量. 对(4)式两边同时积分

$\left\{ {\begin{aligned}&{\frac{{{\rm{d}}x}}{{{\rm{d}}s}} \!=\! \int \!{\frac{1}{{{n_{j,k}}}}\!\left(\! {\frac{{\partial {n_{j,k}}}}{{\partial x}}\! -\! \frac{{{\rm{d}}{n_{j,k}}}}{{{\rm{d}}{L_{x,i}}}}\frac{{{\rm{d}}x}}{{{\rm{d}}{L_{x,i}}}}} \!\right){\rm{d}}{L_{x,i}}\!=\!\cos \Delta {\theta _{x,i}},} }\\&{\frac{{{\rm{d}}y}}{{{\rm{d}}s}} \!=\! \int \!{\frac{1}{{{n_{j,k}}}}\!\left(\! {\frac{{\partial {n_{j,k}}}}{{\partial y}} \!-\! \frac{{{\rm{d}}{n_{j,k}}}}{{{\rm{d}}{L_{y,i}}}}\frac{{{\rm{d}}y}}{{{\rm{d}}{L_{y,i}}}}} \!\right){\rm{d}}{L_{y,i}}\!=\!\cos \Delta {\theta _{y,i}},} }\end{aligned}} \right.$

其中, $\Delta {\theta _{x, i}} = {\theta _{x, i + 1}} - {\theta _{x, i}}$和$\Delta {\theta _{y, i}} = {\theta _{y, i + 1}} - {\theta _{y, i}}$分别表示计算后得出的在x和y方向的偏转角的改变量. (5)式中的$\partial {n_{j, k}}/\partial x = ({n_{j + 1, k}} - {n_{j - 1, k}})/$$2{h_{x, i}}$和$\partial {n_{j, k}}/\partial y = ({n_{j, k + 1}} - {n_{j, k - 1}})/2{h_{y, i}}$采用中心差分来计算, 其中h_{x, i}和h_{y, i}表示x和y方向相邻元胞的间距. 对于CA算法, dx (dy)和Lx (Ly)之间的位置关系只有平行、垂直和成45°角, 因此(5)式中的${\rm{d}}x/{\rm{d}}{L_{x, i}}$和${\rm{d}}y/{\rm{d}}{L_{y, i}}$利用三角关系来计算; ${\rm{d}}{n_{j, k}}/{\rm{d}}{L_{x, i}}$和${\rm{d}}{n_{j, k}}/{\rm{d}}{L_{y, i}}$分别由${\rm{d}}{n_{j, k}}/{\rm{d}}{L_{i, x}} = ({n_{j + 1, k}} - {n_{j, k}})/{h_{x, i}}$和${\rm{d}}{n_{j, k}}/{\rm{d}}{L_{i, y}} = ({n_{j, k + 1}}- {n_{j, k}})/ {h_{y, i}}$ 计算得到.
4)根据上述变换规则, 最终的光线轨迹可由${{r}_N} = {{r}_0} + \displaystyle\sum\nolimits_{i = 1}^N {\Delta {{r}_i}}$得到.

2

-->

3.1.超声速混合层流场

本文基于NPLS技术获得高分辨率超声速混合层和边界层二维密度场分布. NPLS系统主要由光源、示踪粒子发生器、CCD、风洞、同步控制器和数据采集模块等组成^[26]. 超声速混合层流场来源于超声速混合层风洞中隔板上下马赫数不同的超声速气流混合^[27]; 超声速边界层流场来源于超声速风洞靠管壁的壁面薄层. 利用NPLS技术获得超声速剪切层的密度场分布的关键在于使示踪粒子的散射光准确地描述密度场. 风洞及NPLS系统的主要参数、完整结构及实物图易仕和等^[28]和Tian等^[29]已有详细的论述, 这里不再赘述.
实验测量的超声速混合层流场的对流马赫数为0.5, 其中两股来流的马赫数分别为3.509和1.400, 对应的来流速度分别是654.7和421.1 m/s, 属于中等可压缩流场. NPLS图像的像素尺寸为1431 × 281, 单像素分辨率h = 0.16 mm, 入射光波长λ = 1064 nm, 对应的K_GD = 2.195 × 10^–4 m³/kg, 跨帧间隔为15 μs, 样本总数为50帧. 流向为x正向, 光线传输方向为y正向. 实验测量的超声速边界层流场的马赫数为3.0, 对应的来流速度为622.5 m/s. NPLS图像的像素尺寸为2048 × 1436, 单像素分辨率为h = 0.013 mm, 入射光波长λ = 0.5 μm, 对应的K_GD = 2.2 × 10^–4 m³/kg, 跨帧间隔为15 μs, 样本总数为30帧. 流向为x正向, 光线传输方向为y正向. NPLS技术拍摄的某时刻超声速混合层和边界层图像如图1所示.
图 1 NPLS获得的超声速剪切层图像　(a)混合层; (b)边界层
Figure1. The flow visualization results of supersonic shear layers obtained by NPLS: (a) Supersonic mixing layer; (b) supersonic boundary layer.

2 -->

3.2.超声速光学头罩绕流流场

在本文中, 依据模型尺寸参数, 采用基于SST两方程湍流模型的DES模拟来得到超声速光学头罩绕流流场密度场分布^[30]. DES模拟在近物面采用雷诺平均方法(Reynolds averaged Navier-Stokes, RANS), 在其他区域采用大涡模拟方法(large eddy simulation, LES), 兼具前者计算量小的优点和后者能分离湍流流动的优势, DES湍流模型方程和RANS控制方程采用全耦合求解^[31].图2所示为无冷却喷流光学头罩绕流流场归一化密度分布, 相关参数取值为海拔高度10 km, 攻角为0°, 马赫数为6.0. 对密度做了无量纲处理, 假设ρ表示真实密度, 则图中的无量纲密度值可由ρ/0.41270得到. 其中仿真参数为光波波长1 μm, 初始位置位于均匀来流区域. 在仿真计算时, 设定光线传输方向为z方向且垂直于光学窗口入射. 计算时, 沿着平行和垂直于光学窗口方向取二维截面, 尺寸为600 mm × 50 mm.
图 2 光学头罩绕流流场密度场
Figure2. The density field distribution of supersonic flow field surrounding the optical dome obtained by DES.

在混合层流场密度分布变化较大的区域(图1(a)中红框区域), 使用CA和NSRE算法获得的光线路径如图3所示. CA算法和NSRE算法计算得到的出射该区域后的偏角分别为31.1和30.7 μrad.
图 3 CA算法与NSER算法在混合层得的光线路径
Figure3. Beam paths obtained by CA and NSRE in mixing layer.

在气动光学中, 折射率n和密度ρ之间的关系可表示为

$n = 1 + {K_{{\rm{GD}}}}\rho ,$

其中, ${K_{{\rm{GD}}}}$为Gladstone-Dale系数. 对光线沿着路径r上的折射率积分得到OPL

${\rm{OPL}} = \int\nolimits_r {n{\rm{d}}r} ,$

由表达式

${\rm{OPD}} = {\rm{OPL}} - \left\langle {{\rm{OPL}}} \right\rangle, $

可得到光在穿过流场后的光程差. (8)式中$\left\langle {{\rm{OPL}}} \right\rangle $代表空间平均. OPD的计算结果如图4所示.
图 4 CA算法与NSRE算法计算得到的OPD　(a)混合层; (b)边界层; (c)含激波的超声速光学头罩二维剖面流场
Figure4. The OPD results calculated by CA and NSRE: (a) Supersonic mixing layer; (b) supersonic boundary layer; (c) supersonic flow field surrounding the optical dome.

根据统计学中均方根误差的定义

${\rm{OP}}{{\rm{D}}_{{\rm{rms}}}} = \sqrt {\frac{1}{M}\sum\limits_{\tau = 1}^M {{{\left( {{\rm{OP}}{{\rm{D}}_{{\rm{C}}{{\rm{A}}_\tau }}} - {\rm{OP}}{{\rm{D}}_{{\rm{NSR}}{{\rm{E}}_\tau }}}} \right)}^2}} } ,$

其中M表示取样点总个数, $\tau $表示取样点, ${\rm{OP}}{{\rm{D}}_{{\rm{NSRE}}}}$和${\rm{OP}}{{\rm{D}}_{{\rm{CA}}}}$分别表示NSRE算法和CA算法计算的OPD. 超声速混合层、超声速边界层以及超声速光学头罩二维剖面流场的OPD_rms分别为0.0086 λ₁, 0.0014 λ₂和0.0146 λ₃, 其中λ₁, λ₂和λ₃分别对应超声速混合层、超声速边界层以及超声速光学头罩二维剖面流场的波长. 结果表明: 对于超声速混合层、超声速边界层以及超声速光学头罩二维剖面流场, CA计算得到的OPD结果与NSRE算法计算结果具有较好的一致性.
上述结果显示了CA算法在气动光学流场计算中的高精度以及适用性. 然而衡量一种算法的优劣除了其计算结果的正确性, 计算效率也是一个重要的评价指标. tic和toc函数是MATLAB中的内置函数, tic用来保存当前时间, toc用来记录程序完成时间. 两个函数的配套使用可以计算算法的实际运行时间, 从而评价算法的效率. 对于本文的光线追迹, 计算光线传输一个网格的程序运行时间, 混合层、边界层和高速绕流流场的程序执行时间分别如表1、表2和表3所列. 表中数据均保留三位小数. 表中, ${t_1}-{t_5}$表示每种流场采用CA或NSRE算法的计算时间, ${t_{\rm{a}}}$表示对计算时间求平均值, K = 5表示计算总次数, $\varsigma $表示计算次数, ${t_{\rm{a}}}$可表示为

方法	t1/s	t2/s	t3/s	t4/s	t5/s	ta/s
CA	2.032	2.041	2.157	2.085	2.077	2.078
NSRE	8.463	8.470	8.522	8.491	8.511	8.491

表1CA与NSRE算法计算混合层流场的程序执行时间
Table1.The program running time of CA and NSRE in mixing layer.

方法	t1/s	t2/s	t3/s	t4/s	t5/s	ta/s
CA	1.966	2.002	1.983	1.964	1.972	1.977
NSRE	7.241	7.250	7.143	7.250	7.219	7.220

表2CA与NSRE算法计算边界层流场的程序执行时间
Table2.The program running time of CA and NSRE in boundary layer.

方法	t1/s	t2/s	t3/s	t4/s	t5/s	ta/s
CA	2.831	2.884	2.827	2.873	2.840	2.851
NSRE	11.375	11.525	11.401	11.380	11.397	11.416

表3CA与NSRE算法计算高速绕流流场的程序执行时间
Table3.The program running time of CA and NSRE in supersonic flow field surrounding the optical dome.

${t_{\rm{a}}} = \frac{1}{K}\sum\limits_{\varsigma = 1}^K {{t_\varsigma }} .$

表1—表3中的计算时间为两种算法程序在配置为Inter(R) Core(TM) i9-9900k CPU @ 3.6 GHz(内存32 GHz)的电脑上运行得到的. 由表1—表3可见, NSRE算法的程序平均执行时间约为CA算法的4倍. 从计算模拟时间来看, 在相同的条件下, CA算法的计算效率要高于NSRE算法. 从算法结构上来看, 在每一次迭代周期, CA算法的加法和乘法次数分别为11次和7次, NSRE算法的加法和乘法次数分别为39次和32次, NSRE算法的加法和乘法次数约为CA算法的4倍, 这与表1—表3中的时间之间的关系符合得很好.

本文研究并对比了CA光线追迹算法和NSRE算法计算二维超声速气动光学流场OPD的结果. 计算结果显示, 对于二维超声速NPLS流场和二维超声速绕流流场, CA算法的OPD计算结果与NSRE算法计算结果符合较好, 具有很好的计算精度; 从算法效率角度来看, CA的执行时间约为NSRE算法的1/4, 具有比NSRE更高的效率. 研究结果表明, CA光线追迹算法对于二维流场的气动光学计算是适用的, 这为离散非均匀流场的气动光学计算提供了新方案.
感谢国防科技大学空天科学学院易仕和教授、中国空气动力研究与发展中心陈勇研究员在NPLS实验数据及流场CFD数据等方面提供的支持与帮助.