等截面圆柱形肋导热数值模拟研究

2020-02-01刘天昀刘文波

现代计算机 2020年35期

刘天昀，刘文波

（1.华北电力大学能源动力与机械工程学院，保定071003；2.广东生态工程职业学院信息工程系，广州510520）

0 引言

肋片作为二次传热面，扩大传热面积，并促进流体介质扰动来强化换热[1]。在各大领域中，肋片得到广泛的应用。无论是在空调的冷凝器，还是在汽车发动机的散热装置、计算机CPU的散热器中，都能看到它的身影。本文通过数值模拟的方法分别研究了等截面圆柱形肋片在稳态和非稳态情况下，肋片内部温度分布，通过肋片的传热量以及肋片效率。

1 基于MATLAB稳态导热问题数值模拟

1.1 计算区域的离散

本课题所研究的等截面圆柱形肋高为0.06m，横截面直径为0.01m。材料的比热为900J/（kg·K），密度为2700kg/m3，导热系数大小为50W/（m·K）。其中，肋根温度为100℃，周围空气温度为20℃，不考虑肋片与周围可见物体的辐射传热。模型如图1所示。

根据等截面圆肋特点，建立沿径向以及沿轴向二维模型。对该区域离散时，在r方向共有M个节点，空间步长取dr；在长度方向有N个节点，空间步长取dz。

1.2 建立节点方程

本课题根据控制容积热平衡法分别对底面节点、肋端面节点、内部节点以及边界节点建立离散方程。

图1 等截面圆柱形肋模型

沿中心轴节点：

肋端节点：

圆周与肋端交界处节点：

内部节点：

肋根节点：

圆周处节点：

1.3 MATLAB计算结果分析

（1）网格无关性验证

表1 节点数与散热量的关系

通过上表可见，当沿高度方向节点数20以及沿半径方向节点数为20，相对误差为0.022%，因此本次课题将采用高度方向节点数20以及半径方向节点数20进行相关计算。

（2）温度场分析

输入相关节点数后，可得到等面积圆柱肋片温度场。如图所示，对温度场分析可知肋片温度主要沿肋高方向变化，沿半径温度几乎不变，肋端温度为86.48℃。

图2 圆柱体温度分布

图3 肋端温度分布

（3）肋片散热量分析

肋片散热量是肋片边界与周围流体对流换热散热量。根据肋片特点，将肋片分为三部分进行计算散热量。第一部分，肋片圆周但不含与肋端交界边与周围流体对流散热量；第二部分，肋端但不含与肋片圆周交界处与周围流体对流散热量；第三部分，肋片圆周处与肋端交界处离散区域与周围流体对流散热量[1]。根据分块计算思想，相应程序如图4所示，其中Q1代表第二部分散热量，Q2代表第一部分散热量，q（M,N）代表第三部分散热量，计算结果为6.6893W。

图4 散热量计算程序

（4）肋片效率

在基础的散热表面上敷设肋片，主要是想通过增加面积来增大对流传热的热流量。但是肋片表面的温度是随肋片的高度变化的，所以肋片的形状和尺寸（特别是高度）对增强换热的效果有重要的影响。在计算肋片散热量时，为简化计算提出肋片效率ηf一概念，其定义式为：

ηf=（肋片实际散热量）/（假设整个肋片处于肋根温度下的散热量）

根据数值求解散热量，ηf=85.21%，本课题采用数值求解方法采用二维模型，而理论求解采用一维模型，这是产生误差的主要原因[2]。

2 基于MATLAB非稳态导热问题数值模拟

本课题求解非稳态导热问题，计算区域离散方法以及建立节点方程方法与稳态时相似，在此不再赘述，只给出相应结果。

2.1 建立节点方程

本课题采用显示格式建立节点方程。

沿中心轴节点：

肋端节点：

圆周与肋端交界处节点：

内部节点：

肋根节点：

圆周处节点：

2.2 肋片温度

本课题采用显示格式进行相关计算，通过初始化，将所有节点温度都设为20℃。该计算中，实际物理所经过的时间等于时间步长数乘以时间步长，同时为了避免迭代结果离散情况发生，时间步长需满足稳定性条件，本课题根据肋片温度以及几何条件，设定最大时间步长为0.004s，迭代次数10000步。

通过对迭代次数为1000、4000、8000以及10000温度分布图分析，肋片温度变化特点：温度变化位置从肋根开始逐渐至肋端处；随着时间延长，肋片温度逐渐不明显。

图5 1000步迭代

图6 4000步迭代

图7 8000步迭代

图8 10000步迭代

从迭代时间步长以及实际物理经过时间分析，显示格式计算非稳态问题时需考虑稳定性条件，在某些情况下出现时间步长过短，迭代次数过长的情况发生。因此，隐式格式计算非稳态情况适用范围更广。

2.3 肋片散热量

非稳态肋片散热量计算方法与稳态情况相似。肋片散热量是肋片边界与周围流体对流换热散热量。其散热量变化如表2所示，散热量随时间变化逐渐层增加，但由于肋片整体温度逐渐增加，与周围流体换热不断加强，散热量随时间变化斜率逐渐增大。

