© 2025 Journal of Northwestern Polytechnical University. All rights reserved.

导弹结构一般由天线罩、翼舵面以及仪器舱、战斗部舱等舱段组成, 其中舱段结构主要为等直径舱段。对于高速、高超声速导弹, 采用石英灯辐射加热的地面热试验考核是验证高温环境下其舱段结构强度的重要手段, 也是开展飞行试验前的关键环节。石英灯阵列一般采用随形排布, 等直径舱段的石英灯阵列呈圆柱形排布, 设计参数包括石英灯长度、灯管额定功率、灯管环向排布密度、灯管距舱段表面高度以及灯管相对舱段端面的轴向距离等。石英灯阵列排布会直接影响舱段表面热环境的模拟质量, 进而影响试验结果的参考价值, 因此需要计算石英灯阵列在舱段表面产生的辐射热流场, 评估加热阵列设计的合理性。

石英灯辐射热流计算方法种类较多, 孔凡金等[1]比较了包括区域法(zone method)、离散坐标法(discrete ordinate method)、蒙特卡洛法(Monte Carlo method)等在内的常见计算方法以及Fluent、RadTherm、Sinda/Fluint等常用计算软件, 分析了不同热流场计算模型的特点, 指出应结合具体的计算要求和精度目标进行选择, 完成对石英灯辐射热流场的模拟与评价。

蒙特卡洛法能够对光线在石英玻璃壁面以及反射屏作用下的折/反射路径变化进行仿真, 计算精度相对较高, 目前广泛应用于辐射加热场景的热流计算, 文献中多以单灯和平板形灯阵作为计算对象[2-5], 缺点是需要包括光学特性参数在内相对较多的输入参数, 高精度的计算结果依赖于大数量的样本, 计算时间也相对较长, 通常达到小时级和天级, 而热试验设计方案往往需要在分钟级的短时间内完成迭代, 因此蒙特卡洛法难以满足导弹舱段热试验设计中的快速迭代要求。

解析法一般通过简化研究对象, 使其符合兰贝特定律等传热学基本定律所描述的应用条件, 优势是计算过程基本采用显式计算, 容易收敛, 因此计算速度相对较快, 计算时间相对较少, 可以根据用户需求实现自定义功能拓展, 不足则是忽略了传播过程中的光线折/反射效应。在石英灯阵列的快速设计迭代中可不考虑光线折/反射效应, 因此解析法在等直径规则形状的舱段热试验设计中能够快速给出热流分布结果, 可在分钟级的短时间内完成设计方案多次迭代, 能够显著提高设计效率、降低试验风险, 同时便于开发参数化计算工具, 与有限元软件相比不需要舱段材料参数、节省模型前处理时间、计算过程更易收敛, 也没有部分软件中设定的辐射传热网格数量限制, 适合大尺度舱段结构的辐射热流场计算。李昕昕等[6]提出“等效面光源”概念, 融合射线追踪法与解析法特点, 完成了考虑反射屏作用下的热流场解析计算, 但该方法目前只适用于平面形状的辐照面。张肖肖等[7]基于解析法提出适用于3种典型形状石英灯阵列的热流快速计算方法, 将灯丝半圆柱面简化为平面微元, 虽然显著提高了计算速度, 但计算精度相对较低。为了提高计算精度, 又提出了针对平板形阵列灯丝半圆柱面多网格划分时的热流密度解析计算方法[8]，该方法可以用于导弹弹翼及舵面结构平板形石英灯阵列的热流预测, 但不能用于等直径舱段的热流计算, 灯丝半圆柱面采用多网格划分时，舱段圆柱曲面表面上的照射范围判断要更加复杂。

本文拟建立灯丝半圆柱面采用多网格划分时等直径舱段表面热流密度解析计算方法, 推导圆柱形阵列中单根石英灯照射范围判据, 并对比有限元软件热流计算结果, 对该解析计算方法进行验证。

1 多网格解析计算理论推导

研究思路如图 1所示, 首先根据圆柱形阵列辐射加热下的等直径舱段物理模型, 建立单根灯管灯丝半圆柱面照射下的舱段圆柱表面几何模型, 然后结合解析几何理论, 推导出灯丝半圆柱面采用多网格划分时，单个微元在舱段圆柱表面的照射范围判据, 对舱段表面位于照射范围内的接收点, 按照兰贝特定律计算单个微元对该接收点辐射产生的热流密度, 进而得到单根灯管作用下的热流场, 将各支灯管在舱段表面产生的热流场叠加便可得到圆柱形石英灯阵列在等直径舱段表面的辐射热流场。

