石英灯加热固体发动机壳体辐射与温度场计算方法及系统

技术领域

本发明属于热辐射计算及传热学计算领域，主要涉及针对在辐射加热和大空间自然对流环境下的三维圆柱表面辐射热流密度及非稳态瞬态温度场的仿真计算，特别涉及一种石英灯加热固体发动机壳体辐射与温度场计算方法及系统。

背景技术

固体火箭发动机(SRM)等高速飞行器及其部件在大气中超音速飞行时，大气摩擦及阻碍等因素会使得发动机壳体表面(以下简称为“表面”)产生强烈的气动加热效应，导致表面产生极高的温度与温升速率，影响飞行器的安全性能。随着国防工业与航空航天技术的发展，需要对大量不同型号与结构的高速飞行器进行地面结构热试验。

石英灯辐射加热法凭借热惯性小、成本低、便于控制功率与调整结构等优点，在高速飞行器地面结构热试验中，获得了广泛的应用。Turner等基于蒙特卡洛法(Monte Carlomethod)模拟预测了带反射壁面的石英灯管的辐射热流分布(Turner T L,Ash RL.Numerical and experimental analyses of the radiant heat flux produced byquartz heating systems[R].NASA-TP-3387，1994.)，提供了较为完整的模拟计算与实验思路。但是，该研究只能提供同一平面等间距排列的石英灯组加热器在一水平面上的辐射热流密度。朱言旦等提出了石英灯阵在加热平面上的净辐射热流，考虑了辐射散热和对流的传热过程(朱言旦,曾磊,董威,等.石英灯阵热流分布规律计算与试验研究[J].宇航学报,2017,38(10).)，但是仍然没有考虑石英灯阵加热圆柱面等曲面的情况，也没有对辐射散热与对流效应深入研究。另外，利用Fluent等商用软件模拟计算表面温度时，难以准确模拟石英灯辐射热加载的边界条件，也难以实时输出瞬态温度场数据。截止目前，针对多层曲面排列的石英灯阵对圆柱面的蒙特卡洛光子追踪模拟计算方法缺少研究。将蒙特卡洛法模拟计算辐射热流密度与辐射散热、自然对流和导热等多种传热方式耦合，来计算表面瞬态温度场的研究尚未见报道。

深度强化学习(Deep Reinforcement Learning)技术能够实现将高维的表面温度数据矩阵导入到神经网络算法中进行自动训练、更新模型，对石英灯加热功率等参数进行控制，从而实现表面快速、均匀升温。但是，进行一次加热试验能够获得的表面温度场训练数据量有限，待表面冷却后再次加热也需要花费一定时间；为了得到大量训练数据与经验，还需要耗费大量的电能。为了克服这些缺点，可以考虑采用计算机模拟仿真计算得到的表面温度场结果对深度强化学习模型进行提前训练。

发明内容

本发明的目的在于提供一种石英灯加热固体发动机壳体辐射与温度场计算方法及系统，本发明将多种传热过程相耦合的三维柱体表面温度场的建模及仿真计算方法，以优化加热器结构，预测三维柱体表面温度分布情况与变化趋势，并为深度强化学习模型提供大量可靠的训练数据，模拟训练环境。

为达到上述目的，本发明采用了以下技术方案：

一种石英灯加热固体发动机壳体辐射与温度场计算方法，包括以下步骤：

针对多层环状排列的石英灯阵加热器及壳体试件组合，建立石英灯阵辐射物理模型，将物理模型抽象成数学模型；

采用蒙特卡洛光子追踪法模拟计算到达试件表面每个网格的光子总数；根据单个光子的平均功率，通过辐射场模拟计算获得多层石英灯阵加热器对壳体表面的辐射加热热流密度分布；

将辐射加热热流密度分布结果作为边界条件，与多种传热过程相耦合，通过交替方向隐式方法计算三维柱体内部及表面的瞬态温度分布数据。

作为本发明的进一步改进，针对多层环状排列的石英灯阵加热器及壳体试件组合，建立石英灯阵辐射物理模型，具体包括：

确定光子传播过程中的绝对坐标系和旋转坐标系及其变化关系；以壳体底面或者某一参考灯管底面的中心为坐标原点；

确定光子在界面间传播、入射、反射和折射的几何处理方程，用数学方程完整地描绘光子从发射到到达壳体表面的运行历程；

