金字塔凹表面构成通道内超临界CO₂流动二次流特性研究

Characterization of Secondary Vortices in Supercritical CO₂ Flow in a Channel Formed by a Concave Surface of a Pyramid

摘要: 超临界二氧化碳因自身特殊的物性广泛应用于各种能源系统，在高效性、环保性、可持续性等方面均具有明显优势，与常规流体相比，其可以打破传统稳态流动的传热限制，使传热更加均匀和高效。而印刷电路板式换热器因较大的传热面积密度和易于模块化的特性常作为超临界二氧化碳布雷顿循环发电系统的典型回热器被使用，故对以超临界二氧化碳为工质的印刷电路板式换热器进行研究很有必要。于是本文通过一种冲压的手段在流道表面形成金字塔结构壁面来强化换热，即通过二次流强化传热技术有效提高各项性能指标。而通道中诱导的二次流又受限于金字塔单元结构参数的变化，故通过仿真的手段重点探讨了金字塔单元结构参数的不同排列组合对二次流的影响。此外，还简要分析了超临界二氧化碳物性与二次流的关系。

Abstract: Supercritical carbon dioxide is widely used in various energy systems due to its special physical properties, and has obvious advantages in terms of high efficiency, environmental protection, sustainability, etc. Compared with conventional fluids, it can break the heat transfer limitations of traditional steady state flow and make heat transfer more uniform and efficient. The printed circuit plate heat exchanger is often used as a typical return heat exchanger in supercritical carbon dioxide Brayton cycle power generation system due to its large heat transfer area density and easy modularization, so it is necessary to study the printed circuit plate heat exchanger with supercritical carbon dioxide as the work material. Thus, in this paper, the heat transfer is enhanced by forming a pyramid structure wall on the surface of the flow channel by means of stamping, i.e., the performance indexes are effectively improved by the secondary flow-enhanced heat transfer technology. The secondary flow induced in the channel is limited by the variation of the structural parameters of the pyramid unit, so the influence of different combinations of the structural parameters of the pyramid unit on the secondary flow is emphasized by means of simulation. In addition, the relationship between the physical properties of supercritical carbon dioxide and the secondary flow is also briefly analyzed.

1. 引言

近几年高效紧凑节能式换热器快速发展，在此背景下，印刷电路板式换热器(printed circuit heat exchanger，简称PCHE)因较高的传热能力和易于模块化的特性被世界所青睐，一经问世，便受到极大的关注。尤其因CO₂在临界温度附近的可压缩性较大，当将其应用在布雷顿循环发电系统中时会使压缩机的消耗功率显著降低，进而循环效率提高，这便是以S-CO₂为工质的PCHE用于动力循环发电系统的优势所在。Meshram [1]等人对在S-CO₂布雷顿循环中用作回热器的PCHE进行了分析，他们研究并比较了平直流道和Z型流道PCHE的性能。结果发现，在Z型流道PCHE中使用较大的弯曲角度和较小的节距会表现出更好的热工水力性能。Rao [2]等人开发了一种用于再压缩布雷顿循环的热经济性目标函数，该布雷顿循环具有带S形翅片的PCHE，以便其对多用途目标进行优化处理。Park [3]等人通过实验研究了直通道PCHE中CO₂的传热和流体流动特性，结果表明，使用入口和出口数据的平均值是不合适的，因为在临界点附近CO₂的性质发生了显著变化。Dai [4]等人在S-CO₂布雷顿循环的研究中，对以水为工质的锯齿型通道PCHE进行了实验分析，发现与层流状态下的直通道PCHE相比，锯齿型通道在非常低的雷诺数下会发生瞬态流的转变。在有关S-CO₂回热式布雷顿循环PCHE热工水力分析的文献中，多数研究集中在主流结构中的热工水力设计，而PCHE的结构设计与优化引起的关注较少，故评估新型结构的设计与优化改进具有重要的现实意义。

