一种模拟高速流场中星光传输的气动光学效应的方法与流程

本发明属于气动光学领域，涉及光线在飞行器高速绕流场内的传输模拟方法，具体是一种模拟高速流场中星光传输的气动光学效应的方法。

背景技术：

随着高速和高超声速飞行器的快速发展，高精度的光学导航和制导技术成为关注的热点，其中天文导航和红外制导是典型的代表。但在高速流场中光学设备在接收入射的光线时会受到气动光学效应的干扰，使所成的像产生模糊效果。在地面试验中模拟这种气动光学效应并研究相应的图像补偿和复原方法是提升高速飞行器导航和制导精度的有效途径。

目前气动光学效应的精确模拟需要飞行试验、风洞试验或高精度的计算流体力学(cfd)仿真，这些方法都需要高昂的试验成本。因此，开发一种简化的计算模型来模拟气动光学效应可以降低成本。

湍流是引起气动光学效应的主要因素。目前的流体力学理论和实验均表明湍流是由不同尺度的涡构成的。因此，模拟湍流的涡结构是模拟气动光学效应的关键。用人工放置的气体涡结构来模拟湍流中的涡结构是一种降低气动光学效应模拟复杂度的方法。

美国学者最早提出了一种用随机放置的气体球涡来模拟湍流密度脉动的方法，见文献[1]：一种气动光学效应的测试与诊断仿真技术,作者：jamestrolinger，davidweber，williamrose,公开日期：2002年1月；国内学者随后也提出了相似的模拟方法。这些方法的流程是：首先设定气体球涡的尺寸、内外密度差和人工流场中的放置位置等参数；然后通过光线追迹方法计算光束穿过人工球涡流场后的光程差；最后通过分析光程差来评估气动光学效应的强弱。

上述文献中的涡结构模拟技术均采用的是放置气体球涡的方法，球体的空间尺度是各项同性的，不能反映不同入射角度的光线在气动光学效应上的差别，而大量实验均表明气动光学效应对入射光角度具有明显的依赖性，见文献[2]：气动光学效应的物理原理与测量：回顾与展望，作者：ericj.jumper，stanislavgordeyev，发行日期：2017年。此外，现有的涡结构模拟技术中人工涡结构的参数设定均依赖于经验，没有规范的设计流程，可操作性不强。

技术实现要素：

本发明针对上述两个问题，提出了一种用气体椭球涡模拟气动光学效应的方法，不仅可以模拟气动光学效应对入射光角度的依赖性，而且建立了一种依赖于理论结果的气体椭球涡参数设计方法，可操作性更强；具体是一种模拟高速流场中星光传输的气动光学效应的方法。

具体步骤如下：

步骤一、设定待模拟的湍流边界层流场的区域范围；

定义待模拟的湍流边界层流场所在区域的光学窗坐标系oxyz，坐标原点o位于光学窗上表面的中心，x轴正向为流向，z轴垂直于光学窗表面指向飞行器壳体外侧，y轴与x轴和z轴构成右手坐标系；区域范围包括长度、宽度和边界层的厚度。

步骤二、设定气体椭球的控制参数；

控制参数包括椭球的数量knum、位置、尺度、倾角、扁率和内外密度。

knum根据经验取值；在其他参数不变的情况下，增加knum会使畸变波面上波峰/波谷的数量增加；

椭球的位置控制气体椭球在区域内的位置分布，分布规则为：涡的数量随离开壁面的距离增加；

椭球的尺度用气体椭球的长轴长度来表示；用顶点在近壁面高度的二次函数表示椭球尺度λ随高度zⁿ的变化趋势，表达式定义如下:

其中，klen是椭球尺度的控制参数，zⁿ表示归一化后的椭球球心离开壁面的距离，即归一化高度；是最小尺度气体椭球出现的归一化高度，即二次函数的顶点，λmin是最小的椭球尺度。

气体椭球的倾角θ定义为椭球体的长轴与光学窗坐标系oxyz的oxy平面的夹角；

气体椭球的扁率α定义为α＝(a-b)/a；a为气体椭球的长轴长度，b为气体椭球的短轴长度；

气体椭球内外密度是均匀的，仅在椭球面上存在密度的突变。