在选定的坐标系下，描述光子的发射坐标和发射方向，计算一个光子从灯丝发射，经历反射、折射、吸收和再发射过程后到达表面的交点坐标。

作为本发明的进一步改进，采用蒙特卡洛光子追踪法模拟计算到达试件表面每个网格的光子总数，需要对于壳体表面的网格化处理，具体为：

针对壳体表面的几何形状划分均匀网格，并按照网格数定义适当大小的总光子数，网格尺寸为0.1mm×0.1mm；模拟计算中设定的总光子数应为网格总数的300倍以上。

作为本发明的进一步改进，采用蒙特卡洛光子追踪法模拟计算到达试件表面每个网格的光子总数，具体包括：

循环计算大量光子的运行轨迹，统计到达表面每一个网格的光子总数，要求光子总数符合统计学要求；判断方法是：当成倍增加总光子数时，每个网格上光子数也近似呈相同倍数增长，且采用随机数多次模拟计算的结果近似相等。

作为本发明的进一步改进，通过辐射场模拟计算获得多层石英灯阵加热器对壳体表面的辐射加热热流密度分布，具体包括：

针对于辐射波段较窄的石英灯：

假设每个光子携带的能量相等；获得了落入壳体表面每个网格的局部光子数后，按照单个光子的平均功率以及壳体表面的吸收率，将光子数转化为壳体表面吸收的能量值；结合网格微元的面积，获得多层石英灯阵加热器对一定位置的壳体表面的局部辐射加热热流密度；

针对多层环状排列的尺寸和规格相同、排列相似的石英灯管：

假设每一根石英灯产生的辐射分布情况相同，依据其余灯管的相对位置，将获得的一根灯管的热流分布场进行旋转、平移并叠加，获得完整的辐射场结果矩阵；针对不同发热功率的石英灯层，将光子数转化为能量的步骤中，改变石英灯功率的参数，获得不同加热功率下多层石英灯阵对壳体表面的辐射加热热流密度。

作为本发明的进一步改进，所述通过交替方向隐式方法计算三维柱体内部及表面的瞬态温度分布数据，具体包括以下步骤：

对被加热的壳体进行建模；

在柱坐标下，根据控制容积法产生控制容积，根据外节点法产生节点，以节点温度代替整个控制容积的温度；确定控制容积的几何尺寸；

在圆柱坐标下，控制容积的周向尺寸是节点半径的函数；明确壳体所处的热环境，确定各边界条件的计算方法；其中，辐射加热热流密度的边界条件已经获得，辐射散热的边界条件由下式确定：

q

表面局部传热系数h

其中，Nu

通过求解局部努塞尔数Nu

确定了边界条件后，针对壳体内部的导热，需要将柱坐标下的三维导热偏微分方程离散化；将离散化后的导热方程的非稳态项和三个维度的扩散项整理成一系列三元一次方程，并确定式中各项系数的表达方式，利用交替方向隐格式分别在三个维度上更新计算了在每一个时间步长内温度的变化；分别计算三个方向上的中间温度时，分别对三个维度的中间温度进行求解后，获得经历了一个时间步长后壳体的三维温度分布。

作为本发明的进一步改进，所述瞬态温度分布数据在仿真计算辐射加热条件下壳体的温度场时，采用Brian-ADI格式，在求解ADI格式产生的离散方程组时，针对某一维度且始末端不重合的中间温度，需采用三对角阵算法；针对具有周期性的边界条件且计算区域的始末端重合的中间温度时，采用循环三对角阵算法。

一种石英灯辐射加热条件下火箭发动机壳体表面辐射场以及温度场的计算系统，包括：

辐射模型建立模块，用于针对多层环状排列的石英灯阵加热器及壳体试件组合，建立石英灯阵辐射物理模型，将物理模型抽象成数学模型；

辐射场计算模块，用于采用蒙特卡洛光子追踪法模拟计算到达试件表面每个网格的光子总数；根据单个光子的平均功率，通过辐射场模拟计算获得多层石英灯阵加热器对壳体表面的辐射加热热流密度分布；

温度场计算模块，用于将辐射加热热流密度分布结果作为边界条件，与多种传热过程相耦合，通过交替方向隐式方法计算三维柱体内部及表面的瞬态温度分布数据。

一种电子设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述石英灯加热固体发动机壳体辐射与温度场计算方法的步骤。