同时本文研究的二次流技术作为强化传热的主要机理也受到学者们的广泛探索。武永和[5]等在层流工况下对管内插不同涡发生器进行了数值研究，结果表明，管内二次流强度和平均努塞尔数的变化具有一致性。Lin [6]等优化了一种新的扭带涡发生器内插件，发现其在圆管内的强化行为与二次流强度有很大的关系。Habibi [7]等在层流流动下对扭曲管中的流动与换热进行了数值研究，分析了管壁温度、雷诺数的不同与二次流动拓扑结构和传热的内在一致性，以上文献均表明二次流强度与换热有相同的分布现象。故本文基于二次流强化传热机理提出在通道表面冲压凹凸金字塔单元的新颖设计方案，且重点研究了金字塔单元结构参数与二次流强度的内在联系，以期对新型结构PCHE的优化设计提供一定的指导意义。

2. 模型

2.1. 物理模型

如图1所示为金字塔结构壁面PCHE的计算域示意图。由于PCHE在垂直于主流方向的两侧具有周期性结构，故选取其中一段金字塔结构壁面通道作为计算区域，同时为了减少入口效应和消除出口端回流影响，在金字塔结构段的开始与结尾处分别添加10 mm和20 mm的延伸段，而换热区是90.6 mm的核心换热区域。

Figure 1. Schematic of the PCHE computational domain for the walls of a pyramid structure

图1. 金字塔结构壁面PCHE计算域示意图

图2为金字塔结构壁面PCHE的结构参数示意图，通过定义相对长度L/D_i、相对宽度W/C、相对高度H/D_i三个无量纲参数来表征金字塔单元的几何结构。图2(a)中D_i为通道内壁面冲压金字塔单元后通过圆心的两个金字塔塔顶或塔底之间的距离；图2(b)中C是金字塔粗糙元在X方向的周向长度，即以D_i为直径的所在圆周长；而L、W、H分别是金字塔单元的长度、宽度和高度，即图中A、B两点间的水平长度L，B、D两点间的水平宽度W和点S到直线BD之间的垂直高度H。金字塔单元的宽度W由周向金字塔单元个数决定；金字塔单元的长度L由参数相对长度L/D_i和D_i本身决定；而对于金字塔单元高度H也由参数相对高度H/D_i和D_i决定。表1给出了金字塔结构壁面PCHE的主要几何参数，其中完整计算模型的宽度W_s为2.5 mm，高度H_s为1.75 mm，管道直径D是1.51 mm。

2.2. 数学模型

本文假设PCHE通道中的流动与传热过程为无内热源的稳态湍流流动，忽略和环境间的热交换，又因工质为S-CO₂，故考虑了重力对流体的影响，其在三维直角坐标系中的控制方程如下[8]：

连续性方程：

$\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0$ (1)

动量方程：

$\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}=-\frac{\partial P}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\mu +{\mu }_{\text{t}}\right)\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}\left(\mu +{\mu }_{\text{t}}\right)\frac{\partial u_k}{\partial x_k}\right]+\rho g_i$ (2)

能量方程：

$\frac{\partial \left(\rho u_i{c}_{p}T\right)}{\partial x_i}=\frac{\partial }{\partial x_i}\left[\left(\frac{\mu }{Pr}+\frac{{\mu }_{\text{t}}}{P{r}_{\text{t}}}\right)\frac{\partial \left(c_{p}T\right)}{\partial x_i}\right]$ (3)

式中，μ_t为湍流黏度，由密度ρ、特征速度V_t、特征长度l_t和c_μ所组成，其中c_μ是常数；Pr和Pr_t是普朗特数和湍流普朗特数。而剪切应力输运模型(Shear Stress Transport k-ω)因对变物性工质较高的适配性被选择为湍流模型[9]。

Figure 2. Schematic of PCHE 3D structural parameters of the wall of the pyramid structure. (a) 2D fluid domain; (b) Pyramid roughage

图2. 金字塔结构壁面PCHE三维结构参数示意图。(a) 二维流体域；(b) 金字塔粗糙元

Table 1. Main geometrical parameters of the wall PCHE of the pyramid structure

表1. 金字塔结构壁面PCHE的主要几何参数

