单流体方法处理管道流固耦合传热

李一术，朱 波，刘家磊，房军龄

(1．海军工程大学 核能科学与工程系，湖北 武汉430033;2．中国人民解放军92730 部队，海南 三亚572000)

0 引 言

大型舰船为了满足动力需求，通常以核动力作为其动力输出来源。在特定的舱室中一般都包含大量管路，管道的主要作用是输送流体和热量，它的连接使各设备之间形成完整回路。在对舱室温度场的计算分析时，各管路与舱室的传热对舱室内温度场的分布有重要的影响。计算管道传热，通常需要进行流固耦合计算。目前有2 种流固耦合传热方法［1］:一种是强耦合方法，即整场离散整场求解方法;另一种是弱耦合方法，即分区求解边界耦合方法。对于计算模型较大时，一般采用弱耦合计算方法。通常舱室内的设备与管路众多，弱耦合方式则要求设定相应数量的流固交界面，提高了人员工作量。同时相比于普通的单流体流动计算，流固耦合计算需要更多的时间，对计算机的性能也提出了挑战。

针对管道流固耦合传热问题，Al －Nimr 和El －shaarawi［2］的解析方法解决了在2 个半无限大平行平板间或者在一个半无限大有限壁厚有限热容的管道中，不可压缩的热充分发展中的稳定一维的瞬时强制对流耦合热交换问题。利用解析法求解稳态条件下管道温度分布，夏文庆等［3］得出其计算结果与实验结果吻合得很好。张强等［1］以整场离散整场求解为基准，认为在数值模拟中松弛因子取0.75、收敛容差取0.1 可以提高计算精度和求解速度。

为减少计算量，化简数值模拟中的设置，尝试采用单流体方法处理流固耦合传热问题。结合流固耦合传热计算中解析法与数值模拟各自的特点，评价得出:流固耦合数值模拟能够作为管道传热计算的基准解。以此基准，验证单流体模拟计算结果。通过对模拟参数的改进，得到了单流体计算方式下较为理想的计算结果。说明单流体传热计算方法可行，可以应用到复杂的模型计算中。

1 流固耦合数值模拟的验证

1.1 物理模型

假设有一段长l = 3 m 的圆管道，管道内半径ri= 1.5 cm，外半径ro= 2 cm，管道材料的导热系数k2= 0.07 W/ (m·K)。管道内部流动的是空气，导热系数为k1，入口温度To= 20℃，速度V =5 m/s。环境温度Ti= 50℃，环境与圆管外壁的表面换热系数为ho。为了计算的简便，假定ho非常大，使得圆管外壁温度与环境温度趋于一致，这样省略了圆管外壁与环境之间的流固耦合传热计算。其物理模型如图1 所示。

图1 物理模型Fig.1 Physical model

1.2 流固耦合数值模拟

1.2.1 基本的物理守恒定律建立控制方程［4］

连续方程:

动量方程(Navier-Stokes 方程):

其他变量的输运方程:

式中:ρ 为流体密度;t 为时间;φ 为通用变量;Γ 为广义扩散系数;S 为广义源项。

1.2.2 网格划分

根据模型自身特点，利用Ansys－Icem 网格划分软件，对管道的流体域及固体域都采用结构网格划分。通过计算结果对网格敏感性分析，最终确定流体域网格数为3.1×105，固体域网格数为2.9 ×105。其网格划分如图2 所示。

图2 管道网格划分Fig.2 Mesh division of the tube

1.2.3 边界条件

本文所用的流体力学计算软件为Ansys-CFX，该软件具有先进的全隐式耦合多网格线性求解器，且收敛速度较快，具有友好的人机交互界面。为了进行模拟，需要设定计算参数与计算条件。入口边界条件采用速度入口，出口边界条件采用压力出口。流入的空气设为理想气体，设置空气流动模型为k - ε 模型。流体域与固体域网格采用GGI 连接模式对理想气体与管壁进行流固耦合传热计算。

1.2.4 模拟结果

利用可视化软件分析模拟结果，得到如图3 所示结果。图3 (a)为流体边界层温度场分布云图，图3 (b)为对管内流体主流温度随轴向坐标变化提取的数值分布。