图1 等直径舱段表面热流密度解析计算方法研究思路

1.1 圆柱外加热阵列

圆柱外加热阵列多网格照射范围计算模型如图 2所示。可以看到, 对于每一个灯丝微元中心E而言, 存在2个照射范围:

1) 从E点向舱段圆柱面做2条切线EF、EG, 2个切点间的圆弧为该微元在舱段表面理论上所能照射的最大范围;

2) 延长微元所在直线BE, 若BE与圆柱面相交于C, 则圆弧 为该微元实际在舱段表面所能照射的范围; 若BE与圆柱面不相交, 则圆弧 为该微元实际在舱段表面所能照射的范围。

因此在已知某根石英灯在阵列中的空间位置(即图 2中OD相对x轴的逆时针环向角度∠xOD)时, 要确定某个微元在舱段圆柱面上的照射范围, 首先需要确定BE是否与圆柱面相交, 若相交则需要确定以下3个角的大小, 分别是:

1) 圆柱面圆心O和微元中心E的连线OE, 与圆柱面圆心O和灯丝圆心D的连线OD之间的夹角, 即图 2中的θ;

2) 圆柱面圆心O和交点C的连线OC, 与圆柱面圆心O和灯丝圆心D的连线OD之间的夹角, 即图 2中的θ₁;

3) 圆柱面圆心O和微元中心E的连线OE, 与圆柱面圆心O和切点F的连线OF之间的夹角, 即图 2中的θ₂。

以上3个角度描述了圆弧对应的圆心角∠COF的大小, 同时, 描述了照射范围在灯丝空间位置OD连线左右两侧的分布, 对于图 2中微元中心E, 其在舱段表面上的照射范围为从比D点环向角度小|θ₁|到比D点环向角度大|θ|+θ₂的圆弧段, 照射范围大小为|θ|+|θ₁|+θ₂; 相应地, 对于与E点对称的微元中心K, 其在舱段表面上的照射范围为从比D点环向角度小|θ|+θ₂到比D点环向角度大|θ₁|的圆弧段。而对于所在微元平面与圆柱面并不相交的微元中心L, 其在舱段表面上的照射范围为从比D点环向角度小θ₂-|θ|到比D点环向角度大|θ|+θ₂的圆弧段, 照射范围大小为2θ₂; 相应地, 对于与L点对称的微元中心M, 其在舱段表面上的照射范围为从比D点环向角度小|θ|+θ₂到比D点环向角度大θ₂-|θ|的圆弧段。

如图 2所示, 考虑到上述3个角度θ, θ₁, θ₂对于同一阵列中不同环向位置灯管相同顺序的微元均相等, 因此采用图 2中所示位置(∠xOD为90°)对上述角度值进行计算。灯丝圆心D坐标为(0, R+d, 0), 第i个微元中心E的坐标可写为{rcos(π/2n)sin[(i-0.5)·π/n-π/2], R+d-rcos(π/2n)cos[(i-0.5)·π/n-π/2], 0}(此时可保证E点位于与舱段圆柱面相对的灯丝半圆柱面上), 其中r代表灯丝半圆柱面半径, R为舱段圆柱面半径, d为灯丝中心距离圆柱面的高度, n代表灯丝半圆柱面划分网格的数量(图 2中n=4), π/2n代表α的大小, (i-0.5)·π/n-π/2代表∠ODE的大小及符号(为正时表示DE在OD的逆时针方向)。已知E点坐标后, 可以计算得到θ₂的大小。

(1)

θ的大小可以通过对ΔOED运用正弦定理或余弦定理求得, 采用正弦定理时

(2)

(2) 式代表了∠EOD的大小及符号(为正时表示OE在OD的逆时针方向)。

在ΔCOI中, OI线段的长度以及∠OIC、∠OCI(钝角)的大小分别为

(3)

(4)

(5)

计算得到l_OI以及∠OCI之后, 可以得到θ₁

(6)

式代表了∠COD的大小及符号(为正表示OC在OD的逆时针方向)。

判断直线BE是否与圆柱面相交, 则需要比较∠OEG与∠OEI的相对大小, 若∠OEG>∠OEI, 则BE与圆柱面相交; 若∠OEG < ∠OEI, 则BE与圆柱面不相交。考虑到