Case	Ws (mm)	Hs (mm)	D (mm)	L/Di	W/C	H/Di
1	2.5	1.75	1.51	0.26	0.125	0.035
2	2.5	1.75	1.51	0.26	0.125	0.055
3	2.5	1.75	1.51	0.26	0.125	0.065
4	2.5	1.75	1.51	0.26	0.25	0.035
5	2.5	1.75	1.51	0.26	0.083	0.035
6	2.5	1.75	1.51	0.52	0.125	0.035
7	2.5	1.75	1.51	0.65	0.125	0.035

边界条件：

入口处(inlet)：

${\dot{m}}_z=const,T=const$ (4)

出口处(outlet)：

$\frac{\partial u}{\partial x}=\frac{\partial v}{\partial x}=\frac{\partial w}{\partial x}=\frac{\partial T}{\partial x}=0$ (5)

计算域前后壁面及入、出口延伸段的上下壁面采用绝热边界条件，金字塔结构核心换热区域上下壁面速度采用无滑移边界条件且热边界条件采用等热流(UHF)。

在换热通道内的流固耦合接触面上，速度采用无滑移边界条件，热边界条件通过下式定义：

${T|}_f={T|}_s,{\lambda }_f{\frac{\partial T}{\partial n}|}_f={\lambda }_s{\frac{\partial T}{\partial n}|}_s$ (6)

PCHE在X方向上具有周期性，故将计算域左右两侧定义为周期性壁面边界条件[10]：

$u_l=u_r,v_l=v_r,w_l=w_r,T_l=T_r$ (7)

$\frac{\partial u}{\partial x}=\frac{\partial v}{\partial x}=\frac{\partial w}{\partial x}=\frac{\partial T}{\partial x}=0$ (8)

参数定义：

雷诺数定义如下：

$Re=\frac{\rho u_{in}{D}_c}{\mu }$ (9)

T_b是主流方向横截面上的流体质量加权平均温度，定义如下：

$T_{\text{b}}=\frac{{\displaystyle {\int }_A\rho }u{c}_{p}T\text{d}A}{{\displaystyle {\int }_A\rho }u{c}_p\text{d}A}$ (10)

主流方向为Z向的涡量分量定义如下：

${\Omega }_{\text{z}}=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$ (11)

以参数Se来描述通道中的二次流强度，对于Z向为流动方向的平均二次流强度和局部二次流强度，本文采用下式定义[6]：

$Se=\frac{\rho D^2\left|{\Omega }_{\text{z}}\right|}{\mu }$ (12)

$S{e}_{\text{m}}=\frac{1}{V}{\displaystyle {\iiint }_{V}Se\text{d}V},S{e}_{\text{s}}\left(x\right)={\displaystyle {\iint }_{A_X}Se\text{d}{A}_X}/{\displaystyle {\iint }_{A_X}\text{d}{A}_X}$ (13)

对于涡结构的识别，本文采用Liu [11]等人提出的Omega准则，该准则将非均匀流场中任意一点的速度梯度张量分解为对称张量A和反对称张量B。董祥瑞，王义乾[12]等提出ε的近似表达式，是为了防止分母出现极小值时出现较大误差。

Omega (Ω)涡识别准则表达式为：

$\Omega =\frac{{\Vert B\Vert }^2}{{\Vert A\Vert }^2+{\Vert B\Vert }^2+\epsilon }$ (14)

2.3. 网格划分及数值验证

通过Fortran语言对金字塔结构壁面PCHE作网格划分，为了便于展示与绘制，创建的子区域均沿着Y轴左右对称分布，如图3(a)、图3(b)所示。为了提高网格质量，特对图3(a)箭头指示部分进行“Y”型切分。图3(c)是三维流体域网格，此部分通过双边界法进行网格生成。图3(d)是固体域三维网格，在该单元中，较困难的是箭头指示部分的建立，因为这是金字塔结构壁面在流动方向的延伸部分。最后将流体域和固体域的网格进行拼接形成计算域网格系统，如图3(e)所示。