步骤三、在设定的区域范围内按位置控制参数定义的概率密度随机生成若干椭球体，再按照气体椭球的其余控制参数给每个椭球面所包含的空间范围赋予密度值，形成人工密度场；

具体步骤如下：

步骤301、在气体椭球数量区间内随机选取一个整数，作为设定的区域范围内要放置总椭球体数量；

步骤302、按照气体椭球位置坐标的概率分布生成的随机数，给每个椭球分配一个球心位置坐标；

气体椭球位置x和y坐标的概率分布是设定区域范围内的均匀分布，z坐标的概率分布满足下式:

p(zⁿ)表示zⁿ位置放置椭球体的概率；k是概率归一化系数；表示放置椭球体最多的高度；σ是分布的标准差，用于控制气体椭球在z方向分布的集中程度。

步骤303、根据椭球体的球心位置坐标、长轴大小和倾角计算每个椭球面所包含的空间范围；

具体为：设中心坐标为原点的正椭球面方程为

则满足不等式的坐标点均包含在正椭球体内。

对于中心坐标为[x0y0z0]且倾角为θ的斜椭球，采用坐标旋转的办法来选取其内部的点，设斜椭球面上的点的坐标为[x′y′z′]，根据坐标旋转原理，[x′y′z′]与正椭球面上的点[xyz]满足如下关系：

将式(4)右端项带入方程式(3)中，分别替换[xyz]即可得到斜椭球面方程，同理选取斜椭球体内的所有坐标点。

步骤304、根据气体椭球内的密度表达式给每个椭球面所包含的空间坐标集合赋予密度值。

气体椭球内的密度ρin表示为

ρin＝(-1)ⁿkρf(m,z)λρ∞+ρ∞(5)

其中，(-1)ⁿ是取1或-1的随机数，表示气体椭球内的密度可能高于或低于来流密度。kρ是气体椭球的密度控制参数，f(m,z)是马赫数影响因子，ρ∞为气体椭球外的密度。

步骤四、使用光线追迹方法计算某一角度的入射光穿过人工密度场后的光程差，得到一张畸变波前快照。

光程差定义为光线的传播距离s与传播路径上的折射率n的乘积。对于整个光路，光程差用每段路径的光程累加得到，即

δsi表示第i段路径的几何长度，传播距离共有k段。ni表示第i段路径的折射率，且ni与第i段路径上的密度ρi的关系为ni＝1+kgdρi，其中kgd为光线的gladstone-dale系数，只与光线的波长有关。

步骤五、重复步骤三和步骤四共n次，得到n张畸变波前快照，将得到的n张畸变波前快照求平均得到时均波前畸变，即得到模拟的光线穿过湍流边界层后的气动光学效应。

n是快照叠加的帧数。

本发明的优点在于：

1)、一种模拟高速流场中星光传输的气动光学效应的方法，用气体椭球来模拟湍流边界层中涡结构，能够模拟气动光学效应对入射光角度的依赖性，并且相比物理实验和计算流体力学仿真实验手段，成本更低，所需时间更少。

2)、一种模拟高速流场中星光传输的气动光学效应的方法，气体椭球的控制参数是基于目前的气动光学研究的理论成果，具有公式化的显示表达形式，对经验依赖性较小，方法的可操作性较强。

附图说明

图1为本发明人工流场中放置气体椭球的示意图；

图2为本发明一种模拟高速流场中星光传输的气动光学效应的方法的流程图；

图3为本发明气体椭球尺度klen与雷诺数rex的幂函数关系；

图4为本发明气体椭球的示意图

图5(a)为本发明叠加快照的数量n与平均畸变波前均方根opdrms间的关系。

图5(b)为本发明每增加一张快照所造成的平均光程差的变化量总和；

图6为本发明不同角度入射光的归一化波前畸变对比示意图；

图7为本发明模拟得到的垂直入射平行光的畸变波前(光程差)曲面。

图8为本发明模拟得到的垂直入射平行光成像后的像面灰度分布。