图3 数值模拟计算结果Fig.3 Numerical simulation results

1.3 解析法对流固耦合数值模拟结果的验证

夏文庆等通过实验的方法证明，解析法在稳态条件下可以较好的预测实际结果，因此这里用解析法代替实验来验证流固耦合数值模拟的准确性。

解析法中的能量方程为［2］:

对于管内流体:

对于管壁:

式中:ρ1为管内流体密度;ρ2为管壁固体材料的密度;c1为流体定压比热容;c2为固体比热容;ri为管道内半径;ro为管道外半径;V 为流体轴向平均速度;T1为流体平均温度;T2为固体平均温度;Ta为管道外环境温度;hi为管道内壁对流换热系数;ho为管道外壁对流换热系数;t 为时间;qd为热流密度;z 为轴向坐标;k2为固体的导热系数。

将式(4)～式(6)的方程无量纲化，得:

式(7)～式(9)中各参数定义:

式(7)～式(9)受到以下初始条件的限制:

θi(0，Z)=0，i = 1，2，

θi(τ，0)=θoi= 1，i = 1，2，

q(0，Z)=0。

当Z →∞时，θi有限。

当管道内的流体传热处于稳态时，将无量纲方程转化为特征方程，解出其特征值。利用初始条件的限制，对特征值进行根的判别，最终得出流体的温度分布如下:

式中:

这里:

气流在管内的温度变化幅度在20℃～50℃之间，通过查看空气的热物理性质表，发现其物性参数变化很小，基本可以看作是常物性。ρ1=1.165 kg/m3，ρ2= 250 kg/m3，c1= 1 000 J/ (kg·K)，c2= 1 000 J/ (kg·K)，k1= 0.028 W/ (m·K)。在解析法的计算中，气体流动的未知参数为ho与hi。经计算，物理模型中的气流Re=8.8 ×104，属于过度流状态。因此，对于加热型的流体，使用最广泛的为关联式迪图斯－贝尔特公式:

分布函数做出的二维曲线如图4 所示，图4 中同时还加入了之前数值模拟的计算结果，以便两者进行比较。在管道流体初始发展阶段，2 种方法计算得出的流体温度有一定的偏差。沿着管道轴向发展，偏差逐渐变小，最终几乎趋于一致。流体从进入管口开始，需经历一段距离，管断面流速分布和流动状态才能达到定型，初始偏差较大的原因，可能是由于解析解在这个问题上考虑有所欠缺，所以导致偏差较大。就解析解作为基准解而言，数值模拟计算的流体整体温度偏差范围较小，因此，可以很大程度上信赖其计算结果。

图4 数值模拟结果与解析解对比Fig.4 The results of numerical simulation compared with analytical solution

2 单流体处理方法

分析管道的传热，当其长度远大于管道壁厚时，壁内温度可近似看成一维稳态温度场，所以圆管壁面导热公式为［5］:

式中:q 为热流密度;λ 为管道材料导热系数;r1为圆柱体内半径;r2为圆柱体外半径;tw1为圆柱体内壁温;tw2为圆柱体外壁温。

在对流体的传热计算中，CFX 允许流体域外围有一层假想固体域［6］。如图5 所示，内部物理流体边界至最外部物理边界之间有一段厚度的墙壁，但其没有在实际模型中进行模拟。通过模型边界的热流量用下式计算:

式中:To为指定的外部边界温度;Tw为流体域边界温度;qw为从外界流入流体域的热流密度;hc为流体域边界的表面传热系数，其不同于表面换热系数，因为它不涉及流体的温度。

图5 假想边界模型Fig.5 Imaginary boundary model

在处理管道内流固耦合传热计算时，单流体传热计算方法为:假想有一层特定厚度的管壁，实际省略对管壁的建模和流固热耦合的计算。在单流体模拟计算中，To为已知条件，Tw由迭代计算而得出，而hc是需要初始设定的未知参数。根据一维圆柱筒壁面传热的近似，可得hc为:

利用式(15)圆管道的具体尺寸，换算得到在单流体传热计算中需要设定的hc值。

3 单流体模拟计算结果及改进

