基于离散坐标法的辐射模型研究

2022-12-02魏威张耘隆刘汉宇唐伟闫指江

强度与环境 2022年5期

魏威张耘隆刘汉宇唐伟闫指江

（北京宇航系统工程研究所，北京，100076）

0 引言

在工业燃烧炉、航空发动机燃烧室、火箭底部，辐射是一种非常重要的热传递方式[1-3]。在液体火箭底部附近，高温的发动机喷流通过辐射和对流加热火箭的底部防热板，从而可能造成火箭底部温度过高而产生问题，新一代液体运载火箭发动机喷流在底部防热板引起的热流主要为辐射热流。辐射换热可以显著地影响流场的温度，从而影响化学动力学和火焰结构，甚至导致火焰的熄火。此外，近些年来大量的研究工作致力于湍流和辐射的相互作用（TRI）。因此气体辐射的有效计算对于工程中设备的结构设计和校核具有十分重要的意义。气体辐射具有强烈的光谱选择性和容积性（或延程性），其吸收系数随着波长剧烈的变化，因此这成为气体辐射计算的难点，随之而来的有两个问题：第一气体辐射模型问题，第二气体辐射的计算量大为增加。因此在气体辐射中合理地建立辐射模型，对于气体辐射的非灰特性的模化和有效的降低计算量有着重要的意义。

现有的气体辐射的计算方法可以大致分为三类：（1）逐线计算法；（2）谱带模型法；（3）全局模型。逐线计算法是最精确的方法，但其计算量非常大，因此经常用作新模型的校正基准。谱带模型法包括窄谱带模型和宽谱带模型法。统计窄谱带模型（SNB）能够给出高温气体的辐射特性，但是其仅给出窄带内的穿透率，而很多辐射求解方法要求给出气体辐射的吸收系数，因此统计窄谱带模型很难和离散坐标法、有限体积法相结合使用。另外窄谱带模型需要大量的谱带，其计算量相当可观。宽谱带模型可以看做是窄谱带模型的简化，其计算量较窄谱带少，但它给出的是谱带吸收率，因此也不易和离散坐标法、有限体积法等结合使用。此外宽谱带模型法要求光程长度，因此其给出的辐射参数将和光程长度有关。相关k分布法（CK）[4-5]假设谱带足够窄，可以将吸收系数重新排列为一个光滑的单调变化的函数，然后用辐射传递方程对吸收系数的积分取代对波数的积分，因此CK模型更容易应用到辐射传递方程的各种解法中。CK模型需要大量的谱数据，计算量也很大。全局气体辐射模型中主要包括全谱相关k分布法（FSCK）和加权气体和模型（WSGGM）。FSCK模型采用Planck加权k分布函数[6]，在整个谱内采用CK方法，因此相比CK模型减小了大量的计算量。F. Liu[7]等发展了基于Malkmus统计窄谱带模型的FSCK，或称为SNB-FSCK。由于Malkmus统计窄谱带模型的k分布函数和累计k分布函数具有分析表达式，这使得SNB-FSCK计算上更简便。尽管如此，如果和CFD结合，SNB-FSCK方法的计算量还是偏大，不利于工程中应用。加权灰气体和模型使用与温度无关的吸收系数和温度相关的加权因子来拟合气体的发射率，能够通过4-5种灰气体精确地辐射的非灰特性，非常适合与CFD的结合。本文基于离散坐标法求解辐射传递方程，通过和SNB计算结果的对比，比较了几种WSGGM模型的精度。此外还比较了不同离散坐标方式对计算精度的影响，以及计算量的大小。

1 控制方程及数值方法

1.1 控制方程

本文采用离散坐标法对辐射传递方程进行求解，气体辐射模型为加权灰气体和模型，不考虑气体辐射散射，在加权灰气体和模型下辐射传递方程为：

边界条件为：

辐射热源项可以表达为：

壁面热流可通过下式进行计算

加权灰气体和模型是假设几种灰气体具有不依赖于温度的吸收系数和依赖于温度的加权因子式(6)，模型参数通过拟合气体的发射率式(5)得到

1.2 数值求解方法

对于离散坐标体系下的辐射传递方程，其空间离散格式有阶梯格式、指数格式、菱形格式，菱形格式虽然精度高，但是随着网格的加密会出现负的辐射强度值，导致物理上的不真实。指数格式的计算量较菱形格式和阶梯格式大。阶梯格式始终能够给出物理上真实的解，因此本文中选用阶梯格式进行辐射传递方程的空间离散。离散后的方程采用CGSTAB方法对所有节点直接求解。

