吕春静 韩一平

1) (西安电子科技大学通信工程学院,西安 710071)

2) (西安电子科技大学物理与光电工程学院,西安 710071)

为了研究高斯光束在湍流等离子体鞘套中的传输特性,根据广义惠更斯-菲涅耳原理,采用基于快速傅里叶变换的功率谱反演法,用多随机相位屏来模拟湍流带来的影响.根据超声速飞行器绕流等离子体流场厚度在厘米级别的特点,光束在两个相位屏之间的传输过程中采用菲涅耳衍射积分的两次快速傅里叶变换算法(double fast Fourier transform algorithm),利用多随机相位屏模拟等离子体鞘套湍流对光束传输产生的影响,解决了多随机相位屏模拟湍流研究中的超短距离传输问题.当飞行高度为45 km,飞行速度为18马赫时,通过对超声速飞行器绕流等离子体流场的统计分析,发现在此飞行条件下折射率起伏方差的强度范围10-11—10-14.对高斯光束在湍流等离子体流场中的传输特性进行了数值仿真.结果表明：在等离子体鞘套湍流中折射率起伏强度、波长、传输距离等都是影响高斯光束质量的重要因素.折射率方差越大,传输距离越长,光斑弥散越严重,光强起伏越大,光强减弱也越明显.光束的波长越长,高斯光束抑制湍流的能力越强,光斑弥散程度越小,光强起伏也越小.

1 引 言

当飞行器以高超声速在大气中飞行时,其动能会迅速转化为大量的热能导致飞行器周围的空气产生离解和电离,最终形成等离子体包覆的流场[1-4].等离子体的存在对于电磁波信号有很大的衰减作用,甚至导致信号中断形成“黑障效应”[5,6].目前已经提出几种有希望解决这一难题的几种方法[7-11],例如通过提高电磁波频率来减小信号衰减,或者在鞘套层中加入亲电材料,这两种方法的目的都是为了提高电磁波频率与等离子体频率的比值.同时有研究表明等离子体鞘套层是不稳定的[12,13],所以湍流也是影响“黑障效应”的重要因素.近些年,国内外也早对湍流等离子体流场进行了较多的研究,例如Shi等[14]为了研究等离子体鞘套在无线信道上的动态效应,分析飞行条件变化(例如攻角、小尺度湍流)对介质传输性能的影响,提出了一种新的自适应马尔可夫信道建模方法,利用该方法建立的信道的一阶统计特性与无线电波的可变衰减具有良好的一致性,而且该方法具有良好的自适应能力.He等[15]从通信角度采用透射矩阵法研究了再入等离子体鞘套对电磁波传播的动态影响,所建立的模型中包含空间湍流和等离子体电子密度的变化,并对由等离子体的非线性特性引起的幅度和相位变化进行了数值仿真,发现幅度变化服从对数正态分布,相位变化服从瑞利分布.Yao等[16]针对等离子体流场的不稳定性,首先在时-空-频域揭示了电子密度的波动的规律,并基于这些规律提出了时变电子密度分层模型,利用Monte Carlo准静态电磁波数值方法计算了动态等离子鞘套层对Ka波段电磁波的影响.Lin和Sproul[17]重点研究了电子密度波动对电磁波传播的影响,利用Navier-Stokes模型、三维抛物化Navier-Stokes模型、耦合边界层和非黏性流动代码估算了截头锥体的电子密度的均值和方差.Li等[18]基于在高超声速湍流实验中测量的分维数,建立了分形模型和等离子鞘套湍流功率谱,为研究等离子体鞘套湍流中的光束传输特性提供了重要参考.

在采用随机相位屏研究湍流问题的文献中,两个随机相位屏之间的传输过程大多采用单次快速傅里叶变换(single fast Fourier transform,S-FFT)算法[19-21],但是这种菲涅耳衍射积分的S-FFT算法适用于长距离传输,例如大气湍流[22-24]和海洋湍流[25-27].飞行器周围等离子体包覆的流场的厚度很小,在厘米或毫米级别[28,29].对于厘米和毫米级别的短距离传输问题,S-FFT算法会出现欠采样的问题,最终将会导致信号失真.针对等离子体鞘套厚度较小这一特点,采用菲涅耳衍射积分的两次快速傅里叶变换(double fast Fourier transform,D-FFT)算法研究光束在两个相位屏之间的传输过程.本文第2部分首先在高度为45 km、飞行速度为18马赫和0攻角的飞行条件,数值模拟了钝头体飞行目标的高超音速等离子体鞘套湍流流场.然后利用Gladstone-Dale公式根据等离子体密度计算折射率,并对折射率的方差进行了估算;第3部分对比研究了菲涅耳衍射积分的D-FFT算法和S-FFT算法,说明D-FFT算法在厘米级别的短距离传输问题中的显著优势,并利用多随机相位屏模拟了高斯光束在等离子体鞘套湍流中的传输特性.