流固耦合的数值模拟与解析法都能够得到相对准确的结果，但解析法应用的条件更为苛刻，这是因为只有在流体温度变化范围小时，才可以利用其特点取常物性参数。对于数值模拟来说，应用范围更广，而且利用计算机强大计算能力，减轻了人员的工作量。

为了验证单流体假设计算的合理性，做多组数值模拟方案进行分析。在相同物理条件下，参考流固耦合方式的数值模拟结果，验证单流体假设计算的可行性。选取多组不同尺寸管道，其尺寸及相应尺寸下计算所得hc如表1 所示。环境温度Ti提升为100℃，其他物理条件不改变。

表1 管道尺寸及计算所得hcTab.1 Tube size and calculation value of hc

对于加热型管道，流入流体域的热流量是影响管内流体温度的关键因素。为了方便多组计算的比较，直接选用流经管壁的热流量进行对比分析。经过数值模拟，分别使用流固耦合方法和单流体方法计算出的总热流量如图6 所示。可以看出，使用单流体方法得出的总热流量都普遍偏小。这可能是由于管道壁面导热近似处理存在误差所致。设管道外半径与管道内半径的比值为x，流固耦合下的总热流量与单流体得出的总热流量比值为y。仔细分析单流体方式所带来的管壁热流量偏差，可以发现:y 与x 之间存在一定的关系。其关系如图7 所示。

从图7 观察发现，两比值之间近似地服从线性关系，通过线性拟合的方式，得出2 个比值之间的数值关系如下:

图6 两种数值模拟方法所得流经壁面热流量Fig.6 The heat flux through the wall with two kinds of numerical simulation method

图7 x 与y 的关系Fig.7 Relationship of x and y

单流体计算的热流量由式(14)所决定，To为定值，要想改变qw，只能对hc进行修改。当调整hc时，Tw的值也会随着边界条件的改变而改变，那么qw的变化不能够精确掌控。经测试，hc在小范围变化时，Tw在现有物理条件下，相比原有温度变化幅度很小。因此小幅度的改动hc值，Tw的变化对于To与Tw之间的差值影响很小。这样，可以依照管道半径比值与热流量比值之间的关系来修改hc，以达到修正qw的目的。修改后的hc变为h′c，其表达式如下:

利用h′c重新计算各方案通过管壁的热流量，并再次与之前计算的流固耦合下热流量进行对比，得到如图8 所示的结果。可以看到，利用这种方法很好的解决了单流体传热计算带来的管壁热流量偏差，使其计算结果与基准结果更加接近。这也说明单流体方法处理管道流固耦合传热有效。

图8 改进后两种数值模拟所得流经壁面热流量Fig.8 The heat flux through the wall with two kinds of improved numerical simulation method

4 结 语

对于管道内流固耦合传热这类实际问题，解析法与流固耦合数值模拟法都能够得到较为准确的计算结果。解析法通常只能针对简单模型抽象处理，当实际问题较为复杂时，如求解舰船舱室内温度场分布问题，流固耦合数值模拟更有其实用价值。

在流固耦合数值模拟的基础上，改进模拟思路，提出了单流体数值模拟方法。经验证，其同样可行。应用单流体数值模拟方法，不仅令设置更加简单，而且减少了计算量，降低了对计算机硬件的要求。所以，在传热模型复杂时，单流体数值模拟方法更为有效。

［1］张强，刘巨保，赵晓荣.换热管流固耦合传热的分区求解数值算法［J］. 科学技术与工程，2010，10(28):1671-1815.

［2］Al-Nimr M A，El-Shaarawi M. Analytical solutions for the transient conjugated heat transfer in parallel plate an circular ducts［J］.Int Comm Heat Mass Transfer，1992，19:869 -878.

［3］夏文庆，明晓，朱春玲.基于状态空间的流固耦合传热的分析［J］.南京航空航天大学学报，2004，36(1).

［4］王福军. 计算流体动力学分析—CFD 软件原理与应用［M］.北京:清华大学出版社，2004.

［5］章熙民，任泽霈，梅飞鸣.传热学［M］.北京:中国建筑工业出版社，2007.

［6］ANSYS CFX - Pre User′ s Guide［M］. Pennsylvania:ANSYSY，2009.