一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现所述石英灯加热固体发动机壳体辐射与温度场计算方法的步骤。

本发明的有益效果体现在：

本发明在数值传热学的基础上，提出了利用蒙特卡洛光子追踪法模拟石英灯阵对表面每一微小网格的辐射加热热流密度，并与自然对流、辐射散热和壳体内部导热等传热过程相耦合，计算了变功率辐射加热下表面瞬态温度场分布，实现了瞬态温度场分布数据的输出和深度强化学习算法进行自动模拟训练的目标。本方法具有计算速度快、模拟结果可靠、模型适应性强等特点。计算获得的瞬态温度场数据，可以供后续的数据分析以及导入深度强化学习的DQN神经网络算法进行模拟数据预训练。

进一步，对于辐射加热条件下的壳体表面，将采用蒙特卡洛法模拟获得的辐射分布数据，转化为热流密度分布；将热流密度分布结果矩阵应用到温度场模拟计算上，作为边界条件之一，并考虑了自然对流、辐射散热和壳体内部导热等复杂传热过程，以模拟真实的热加载环境，获得可靠的温度场数据。

附图说明

图1是采用蒙特卡洛光子追踪法模拟计算的流程图；

图2是光子在不同坐标系下的发射点及发射方向示意图；

图3是温度场耦合计算时柱坐标下的控制容积立体示意图和俯视图；

图4是壳体瞬态温度场耦合计算的流程图；

图5是壳体纵向剖面网格、节点及边界条件示意图。

图6为石英灯加热固体发动机壳体辐射与温度场计算系统结构示意图；

图7为本发明优选实施例电子设备结构示意图。

具体实施方式

下面将参考附图并结合实施例来详细说明本发明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

以下详细说明均是示例性的说明，旨在对本发明提供进一步的详细说明。除非另有指明，本发明所采用的所有技术术语与本申请所属领域的一般技术人员的通常理解的含义相同。本发明所使用的术语仅是为了描述具体实施方式，而并非意图限制根据本发明的示例性实施方式。

首先，针对多层环状排列的石英灯阵加热器及壳体试件组合，建立石英灯阵辐射物理模型。将实际物理模型抽象成数学模型。随后，采用蒙特卡洛光子追踪法模拟计算到达试件表面每个网格的光子总数。

其中，蒙特卡洛法中的几何处理代码被重构为由CUDA封装的代码，实现GPU加速；计算总光子数应当满足使得到达每个网格的光子数符合统计学要求。获得了落入壳体表面每个网格的局部光子数后，按照单个光子的平均功率，将光子数转化为壳体表面吸收的能量值，从而获得多层石英灯阵加热器对壳体表面的辐射加热热流密度。

得到表面辐射加热热流密度分布后，将其结果作为边界条件，与对流、辐射散热和壳体内部导热等复杂传热过程相耦合，通过交替方向隐式方法(ADI)计算三维柱体内部及表面的瞬态温度分布。其中，为了满足较小网格和较大时间步长的条件，需采用无条件稳定的Brian-ADI格式。在求解ADI格式产生的离散方程组时，针对具有环形特征的壳体，可以采用三对角阵算法(TDMA)和循环三对角阵算法(CTDMA)计算每一个维度下每一个时间步长中的中间温度。设定需要仿真计算的总时间后，按时间步长推进并循环计算，获得瞬态温度分布结果。

计算完成后，导出并保存瞬态温度场数据，以供后续的数据分析以及导入深度强化学习的DQN神经网络算法进行模拟数据预训练。

为了提高运算速度，计算程序采用C++语言进行编写；另外，对于辐射场分布的计算，蒙特卡洛法中的几何处理代码被重构为由CUDA封装的代码，实现GPU加速，以减少模拟计算中并行运算部分的所需的时间；在瞬态温度场的计算中，为了将每个时间步长计算得到的瞬态温度场结果传输到DQN算法的神经网络中，实时进行模拟数据预训练，在Python中使用ctypes库调用C++程序，将C++计算获得的存有温度场数据的三维数组实时传递到Python中。

下面结合附图和实施例对本发明做详细描述。

如图1的计算流程图所示，利用蒙特卡洛法计算到达表面的石英灯阵的辐射加热热流密度的基本思路是：依据石英灯阵尺寸和排布方式，以及壳体的尺寸与相对位置的物理模型，抽象出数学模型，包括如下步骤：

