任意折射率介质中的光束追迹方法.docxVIP

任意折射率介质中的光束追迹方法 1 网格节点不同时点折射波场模拟 由于存在复杂的结构，如激波、边缘层、剪切层和水流，高速流动场的流场介质密度不均，没有规律性，因此折射不可避免。计算流体力学所计算网格点或者网格单元中心点的密度(折射率)代表了该点处介质的密度(折射率),所有网格点信息的集合则描述了整个流场的折射率分布,这是对无规则非均匀折射率场的一种描述方法。 目前对折射率分布规律已知的介质进行光线追迹的研究已经较为普遍,并且得到了广泛的应用,但是针对任意无规则折射率分布介质的研究则相对较少。这些研究的基本思想是将每个网格单元视为等折射率单元,以两个格子之间的界面作为折射面,根据Snell折射定律进行光线追迹。然而在实际情况中,计算高速流场所采用的网格在大部分情况下是非均匀的,网格尺寸相差较大,而且存在激波等折射率变化剧烈的区域,不能简单认为网格单元内折射率处处相等,而且以一个网格单元作为一次追迹步长也不合理。 本文基于三维空间离散网格,网格节点介质折射率已知,通过实时计算光传播路径上所需各点的折射率和折射率梯度,采用Runge-Kutta法实现光线追迹。在追迹过程中,每一步追迹的起点和终点不一定在网格点或者网格单元分界面上,并且追迹步长根据当地折射率梯度和网格尺寸自适应调节,以提高追迹的精度和效率。 2 限制节点法选取 下式描述了折射率任意分布介质中的光线传播,是研究梯度折射率介质光线传播的基本方程之一: 式中s为光线传播路径上的弧长,r为光线矢径,n为折射率,▽n为折射率梯度。 该方程在大多数情况下是难以求解的,仅对等折射率面为平面、球面和圆柱面等比较特殊的形式有解析解。对于任意非均匀折射率介质中的光线传播路径,基本上没有解析解,一般采用数值方法进行求解。由于Runge-Kutta法在达到相同精度的情况下,花费的运算量较少,因此本文选用RungeKutta法作为基本光线追迹手段。 对(1)式,令光线矢量Tray=ndr/ds=dr/dt,式中dt=ds/n是一个引入的参量,则光线传播方程可写为 式中有3个分量,可用三元一维数列解得,记 R为光线传播路径上的全坐标,T为光线矢量分量,D为(2)式右端项。经过推导可得Runge-Kutta法计算公式如下: 这里矩阵 式中Δt=Δs/n,Δs为空间追迹步长。当已知光线初始条件R=R0(x0,y0,z0),T=T0时,可以逐次进行追迹,在每一步追迹完成后,需要计算该步末端点的折射率和折射率梯度,作为下一次追迹的起始条件,依此往后,直到完成整个光线追迹过程。 3离散空间折射率及其梯度求解 采用Runge-Kutta法对折射率空间任意离散分布的流场介质进行光线追迹时,其追迹精度除了Runge-Kutta法本身精度的限制外,还受到3个因素的影响:1)追迹路径上的折射率的计算精度;2)追迹路径上的折射率梯度的计算精度;3)追迹步长的大小。下面将分别对这3个因素展开讨论。 3.1 距离加权平均插值方法 在离散空间内,已知网格节点上的折射率,需要求解其他空间任意点的折射率,通常采用插值方法。常用的插值方法有距离加权平均插值方法、补丁法、双三次样条插值和三线性插值法,此处选用精度和稳定性均较好的距离加权平均插值方法。 在三维离散空间内,如图1所示网格单元,已知单元的八个顶点的折射率为ni(i=1,2,…,8),要求解单元内部任一点P(x,y,z)的折射率nP。设点P到八个点(xi,yi,zi)的距离为di,距离加权平均插值方法可表示为 3.2 充分利用有限差分方法求解流场 已知空间网格节点的折射率,求解网格内任意点的折射率梯度,需要分两步:第一步求解网格点本身的梯度值,第二步利用网格单元顶点的梯度插值求解网格内任意点的梯度值。其中第二步的插值方法与求解折射率的方法相同,这里重点论述第一步的求解方法。 求解网格节点的梯度需要利用节点周围其它多个点进行求解,常用的方法有Horn梯度算子,Barron梯度算子等,此处介绍计算流体力学中具有二阶精度的有限差分方法求解梯度。 计算流体力学一般采用曲线坐标系(ξ,η,ζ)下的N-S方程求解流场,如图2所示,而光线追迹一般在笛卡儿坐标系(x,y,z)下进行,所以需要求解笛卡儿坐标下的折射率梯度▽n。在曲线坐标系下,首先运用有限差分方法求得折射率梯度▽n的三个分量▽nξ,▽nη,▽nζ,再利用Jacobi坐标变换矩阵得到笛卡儿坐标下的分量▽nx,▽ny,▽nz。 1)求解▽nξ,▽nη,▽nζ 以ξ方向为例,采用如下的差分形式求解网格点(i,j,k)处的折射率的梯度 左边界网格点采用三点向前偏心差分 中间网格点采用中心差分 右边界网格点采用三点向后偏心差分 上面各式中,下标i,j,k与ξ,η,ζ的关系如图2所示。 2)利用Jacobi坐标变换矩阵得到▽nx,▽n