具体实施方案

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种模拟高速流场中星光传输的气动光学效应的方法，用人工放置的气体椭球生成人工密度场来模拟光线穿过湍流边界后受扰效果。用气体椭球涡模拟的人工流场如图1，采用了气体椭球涡模拟气动光学效应的方法，通过按一定规则在人工流场中随机放置气体椭球，使不同角度入射光产生了与入射角有关的模糊效果，达到了模拟湍流边界层气动光学效应的作用，具有计算量小，参数易于确定的优势。

如图2所示，具体步骤如下：

步骤一、设定待模拟的湍流边界层流场的区域范围；

定义待模拟的湍流边界层流场所在区域的光学窗坐标系oxyz，坐标原点o位于光学窗上表面的中心，x轴正向为流向，z轴(表示壁面高度)垂直于光学窗表面指向飞行器壳体外侧，y轴与x轴和z轴构成右手坐标系；模拟湍流边界层流场的区域范围设定为：长(x轴)200mm，宽(y轴)100mm，厚(z轴)30mm。

步骤二、设定气体椭球的控制参数；

控制参数包括椭球的数量knum、位置、尺度、倾角、扁率和内外密度。

(1)气体椭球的数量knum

knum会影响本方法输出的畸变波面的形状，在其他参数不变的情况下，增加knum会使畸变波面上波峰/波谷的数量增加。knum的取值依赖于经验，一般经验值是几百到几千。在设定的流场区域尺寸下，取knum＝1000能够使输出较为合理的畸变波面形状。实际操作时，在模拟每帧波前快照时放置的气体椭球数在[0.9knum，1.1knum]范围内取随机整数，这样能够体现湍流的随机性。

(2)气体椭球的位置

目前的气动光学实验表明：边界层的外区是光学活跃区域，对波前畸变起主要贡献；湍流边界层中大尺度结构对穿过的光束的畸变起主导作用。因此，本发明设定气体椭球的位置控制气体椭球在区域内的位置分布，分布规则为：涡的数量随离开壁面的距离有增加的趋势；用如下表达式定义:

其中，p(zⁿ)表示zⁿ位置放置椭球体的概率。zⁿ＝z/δ表示归一化后的气体椭球的球心离开壁面的距离，即归一化高度；δ表示模拟区域的边界层厚度，表示放置椭球体最多的高度。目前的气动光学试验结果表明壁面高度0.8δ处的流场结构对光学畸变的贡献最大，因此取σ是此分布的标准差，用于控制气体椭球在z方向分布的集中程度，按经验值取σ＝0.5δ。k是概率归一化系数，满足

zⁿi表示第i个气体椭球的球心离开壁面的距离。

(3)气体椭球的尺度

本发明用气体椭球的长轴长度来表示椭球体的特征尺度。目前的气动光学试验结果表明湍流边界层外区的密度关联长度λ(zⁿ)随zⁿ递增，在近壁面区域密度关联长度显著增加。考虑到规则设计的简洁性，本发明用顶点在近壁面高度的二次函数表示椭球体特征尺度λ随高度zⁿ的变化趋势，用如下表达式定义:

其中，是最小尺度气体椭球出现的归一化高度，即二次函数的顶点，取为0.05δ。λmin是最小的椭球体尺度，根据现有技术中对边界层气体球涡尺度下界的设定，可知湍流边界层中最小尺度的涡结构在0.05δ～0.1δ范围内，因此取λmin＝0.05δ。klen是椭球体尺度的控制参数，气动光学实验表明雷诺数rex会影响湍流边界层中的密度相关尺度，且雷诺数越大相关尺度和opdrms越小，因此，本发明通过物理意义明确的雷诺数rex来设定klen。在高速飞行器的典型雷诺数范围rex∈[1×10⁶，9×10⁸]内通过经验观察畸变波面破碎程度的相对大小，给定klen的取值范围建立如图3所示的klen与rex的幂函数关系，表达式如下

klen＝g1(rex)＝242.7(rex)^-0.1883(4)

(4)气体椭球的倾角和扁率

湍流边界层中的结构倾角在试验和流体力学仿真中均得到了证实。这种结构倾角的大小很难从流场中直接辨识，通过测量不同角度入射光线的光程差均方根opdrms的大小能够间接反映结构倾角。椭球体尺度具有空间各向异性特点，适合表现湍流边界层中的结构倾角。