(7)

(8)

因此通过比较∠OHG与∠EID的相对大小, 也可判断BE是否与圆柱面相交。当BE恰好与圆柱面相切时, C与G重合, I与H重合, 此时

(9)

(10)

因此对于第i个微元所在的直线, 若π/2-|(i-0.5)·π/n-π/2| < arcsin(R/l_OI), 则BE与圆柱面相交; 若π/2-(i-0.5)·π/n-π/2 >arcsin(R/l_OI), 则BE与圆柱面不相交。

计算热流时先判断BE是否与圆柱面相交, 然后采用(1)~(2)式和(6)式计算得到θ, θ₁, θ₂, 再结合当前灯丝所在位置D点环向角度∠xOD即可确定某段平面微元在舱段圆柱外表面上的照射范围。对于位于当前微元在舱段表面可照射范围内的接收点J, 其受该微元照射获得的热流密度为

(11)

式中: Q_e为将灯丝沿长度方向N等分，环向采用n个网格进行等分后每一段微元在其半球空间内对外辐射产的总热流, Q_e=πrL·σT⁴/(N·n), 其中r为灯丝半径, L为灯管长度, σ为斯蒂芬玻尔兹曼常数, T为灯丝温度; θ₃为接收点J点处法线与JE连线的夹角; θ₄为微元中心E点处法线与JE连线的夹角。θ₃和θ₄的余弦值可通过向量点积计算得到, 例如将微元中心E点处的法向量记为a(a可写为{rcos(π/2n)sin[(i-0.5)·π/n-π/2+∠DOx-π/2], -rcos(π/2n)cos[(i-0.5)·π/n-π/2+∠DOx-π/2], 0}), 圆柱面上J点处法向量记为b(b可写为(Rcos∠JOx, Rsin∠JOx, 0)), J, E之间向量记为c, 则cosθ₄=|c·a/(|c|·|a|)|, cosθ₃=|c·b/(|c|·|b|)|。

将同一灯管各段平面微元对J点产生的热流密度累加, 得到在单根石英灯作用下J点的热流密度。对于多根石英灯组成的阵列, 可以按照上述过程计算每根石英灯对各接收点产生的热流密度, 累加后便得到不同网格数量划分下圆柱外加热阵列的热流密度场分布情况。

图2 圆柱外加热阵列单根灯管多网格照射范围计算模型

1.2 圆柱内加热阵列

圆柱内加热阵列多网格照射范围计算模型如图 3所示。

可以看到, 对于每一个灯丝微元中心E而言, 其在舱段表面上的照射范围为直线BE与圆柱面产生的2个交点G, F之间形成的圆弧 , 即其在舱段表面上的照射范围为从比D点环向角度小|θ₂|到比D点环向角度大|θ₁|的圆弧段, 照射范围大小为|θ₁|+|θ₂|; 相应地，对于与E点对称的微元中心K, 其在舱段表面上的照射范围为从比D点环向角度小|θ₁|到比D点环向角度大|θ₂|的圆弧段。

在ΔEDI中

(12)

(13)

在ΔOHI中

(14)

在ΔOHG中

(15)

因此, θ₁和θ₂可表示为

(16)

(17)

上述公式分别代表了∠COF与∠GOC的大小及符号(为正时表示OG和OF在OD的逆时针方向)。

计算热流时采用(15)~(17)式得到θ, θ₁, θ₂, 再结合当前灯丝所在位置D点环向角度∠xOD即可确定某段平面微元在舱段圆柱内表面上的照射范围。对于位于当前微元在舱段表面可照射范围内的接收点J, 其受该微元照射获得的热流密度可参照(11)式进行计算, 区别在于对于内加热阵列, 微元中心E处法向量a为{rcos(π/2n)sin[(i-0.5)·π/n+π/2+∠DOx-π/2], -rcos(π/2n)cos[(i-0.5)·π/n+π/2+∠DOx-π/2], 0}, 圆柱面上J点处法向量b为(-Rcos∠JOx, -Rsin∠JOx, 0)。

图3 圆柱内加热阵列单根灯管多网格照射范围计算模型

2 计算结果

2.1 圆柱外加热阵列