表2 散热量随时间变化

由于计算机性能局限，无法得到肋片达到稳态时变化情况，因此对散热量随时间变化做线性预测，如图9所示，通过幂函数预测曲线，当稳态时散热量6.6893W，所需时间约为85.41s，迭代次数约为21步。

散热量与时间关系预测曲线：

图9 散热量随时间变化线性曲线

3 基于商业软件Fluent数值模拟

3.1 几何模型建立

本课题所研究的等截面圆柱形肋高为0.06m，横截面直径为0.01m。

3.2 网格划分

在确定网格数量时，既需要满足计算精度要求，又要避免计算机运算负担过重[3]。以不同的网格划分数量为变量，以端部温度为对比参数进行网格无关性分析。

3.3 网格无关性验证

通过表3可见当网格规模于30万以上时，肋端温度相对误差小于1%，所以取30万以上规模的网格进行数值计算。

表3 网格总数与肋端温度的关系

3.4 物理模型构建

本次课题中，使用Gambit 2.4.6建模.使用直接生成体的方法绘制出高为0.06m，横截面直径为0.01m的圆柱形。在本例中利用体生成的方式直接生成网格，网格类型设置为map，选取合适的interval size，网格总数要30万以上。网格划分完毕后，需要对边界条件进行设定。在本例中，所有壁面设置为wall。但是需要注意的是，肋端、肋根以及侧面的边界条件需要分开定义。建立好后的几何模型如图10所示。

图10 导入Fluent后的网格文件

3.5 数值模拟分析

在Fluent中设置好边界条件后，分别对稳态和非稳态条件下等截面圆柱形肋片导热问题进行迭代计算，然后得到残差曲线图、温度分布云图、等温线图等相关数据。

稳态状况下时，通过做圆柱体纵向切片观察温度分布云图，发现在距离肋根0.04m以内的范围内基本只沿肋高方向变化，径向的温度梯度较小。而在端部附近的区域，径向的温度梯度不可忽略，如图11所示（右边为肋根，左边为肋端）。通过画plot图，可以发现随着与肋根距离的增加，温度不断下降，温度梯度不断减小，最终在肋端处，温度达到最小值359.3K即86.15℃，如图12所示。通过report可以读取通过肋片的热流量为6.668 W。

图11 圆柱形肋温度分布云图

图12 沿肋高方向的温度曲线图

在非稳态状况下，改变General选项中Time的设置，选中Time中的Transient选项，选择合适的时间步长，初始温度设为20℃，其余边界条件保持不变。使用同样的方式做不同时刻下的纵向切片的温度分布云图，观察温度场随时间的变化情况，如图13所示（右边为肋根，左边为肋端）。

通过做不同时刻下的纵向切片的温度分布云图，观察到肋片整体温度随时间而升高。到一定时刻后，温度场不再随时间变化，沿着肋高方向温度不断下降，温度梯度不断减小。

不同时刻下，还可以通过report读取对应时刻通过肋片的热流量。通过读取可得随着时间的增加，通过肋片的热流量不断减小，最终达到稳态时热量为6.668W，达到稳态时间所需时间约为116s左右。

图13 不同时刻下温度分布云图

4 Fluent数值模拟与MATLAB模拟对比分析

本课题分别使用MATLAB以及Fluent等软件对等截面圆柱肋片进行数值模拟工作，为使模拟工作更加精确，对两种方式得到相关结论进行比较并分析误差产生原因。

4.1 稳态情况

温度场分析，如表4所示，通过对比两种方法得出肋端温度以及散热量，两者相差较小。稳态情况下，两者模拟结论差别较小。这主要是因为：两者都是数值解法，基本思想都是把原来在空间与时间域内连续的物理量场用有限个离散点上的值的集合来代替[4]，数值解法的精确度与网格的密度及数量密切相关，网格划分越小，结果越精确。在网格密度都达到一定标准时，使用编程和商业软件进行数值模拟得到的结果基本一致，误差几乎可以忽略不计。

4.2 非稳态情况

非稳态时，通过对达到稳态情况所需时间分析，两种模拟方式差别较大，其中MATLAB模拟所需时间为85.61s而Fluent模拟所需时间为116s。

在非稳态时通过商业软件得到的结果与编程得到的略微存在一定偏差，这主要因为：与编程数值模拟不同的是，Fluent是基于有限元法的商业软件。有限元分析是用较简单的问题代替复杂问题后再求解，它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的近似解，然后推导求解这个域总的满足条件，从而得到问题的解[5]。与编程使用的有限差分法相比，Fluent使用的算法对于不同的算例的求解都有更好的收敛性、稳定性和精度。

使用基于有限差分法的MATLAB编译处理非稳态问题时，保持空间步长不变，增大时间步长，此时傅里叶数会增大，甚至不再满足节点离散方程的稳定性条件，因此其结果将进一步发散，与Fluent相比没那么精确。再者由于计算机性能限制以及程序算法的不足，使用MATLAB方法得到肋片达到稳态所需时间经过数据拟合得到，因此此种方式本就会带来一定误差。