2 高超声速钝头体等离子体鞘套湍流流场计算

2.1 钝头体飞行目标绕流流场数值模拟

利用三维Navier-Stokes (N-S)方程和湍流两方程k-ε模型数值模拟湍流流场.在笛卡尔坐标系中,考虑1989年提出的双温热力学模型[30]和7组分空气化学反应模型(七组分化学反应模型如表1所列),并对三维N-S方程进行适当封闭,相应模型如下[31-33]：

其中,Q和S分别是守恒变量和源项矢量;以x方向为例,E和Ev分别是对流通量矢量和黏性通量矢量;κ和κt分别是热传导系数和湍流热传导系数;u,v和w分别是3个方向的瞬时速度分量;k是湍动能;是振动源项,是每个组分所占的质量分数;τij(ij=x,y,z) 为各个方向的剪应力张量,需要利用湍流模型进行封闭[34];ev,i,hi和ρi分别是每一种组分的振动能量、显焓和密度,下标i代表每一种组分 (i=1,2,3,···,7);ρ,p和Ua分别为混合物的密度、压强和平均速度;Eint,Et和Tint分别是混合物的内部能量、湍流内部能量和内部温度.

基于有限体积法,利用AUSMPW+(the advection upstream splitting method by pressurebased weight function)方法[35]求解N-S方程,其中时间推进采用LU-SGS (lower-upper symmetric Gauss-Seidel relaxation)方法[36].值得注意的是,文中采用的湍流k-ε模型是利用不同的方程来控制湍流能量传输并适合高雷诺数环境的两方程模型,其中控制湍动能k和湍动能耗散率ε的两方程分别为

表1 化学反应模型中的方程Table 1.Reactions considered in chemistry model.

其中,µ,µt,µi分别是等离子体的动力学黏性系数、湍流黏性系数和质量平均速度;常数Cε1=1.44,Cε2=1.92,σk=1.0,σε=1.3;P被定义为ui,uj和um为张量表示法中时均速度,δij是克罗内克尔符(当i=j时,δij=1,否则δij=0 ).

2.2 对流场数据做统计分析,估算等离子体折射率起伏方差

基于Gladstone-Dale公式[37],根据2.1节中计算的等离子体密度求解其折射率,Gladstone-Dale公式为

图1(a)和图1(b)分别为飞行目标的区域划分图和五个区域的折射率统计方差图.图2(a)和图2(b)分别为飞行目标头部和侧部选取的一条传输路径,头部路径三维坐标范围分别为：x=[0.0047 m,0.0170 m],y= [—0.0581 m,—0.0475 m],z= [0.0627 m,0.0512 m],侧部路径三维坐标范围分别为：x= [0.0674 m,0.0788 m],y= [—0.1555 m,—0.1225 m],z= [0.0628 m,0.0527 m].对两条传输路径上的折射率进行方差统计,统计方法为：以传输路径中的每个三维流场坐标点为中心,以超声速湍流中典型的相关长度[38]—L/10为半径(L为等离子体鞘套的厚度)建立球模型,对球内所有数据进行统计,计算每个网格点处的折射率方差,流场边界值采用外推法求解.通过对比图2(a)和图2(b)发现,在接近飞行目标壁面时,头部折射率方差较大,比两侧的折射率方差要高出一个数量级,这是因为等离子鞘套中的湍流状态还会受到激波、高温、高压等因素的影响,在靠近头部的区域内温度、压强较高且起伏较大.由图1和图2的计算结果显示,在高度为45 km、飞行速度为18马赫和0攻角的飞行条件下,飞行目标的绕流流场中折射率起伏方差的强度的为10-11—10-14.

图1 (a)五个区域;(b)五个区域的折射率方差Fig.1.(a) Five regions;(b) refractive index variance in five regions.

图2 折射率方差 (a)头部区域1;(b)侧部区域2Fig.2.Refractive index variance in different regions：(a) Head area 1;(b) side area 2.

3 湍流等离子体鞘套中的随机相位屏模拟

3.1 相邻两屏之间传输算法讨论及仿真实验验证

S-FFT算法和D-FFT算法都是从菲涅耳衍射公式出发,但对于衍射距离d的要求不同.SFFT算法主要适用于衍射距离较大的情况,而DFFT算法可适用于衍射距离较小的问题,下面对两种算法进行理论分析和数值仿真,说明选用DFFT算法的原因.

