图1给出了圆管的物理模型。圆管外径为3 mm,内径为2 mm。进口绝热段和出口绝热段长度均为150 mm,用以保证进口流动充分发展和避免出口效应影响; 中间加热段长度为500 mm,外壁面给定均匀热流密度。进口为质量流速边界,给定质量流速和进口温度,出口为静压边界。进口和出口的圆环壁面设定为绝热边界。通过温度和热流密度相等来实现固壁和流体的耦合机制。重力加速度a施加在y方向,范围为0～6g,g=-9.8 m/s²。取管顶部(Top)和管底部(Bottom)做换热分析。

图1 水平圆管
Fig.1 Horizontal circular tube

连续性方程为

(1)

动量方程为

(2)

能量方程为

(3)

式中:ρ为密度; u为流速; c_p为定压比热容; T为温度; m_e为有效动力黏度; K_e为有效热导率; δ_ij为克罗内克符号。

选取RNG k-ε湍流模型结合增强壁面处理解决湍流换热问题。湍流模型如下

(5)

式中:k为湍动能; ε为耗散率; a_k和a_ε为湍流普朗特数; 常数项C_1ε=1.42,C_2ε=1.68,C_3ε=0.084 5; G_k和G_b均为湍流产生项; R_ε为附加项。

热传导方程

(6)

不锈钢圆管的热导率K为20 W/(m·K)。

正癸烷的临界压力p_cr和临界温度T_cr分别为2.1 MPa和618 K^[21]。超临界压力下其密度、定压比热容、热导率和动力黏度^[21]采用分段线性函数加到Fluent 19。图2给出了密度和定压比热容随温度的变化情况。

图2 超临界条件下热物性随温度的变化情况
Fig.2 Thermo-physical properties variations with temperature under supercritical condition

通过结构化“O型”网格划分计算域。流体域近壁网格需要加密,第1层网格满足y⁺<1,前3层网格满足y⁺≤5。管壁网格尺寸由外部向内部沿径向逐渐减小,共15层。表1制定5种方案,即改变截面网格数量和轴向网格数量考察网格无关性。进口温度为600 K,质量流速为1 200 kg/(m²·s),壁面热流密度为600 kW/m²,压力为3 MPa,重力加速度为2g。研究发现,网格方案为3 200×800时,截面网格数量和轴向网格数量继续增加,对出口温度T_out和出口流速u_out基本没有影响,满足网格无关性要求。管截面网格划分见图3。基于Fluent 19求解控制方程,通过二阶迎风差分格式离散控制方程,压力和流速耦合通过SIMPLEC算法处理,隐式Gauss-Seidel迭代处理。

表1 网格无关性分析
Tab.1 Grid-independence analysis

图3 管截面网格
Fig.3 Mesh configuration in the tube cross section

选择文献[11]常规重力下的实验方案验证湍流模型。水平圆管的管外径为2.20 mm,管内径为1.86 mm,管加热长度为500 mm。RP-3航空煤油进口温度T_in为573 K,压力p为3 MPa,质量流速G为538 kg/(m²·s),管壁热流密度q为320 kW/m²。图4给出管内壁温度沿管轴向(l为局部加热长度)的变化情况。数值结果与实验结果的相对偏差基本落在±5%的范围,说明湍流模型有效且合理。超重力条件下尚无实验数据,增大重力加速度基本不影响湍流模型的有效性^[19]。

图4 数值模型验证
Fig.4 Numerical model validation

图5给出了管内壁温度、管内壁热流密度和换热系数沿管轴向的分布情况。进口温度设定为600 K,质量流速为1 200 kg/(m²·s),外壁热流密度为600 kW/m²。重力加速度范围为0～6g。由图可以看到,运行压力为3 MPa且重力加速度为0g时,从l=0.075 m位置开始出现管壁温度陡增和换热系数剧降,即突出的传热恶化问题。过程中热流密度基本恒定。因为无浮升力作用,管顶部和管底部的换热参数相同。重力加速度为2g时,管周向出现壁温的非均匀分布问题,管顶部壁温高于管底部,这源于超重力产生的浮升力作用。随着重力加速度增加,管壁温度下降,超重力改善了热防护效果; 管顶部和管底部的壁温差增大,说明浮升力作用加剧。管底部热流密度高于管顶部,分别出现了热流密度峰值和谷值,说明浮升力致使周向热流重新分配,更多热流由管顶部传递到管底部吸收利用。重力加速度增加,浮升力增强,热流密度周向差别增大。管壁温度和热流密度的分布状况共同导致了管顶部先恶化后强化过渡到管底部只强化的换热特征,即周向不均匀的换热问题。重力加速度为6g时,运行压力提高,管壁温度、热流密度、换热系数的异常变化减弱,不仅传热恶化受到抑制,周向不均匀换热也减弱。