表 1给出了n=4时圆柱外加热阵列8个典型环向位置处单根灯管不同顺序的平面微元照射范围(环向从0°开始, 每间隔1°取一个接收点, 取至359°), 其中灯丝直径为4 mm, 灯丝长度为310 mm, 舱段圆柱面高350 mm、直径400 mm, 灯丝中心距圆柱面高度为60 mm, 灯丝温度2 500 K。可以看到, 各平面微元的照射范围随着所在灯管环向位置和微元顺序的改变, 按θ, θ₁, θ₂的计算值呈现出规律性变化, 与1.1节中关于照射范围的描述一致, i=1和i=4的微元所在平面与舱段圆柱面相交, i=1的微元照射范围为比灯管环向角度小6.73°(|θ₁|)到比环向角度大39.91°(|θ|+θ₂)的圆弧段; i=4的微元照射范围为比灯管环向角度小39.91°(|θ|+θ₂)到比环向角度大6.73°(|θ₁|)的圆弧段; i=2和i=3的微元所在平面与舱段圆柱面不相交, i=2的微元照射范围为比灯管环向角度小39.10°(θ₂-|θ|)到比环向角度大39.42°(|θ|+θ₂)的圆弧段, i=3的微元照射范围为比灯管环向角度小39.42°(|θ|+θ₂)到比环向角度大39.10°(θ₂-|θ|)的圆弧段。表 2给出了不同网格数量下各微元在舱段外表面的照射范围, 可以看到照射范围的大小主要受到微元所在平面和舱段表面是否相交影响, 总体呈现两端微元照射范围小, 中间微元照射范围大; 位于两端的微元(i=1和i=n)的微元照射范围随着n的增大逐渐减小; n为奇数时, 处于中间位置的微元(i=(n+1)/2)的照射范围也随着n的增大逐渐减小; n为偶数时, 处于中间位置微元(i=n/2和i=n/2+1)的照射范围也随着n的增大逐渐减小。

图 4给出了n=1和n=8时单根石英灯的热流分布云图, 图 5给出了n=1~8时灯管正下方母线ϕ=180°和环向圆周z=175 mm上的热流分布曲线。总体来看:

1) 随着n的增大, 灯丝正下方区域的热流水平呈下降趋势, 这是由于随着网格数的增大, 辐射到灯丝正下方区域的辐射能量不断减小, 更多的辐射能量向周围空间传播;

2) 随着n的增大, 舱段表面上的照射范围没有明显变化, 这是由于外加热阵列中灯管的照射范围是由θ₂主导的, 灯丝直径与舱段直径以及灯丝中心距圆柱面高度相比很小, 因此不同网格数下θ₂变化不大;

3) 当n≥4后, 热流密度变化趋于平缓, 热流分布趋于稳定。

图4 不同网格数量下热流分布云图(圆柱外加热阵列, 单根灯管)

图5 不同网格数量下热流分布曲线(圆柱外加热阵列, 单根灯管)

2.2 圆柱内加热阵列

表 3给出了n=4时圆柱内加热阵列8个典型环向位置处灯管不同顺序的平面微元照射范围, 计算参数同2.1节。可以看到各平面微元的照射范围随着所在灯管环向位置和微元顺序的改变, 按照θ, θ₁, θ₂的计算结果呈现规律性变化, 与1.2节中关于照射范围的描述一致, i=1的微元照射范围为比灯管环向角度小141.41°(|θ₂|)到比环向角度大6.41°(|θ₁|)的圆弧段; i=4的微元照射范围为比灯管环向角度小6.41°(|θ₁|)到比环向角度大141.41°(θ₂|)的圆弧段; i=2的微元照射范围为比灯管环向角度小71.51°(|θ₂)到比环向角度大26.51°(|θ₁|)的圆弧段, i=3的微元照射范围为比灯管环向角度小26.51°(|θ₁|)到比环向角度大71.51°(|θ₂|)的圆弧段。表 4给出了不同网格数量下各微元在舱段内表面的照射范围, 可以看到照射范围的大小与微元顺序相关, 总体呈现两端微元照射范围大, 中间微元照射范围小; 位于两端(i=1和i=n)的微元照射范围随着n的增大逐渐增大; n为奇数时, 处于中间位置的微元(i=(n+1)/2)的照射范围随着n的增大逐渐减小; n为偶数时, 处于中间位置的微元(i=n/2和i=n/2+1)的照射范围也随着n的增大逐渐减小。图 6给出了n=1和n=8时单根石英灯的热流分布云图, 图 7给出了n=1~8时灯管正下方母线ϕ=180°和环向圆周z=175 mm上的热流分布曲线。与外加热阵列不同的是随着n的增大, 舱段表面上的照射范围呈现较为明显的变化, 这是由于内加热阵列中灯管的照射范围是由θ₁和θ₂共同决定, 当网格数增大时, θ₁和θ₂受到直接影响, 变化较为明显。