首先，如图2所示，确定光子传播过程中的绝对坐标系和旋转坐标系及其变化关系，如以被加热件或某一石英灯管的底面的中心等为原点，要求使得几何处理更加简便；

其次，确定光子在界面间传播、入射、反射和折射的几何处理方程，用数学方程完整地描绘光子从发射到到达壳体表面的运行历程；另外，需要针对壳体表面的几何形状，进行均匀化网格划分，要求网格划分均匀，便于计算，并按照网格数定义适当大小的总光子数。为了保证模拟计算精度，网格尺寸应在0.1mm×0.1mm左右；考虑到有部分光子无法到达目标表面，总光子数应为网格总数的300倍以上。

随后，在选定的坐标系下，描述光子的发射坐标和发射方向，计算一个光子从灯丝发射，经历反射、折射、吸收和再发射等复杂过程后到达表面的交点坐标。

由于大量光子的运算相互独立，即每个光子的轨迹追踪采用并行计算，本发明采用了C++编写运算程序，并将串行计算的几何处理代码重构为由CUDA封装的代码，实现GPU加速计算，可以极大地加快运算速度。循环计算大量光子的运行轨迹，统计到达表面每一个网格的光子总数，要求光子总数符合统计学要求。判断满足要求的方法是，当成倍增加总光子数时，每个网格上光子数也近似呈相同倍数增长，且采用随机数多次模拟计算的结果近似相等。

对于辐射波段较窄的石英灯，可以假设每个光子携带的能量相等。获得了落入壳体表面每个网格的局部光子数后，按照单个光子的平均功率以及壳体表面的吸收率，将光子数转化为壳体表面吸收的能量值；结合网格微元的面积，可以获得多层石英灯阵加热器对一定位置的壳体表面的局部辐射加热热流密度。

其中，为了加快模拟计算速度，针对多层环状排列的尺寸和规格相同、排列相似的石英灯管，本发明提出了如下的快速计算方法：假设每一根石英灯产生的辐射分布情况相同，依据其余灯管的相对位置，将上述获得的一根灯管的热流分布场进行旋转、平移并叠加，最终获得完整的辐射场结果矩阵；针对不同发热功率的石英灯层，只需在将光子数转化为能量的步骤中，改变石英灯功率的参数，即可获得不同加热功率下多层石英灯阵对壳体表面的辐射加热热流密度。

得到表面辐射加热热流密度分布后，将其结果作为边界条件，与对流、辐射散热和导热等复杂传热过程相耦合，通过交替方向隐式方法(ADI)计算三维柱体内部及表面的瞬态温度分布。计算流程如图3所示。具体包括以下步骤：

首先，对被加热的壳体进行建模，如对于圆柱壳体，通常以底面的几何中心为坐标原点建立圆柱坐标系；

其次，如图4(a)所示，在柱坐标下，根据控制容积法产生控制容积，根据外节点法产生节点，以节点温度代替整个控制容积的温度；同时，结合图4(a)与(b)，确定控制容积的几何尺寸，要求控制容积的尺寸与上述辐射场仿真计算时确定的表面网格尺寸相近，在圆柱坐标下，控制容积的周向尺寸是节点半径的函数；另外，如图5所示，明确壳体所处的热环境，确定各边界条件的计算方法。其中，辐射加热热流密度的边界条件已经获得，辐射散热的边界条件可由下式确定：

q

为了较为准确地计算表面局部传热系数h

其中，Nu

通过求解局部努塞尔数Nu

确定了边界条件后，针对壳体内部的导热，需要将柱坐标下的三维导热偏微分方程离散化；将离散化后的导热方程的非稳态项和三个维度的扩散项整理成一系列三元一次方程，并确定式中各项系数的表达方式，利用交替方向隐格式分别在三个维度上更新计算了在每一个时间步长内温度的变化；为了满足较小网格和较大时间步长的条件，本发明采用无条件稳定的Brian-ADI格式。对于圆柱坐标系下的壳体采用Brian-ADI方法，分别计算三个方向上的中间温度时，计算某一维度且始末端不重合的中间温度时，如圆柱壳体的轴向与径向方向，可以采用三对角阵算法(TDMA)；计算具有周期性的边界条件且计算区域的始末端重合的中间温度时，如圆柱壳体周向方向，需要采用循环三对角阵算法(CTDMA)。使用该算法分别对三个维度的中间温度进行求解后，即可获经历了一个时间步长后壳体的三维温度分布。