3.1.1 菲涅耳衍射积分的S-FFT算法

在傍轴近似情况下,菲涅耳衍射积分可以表示为

其中,λ为光波波长,kλ=2π/λ为波数,d是观察面到衍射屏的距离,(x0,y0)和(x,y) 分别是光场在衍射面上和观察面上某一点的复振幅分布.对(10)式进行适当变换,并根据傅里叶变换法则,上式可改写为

由于这种算法只需要进行一次傅里叶变换,所以称为S-FFT算法.由(11)式可知其空间频率为,即观察屏上不同的点(x,y) 其实对应不同的空间频率(方向)的光.又因为光传播方向与空间频率的关系为fµ=1/X=cosα/λ,所以fµ=sin(π/2-α)/λ≈1/λ·(x/d) .

为了确定观察面的尺寸,以x方向为例进行讨论说明,y方向的情况与之类似.设衍射面在x方向上的宽度为Lx0,共有Nx个抽样点,则按照抽样定理,在用FFT变换计算时得到的空间频率最高为fµmax=Nx/2Lx0,空间频率为fµmax的这些衍射光对应的传播方向为fµmax=1/X=cosαmax/λ.这些衍射光出现在衍射观察屏上坐标的位置为fµmax=sin(π/2-α)/λ≈1/λ·(Lxmax/2/d) ,即衍射屏上的最大坐标为Lxmax=Nxλd/Lx0.由此可见在使用S-FFT算法时,衍射观察屏的大小不仅是波长的函数,而且还是再现距离和抽样数的函数,当d很小,如果保持取样数不变,则观察屏上的再现结果只是邻近光轴的很小区域,所以该算法主要适用于衍射距离较大的情况.

协议书附图是描述边界线地理位置的法律图件，是协议书的附图，由乡级行政区域界线信息与修测后的边界地形图叠加一起制作形成，与界线协议书具有同等的法律作用。

3.1.2 菲涅耳衍射积分D-FFT算法

同S-FFT算法相类似,在直角坐标系中,傍轴近似下,当衍射距离满足菲涅耳衍射条件时,根据衍射积分表达式(10)可以得到其脉冲响应：

显然,菲涅耳衍射下的脉冲响应是不变的,所以可以用卷积运算来完成,即

由于在空域完成卷积运算需要花费很长的机器时间,所以选择到频域完成.假设菲涅耳衍射的传递函数为(fµ,fv) ,则菲涅耳积分可以分为两次傅里叶变换完成,所以称其为D-FFT算法.可以通过傅里叶变换的性质直接得到菲涅耳衍射传递函数(fµ,fv) 的解析表达式：

综上所述,菲涅耳衍射公式可以这样计算：

D-FFT算法中根据已经确定的传递函数(fµ,fv)的解析表达式,只需要将频域中x方向和y方向的频率fµ和fv代入解析式,就可以得到传递函数的数值.假设实际问题是在宽度为Lx0×Ly0的空间采样,并且有Nx×Ny个抽样点,最高空间频率分别为.所以,fµ的取值范围从到取Nx个点,同样,fv的取值范围从到取Ny个点.

3.2 S-FFT算法和D-FFT算法的数值模拟结果的比较

在衍射计算过程中,应该特别注意欠采样的问题,在保证抽样数和光束波长不变的情况下,改变衍射距离d,比较S-FFT算法和D-FFT算法抑制欠采样发生的能力,说明D-FFT算法在处理短距离衍射积分问题的优势.

图3 在不同衍射距离下S-FFT和D-FFT两种算法的光强分布 (a) S-FFT,d = 100 mm;(b) D-FFT,d = 100 mm;(c) S-FFT,d = 10 mm;(d) D-FFT,d = 10 mmFig.3.Intensity distribution of S-FFT and D-FFT algorithm at different diffraction distances：(a) S-FFT,d = 100 mm;(b) D-FFT,d = 100 mm;(c) S-FFT,d = 10 mm;(d) D-FFT,d = 10 mm.

图3(a)和图3(b)分别为衍射距离为100 mm时利用S-FFT算法和D-FFT算法得到的衍射光强分布,图3(c)和图3(d)分别为衍射距离为10 mm时利用S-FFT算法和D-FFT算法得到的衍射光强分布.通过对比图3(a)和图3(b)发现：当衍射距离为100 mm时,利用S-FFT和D-FFT两种算法得到的衍射光强分布均未出现欠采样现象.通过对比图3(c)和图3(d)可知：随着衍射距离的减小,采用S-FFT算法得到的衍射光强出现了明显的欠采样现象,但是D-FFT算法未出现,可见D-FFT算法的抗欠采样能力要明显高于S-FFT算法.结果表明,D-FFT算法不惧怕计算衍射距离近的问题,即使衍射距离d= 1 mm时,用D-FFT算法依然可以得到很好的结果.超高音速飞行器周围形成的等离子体鞘套的厚度属于厘米数量级.由此可见,用FFT算法数值模拟光束在相邻两屏之间的传输过程时,选用D-FFT算法更合适.