本发明将长轴长度为a、短轴长度为b的椭圆绕其长轴旋转后得到旋转椭球体，再将此椭球体旋转θ角，如图4所示。用此旋转后的椭球体来表示湍流边界层中的密度相关结构。

气体椭球的倾角θ定义为椭球体的长轴与光学窗坐标系oxyz的oxy平面的夹角，取θ＝24.17°，气体椭球的扁率α定义为α＝(a-b)/a；取扁率α＝0.447。

以γ角射出/射入的光线(γ角定义为光线与上游来流方向的夹角)，其opdrms与垂直入射光的opdrms具有相似的比例关系：

其中c(γ)是与光束倾角有关的比例系数，根据实验结果，其上界大概在1.4～2范围内，本发明取1.7，其下界c(γ)约为0.94。这也就是指向下游的光束比指向上游的光束的畸变更严重，这可以解释为指向下游的光束在湍流结构中的光程更长。

本发明将整个流场等效为一个如图4所示的气体椭球，则沿长轴a和沿短轴b方向的光束分别具有最大和最小畸变。由于气体椭球内各处具有相同的密度，所以光程可由几何长度代替，可得

其中，lv的定义见图4。

由上式可解得气体椭球的倾角θ＝24.17°，扁率α＝(a-b)/a＝0.447，对于人工流场区域内的所有气体椭球均采用相同的倾角和扁率设定。

(5)气体椭球内的密度

不同于真实流场中密度是连续变化的，为了便于计算，本发明设定气体椭球内外密度都是均匀的，仅在椭球面上存在密度的突变。

通过将模拟得到的畸变波前的opdrms与现有技术中具有理论公认性的湍流边界层气动光学效应预测模型(以下简称blpm)锚定，来设计气体椭球的密度参数。由于blpm对马赫数m具有依赖性，所以在气体椭球的密度表达式中加入了马赫数影响因子。由于目前的气动光学试验表明大尺度结构对波前畸变起主要贡献，所以将气体椭球的尺度λ也加入到椭球体密度表达式中，则气体椭球内的密度ρin表示为

ρin＝(-1)ⁿkρf(m,z)λρ∞+ρ∞(6)

其中，ρ∞为气体椭球外的密度，设为来流密度。(-1)ⁿ是取1或-1的随机数，表示气体椭球内的密度可能高于或低于来流密度。f(m,z)是马赫数影响因子，表达式为

其中，κ是比热比。f(z/δ)＝uz/δ/u∞是平均速度廓线；uz/δ是高度z/δ处的平均速度，u∞是来流速度。根据流体力学理论，边界层内区和外区的平均速度分别具有线性和对数形式

其中，系数c1,c21和c22可由平均速度的边界条件和其在内外区的连续性条件解得。

kρ是气体椭球的密度控制参数，本发明通过标定kρ将所模拟的畸变波面的opdrms与blpm锚定，标定kρ的步骤如下：

1)在给定的一个雷诺数rex0下，在不同马赫数下拟合得到opdrms关于kρ的函数opdrms(kρ)；

2)由锚定条件(畸变波面的opdrms＝blpm理论计算值)和拟合函数opdrms(kρ)确定kρ在不同马赫数下取值的散点集合{kρ(m,rex0)i}；

3)选取不同的rex0，重复步骤1)和2)，得到不同马赫数和不同雷诺数下kρ取值的散点集合{kρ(m,rex)i}，使用曲面拟合技术得到kρ关于马赫数m和雷诺数rex的函数表达式kρ(m,rex)。

步骤三、在设定的区域范围内按位置控制参数定义的概率密度随机生成若干椭球体，再按照气体椭球的其余控制参数给每个椭球面所包含的空间范围赋予密度值，形成人工密度场；

首先在气体椭球数量区间内随机选取一个整数，作为要放置总椭球体数量；然后按照气体椭球的位置分布概率生成随机数，给每个椭球分配一个球心位置坐标；接着根据椭球体的球心位置、长轴大小、倾角和扁率计算每个椭球面所包含的空间范围；最后根据气体椭球内的密度表达式给每个椭球面所包含的空间坐标集合赋予密度值。

具体步骤如下：