Figure 3. Computational domain grid system. (a) 2D fluid domain; (b) 2D solid domain; (c) Fluid domain; (d) Solid domain; (e) Computational domain

图3. 计算域网格系统。(a) 二维流体域；(b) 二维固体域；(c) 流体域；(d) 固体域；(e) 计算域

本文选取金字塔单元相对长度、相对宽度、相对高度分别为0.26、0.125、0.035条件下的金字塔壁面流道PCHE进行网格无关性验证。验证结果如表2所示，以第四套网格为基准，计算其他几套与其之间的相对误差，随着网格数量的增加，Se_m的最大相对误差在减小，尤其第三套网格的Se_m最大相对误差仅有0.71%，说明此时的数值解已满足网格独立性要求，故综合考虑下选择第三套网格进行数值模拟。同时本文选取与相梦如[13]相同的几何模型进行数值方法的准确性验证，将模拟结果与文献的数据进行对比，结果如图4所示。上下壁面的温度和周向平均对流换热系数的最大相对误差在5%之内，表明本文所采用的数值计算方法是合理的。

Table 2. Grid independence assessment

表2. 网格独立性考核

编号	数量(万)	Sem	ΔSem (%)
1	6,583,289	6639.54	10.94
2	7,501,495	7117.11	4.54
3	9,368,040	7402.38	0.71
4	12,442,587	7455.21	-

Figure 4. Validation of numerical accuracy

图4. 数值准确性验证

3. 结果分析

3.1. 金字塔相对高度对二次流特性影响

较大的相对高度使金字塔塔顶相对于塔底的倾斜程度更高，故增强了对流体的扰动，而流体冲刷金字塔单元诱导产生了更多的涡旋，从而有效提高了Se。如图5(a)所示，当Re相同时，相对高度越大，Se_m越大；当相对高度不变时，Se_m随着Re的增大而增大，这是由于有更多的流体冲刷金字塔单元的缘故。如图5(b)所示，当H/D_i = 0.035时，Se_s在拟临界点前达到局部极大值，且随着相对高度的增大，Se_s在拟临界点之后明显增大，从而使换热性能明显提升。

3.1.1. 不同相对高度下主流方向截面局部二次流云图

如图6所示，当相对高度较小时，金字塔塔顶相对于塔底的倾斜程度较低，对流体诱导产生的涡旋有限，导致对流体的扰动不大，因此Se较小。但随着相对高度的增大，金字塔单元对流体的扰动作用诱导产生了更多的涡旋，使得通道内的Se得到显著提升。

3.1.2. 不同相对高度下主流方向涡结构分布

Ω准则作为一种新型的涡识别方式被提出，其阈值进行了归一化处理，通常在0~1之间变化。当Ω小于0.5时，对称张量占优，此时涡的浓度低，当Ω大于0.5，涡量超过变形，就认为漩涡形成，故本小节选取了Ω = 0.52、Ω = 0.6、Ω = 0.7的等值面去识别漩涡表面，进而显示流场中的涡结构，如图7所示。可以发现，金字塔单元诱导流体产生了涡旋，且金字塔塔顶处的涡旋形状较大，而金字塔塔底处仅诱导了小漩涡，同时随着相对高度的增大，金字塔单元处产生的涡旋形状以及数量均有连续增强，这表明较大的相对高度会给予流体更强的扰动，使管内的Se水平较高。

3.2. 金字塔相对长度对二次流特性影响

在保证相对高度和宽度一定时，即不改变金字塔单元在垂直主流方向和周向上对流体的阻碍作用，轴向上较小数目的金字塔单元使管壁在沿着流动方向上趋于流线型，即金字塔单元在轴向上的倾斜程度较低，此时诱导的二次流强度有限。当相对长度减小，即主流方向的金字塔单元数目增多，此时较大的倾斜程度对流体的扰动诱导产生了较多涡旋，故二次流强度较大。如图8(a)所示，Se_m随着Re的增大而增大，说明当质量流量增大时，通道内所形成的涡旋更加显著。又如图8(b)所示，不同相对长度下的Se_s分布趋势类似，均在拟临界点前达到局部极大值。

