基于对流传热反问题的矩形通道内扰流横肋的形状优化

陈金月，王 坤，闵春华，杨旭光

（河北工业大学 能源与环境工程学院，天津 300401）

0 引言

通道表面增设扰流横肋是无源强化换热的一种常用手段。因其具有操作简单、强化效果较好等优点，在电子器件散热[1]、冷却设备[2]及太阳能空气加热器[3]等领域均得到了广泛应用。扰流横肋的形状和排布方式是影响强化换热的重要因素。已有很多学者对此开展了相关研究工作。表1给出了文献中出现的几种不同形状的扰流横肋结构。

表1 文献中典型横肋扰流横肋汇总Tab.1 Summary of typical structures of transverse ribs

Vanaki 等[9]对矩形、三角形、楔形等不同扰流横肋形状进行了优化和性能对比，结果表明三角形肋能够实现最大综合强化换热效果；Yadav等[10]对太阳能空气加热器内的流动与换热进行了数值分析，研究了湍流状态下的12种不同尺寸的等边三角形剖肋构型，结果表明肋片存在最优的螺距与高度。除肋片形状外，排布方式对加肋通道的性能影响也很大，以下是研究者对不同排布方式的肋片的研究。Vanaki 等[9]对湍流状态下，上底面加热的2 种布置方式（直排式和交错式）进行数值模拟对比，结果表明直排式的Nu值最高，综合强化换热效果最好；Liu 等[11]重点研究了不同截断方式肋片流道的传热及阻力特性，对比研究发现，截断型肋片在不降低换热效果的基础上，可以大大减小流道的压降；Abraham等[12]对湍流状态下V型和W型肋的中心顶点分别指向上游和下游的4种布置方式进行了实验和模拟研究，最终得出建议使用W型排布肋。

以上研究针对不同形状和排布的扰流横肋的性能进行了对比分析，未对横肋的形状进行优化。近年来反演优化的思想越来越受到大家的重视。早期的研究主要关注于导热反问题[13-15]或是反演未知的热边界条件[16-18]，而将其应用于对流换热几何参数反演优化的还较少。在本研究团队之前的研究工作中[19-20]，对于布置方形肋片的二维流道进行了对流换热反问题的研究，分别以流道出口温度（Tout定壁温边界条件）和火积耗散（Eh定热流边界条件）作为目标函数，反演优化了流道内肋片的最佳节距值。然而，这些研究工作仅对传热性能进行了优化，并未考虑阻力特性。本文基于共轭梯度法，以综合考虑传热与阻力特性的综合换热性能作为目标函数，对不同排布方式下矩形通道内扰流横肋的形状进行反演优化。

1 模型及计算方法

1.1 计算模型

图1 a）是三维矩形通道，图1 b）是简化的二维模型。主流道几何尺寸为：高H=20 mm，长L=12.5H，宽W=2H。计算过程中为保证进口流体速度稳定及出口无回流，流道的实际计算长度分别向两端延伸了l，其中l=5H。初始扰流横肋结构为等腰直角三角形横肋h=1/4H，w=1/2H，第一个三角形横肋距流道入口端的距离为11.25H，S为扰流横肋顶点之间的距离，扰流横肋个数用N表示，图中N=4。

图1 计算模型示意图Fig.1 Schematic diagram of model

流道中的流动介质为空气，假设空气的物性参数为常数，通道左侧为入口，流体以恒定流速流入，给定雷诺数为8 000，进口温度Tin=300 K，右侧为出口，出口截面为压力出口。主流道底面为第二类边界条件，q=5 000 W/m2，其他边界采用绝热边界条件，忽略该系统的辐射换热。

1.2 计算方法

压力和速度的耦合项采用SIMPLE 算法，无滑移边界条件，对有限体积法中控制体进行离散化时，对流项采用二阶迎风格式离散，扩散项采用中心差分格式离散，湍流模型选择standardk-ε。正问题的收敛准则是当连续方程、动量方程的残差均小于10-6，能量方程的残差小于10-7时。模型的处理、网格的划分均采用ANSYS-FLUENT软件实现。

流动与换热过程的控制方程如下。

连续性方程为

式中，ρ为流体的密度。

动量守恒方程为

式中，

湍流动能方程为

耗散率方程为

式中：S为应变张量的平均速率之模，其定义式为S=(2SijSij)1/2，Sij=1/2(∂uj/∂xi+∂uj/∂xi)；αk和αε分别为k和ε的有效普朗特数（Prandtl number）的倒数，其值αk=1.0，αε=1.3；C1ε=1.42，C2ε=1.68；应变率，η=Sk/ε，η0=4.38，β=0.012。

能量守恒方程为

式中，λ是导热系数。

1.3 模型与网格验证

对于计算区域进行网格独立性验证时，矩形通道采用结构化网格，扰流元左右20 mm附近采用非结构性网格加密，网格数大于109 942后，流道的Nu及f均不发生大的变化，其结果如图2所示。因此，选择此时的网格数进行数值模拟计算。

图2 网格无关性考核Fig.2 Grid independence test