上述步骤详细地刻画了本发明提出的针对石英灯辐射加热条件下三维柱体表面的温度场仿真计算的具体实施方法。计算结束后，需要保存在辐射加热条件下的瞬态温度场数据，以供后续的数据分析以及导入深度强化学习的神经网络中进行模拟数据预训练使用。为了实现将每个时间步长计算得到的瞬态温度场结果传输到DQN算法的神经网络中，实时进行模拟数据预训练，在Python中使用ctypes库调用了C++程序。即在模拟训练时，调用C++计算瞬态温度场，将计算获得的存有温度场数据的三维数组实时传递到Python中的神经网络框架中。

本发明的特点如下：针对复杂的多层石英灯阵辐射加热条件下的辐射场进行仿真计算时，采用了蒙特卡洛光子追踪法，考虑了光子在石英灯内部的发射、反射、吸收、再发射以及折射等过程；针对三维柱体表面瞬态温度场的变化进行预测时，将辐射场的结果作为辐射加热边界条件，并将辐射加热、自然对流、辐射散热以及壳体内部导热等复杂传热过程相耦合。

总之，本发明提出的基于蒙特卡洛光子追踪法的对复杂石英灯阵加热器产生的辐射场的模拟计算，考虑了光子运行过程中的多种效应，能够较为精确地获得各种复杂石英灯阵排列方式下的壳体表面的辐射分布以及辐射热流密度分布；将上述的结果作为边界条件，本发明提出的温度场耦合计算方法，能够获得不同加热条件下，在复杂传热环境中的圆柱壳体的瞬态三维温度场模拟计算结果。该结果不仅能够指导石英灯阵加热器的优化设计，以达到一定的加热需求，并与试验结果进行对比，改进模型，还能够将准确可靠的模拟数据导入到深度强化学习中的DQN算法的神经网络中进行实时预训练，以提高在缺乏实验数据的条件下模型训练的效率与效果。

如图6所示，本发明的另一目的在于提出一种石英灯加热固体发动机壳体辐射与温度场计算系统，包括：

辐射模型建立模块，用于针对多层环状排列的石英灯阵加热器及壳体试件组合，建立石英灯阵辐射物理模型，将物理模型抽象成数学模型；

辐射场计算模块，用于采用蒙特卡洛光子追踪法模拟计算到达试件表面每个网格的光子总数；根据单个光子的平均功率，通过辐射场模拟计算获得多层石英灯阵加热器对壳体表面的辐射加热热流密度分布；

温度场计算模块，用于将辐射加热热流密度分布结果作为边界条件，与多种传热过程相耦合，通过交替方向隐式方法计算三维柱体内部及表面的瞬态温度分布数据。

如图7所示，本发明第三个目的是提供一种电子设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述石英灯加热固体发动机壳体辐射与温度场计算方法的步骤。

所述石英灯加热固体发动机壳体辐射与温度场计算方法包括以下步骤：

针对多层环状排列的石英灯阵加热器及壳体试件组合，建立石英灯阵辐射物理模型，将物理模型抽象成数学模型；

采用蒙特卡洛光子追踪法模拟计算到达试件表面每个网格的光子总数；根据单个光子的平均功率，通过辐射场模拟计算获得多层石英灯阵加热器对壳体表面的辐射加热热流密度分布；

将辐射加热热流密度分布结果作为边界条件，与多种传热过程相耦合，通过交替方向隐式方法计算三维柱体内部及表面的瞬态温度分布数据。

本发明第四个目的是提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现所述石英灯加热固体发动机壳体辐射与温度场计算方法的步骤。

所述石英灯加热固体发动机壳体辐射与温度场计算方法包括以下步骤：

针对多层环状排列的石英灯阵加热器及壳体试件组合，建立石英灯阵辐射物理模型，将物理模型抽象成数学模型；

采用蒙特卡洛光子追踪法模拟计算到达试件表面每个网格的光子总数；根据单个光子的平均功率，通过辐射场模拟计算获得多层石英灯阵加热器对壳体表面的辐射加热热流密度分布；

将辐射加热热流密度分布结果作为边界条件，与多种传热过程相耦合，通过交替方向隐式方法计算三维柱体内部及表面的瞬态温度分布数据。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。