搜索

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于元胞自动机的气动光学光线追迹算法

雒亮 夏辉 刘俊圣 费家乐 谢文科

基于元胞自动机的气动光学光线追迹算法

雒亮, 夏辉, 刘俊圣, 费家乐, 谢文科
PDF
HTML
导出引用
  • 对于含激波等大密度脉动结构气动光学流场, 沿着直线路径对折射率积分会带来较大的光程误差. 因此, 数值求解光线方程进行光线追迹是必要的. 与数值求解光线方程不同, 元胞自动机(cellular automata, CA)通过给定光线位置和方向变换规则来模拟光线在介质中传输路径. 本文基于已有的实验测量、数值仿真所获得的高超声速流场密度场数据, 分别采用数值求解光线方程法和CA光线追迹算法进行光线追迹, 进而得到光线出射流场后的光程差. 结果表明, CA算法对于二维气动光学流场中光线追迹的适用性, 且较数值求解光线方程方法具有更高的效率.
      通信作者: 谢文科, wenkexiedan@163.com
    • 基金项目: 装备预研领域基金(批准号: 6140415020311)、高能激光技术湖南省重点实验室开放基金(批准号: GNJGJS04)和湖南省光电惯性工程技术研究中心开放基金(批准号: HN-NUDT1908)资助的课题
    [1]

    Gordeyev S, Jumper E 2010 Prog. Aerosp. Sci. 46 8

    [2]

    谢文科, 刘俊圣, 费家乐, 周全, 夏辉, 陈欣, 张盼, 彭一鸣, 于涛 2019 物理学报 68 094202

    Xie W K, Liu J S, Fei J L, Zhou Q, Xia H, Chen X, Zhang P, Peng Y M, Yu T 2019 Acta Phys. Sin. 68 094202

    [3]

    Pond J E, Sutton G W 2006 J. Aircraft 43 3

    [4]

    Ding H L, Yi S H, Zhu Y Z, He L 2017 Appl. Opt. 56 27

    [5]

    于涛, 夏辉, 樊志华, 谢文科, 张盼, 刘俊圣, 陈欣 2018 物理学报 67 134203

    Yu T, Xia H, Fan Z H, Xie W K, Zhang P, Liu J S, Chen X 2018 Acta Phys. Sin. 67 134203

    [6]

    Zhu K C, Li S X, Tang Y, Yu Y, Tang H Q 2012 J. Opt. Soc. Am. A 29 3

    [7]

    Montagnino L 1968 J. Opt. Soc. Am. A 58 12

    [8]

    Chang X F, Wang T, Wan S Z, Yan J, Fu W X 2015 Optik 126 23

    [9]

    Xu L, Xue D T, Lv X Y 2018 Opt. Express 26 1

    [10]

    Tang L P, Tang L M , Wang D, Deng H X, Chen K Q 2018 J. Phys.: Condens. Matter 30 465301

    [11]

    Chen Q, Wang Y 2015 Physica A 432 15

    [12]

    Sun G Q, Jin Z, Song L P, Chakraborty A, Li B L 2011 Ecol. Res. 26 2

    [13]

    Chen C K, Li J, Zhang D 2012 Physica A 391 7

    [14]

    Ahmadpour S S, Mosleh M 2018 J. Supercomput. 74 9

    [15]

    Qin Y, Feng M Y, Lu H C, Cottrell G W 2018 Int. J. Comput. Vision. 126 751

    [16]

    Zhang H, Wei J, Gao X L, Hu J 2019 Int. J. Mod. Phys. C 30 5

    [17]

    朱杨柱, 易仕和, 孔小平, 全鹏程, 陈植, 田立丰 2014 物理学报 63 134701

    Zhu Y Z, Yi S H, Kong X P, Quan P C, Chen Z, Tian L F 2014 Acta Phys. Sin. 63 134701

    [18]

    Yi S H, He L, Zhao Y X, Tian L F, Cheng Z Y 2009 Sci. China, Ser. G 52 12

    [19]

    Zhu J, Li X L, Tang H Q, Zhu K C 2017 Opt. Express 25 17

    [20]

    Yu T, Xia H, Fan Z H, Xie W K, Zhang P, Liu J S, Chen X, Chu X X 2018 Opt. Commun. 436 1

    [21]

    Guo G M, Liu H 2017 Appl. Opt. 56 16

    [22]

    Ji B, Long Y, Long X P, Qian Z D, Zhou J J 2017 J. Hydrodyn. 29 1

    [23]

    Weghorst H, Hooper G, Greenberg D P 1984 Acm. T. Graphic. 3 1

    [24]

    Huang Y, Shi G D, Zhu K Y 2016 J. Quant. Spectrosc. Radiat. 176 24

    [25]

    Jiang H, Ren G, Zheng L, Cheng J X, Huang Z F 2014 Int. J. Mod. Phys. B 28 16

    [26]

    赵玉新 2008 博士学位论文(长沙: 国防科学技术大学)

    Zhao Y X 2008 Ph. D. Dissertation (Changsha: National University of Defense Technology) (in Chinese)

    [27]

    易仕和, 赵玉新, 田立丰, 何霖, 程忠宇 2009 空气动力学学报 27 114

    Yi S H, Zhao Y X, Tian L F, He L, Cheng Z Y 2009 Acta Aerodyn. Sin. 27 114

    [28]

    易仕和, 陈植, 何霖, 武宇, 田立丰 2014 实验流体力学 28 1

    Yi S H, Chen Z, He L, Wu Y, Tian L F 2014 J. Fluid. Mech. 28 1

    [29]

    Tian L F, Yi S H, Zhao Y X, He L, Chen Z Y 2009 Sci. China, Ser. G 52 9

    [30]

    Lyons D C, Peltier L J, Zajaczkowski F J, Paterson E G 2007 J. Fluids Eng.-Trans. ASME 131 11

    [31]

    Usta O, Korkut E 2018 Ocean Eng. 160 15

  • 图 1  NPLS获得的超声速剪切层图像 (a)混合层; (b)边界层

    Fig. 1.  The flow visualization results of supersonic shear layers obtained by NPLS: (a) Supersonic mixing layer; (b) supersonic boundary layer.

    图 2  光学头罩绕流流场密度场

    Fig. 2.  The density field distribution of supersonic flow field surrounding the optical dome obtained by DES.

    图 3  CA算法与NSER算法在混合层得的光线路径

    Fig. 3.  Beam paths obtained by CA and NSRE in mixing layer.

    图 4  CA算法与NSRE算法计算得到的OPD (a)混合层; (b)边界层; (c)含激波的超声速光学头罩二维剖面流场

    Fig. 4.  The OPD results calculated by CA and NSRE: (a) Supersonic mixing layer; (b) supersonic boundary layer; (c) supersonic flow field surrounding the optical dome.

    表 1  CA与NSRE算法计算混合层流场的程序执行时间

    Table 1.  The program running time of CA and NSRE in mixing layer.

    方法t1/st2/st3/st4/st5/sta/s
    CA2.0322.0412.1572.0852.0772.078
    NSRE8.4638.4708.5228.4918.5118.491
    下载: 导出CSV

    表 2  CA与NSRE算法计算边界层流场的程序执行时间

    Table 2.  The program running time of CA and NSRE in boundary layer.

    方法t1/st2/st3/st4/st5/sta/s
    CA1.9662.0021.9831.9641.9721.977
    NSRE7.2417.2507.1437.2507.2197.220
    下载: 导出CSV

    表 3  CA与NSRE算法计算高速绕流流场的程序执行时间

    Table 3.  The program running time of CA and NSRE in supersonic flow field surrounding the optical dome.

    方法t1/st2/st3/st4/st5/sta/s
    CA2.8312.8842.8272.8732.8402.851
    NSRE11.37511.52511.40111.38011.39711.416
    下载: 导出CSV
  • [1]

    Gordeyev S, Jumper E 2010 Prog. Aerosp. Sci. 46 8

    [2]

    谢文科, 刘俊圣, 费家乐, 周全, 夏辉, 陈欣, 张盼, 彭一鸣, 于涛 2019 物理学报 68 094202

    Xie W K, Liu J S, Fei J L, Zhou Q, Xia H, Chen X, Zhang P, Peng Y M, Yu T 2019 Acta Phys. Sin. 68 094202

    [3]

    Pond J E, Sutton G W 2006 J. Aircraft 43 3

    [4]

    Ding H L, Yi S H, Zhu Y Z, He L 2017 Appl. Opt. 56 27

    [5]

    于涛, 夏辉, 樊志华, 谢文科, 张盼, 刘俊圣, 陈欣 2018 物理学报 67 134203

    Yu T, Xia H, Fan Z H, Xie W K, Zhang P, Liu J S, Chen X 2018 Acta Phys. Sin. 67 134203

    [6]

    Zhu K C, Li S X, Tang Y, Yu Y, Tang H Q 2012 J. Opt. Soc. Am. A 29 3

    [7]

    Montagnino L 1968 J. Opt. Soc. Am. A 58 12

    [8]

    Chang X F, Wang T, Wan S Z, Yan J, Fu W X 2015 Optik 126 23

    [9]

    Xu L, Xue D T, Lv X Y 2018 Opt. Express 26 1

    [10]

    Tang L P, Tang L M , Wang D, Deng H X, Chen K Q 2018 J. Phys.: Condens. Matter 30 465301

    [11]

    Chen Q, Wang Y 2015 Physica A 432 15

    [12]

    Sun G Q, Jin Z, Song L P, Chakraborty A, Li B L 2011 Ecol. Res. 26 2

    [13]

    Chen C K, Li J, Zhang D 2012 Physica A 391 7

    [14]

    Ahmadpour S S, Mosleh M 2018 J. Supercomput. 74 9

    [15]

    Qin Y, Feng M Y, Lu H C, Cottrell G W 2018 Int. J. Comput. Vision. 126 751

    [16]

    Zhang H, Wei J, Gao X L, Hu J 2019 Int. J. Mod. Phys. C 30 5

    [17]

    朱杨柱, 易仕和, 孔小平, 全鹏程, 陈植, 田立丰 2014 物理学报 63 134701

    Zhu Y Z, Yi S H, Kong X P, Quan P C, Chen Z, Tian L F 2014 Acta Phys. Sin. 63 134701

    [18]

    Yi S H, He L, Zhao Y X, Tian L F, Cheng Z Y 2009 Sci. China, Ser. G 52 12

    [19]

    Zhu J, Li X L, Tang H Q, Zhu K C 2017 Opt. Express 25 17

    [20]

    Yu T, Xia H, Fan Z H, Xie W K, Zhang P, Liu J S, Chen X, Chu X X 2018 Opt. Commun. 436 1

    [21]

    Guo G M, Liu H 2017 Appl. Opt. 56 16

    [22]

    Ji B, Long Y, Long X P, Qian Z D, Zhou J J 2017 J. Hydrodyn. 29 1

    [23]

    Weghorst H, Hooper G, Greenberg D P 1984 Acm. T. Graphic. 3 1

    [24]

    Huang Y, Shi G D, Zhu K Y 2016 J. Quant. Spectrosc. Radiat. 176 24

    [25]

    Jiang H, Ren G, Zheng L, Cheng J X, Huang Z F 2014 Int. J. Mod. Phys. B 28 16

    [26]

    赵玉新 2008 博士学位论文(长沙: 国防科学技术大学)

    Zhao Y X 2008 Ph. D. Dissertation (Changsha: National University of Defense Technology) (in Chinese)

    [27]

    易仕和, 赵玉新, 田立丰, 何霖, 程忠宇 2009 空气动力学学报 27 114

    Yi S H, Zhao Y X, Tian L F, He L, Cheng Z Y 2009 Acta Aerodyn. Sin. 27 114

    [28]

    易仕和, 陈植, 何霖, 武宇, 田立丰 2014 实验流体力学 28 1

    Yi S H, Chen Z, He L, Wu Y, Tian L F 2014 J. Fluid. Mech. 28 1

    [29]

    Tian L F, Yi S H, Zhao Y X, He L, Chen Z Y 2009 Sci. China, Ser. G 52 9

    [30]

    Lyons D C, Peltier L J, Zajaczkowski F J, Paterson E G 2007 J. Fluids Eng.-Trans. ASME 131 11

    [31]

    Usta O, Korkut E 2018 Ocean Eng. 160 15

  • [1] 丁浩林, 易仕和, 朱杨柱, 赵鑫海, 何霖. 不同光线入射角度下超声速湍流边界层气动光学效应的实验研究. 物理学报, 2017, 66(24): 244201. doi: 10.7498/aps.66.244201
    [2] 朱杨柱, 易仕和, 陈植, 葛勇, 王小虎, 付佳. 带喷流超声速光学头罩流场气动光学畸变试验研究. 物理学报, 2013, 62(8): 084219. doi: 10.7498/aps.62.084219
    [3] 陈灿, 佟亚军, 谢红兰, 肖体乔. Laue弯晶聚焦特性的光线追迹研究. 物理学报, 2012, 61(10): 104102. doi: 10.7498/aps.61.104102
    [4] 吴海英, 张淳民, 赵葆常, 李英才. 改型Wollaston棱镜的光程差及其特性分析. 物理学报, 2009, 58(3): 1642-1647. doi: 10.7498/aps.58.1642
    [5] 张天天, 易仕和, 朱杨柱, 何霖. 基于背景纹影波前传感技术的气动光学波前重构与校正. 物理学报, 2015, 64(8): 084201. doi: 10.7498/aps.64.084201
    [6] 穆廷魁, 张淳民, 赵葆常. 偏振干涉成像光谱仪中Wollaston棱镜光程差及条纹定位面的精确计算与分析. 物理学报, 2009, 58(6): 3877-3886. doi: 10.7498/aps.58.3877
    [7] 张书赫, 梁振, 周金华. 运用四元数分析椭球微粒所受的光阱力. 物理学报, 2017, 66(4): 048701. doi: 10.7498/aps.66.048701
    [8] 永贵, 黄海军, 许岩. 菱形网格的行人疏散元胞自动机模型. 物理学报, 2013, 62(1): 010506. doi: 10.7498/aps.62.010506
    [9] 梁经韵, 张莉莉, 栾悉道, 郭金林, 老松杨, 谢毓湘. 多路段元胞自动机交通流模型. 物理学报, 2017, 66(19): 194501. doi: 10.7498/aps.66.194501
    [10] 花 伟, 林柏梁. 考虑行车状态的一维元胞自动机交通流模型. 物理学报, 2005, 54(6): 2595-2599. doi: 10.7498/aps.54.2595
    [11] 牟勇飚, 钟诚文. 基于安全驾驶的元胞自动机交通流模型. 物理学报, 2005, 54(12): 5597-5601. doi: 10.7498/aps.54.5597
    [12] 吴可非, 孔令江, 刘慕仁. 双车道元胞自动机NS和WWH交通流混合模型的研究. 物理学报, 2006, 55(12): 6275-6280. doi: 10.7498/aps.55.6275
    [13] 郭四玲, 韦艳芳, 薛 郁. 元胞自动机交通流模型的相变特性研究. 物理学报, 2006, 55(7): 3336-3342. doi: 10.7498/aps.55.3336
    [14] 张文铸, 袁 坚, 俞 哲, 徐赞新, 山秀明. 基于元胞自动机的无线传感网络整体行为研究. 物理学报, 2008, 57(11): 6896-6900. doi: 10.7498/aps.57.6896
    [15] 岳 昊, 邵春福, 陈晓明, 郝合瑞. 基于元胞自动机的对向行人交通流仿真研究. 物理学报, 2008, 57(11): 6901-6908. doi: 10.7498/aps.57.6901
    [16] 梅超群, 黄海军, 唐铁桥. 高速公路入匝控制的一个元胞自动机模型. 物理学报, 2008, 57(8): 4786-4793. doi: 10.7498/aps.57.4786
    [17] 梅超群, 黄海军, 唐铁桥. 城市快速路系统的元胞自动机模型与分析. 物理学报, 2009, 58(5): 3014-3021. doi: 10.7498/aps.58.3014
    [18] 李庆定, 董力耘, 戴世强. 公交车停靠诱发交通瓶颈的元胞自动机模拟. 物理学报, 2009, 58(11): 7584-7590. doi: 10.7498/aps.58.7584
    [19] 单博炜, 林鑫, 魏雷, 黄卫东. 纯物质枝晶凝固的元胞自动机模型. 物理学报, 2009, 58(2): 1132-1138. doi: 10.7498/aps.58.1132
    [20] 田欢欢, 薛郁, 康三军, 梁玉娟. 元胞自动机混合交通流模型的能耗研究. 物理学报, 2009, 58(7): 4506-4513. doi: 10.7498/aps.58.4506
  • 引用本文:
    Citation:
计量
  • 文章访问数:  390
  • PDF下载量:  51
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-04-10
  • 修回日期:  2020-06-16
  • 上网日期:  2020-06-17
  • 刊出日期:  2020-10-05

基于元胞自动机的气动光学光线追迹算法

  • 1. 中南大学物理与电子学院, 长沙 410083
  • 2. 武汉大学高等研究院, 武汉 430072
  • 通信作者: 谢文科, wenkexiedan@163.com
    基金项目: 装备预研领域基金(批准号: 6140415020311)、高能激光技术湖南省重点实验室开放基金(批准号: GNJGJS04)和湖南省光电惯性工程技术研究中心开放基金(批准号: HN-NUDT1908)资助的课题

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

English Abstract

    • 气动光学是研究稳态或非稳态气体流场与在其中传输的光束相互作用及传输规律的学科[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)和折射率nj, k来表示. 光子在元胞间移动遵循一定位置变换规则和方向变换规则. 位置变换规则用来获得光线矢量矢端点的位置信息和判断本次迭代后是否跨越网格, 如果跨越网格, 则需要根据方向变换规则计算光线偏转角, 反之则保持方向进入下一次位置变换. 可见, CA追迹算法与差分数值求解光线方程追迹法最大的不同在于: 将光线位置的改变与光线方向的改变拆分为两个相互独立的过程. 避免了光线位置改变和方向改变的同时计算, 因此编程更容易实现且具有较高的效率. CA光线追迹算法的具体步骤如下:

      1)给定初始光线矢量${{r}_0} = (0, 0)$. ${{V}_1}$表示初始元胞(j0, k0)处归一化光速, $\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)式可写成xy方向的分量的形式:

      $\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}}$分别表示Lxy方向的分量. 对(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}}$分别表示计算后得出的在xy方向的偏转角的改变量. (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}}$采用中心差分来计算, 其中hx, ihy, i表示xy方向相邻元胞的间距. 对于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}}$得到.

    • 本文基于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, 对应的KGD = 2.195 × 10–4 m3/kg, 跨帧间隔为15 μs, 样本总数为50帧. 流向为x正向, 光线传输方向为y正向. 实验测量的超声速边界层流场的马赫数为3.0, 对应的来流速度为622.5 m/s. NPLS图像的像素尺寸为2048 × 1436, 单像素分辨率为h = 0.013 mm, 入射光波长λ = 0.5 μm, 对应的KGD = 2.2 × 10–4 m3/kg, 跨帧间隔为15 μs, 样本总数为30帧. 流向为x正向, 光线传输方向为y正向. NPLS技术拍摄的某时刻超声速混合层和边界层图像如图1所示.

      图  1  NPLS获得的超声速剪切层图像 (a)混合层; (b)边界层

      Figure 1.  The flow visualization results of supersonic shear layers obtained by NPLS: (a) Supersonic mixing layer; (b) supersonic boundary layer.

    • 在本文中, 依据模型尺寸参数, 采用基于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  光学头罩绕流流场密度场

      Figure 2.  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算法在混合层得的光线路径

      Figure 3.  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)含激波的超声速光学头罩二维剖面流场

      Figure 4.  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. 超声速混合层、超声速边界层以及超声速光学头罩二维剖面流场的OPDrms分别为0.0086 λ1, 0.0014 λ2和0.0146 λ3, 其中λ1, λ2λ3分别对应超声速混合层、超声速边界层以及超声速光学头罩二维剖面流场的波长. 结果表明: 对于超声速混合层、超声速边界层以及超声速光学头罩二维剖面流场, 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/st2/st3/st4/st5/sta/s
      CA2.0322.0412.1572.0852.0772.078
      NSRE8.4638.4708.5228.4918.5118.491

      表 1  CA与NSRE算法计算混合层流场的程序执行时间

      Table 1.  The program running time of CA and NSRE in mixing layer.

      方法t1/st2/st3/st4/st5/sta/s
      CA1.9662.0021.9831.9641.9721.977
      NSRE7.2417.2507.1437.2507.2197.220

      表 2  CA与NSRE算法计算边界层流场的程序执行时间

      Table 2.  The program running time of CA and NSRE in boundary layer.

      方法t1/st2/st3/st4/st5/sta/s
      CA2.8312.8842.8272.8732.8402.851
      NSRE11.37511.52511.40111.38011.39711.416

      表 3  CA与NSRE算法计算高速绕流流场的程序执行时间

      Table 3.  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数据等方面提供的支持与帮助.

参考文献 (31)

目录

    /

    返回文章
    返回