步骤301、在气体椭球数量区间内随机选取一个整数，作为设定的区域范围内要放置总椭球体数量；

步骤302、按照气体椭球位置坐标的概率分布生成的随机数，给每个椭球分配一个球心位置坐标；

气体椭球位置x和y坐标的概率分布是设定区域范围内的均匀分布，z坐标的概率分布满足式(1)。

步骤303、根据椭球体的球心位置坐标、长轴大小和倾角计算每个椭球面所包含的空间范围；

具体为：设中心坐标为原点的正椭球面方程为

则满足不等式的坐标点均包含在正椭球体内。

对于中心坐标为[x0y0z0]且倾角为θ的斜椭球，采用坐标旋转的办法来选取其内部的点，设斜椭球面上的点的坐标为[x′y′z′]，根据坐标旋转原理，[x′y′z′]与正椭球面上的点[xyz]满足如下关系：

将式(10)右端项带入方程式(9)中，分别替换[xyz]即可得到斜椭球面方程，同理选取斜椭球体内的所有坐标点。

步骤304、根据气体椭球内的密度表达式(6)给每个椭球面所包含的空间坐标集合赋予密度值。

步骤四、使用光线追迹方法计算某一角度的入射光穿过人工密度场后的光程差，得到一张畸变波前快照。

光程差定义为光线的传播距离s与传播路径上的折射率n的乘积。对于整个光路，光程差用每段路径的光程累加得到，即

δsi表示第i段路径的几何长度，传播距离共有k段。ni表示第i段路径的折射率，且ni与第i段路径上的密度ρi的关系为ni＝1+kgdρi，其中kgd为光线的gladstone-dale系数，只与光线的波长有关。

步骤五、重复步骤三和步骤四共n次，得到n张畸变波前快照，将得到的n张畸变波前快照求平均得到时均波前畸变，即得到模拟的光线穿过湍流边界层后的气动光学效应。

n是快照叠加的帧数，属于本方法的控制参数。

用涡结构模拟技术生成的一个边界层流场模拟得到的畸变波前代表了一个瞬时的结果，将一定数量的瞬态畸变波前快照平均后可以模拟曝光时间内的平均波前。为了确定合适的叠加快照数量n，考虑两方面的因素：一方面是平均畸变波前随n变化的稳定度，当平均畸变波前随n的增加变化不大时认为输出稳定。另一方面是计算量，叠加的快照数越少耗时越短。

平均畸变波前随n的变化如图5所示，图5(a)中平均畸变波前均方根随n增大而减小，但减小的速度在n超过20后迅速降低并逐渐趋于0，图5(b)是每叠加一张快照所造成的平均波前光程差的变化量总和，可以看出在n超过25后，再增加快照数量对波面光程差的影响很小，说明仿真器的输出达到了稳定。综合考虑输出稳定性和计算量，本方法将叠加的快照数量n取为25。

仿真条件如表1所示：

表1

本发明对不同角度入射光束气动光学效应(用畸变波前的光程差均方根opdrms来衡量)的仿真结果如图6所示。图6中横轴表示出射/入射光线与上游的夹角，纵轴表示对垂直入射光归一化后的opdrms，本发明的模拟方法与现有技术中通过计算流体力学(cfd)仿真得到的实验结果进行了对比，simulator表示本发明的模拟方法；white&visbal(2012)和wang&wang(2012)均为现有技术中通过计算流体力学(cfd)仿真得到的实验结果；可以看出本发明的模拟方法与文献数据十分接近，说明本发明的气动光学模拟方法能够有效模拟湍流边界层中气动光学效应的各向异性特征。

本发明对垂直入射平行光的气动光学效应的畸变波前和成像效果的仿真结果如图7和图8所示。图7与现有技术中通过计算流体力学方法得到的畸变波前形状具有很高的相似度，但在表1中所列的硬件仿真条件下，本发明的模拟方法耗时约2400～2600秒，而计算流体力学方法通常需要在大型服务器上运行几十个小时才能获取相似的仿真结果。图8中使用本发明的方法模拟天文导航中星光受到气动光学效应扰动后成像模糊的效果，可以看出像面星点能量扩散的效果被真实地反映了出来。