在进行优化之前，为了验证模型的准确性，将扰流横肋周围的局部努塞尔数数值结果与Liou等[21]的实验结果进行了比较，其结果如图3所示。可以看出，模拟结果与实验结果具有很好的吻合度，验证了数值模型的可靠性。

图3 模型验证Fig.3 Model validation

2 反演优化方法

2.1 优化问题的目标函数

上述的扰流横肋几何边界反问题的目的是找到扰流元几何边界r（θ），使目标函数值最大，目标函数可以定义为

式中：a为系统传热能力与系统阻力的权重系数，其大小反映了优化过程中系统换热能力及减阻效果的侧重程度；Nu0和f0分别代表光滑通道的努塞尔数和摩擦系数，由Gnielinski和Filonenko公式定义：

2.2 扰流横肋边界几何形状的拟合

本文采用三次Spline函数将离散的数据点连成光滑曲线，与其他方法相比，此方法操作性强，光滑性好。设有m个边界点(θi,ri)(i=1,2,…)，ri(i=1,2,3,…)为极半径，ri值可变；θi为极角，(0＜i≤m)，下角标i表示沿顺时针方向的离散点序号。由插值函数的边界条件可得

则可以唯一确定一个φΔ(θ)表示成如下的分段形式：

式中，Sk(θ)(1 ≤k≤m)是1个三次多项式，且满足插值条件

由式（13）可得

其系数的确定过程为：

对式（14）两边同时求导得，

利用αk(θ)和βk(θ)的性质，式（15）等价于

式中，hk=θk+1-θk。

由边界条件公式（11）可导出2个方程：

由式（17）和式（18）得矩阵形式的线性方程组为

由此可得到m0,m1,···,mm-1,mm。

2.3 简化共轭梯度法求解反问题的计算流程

本文所用共轭梯度法构造一组搜索方向，求出目标函数的极值点。其求解的计算流程如下。

每个元素ri的更新方式为

式中：βi为搜索步长，取搜索步长为定值，其值为0.0002mm[19]；di为搜索方向，其计算式（21）为

反问题求解的计算过程如下：

1）采用数值模拟的方法，求解整个计算区域的速度和温度分布，然后求得目标函数值；

2）根据正问题的结果和收敛准则ε，（ε是很小的正数），此时令k=1；

3）给研究变量一个小的扰动Δr，根据数值计算结果，按照式（8）和式（23）计算对应的目标函数J[r(θ)]及其梯度，判断梯度值是否满足收敛准则，若≤ε，则迭代终止，否则继续；

4）按照式（22）和式（21）分别计算共轭系数γi和搜索方向di；

5）根据式（20）计算更新ri的值，并令n=n+1，返回步骤3。

3 结果与讨论

3.1 不同目标函数的对比

以具有单个扰流横肋的通道为例，初始扰流横肋形状为图1所示结构，此次迭代为未拟合边界所得结果，离散边界点为11，其他边界条件不变。Fan等[22]曾提出一种传热综合效果评价图，给出了等流量约束条件下、等压降约束条件下、等泵功约束条件下分别为a=1、1/2、1/3。根据式（8），目标函数分别为。

图4和图5分别为不同a值反演优化后所对应的优化结果。由图可知，随着a值逐渐增大，即随着系统阻力权重的逐渐增大，优化后扰流横肋结构的尖角逐渐平缓甚至消失。这是由于尖角的存在可以提高传热性能，但同时也会增大流动阻力。当a=1/3时，目标函数值达到1.869，但增量最小，仅为3.32%。可根据不同的工程要求来选取不同的a值。

图4 不同a 值下的优化形状Fig.4 Optimal shape with different a

图5 不同a 值下优化结果Fig.5 Optimization results under different a

3.2 不同排布方式的对比

以具有2 个扰流元的通道为例（N=2），初始扰流元形状为图1 所示结构，此时目标函数为，其他边界条件不变。图6为扰流元不同排布方式的结构示意图，之后分别对4 种布置方式下扰流元的间距S进行讨论，得出底排随着S增大J值一直在减小，最终取S=10 mm；叉排后置S增大的过程中，出现J值先增大后变小，当S=65 mm 时，此时J值最大；叉排前置随着S增大J值一直增大，当S=120 mm 时，此时扰流元已经位于加热壁面的边缘，因此取此时的位置对扰流元进行反演优化，以下均采用四种间距下的初始结构进行优化。

图6 不同排布方式的结构（N=2）Fig.6 Schematic diagram of different arrangement（N=2）

图7中4种扰流元布置方式，可以看出布置方式为底排的扰流元优化后所得的目标函数最大；顺排布置的扰流元优化后所得的目标函数与底排布置的扰流元优化结果相比稍小；叉排布置的2 种扰流元优化后所得的目标函数远小于底排布置的结构。

图7 不同排布方式优化性能随迭代步数的变化（N=2）Fig.7 Optimization of the performance-changing process with different different arrangement（N=2）

3.3 不同扰流横肋个数的对比

3.3.1 边界离散点数目的影响