利用谱反演法生成等离子体鞘套湍流随机相位屏,首先需要在频域生成一个均值为1,方差为0的N×N复随机矩阵H(kx,ky) ,利用等离子体鞘套湍流功率谱[18]对其进行滤波,然后再进行傅里叶逆变换便可以得到等离子体鞘套湍流的随机相位屏Φ(kx,ky) .等离子体鞘套湍流功率谱为

其中为折射率方差,L0为湍流外尺度,kx和ky分别为不同方向的频率.由折射率频谱可得到相位频谱为Φ(kx,ky)=2πk2Φn(kx,ky) ,则相位频谱的方差为σ2(kx,ky)=(2πd/NΔx)2Φ(kx,ky) ,其中d是每两个随机相位屏的间距,N为采样数,Δx是栅格间距.通过傅里叶变化便可得到空域中的随机相位屏[39]：Φ(x,y)=F[H(kx,ky)σ(kx,ky)] .

利用菲涅耳衍射积分的D-FFT算法求解高斯光束在相邻两屏之间的传输过程,利用功率谱反演法生成的多随机相位屏模拟等离子鞘套湍流带来的影响,多随机相位屏示意图如图4所示.

图4 多随机相位屏示意图Fig.4.Schematic diagram of multi-random phase screen.

图5 不同传输条件时的高斯光束光强数值仿真 (a)λ = 1.55 × 10-6 m,ω0 = 40 mm,d = 13 mm;(b)ω0 = 40 mm,d =13 mm, = 10-12;(c)λ = 1.55 × 10-6 m,ω0 = 40 mm, = 10-12Fig.5.Numerical simulation of intensity of Gaussian beam in different propagation conditions：(a)λ = 1.55 × 10-6 m,ω0 = 40 mm,d = 13 mm;(b)ω0 = 40 mm,d = 13 mm, = 10-12;(c)λ = 1.55 × 10-6 m,ω0 = 40 mm, = 10-12.

3.3 等离子体鞘套湍流中影响高斯光束质量的因子分析

图5的传输参数为：束腰半径ω0= 40 mm,相位屏间隔d= 0.013 m,相位屏采样数为 1024×1024.分析折射率方差、波长和传输距离对高斯光束质量的影响.

图5(a)是高斯光束在折射率起伏方差分别为=10-13,10-12,10-11时通过10个随机相位屏后系统归一化的三维(3D)和二维(2D)光强数值仿真结果.由图5(a)可知,在湍流折射率方差不断增大的过程中,光斑分裂成形状不规则的若干小光斑,光斑弥散程度加剧,而且光强起伏逐渐增大,强度逐渐减小.

图5(b)是高斯光束波长分别为λ=0.6328×10-6,1.05×10-6,1.55×10-6m时通过10个随机相位屏后系统归一化的3D和2D光强数值仿真结果.通过对比分析可知,随着波长的不断增加,高斯光束抵抗湍流的能力逐渐增强,具体表现在光斑弥散程度减小、光强起伏减弱上,所以在等离子体鞘套通信系统中选用波长较长的激光光束作为信息载体比较好.

图5(c)是高斯光束沿z轴传输时距离分别为6.5,13,19.5 cm,即分别通过5,10,15个随机相位屏后系统归一化的3D和2D光强数值仿真结果.结果反映出：传输距离越大,光斑弥散越严重,光强起伏也越明显.

4 结 论

首先利用菲涅耳衍射积分的D-FFT算法解决了湍流问题中的短距离传输问题,然后利用多随机相位屏数值模拟了高斯光束在等离子体鞘套湍流中的传输特性.研究表明折射率方差越大,高斯光束的光斑弥散程度越大,光强越弱,光强起伏越明显;随着波长的不断增加,高斯光束抵抗湍流的能力逐渐增强,所以在等离子体鞘套通信系统中选用波长较长的激光光束作为信息载体比较好;传输距离越大,光斑弥散越严重,光强起伏也越明显.综上所述,高斯光束在等离子体鞘套湍流中的传输质量与外界因素折射率起伏方差、传输距离及自身因素波长、束腰半径等密切相关.本文对研究激光在湍流等离子体鞘套中的传输问题具有一定的意义.