Figure 5. Comparison of PCHE secondary flow intensities at the wall of a pyramid structure at different relative heights. (a) Re-Se_m; (b) Z-Se_s

图5. 不同相对高度下金字塔结构壁面PCHE二次流强度比较。(a) Re-Se_m；(b) Z-Se_s

Figure 6. Localized secondary flow intensity cloud along the flow direction section at different relative heights. (a) H/D_i = 0.035; (b) H/D_i = 0.055; (c) H/D_i = 0.065

图6. 不同相对高度下沿流动方向截面局部二次流强度云图。(a) H/D_i = 0.035；(b) H/D_i = 0.055；(c) H/D_i = 0.065

Figure 7. Distribution of vortex structure at different relative heights. (a) H/D_i = 0.035; (b) H/D_i = 0.055; (c) H/D_i = 0.065

图7. 不同相对高度下的涡结构分布。(a) H/D_i = 0.035；(b) H/D_i = 0.055；(c) H/D_i = 0.065

Figure 8. Comparison of PCHE secondary flow intensities at the walls of pyramidal structures with different relative lengths. (a) Re-Se_m; (b) Z-Se_s

图8. 不同相对长度下金字塔结构壁面PCHE二次流强度比较。(a) Re-Se_m；(b) Z-Se_s

3.2.1. 不同相对长度下主流方向截面局部二次流云图

如图9所示，不同相对长度下的Se_s分布趋势类似，壁面处的Se稍大，这是由于金字塔单元的存在增大了此处对流体的扰动。当相对长度越小时，沿着流动方向的金字塔数目增多，流体冲刷凹凸的金字塔单元受到的扰动增大，涡旋数量增加，这使得通道内诱导产生的Se剧烈。然而随着相对长度的增大，由于金字塔数目在轴向的减少，导致Se受限。

3.2.2. 不同相对长度下主流方向涡结构分布

本小节选取了Ω = 0.52、Ω = 0.55、Ω = 0.6的阈值进行对比分析，如图10所示。在相对长度为0.26时，涡旋形状小、数量多，且金字塔塔底处出现了小漩涡；而相对长度为0.52时，金字塔塔底处虽未出现小漩涡，但金字塔塔顶处的漩涡形状明显增大，故二者所诱导的Se相差无几。当相对长度为0.65时，漩涡的数量较少，Se也较低。

Figure 9. Localized secondary flow intensity cloud along the flow direction section for different relative lengths. (a) L/D_i = 0.26; (b) L/D_i = 0.52; (c) L/D_i = 0.65

图9. 不同相对长度下沿流动方向截面局部二次流强度云图。(a) L/D_i = 0.26；(b) L/D_i = 0.52；(c) L/D_i = 0.65

Figure 10. Distribution of vortex structures at different relative lengths. (a) L/D_i = 0.26; (b) L/D_i = 0.52; (c) L/D_i = 0.65

图10. 不同相对长度下的涡结构分布。(a) L/D_i = 0.26；(b) L/D_i = 0.52；(c) L/D_i = 0.65

3.3. 金字塔相对宽度对二次流特性影响

随着相对宽度的增加，换热器沿着周向越接近于流线型，使诱导的Se有限。而当相对宽度减小时，过多的金字塔单元对流体造成的扰动较大，流体冲刷金字塔单元诱导了较多的漩涡，使Se提升明显。正如图11(a)中当Re相同时，相对宽度越大，Se_m就越小。同时还发现管内Re在较小时，出现了大相对宽度条件下的Se小于平直通道的现象，随着Re的增大，这种现象才被消退。如图11(b)所示，不同相对宽度下的Se_s分布趋势类似。随着相对宽度的增大，导致在管壁周向的金字塔单元数量减少，故金字塔单元所诱导的涡旋减少，使Se_s有限。

Figure 11 Comparison of PCHE secondary flow intensities at the walls of pyramidal structures with different relative widths. (a) Re-Se_m; (b) Z-Se_s

图11. 不同相对宽度下金字塔结构壁面PCHE二次流强度比较。(a) Re-Se_m；(b) Z-Se_s

