本研究由成功大学机械工程学系的曾弘毅（Hung-Yi Tseng）、洪翊斌（Yi-Bin Hong）、吴志阳（Chih-Yang Wu）、洪道驰（Dao-Chi Hong）与陈彦辰（Yen-Chen Chen）共同完成，发表于 Applied Mathematical Modelling 期刊，文章于2017年3月11日在线发表。该研究属于辐射传热（radiative transfer）与计算光学领域，核心目标是发展一种能够在折射率（refractive index）仅于离散点上已知的介质中进行辐射传输模拟的蒙特卡洛方法（Monte Carlo method）。
在自然界与工程应用中，许多介质的折射率随空间位置连续变化，例如激光引起的大气热晕（thermal blooming）、生物组织的光学成像以及大气遥感等场景。当折射率连续变化时，辐射能量不再沿直线传播，而是沿弯曲路径行进。传统辐射传输方程（radiative transfer equation, RTE）的数值求解方法，如离散坐标法（discrete ordinates method, DOM）和积分方程法（integral equation, IE），通常要求折射率分布具有显式解析表达式。然而在实际测量或耦合计算中，折射率常常只能通过离散采样点获得。以往已有学者采用多项式插值结合蒙特卡洛法处理离散折射率分布，但相关工作较少且方法适用范围有限。因此，本文旨在开发一种通用的、基于离散射线追踪（discrete ray tracing）的蒙特卡洛方法，用于处理折射率仅在一组任意离散点上给定的辐射传输问题，并为后续复杂几何与散射条件下的模拟提供工具。
为了实现这一目标，作者建立了一套完整的数值模拟流程。首先，研究采用了几何光学中的射线方程（ray equation）描述弯曲的光线轨迹。为便于数值积分，作者引入了Sharma等人提出的变量变换，将射线方程改写为关于新变量t的一阶常微分方程组，其中包括位置矢量、光学射线矢量（optical ray vector）以及路径长度（path length）的演化方程。随后，作者选择Runge-Kutta Dormand-Prince（RKDP）格式对该方程组进行逐步积分。RKDP方法能够在每一步同时获得四阶与五阶精度的解，并利用两者的差异估计局部误差，从而自动调整步长，以保证轨迹计算精度与效率的平衡。在射线追踪过程中，每一个Runge-Kutta节点上都需要知道折射率的值与梯度，而这些信息必须从离散采样点中重建。
针对一维问题，作者采用三次样条插值（cubic spline interpolation）方法。该方法在每个区间内使用三次多项式逼近折射率分布，并要求节点处函数值、一阶导数与二阶导数连续，同时在边界处令二阶导数为零，从而唯一确定插值多项式。对于二维问题，作者采用移动最小二乘（moving least square, MLS）方法。MLS方法并非严格通过所有数据点，而是通过局部加权最小二乘拟合来逼近离散折射率分布，更适合处理非结构采样点或含有测量误差的数据。在MLS方法中，作者比较了二次基函数与三次基函数的效果，并研究了采样数据点数量对拟合精度的影响。研究还引入了kd-tree（k-dimensional tree）数据结构加速最近邻搜索，以提高局部拟合的计算效率。MLS权重函数使用了Most与Bucher推荐的带正则化参数的函数形式，正则化参数ε_w取为10⁻⁵。
在完成折射率重建后，作者将其嵌入完整的蒙特卡洛模拟框架中。对于一维平面介质中的辐射平衡（radiative equilibrium）问题，介质被划分为若干子体积，从热边界发射的大量光子束（photon bundles）携带由边界温度与折射率决定的能量。光子束在介质中经历吸收、散射等随机过程，其轨迹由数值射线追踪确定。辐射平衡条件下，每个子体积吸收的能量等于其发射的能量，由此可迭代确定介质内部温度分布。对于给定温度场的二维介质，作者从介质体积内发射光子束，并统计到达各个边界面元的能量，从而计算边界辐射通量。模拟中考虑了线性各向异性散射相函数（scattering phase function），边界假设为黑体表面。
数值验证结果表明，该方法具有较高的精度与可靠性。在一维平面介质情形中，作者比较了使用不同离散点数（n_k+1从6到26）得到的无量纲辐射通量。结果显示，当离散点数达到16或21以上时，结果趋于稳定，与基于积分方程法的基准解相比，偏差不超过0.0003至0.0004。对于不同折射率分布（包括m=-1、1、3的参数化形式）、不同散射反照率（scattering albedo）以及不同光学厚度（optical thickness）的情况，本方法均能取得一致的高精度结果。温度分布与DOM方法的结果也吻合良好。
在二维射线追踪精度测试中，作者选取了具有解析射线路径解的折射率分布，并计算了平均位置误差与平均路径长度误差。结果表明，采用三次基函数的MLS方法明显优于二次基函数，且采样数据点数量存在一个最优值。当采用三次基函数并使用144个采样点时，位置误差与路径长度误差均显著降低，且数值射线追踪仅在求解射线方程方面就表现出很高的效率。然而，从随机离散点重建折射率的过程占用了较多计算时间，说明MLS拟合是二维方法中主要的计算瓶颈。
在二维梯形介质辐射传输模拟中，作者以折射率线性分布n(x,y)=1+2(x+y)/h为例，比较了不同光子束数量（10⁷与10⁸）下的边界辐射通量。结果显示，10⁸个光子束已足够保证结果稳定。与基于显式折射率分布求解积分方程得到的基准解相比，本方法的模拟结果几乎完全吻合。散射反照率越低，介质发射功率越大，边界通量也随之增大。由于梯形介质几何形状的影响，底部边界通量在靠近角点处下降，最大值出现在中间偏右的位置。
作者还模拟了高宽比很大的矩形介质（h/w=10）中的辐射平衡问题。折射率仅在一组随机点上给定。模拟结果表明，在距左右冷边界较远的中心区域，温度分布与一维平面介质的结果几乎一致；而在靠近侧边界的区域，由于能量向冷侧壁逃逸，温度显著降低。顶部边界通量也呈现类似趋势。该结果进一步验证了MLS近似折射率分布在随机离散点条件下的有效性。
此外，研究还模拟了具有保守散射（ω=1）的方形介质，折射率分布采用径向对称形式。模拟结果与积分方程解吻合。在最后一个算例中，作者考虑了折射率具有高斯型径向分布的方形介质，温度场针对不同α参数（α=2、8、32）进行了模拟，并与DOM结果进行了对比，验证了方法对不同折射率梯度分布形态的适应能力。
本研究的结论指出，所发展的基于离散射线追踪的蒙特卡洛方法能够有效处理折射率仅在离散点给定的辐射传输问题。对于一维问题，三次样条插值能够提供足够精度的折射率重建；对于二维问题，MLS方法结合RKDP射线追踪能够在随机采样点条件下获得高精度结果。文中所有算例均表明，采用10⁸个光子束以及适当的离散折射率重建策略时，该方法的结果与基于显式折射率分布的积分方程法或离散坐标法解高度一致。
本文的亮点在于：第一，提出了一种将RKDP自适应步长数值积分与离散折射率重建相结合的蒙特卡洛方法，适用于任意几何形状和复杂折射率分布；第二，系统比较了MLS方法中基函数阶数与采样点数量对射线追踪精度的影响，为实际应用提供了参数选择依据；第三，验证了方法在一维辐射平衡、二维辐射传输以及具有强折射率梯度介质中的广泛适用性。该研究为折射率仅能通过离散测量或网格计算获得的实际工程问题提供了有效且通用的辐射传热分析工具，具有重要的科学意义与应用价值。