图5 换热参数的轴向分布情况
Fig.5 Axial distribution of heat transfer parameters

图6给出了管截面温度和无量纲质量流速的分布情况。无量纲质量流速定义为π_G=ρu/G。管截面选取l=0.2 m位置。由图可以看到,运行压力为3 MPa且重力加速度为0g时,流体域和固体域温度等值线均是规则的圆形。质量流速也呈周向均匀分布的特征,近壁面质量流速相对主流质量流速较小,即紧贴壁面的局部流量较小,出现流动性能较差(π_G<0.6)的流体层,称为拟膜态层。拟膜态层阻隔了热流的有效吸收,引起传热恶化现象。沿管轴向,首先近壁燃料密度急剧减小,拟膜态层不断增厚,传热恶化加剧; 随后近壁流速迅速上升,拟膜态层逐渐减薄,换热恢复。

图6 温度和无量纲质量流速的分布情况
Fig.6 Distribution of temperature and dimensionless mass flow rate

通过图6还可以观察到,重力加速度为2g时,固体域和流体域温度沿重力方向均出现了非规则圆环的异常分层问题。固体域温度等值线向下偏移,说明固壁热传导受到影响,出现周向非均匀的热传导过程,热流密度由管顶部转移到管底部。流体域温度中心向下偏移,温度等值线呈心形。流体温度异常分层导致密度的异常分层,产生周向密度梯度,出现周向不平衡动能。拟膜态层沿管周向出现非均匀现象,管顶部较厚,沿圆周方向逐渐减薄。重力加速度提高减小了拟膜态层厚度; 而流体温度异常分层减弱,整体趋于低温区,密度对温度的变化更加敏感,浮升力作用增强。提高运行压力,流体温度和质量流速等值线逐渐过渡为规则圆环,拟膜态层弱化,浮升力作用和拟膜态机制减弱。

图7给出了二次流动能的分布情况。二次流动能定义为k_sec=(u_x)²+(u_y)²。由图可以看到,相比无重力工况,周向不平衡动能驱动下管截面产生二次流问题。管中心冷流体向下流动冷却底壁面,对换热有强化作用; 热流体在管顶部形成驻点,流动性能最差,导致换热减弱。二次流动能最大值位于两侧接近壁面,说明该位置周向密度差别最大,对二次流的驱动作用最强。二次流对周向非均匀换热和固体温度异常分层发挥了作用。重力加速度提高使二次流动能增大,运行压力增大则使二次流动能受到严重削弱。

图7 二次流动能的分布情况
Fig.7 Distribution of secondary flow kinetic energy

图8给出了二次流涡量Jⁿ沿流动方向的变化情况。二次流涡量J^n[22]的定义为

(7)

式中A为流通面积。

图8的结果表明,二次流涡量沿流动方向先上升后下降,其峰值与管顶部和管底部最大壁温差值出现的位置一致。随着重力加速度提高,浮升力作用增强,二次流涡量增大,二次流增强; 同时,二次流涡量峰值向进口方向前移,说明重力加速度变化对壁温陡增阶段二次流的影响更为显著。运行压力提高,二次流涡量显著减小,二次流减弱,其峰值移向出口方向。

图8 二次流涡量的轴向分布情况
Fig.8 Axial distribution of secondary flow vorticity

图9给出了不同热质比(q/G)和重力加速度下管顶部和管底部最大内壁温度差随压力的变化情况。由图可以看到,运行压力提高造成最大壁温差逐渐减小,而重力加速度增加和热质比增大均使周向壁温差扩大。

图9 最大内壁温差随压力的变化情况
Fig.9 Variation of maximum inner-wall temperature difference with pressure

基于数值结果,建立了最大内壁温差的预测准则为

ΔT_wi,max=(238.25(q/G)^2.12(a/g)^0.49)/((p/p_cr)^2.40)(8)

图 10给出了最大壁温差数值结果与预测数据的比较情况。可以发现,两者误差满足±15%的要求,预测公式具有较高的精度。其适用的运行范围:0.33 kJ/kg≤q/G≤0.50 kJ/kg、3 MPa≤p≤6 MPa、2g≤a≤6g。

图 10 最大内壁温差预测数据与数值结果对比情况
Fig.10 Comparison between prediction data and numerical results of maximum inner-wall temperature difference