以具有单个扰流元的通道为例，初始扰流元形状为图1 所示结构，此时目标函数为，其他边界条件不变。分别对离散边界点个数为m=3、7、11 的3 种情况进行反演。表2 给出了拟合前后取不同离散边界点优化结果的速度与压力分布。可以看出，拟合前后大致形状基本类似。

表2 拟合前后云图对比Tab.2 Contour plot before and after fitting

图8 给出了拟合后不同离散边界点所得的目标函数值，m=3时，J值达到1.866，此时的目标函数值远高于多个点所取得的值，分别比m=7 时增加0.865%，比m=11 时增加0.92%。因此，选用m=3的反演结果。对于拟合，并非所选点越多，拟合阶次越高越能达到我们的要求，传统的高阶曲线拟合方法，如拉格朗日法会在拟合曲线的两端出现大幅震荡的“龙格”现象，造成求解区域的严重失真，反演过程无法实现。

图8 不同点数下的优化结果Fig.8 Optimization results under different m

图9 为不同初始扰流横肋形状示意图，分别为：等腰直角三角形横肋、正五边形、半圆形结构，主要尺寸w，h与等腰直角三角形相同，3个初始形状在离散边界点为3时进行反演优化。图10为离散边界点为3时不同初始形状下的优化结果，拟合优化后形状基本类似，如半个水滴形，最高点像左偏移，文献[23]中曾模拟了扰流元柱的流动和换热，优化得到插排布置的扰流元柱为水滴形，与本文离散边界点为3时反演优化得到的结果类似。

图9 不同初始扰流元形状Fig.9 Different initial rib shapes

图10 不同初始形状下的优化结果Fig.10 Optimization results under different initial shapes

3.3.2 扰流横肋个数的影响

取图1所示的三角形为初始结构，边界离散点为3，目标函数为，其他边界条件不变。图11为不同扰流横肋个数拟合后的目标函数值的比较，由图11可知，不同扰流元个数拟合所得结果，N=1与N=3 所得目标函数结果相差不大，能达到1.866，N=1 时，减少2.84%的换热，降低了18.47%的阻力，得到2.89%的目标函数的提升；N=3时，增加了2.6%的换热，降低了2.12%的阻力，得到3.33%的目标函数的提升。

图11 不同扰流元个数下的优化结果Fig.11 Optimization results under different number of ribs

图12 为不同个数扰流横肋的速度等值线图，可以看出，N=1 时扰流元的形状为顶点向左偏移的半个水滴形；N=2 时，为2 个向内凹的三角形结构；N=3 时，前2 个为向内凹的三角形结构，第3 个结构为顶点较低且向左偏移的半个水滴形；N=4时，前2个同为向内凹的三角形结构，第3个为较低的三角形结构，第4 个为顶点较低且向左偏移的半个水滴形，N=5 时，此时目标函数值较小，不与考虑。这是因为N=3 时的结构是由N=1 和N=2 结构的组合，初始形状需对来流的扰动尽可能大，所以开始为内凹三角形，最后1个需要减小由于三角形的强扰动引起的阻力，所以结构会较平缓，最终会出现这种形状的多个扰流元的组合结构。

图12 不同扰流元个数下的速度云图Fig.12 Velocity contour with different numbers of rib transverse

由于N=1与N=3所得目标函数结果相差不大，所以我们观察在流动过程中，如图13所示，不同扰流元个数沿x轴局部Nux和fx变化，当气体流到扰流元之前，矩形通道内气体的局部阻力损失相同，Re一定时，流体刚进入通道时，由于入口段的边界层较薄，局部Nux数较高，随着流动的进行，入口效应逐渐减弱，Nux数快速降低，当流体流经扰流元结构时，在扰流元尾部，在压差的作用下，产生了二次流，加强了扰动和混合，Nux数有所增加，之后Nux会趋于稳定，N=1与N=3的Nux数基本相似，但沿着通道到达扰流元尾部时，N=3时的局部阻力值远小于N=1时的局部阻力，在二次流完全结束之后，两者的局部阻力值基本相同。因此，最终选用N=3时的优化结构。

图13 不同扰流横肋个数沿x 轴局部Nux 和fx 变化Fig.13 The number of different rib transverse varies locally along the X-axis in Nux and fx

4 结论

本文以综合换热性能最佳为目标，采用共轭梯度法对矩形通道内扰流横肋的最佳形状进行了反演优化，对优化目标函数的选取、扰流横肋布置方式及数目的影响进行了讨论分析。得到如下结论：

1）采用不同的综合换热性能评价因子作为优化目标，反演优化得到不同的最优肋片形状，在横肋的肋基面积给定的情况下，随着优化目标中阻力特性权重的增大，最优横肋形状的尖角逐渐平缓甚至消失；

2）扰流横肋底部顺列的布置方式优于上下错列的布置方式，能够获得更高的综合换热性能；

3）在本文研究工况下，扰流横肋数目为3时，能够取得最佳的强化换热效果，且沿着流动方向3个肋片的最优形状逐渐从向内微凹的三角形变成向外微凸的水滴状。