3.3.1. 不同相对宽度下主流方向截面局部二次流云图

如图12所示，当相对宽度越小时，表明在周向的金字塔单元宽度变窄，这也同时说明在沿着流动方向上的金字塔数目增多，过多的金字塔单元在冲刷流体时造成较强扰动，故使得通道内产生的Se较大。且由于金字塔单元塔顶和塔底对流体的涡旋作用，使得在此处的Se明显增大。而在相对宽度为0.25时，周向金字塔数量显著降低，此时对流体造成扰动的因素减少，故此时通道内诱导产生的Se水平较低。

Figure 12. Localized secondary flow intensity cloud along the flow direction cross-section for different relative widths. (a) W/C = 0.125; (b) W/C = 0.25; (c) W/C = 0.083

图12. 不同相对宽度下沿流动方向截面局部二次流强度云图。(a) W/C = 0.125；(b) W/C = 0.25；(c) W/C = 0.083

3.3.2. 不同相对宽度下主流方向涡结构分布

本小节选取了Ω = 0.52、Ω = 0.57、Ω = 0.62的阈值进行对比分析，如图13所示。当相对宽度为0.083时，由于周向的金字塔数量增多，使金字塔单元处的涡旋数量也明显增加，且金字塔塔顶和塔底处均出现了大小不一的涡旋，故此时的Se很高。然而随着相对宽度的增大，金字塔数量的减少导致涡旋的数量减少，尤其在相对宽度为0.25时，涡旋的数量及强度明显降低，诱导的Se水平也有限。

Figure 13. Distribution of vortex structures at different relative widths. (a) W/C = 0.125; (b) W/C = 0.25; (c) W/C = 0.083

图13. 不同相对宽度下的涡结构分布。(a) W/C = 0.125；(b) W/C = 0.25；(c) W/C = 0.083

3.4. S-CO₂物性对二次流特性影响

不同热流密度、进口温度下的Se_s分布趋势类似，如图14所示。在局部极大值之前，Se_s随着热流密度、进口温度的增大而增大；而在局部极大值之后，Se_s随着热流密度、进口温度的增大而减小，这种现象是由于S-CO₂的定压比热容发生了改变导致的。如图15(a)所示，在拟临界压力不变的情况下，随着进口温度的增大，不同流体温度下S-CO₂的拟临界点向入口端前移。这一点在图15(b)中也有所展示，定压比热容随着进口温度的增大而前移，同时均在拟临界点附近达到最大值且峰值基本保持不变，且S-CO₂的物性定压比热容在峰值之前随着进口温度的增大而增大，而在峰值之后随着进口温度的增大而减小，此分布可以很好的解释图14(b)中Se_s的变化。

Figure 14. Localized secondary flow distribution for different heat flux and inlet temperatures. (a) q; (b) T_in

图14. 不同热流密度和进口温度下的Se_s分布。(a) 热流密度；(b) 进口温度

Figure 15. Distribution of fluid temperature and constant pressure specific heat capacity at different inlet temperatures. (a) T_b; (b) C_p

图15. 不同进口温度下的流体温度和定压比热容分布。(a) 流体温度；(b) 定压比热容

4. 结论

本文研究了结构参数和物性变化对金字塔结构壁面PCHE通道内Se的影响，得到的主要结论如下：

(1) 当相对高度增大，塔顶相对于塔底的倾斜程度更高，流道内诱导产生较多涡旋，使通道内的Se明显提升。

(2) 当相对长度和相对宽度减小，较多的金字塔数目给予流体更强的扰动，Se水平也得到较大提升。

(3) S-CO₂的定压比热容在拟临界点附近达到峰值，且定压比热容分别在峰值之前、之后随着进口温度的增大而增大、减小，此变化解释了不同进口温度下Se_s的分布情况。

致谢

感谢甘肃省高等学校产业支撑计划项目(2023CYZC-37)的资助。

基金项目

甘肃省高等学校产业支撑计划项目(2023CYZC-37)。

NOTES

^*通讯作者。