1.3 计算模型

本文采用Liu[8]采用SNB辐射模型计算的模型算例如图1所示，三维的计算域为2m×2m×4m，x方向长度为2m，y方向长度为2m，z方向的长度为4m，四个壁面为黑体冷壁面，温度为300K。计算域中包含的混合气体体积分数组成为10% CO2, 20% H2O 和70% N2，气体的压力为1atm。模型中气体的温度是不均匀的，计算域中温度分布为T=(Tc-Te)f(r/R)+Te，其中Tc为计算域中沿z方向的中心线上温度分布，Te为z=4m处的温度，中心线上的温度从z=0m处的400K线性增加到z=0.375 m处的1800K，然后线性减小到z=4m处的800K。长方体区域中有一个R=1m的同心圆柱形区域，圆柱形区域内的温度由函数f(r/R)=1-3(r/R)2+2(r/R)3描述，其中r为响应点到中心线的距离，在这个区域之外气体的温度为800K。本文中网格划分为41×41×80，分别为x、y、z方向的网格数。

图1 辐射计算域Fig.1 Radial region for calculation

2 结果及讨论

2.1 不同辐射模型的结果比较

本文比较了三种WSGGM模型，分别为Smith(4GG)[9]，Perry(5GG)[10]和Truelove(4GG)[11]，坐标离散方式采用S4离散，并和文献[5]的统计窄谱带模型法（SNB）方法进行了比较，SNB模型在计算量上很大，精度非常高，因此常用作其它模型的基准。图2（a）和（b）分别给出了中心线（0m，0m，z）上的辐射热源和壁面（1m，1m，z）上的热流分布。可以看到不同的辐射模型无论对辐射热源，还是对壁面热流结果都有非常大的影响。

图2 三种WSGGM模型计算的结果Fig.2 Results of three WSGGM models

从图 2（a）可以看到 Smith(4GG)和Truelove(4GG)计算的中心点热源在最低点和SNB方法的计算结果比较接近，但在z大于0.375 m的大区域内都低估了辐射热源。Perry(5GG)在中心热源的最低点误差较大，但在最低点之外的大部分区域和SNB的计算结果吻合的很好。从图2(b)可以看到Smith(4GG)在壁面热流最高点附近较大的区域内高估了辐射热流。Truelove(4GG)在整个范围内都低估了辐射热流。Perry(5GG)在除热流最高点附近的大部分区域略低于SNB的结果，但整体上Perry(5GG)的结果较好。

2.2 不同坐标离散方式对结果的影响

在离散坐标法中，非常重要的是离散方向的选取，不同的离散方向和权重，产生了多种离散坐标方式，这里对比了六种坐标离散方式对此算例结果的影响，辐射模型为Perry(5GG)。图3（a）和（b）分别给出了中心线（0m，0m，z）上的辐射热源和壁面（1m，1m，z）上的热流分布。

从结果可以看到不同的坐标离散方式对辐射热源影响非常小，但是对辐射热流的计算产生了非常大的影响。从图3（b）可以看到S2对称和非对称的离散方式都低估了辐射热流，并且误差非常大。层对称奇阶格式的S6 LSO和S8 LSO以及层对称偶阶格式S4 LSE的计算结果和SNB的结果吻合的比较好，层对称奇阶格式S4 LSO在热流最高点附近高估了辐射热流，并且具有较大误差。

图3 不同坐标离散方式对结果的影响Fig.3 Effect of different coordinate discrete methods

2.3 不同坐标离散方式对计算效率的影响

上一小节研究了不同的坐标离散方式对加权灰气体和模型计算精度的影响，本文同时考察了不同的离散方式的计算效率。表1中给出了在采用辐射模型Perry(5GG)时不同的坐标离散方式的时间消耗，可以看出随着离散方向的增加，辐射计算的时间消耗逐渐增加，相同离散方向数的两种S2以及两种S4的时间消耗基本相同。在辐射模型确定时，S4离散的计算时间是S2离散计算时间的3.5倍左右，S6计算时间是S2的10倍左右，S8离散计算时间是S2离散的23倍左右。表1的结果表明，离散方式对辐射计算时间的影响主要和离散方向数有关。这是由于当离散方向数目增加时，求解方程的数目也相应增加，求解的时间消耗也相应增加。在流场与辐射耦合计算时，可以采用相对低精度的坐标离散方式来显著的减小辐射计算部分的计算量，从而加速辐射计算速度。