图6 不同网格数量下热流分布云图(圆柱内加热阵列, 单根灯管)

图7 不同网格数量下热流分布曲线(圆柱内加热阵列, 单根灯管)

3 有限元计算结果对比

采用2.1节中的石英灯外加热阵列与舱段尺寸参数, 在有限元软件中建立单灯内加热与外加热模型(见图 8), 灯丝半圆柱面环向划分网格数量n=8, 轴向网格数量为100, 舱段表面环向网格数量为120, 轴向网格数量为30, 舱段温度设置为1 K, 环境温度为0 K。采用半立方体(hemicube)算法开展面面辐射传热计算。根据文献[8], 为保证有限元计算过程中灯丝对外辐射总能量与解析法一致, 需要将有限元计算模型中的灯丝温度T调整为(18)式中的T_e。

(18)

计算得到的T_e为2 504 K。外加热热流分布云图如图 9所示。提取灯管正下方母线和环向圆周z=175 mm上节点的热流计算结果, 与解析法计算结果进行对比, 环向和轴向对比结果如图 10~11所示。可以看到, 计算结果的一致性良好, 外加热时轴向热流结果相对误差小于1%;内加热时中间位置节点轴向热流结果相对误差小于1%, 靠近两端位置节点处相对误差小于3%。

图8 有限元计算模型(外加热)

图9 热流分布云图(外加热)

图10 热流计算结果对比(外加热)

图11 热流计算结果对比(内加热)

4 结论

针对导弹等直径舱段圆柱形石英灯阵列辐射热流的快速计算, 提出了灯丝半圆柱面采用多网格划分时的辐射热流解析计算方法, 推导了单根石英灯的照射范围判据, 分析了不同网格数量下的热流密度分布差异, 并与有限元软件计算结果进行对比, 得到了以下结论:

1) 对于外加热阵列, 若微元所在平面与舱段表面相交, 则照射范围大小为|θ|+|θ₁|+θ₂, 不相交时微元照射范围大小为2θ₂; 对于内加热阵列, 微元照射范围大小为|θ₁|+|θ₂|;

2) 对于圆柱形加热阵列, 随着n的增大, 灯丝正下方区域的峰值热流水平呈下降趋势, 为保证石英灯辐射热流计算精度, 建议对灯丝半圆柱面进行网格划分的数量不小于4;

3) 同一阵列中不同环向位置灯管相同顺序的微元产生的照射范围大小相等, 各微元照射范围的变化随微元左右位置呈现对称性;

4) 灯丝外加热时总体呈现两端微元照射范围小、中间微元照射范围大, 内加热则与之相反; 随着n的增大, 灯丝半圆柱面上各微元的照射范围变化趋势与微元顺序相关, 中间微元照射范围逐渐较小, 两端微元照射范围在外加热时逐渐减小、在内加热时逐渐增大;

5) 单根灯管外加热和内加热的解析法计算结果与有限元软件一致性良好, 验证了本文提出的辐射热流解析计算方法的正确性, 可以有效缩短大规模阵列的热流场计算时间, 支撑导弹等直径舱段的热试验快速设计迭代。

References

All Tables

All Figures

	图1 等直径舱段表面热流密度解析计算方法研究思路
In the text

	图2 圆柱外加热阵列单根灯管多网格照射范围计算模型
In the text

	图3 圆柱内加热阵列单根灯管多网格照射范围计算模型
In the text

	图4 不同网格数量下热流分布云图(圆柱外加热阵列, 单根灯管)
In the text

	图5 不同网格数量下热流分布曲线(圆柱外加热阵列, 单根灯管)
In the text

	图6 不同网格数量下热流分布云图(圆柱内加热阵列, 单根灯管)
In the text

	图7 不同网格数量下热流分布曲线(圆柱内加热阵列, 单根灯管)
In the text

	图8 有限元计算模型(外加热)
In the text

	图9 热流分布云图(外加热)
In the text

	图10 热流计算结果对比(外加热)
In the text

	图11 热流计算结果对比(内加热)
In the text