耦合氢正仲催化转化过程的低温板翅式换热器的校核方法

1.本发明涉及耦合氢正仲催化转化的低温板翅式换热器技术领域，尤其涉及一种耦合氢正仲催化转化过程的低温板翅式换热器的校核方法。


背景技术：

2.氢是双原子分子，两个氢原子核是绕轴自转的。根据两个核自旋的相对方向，氢分子可分为正氢(ortho-h2)和仲氢(para-h2)，简写为o-h2和p-h2，通常的氢是这两种形式氢分子的混合物，正仲氢之间的平衡百分比仅与温度有关，室温以上的温度时，一般称为正常氢，含正氢75％，仲氢25％。平衡氢中仲氢浓度百分比会随着温度的降低而降低，当温度降低氢气液化时，正氢会自发转换为仲氢，并且释放出热量，引起储存的液氢大量汽化，甚至使得储存第一天的蒸发总量达到了总储存量的20％以上，即使把液氢储存在一个理想绝热的容器中，它同样也会发生汽化，在开始的24h内，液氢大约要蒸发损失18％，100h后损失将超过40％，这会导致液氢的损失，同时使得罐内压力上升，因此在设计液氢储罐时要留有一定的裕量空间，这导致氢生产、储存和使用的成本增加。在自然状态下，正仲转化的速率非常缓慢，因此在国外成熟的氢液化设备中，一般都采用了多级催化，在氢液化的降温过程中加速将正氢转化为接近平衡浓度的仲氢，得到仲氢含量95％以上的液氢产品，以减少正仲氢转化引起的液氢蒸发损失。
3.板翅式换热器换热器是一种高效、紧凑、轻巧的换热设备，其中的翅片结构是一种理想的催化剂载体，可将其应用于氢正仲连续转化过程，连续转化是将氢正仲转化热及时通过换热器内冷热流体热交换过程移除，使得氢在转化过程中始终保持对应温度下的接近平衡氢的浓度，是公认的耗能最小的转化方式。板翅式换热器流动换热与氢正仲催化转化一体化整合的优势为(1)氢气的冷却换热和正仲转化反应同时进行，不仅加速氢的正仲转化速率，且可有效提高低温冷能的利用效率；(2)填充至换热器通道内的催化剂颗粒使得高压侧氢气流体换热明显增强，可有效减少换热面积和缩小换热设备体积；(3)换热器与正仲转化器合二为一，略去传统的正仲转化反应器及其配套冷却设备，提高了氢液化工艺系统的紧凑性。
4.工业中采用的板翅式换热器一般超过两股流体，最多可能会发生20多股流体的流动换热过程，翅片通道数量达到几百层，若其中再集成氢正仲催化转化过程，则整个流动换热过程非常复杂，常规的cfd数值模拟过程显得力不从心，而主流的换热器设计软件如htri、aspen muse和apsen edr还未在换热器设计校核模型中集成氢正仲催化转化过程。总体来看，构建一种耦合氢正仲催化转化的板翅式换热器设计校核模型，并完善为一种系统的方法，编写成通用的计算程序，对促进氢能源高效利用，保障国家能源安全、培育新的经济增长点具有重要的战略意义。


技术实现要素：

5.为了克服上述现有技术的缺点，本发明的目的在于提供一种耦合氢正仲催化转化
过程的低温板翅式换热器的校核方法，在换热器计算中引入氢正仲催化转化，计算结果包括板翅式换热器整体的温度场、压力场和仲氢浓度场分布，促进氢能源高效利用。
6.为了达到上述目的，本发明采取的技术方案为：
7.一种耦合氢正仲催化转化过程的低温板翅式换热器的校核方法，包括以下步骤：
8.1)给定参数：
9.选取多股流的板翅式换热器，确定板翅式换热器整体结构参数、工况参数以及各股流体的翅片结构参数；
10.整体结构参数包括换热器的排布方式、换热器芯体的长l、宽w；
11.工况参数包括多股流的流体、流体通道层、质量流量、入口温度、入口压力；
12.翅片的结构参数为：选择平直翅片时，参数为翅高hi、翅距si、翅厚ti；选择打孔翅片时，参数为翅高hi、翅距si、翅厚ti、孔隙率φi；选择锯齿翅片时，参数为翅高hi、翅距si、翅厚ti、节距li；
13.流体通道中填充的氢正仲转化催化剂颗粒的参数是：空隙率ε是0.35，空隙率的定义是填料层内颗粒间空隙体积占总体积的百分比，颗粒直径是0.371mm，催化剂颗粒的导热系数是0.58w
·
m-1
·
k-1
；
14.2)模型简化及初始化：
15.假设流体通道中流体沿着z方向均匀分布，将三维模型简化至二维平面内，另外，在计算中只考虑稳态，认为换热器与周围环境绝热，在流体计算中忽略了导热项,模型的构建中考虑了流体域、翅片和隔板固体域，以及催化剂颗粒固体域，在发生氢正仲催化反应的通道中，采用假设：氢气是由正氢和仲氢组成的单相可压缩气体，通道内的催化剂填料均匀，颗粒近似为球体且直径一致；
16.一共有n+1层隔板和n层翅片通道，在隔板中沿着流动方向布置n个节点，而在流体通道中沿着流动方向布置n+1个节点，由隔板中的每个隔板温度节点t
w,i,k
定义一个隔板单元，而由流体通道中的每两个流体温度节点t
i,k
和t
i,k+1
定义一个流体单元，流体温度节点还用来储存压力数据p
i,k
和p
i,k+1
，若流体通道中还发生氢的正仲催化转化，则在相应的温度节点处储存仲氢浓度pa
i,k
和pa
i,k+1
，在一个板翅式换热器中一共设置了n
×
(n+1)个隔板单元和n
×
n个流体单元；
17.整个计算过程是迭代，在第一个外迭代轮次的计算之前，需要给定一个初始化的温度场、压力场和仲氢浓度场；每层流体通道中的流体温度节点初始化为步骤1)中给定的入口温度，温度节点中储存的压力值初始化为步骤1)中给定的入口压力，在发生氢正仲催化转化的流体通道中，温度节点中还储存了仲氢浓度值，初始化为步骤1)中给定的入口仲氢浓度；
18.3)流体物性拟合：
19.所需要的流体物性包括比热容、导热系数、动力粘度、密度和普朗特数，物性数据从nist数据库中提取，把流体的密度与温度和压力的关系拟合在样条曲面上，其他类型流体物性则分别与温度拟合在样条曲线上，对于固体的导热系数，拟合成导热系数-温度样条曲线，最后将所有拟合好的样条曲线和样条曲面的参数加载在内存中；
20.在发生氢正仲催化转化的翅片通道中，正氢和仲氢的物性需分开拟合，若需要计算某仲氢浓度下氢流体的物性时，根据拟合好的正氢和仲氢的物性插值计算；
21.4)计算板翅式换热器流体通道中的压力场分布：
22.先计算得到流体通道中每个流体单元的压力降，根据流体通道中是否发生氢的正仲催化转化，将压力降的计算方式分为两种：
[0023][0024]
式(1)是计算由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元中的压力降，若该层流体通道中未填充催化剂颗粒，即未发生氢正仲催化转化，则采用第一种方法，通过该流体单元中的摩擦因子f
i,k
来计算，另外gi是每层翅片通道中流体的质量流速，δl
x
是流体单元沿着流动方向的长度，di该层流体通道中翅片结构的当量直径，ρ
i,k+1
和ρ
i,k
由流体温度节点t
i,k+1
和t
i,k
以及储存在其中的压力值p
i,k+1
和p
i,k
所确定的密度，是是流体单元的密度，f
i,k
和根据该流体单元中的定性温度和定性压力来确定：
[0025][0026]
在流体通道中布置平直翅片、打孔翅片或者锯齿翅片，对锯齿翅片采用如下的方法来计算摩擦因子：
[0027][0028]
式中的α,θ,γ分别是：
[0029][0030]
对于打孔翅片采用如下的方法计算摩擦因子：
[0031]
lnf
i,k
＝28.78906-12.31399lnre
i,k
+1.565191(lnre
i,k
)
2-0.06736098(lnre
i,k
)3ꢀꢀ
(5)
[0032]
对于平直翅片采用如下的方法计算摩擦因子：
[0033][0034]
re
i,k
是雷诺数：
[0035][0036]
是由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元的动力黏度，根据其中的定性温度和定性压力来确定；
[0037]
若该层流体通道中填充了氢正仲催化转化的催化剂颗粒，则需要采用式(1)中的第二种方法来计算流体单元中的压力降，式(1)中的是该流体单元中的流体表观速度，即不考虑流体通道中填充的催化剂颗粒所计算得到的流体单元中的平均速度，α是催化剂
颗粒的渗透率，c1是催化剂颗粒的惯性阻力系数，两个参数的计算方法分别是：
[0038][0039][0040]
式中d
p
是催化剂颗粒的直径，ε是空隙率；
[0041]
根据式(1)计算得到每个流体单元的压力降，在给定每层流体通道入口压力的情况下，即沿着流动方向逐点更新储存在每个温度节点中的压力值；将更新得到的压力场分布与旧的压力场分布进行比较，计算压力场残差，若残差不满足收敛性条件，则重复步骤2)，即在旧压力场的基础上继续更新压力场，直至压力场残差满足收敛性条件，进入步骤5)，在后续步骤的计算中采用了此步更新得到的压力场分布作为各个流体单元的定性压力；
[0042]
5)板翅式换热器隔板温度场分布计算：
[0043]
当隔板相邻的流体通道中填充了催化剂颗粒，此时流体向隔板传导的热量分为三部分，一部分是流体热量通过翅片表面传至翅片固体域，然后热量经翅片固体域传导至隔板固体域中；另一部分是流体热量通过催化剂表面传至催化剂固体域，然后通过催化剂固体域传导至隔板固体域中；第三部分是流体的热量直接通过隔板表面传至隔板固体域中。若隔板相邻的流体通道中未填充催化剂，则流体向隔板传导的热量仅包含第一部分和第三部分热量；
[0044]
对于翅片和催化剂颗粒固体域，为了简化计算，不考虑沿着轴向即+x或者-x方向的导热，那么由温度节点t
i,k
和t
i,k+1
所定义的流体单元中的翅片和催化剂颗粒固体域的一维导热方程是：
[0045][0046]
式(10)求得翅片和催化剂固体域的温度t
s,i,k
沿着y方向的温度分布，λ
s,i,k
和λ
p,i,k
分别是翅片和催化剂颗粒固体域的导热系数，是一维导热方程中的换热源项，包括了翅片表面与流体之间的换热，以及催化剂颗粒表面与流体之间的换热：
[0047][0048]
式(11)中λ
s,i,k
和λ
p,i,k
分别由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元中翅片和催化剂颗粒固体域的导热系数，α
p
是催化剂颗粒的比表面积，在假设催化剂颗粒大小均匀且均为正球体的情况下，计算方法是：
[0049]
α
p
＝6(1-ε)/d
p
ꢀꢀ
(12)
[0050]hc,i,k
是流体通道中流体和固体域表面的对流换热系数，当流体通道中未填充催化
剂时，h
c,i,k
的计算方法是：
[0051][0052]
式中pr
i,k
是由温度节点t
i,k
和t
i,k+1
所定义的流体单元的普朗特数，λ
fl,i.k
是该流体单元的导热系数，j
i,k
该流体单元中用来评价翅片表面换热能力的评价指标，对锯齿型翅片，采用如下的j因子关联式：
[0053][0054]
对于打孔翅片，采用关联式：
[0055]
lnj
i,k
＝34.57583-15.92678lnre
i,k
+2.137607(lnre
i,k
)
2-0.09544151
×
(lnre
i,k
)3ꢀꢀ
(15)
[0056]
对于平直翅片，采用关联式：
[0057][0058]
若流体通道中填充了催化剂颗粒，则流体单元中固体域表面的对流换热系数h
c,i,k
的计算方法是：
[0059][0060]
式中re
p
是取颗粒直径为特征尺寸，以及流体单元中表观速度所计算得到的雷诺数；
[0061]
根据式(10)、(11)推导得到由温度节点t
i,k
和t
i,k+1
所定义的流体单元中固体域沿着y方向的温度分布：
[0062][0063]
式中，sh和ch分别是双曲余弦和双曲正弦，θ
i,k
、θ
i-i,k
和θ
i-(i+1),k
为翅片、翅片顶部和根部的过余温度；m
e,i,k
、θ
i,k
，θ
i-i,k
和θ
i-(i+1),j
表示如下：
[0064][0065]
当得到流体单元中包含的固体域中的温度分布，通过傅里叶导热定律求得流体单元包含的固体域与隔板固体域的导热量；对于一个隔板单元，来源于相邻流体单元的热量还有第四部分热量，是该隔板单元与相邻隔板单元的导热量，隔板单元中的这四部分热量满足能量守恒，经推导得到计算隔板温度节点的代数方程组：
[0066][0067][0068]ai-1,k+1
＝a
i-1,k
[0069][0070]ai,k+1
＝a
i,k
[0071][0072][0073][0074]
[0075]aw,i,k
＝a
i-1,k
+a
i-1,k+1
+a
i,k
+a
i,k+1
+a
w,i-1,k
+a
w,i+1,k
+a
w,i,k-1
+a
w,i,k+1
ꢀꢀ
(21)
[0076]
式(21)中t
sp
是板翅式换热器的隔板厚度，λw是隔板的导热系数，式中在括号外标注了下标(i,k),则括号内的参数的下标标注分为三种情况，若是参数δl
x
,w,t
sp
,ε,α
p
，则不需要有下标，若是参数s,t,h，则下标为i，若是参数λw,λs,hc,me，则下标为(i,k)；若中括号外的下标是(i-1,k)或者(i,k-1)，则中括号内的参数按照同样的思路标注，总之括号内的参数标注与文中其它公式中出现的该参数的标注方法相同；λ
w,i,k
是隔板温度节点t
w,i,k
和t
w,i,k+1
的界面上的导热系数，若k＝n,则直接根据隔板温度节点t
w,i,n
来确定λ
w,i,n
,λ
w,i,k-1
是隔板温度节点t
w,i,k-1
和t
w,i,k
的界面上的导热系数，若k＝1，则直接根据隔板温度节点t
w,i,1
来确定λ
w,i,1
，另外，式(21)的中m
e,i,k
的计算要根据式(19)分为两种情况；
[0077]
对每一个由隔板温度节点所定义的隔板单元，都要求解形如式(20)所示的代数方程，一共需要求解n
×
(n+1)个方程，为了加快隔板温度节点的计算速度，采用结合了高斯-赛德尔(gauss-seidel)法的线迭代法(line iteration)，把隔板固体域求解分解为平行y轴的多条线，每条线只包含一层温度节点，在同一条线上各节点的值是用代数方程的直接解法来获得，aa'线上的隔板温度节点所构成的块，同一块内各节点的值是以隐含的方式联系着，但是沿着+x方向从一条线向另一条的推进是用迭代的方式进行：
[0078][0079]
式中t
p
代表了所要求的隔板温度节点，tn,tw,ts,te(north,west,south,east)是相邻的隔板温度节点,n0理解为外迭代轮次的标识符，b代表了不能归类为相邻隔板温度节点的影响的换热源项，由于扫描方向是沿着+x方向，即从左至右，因此tn,tw,ts采用了已经更新的隔板温度(图4所示)，而te采用了前一个外迭代轮次更新得到的温度值，对于一条线上的每一个隔板温度节点，都要求解形如式(22)的代数方程组，该条线上一共包含了n+1方程，采用基于gauss消元法的thomas算法对方程组进行联立求解；
[0080]
当从左至右更新了所有隔板温度节点，进入步骤6)；
[0081]
6)流体通道中温度分布及仲氢浓度分布计算：
[0082]
若编号为i的流体通道中发生氢正仲催化转化过程，则温度节点t
i,k
和t
i,k+1
所定义的流体单元中的氢正仲催化转化热是：
[0083][0084]
式中s
o-p,i,k
是由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元中的氢的正仲转化热源项，m是氢分子的摩尔质量，δh
i,k
是该流体单元中氢正仲催化反应的反应热，反应热是温度的函数：
[0085][0086]
单位体积内氢正仲转化的速率根据elovich模型计算：
[0087]
[0088]
式中，pc是氢的临界压力，是由温度节点t
i,k
和t
i,k+1
所定义的流体单元中的仲氢浓度，为elovich计算模型的反应速率常数；a、b1、c和d为实验数据拟合系数，取值分别是1.0924,59.7mol
·
m-3
·
s-1
，-253.9mol
·
m-3
·
s-1
，-11.6mol
·
m-3
·
s-1
；为平衡氢中仲氢含量，计算方法是：
[0089][0090]
式中tc是氢的临界温度；
[0091]
流体通道中温度节点的计算分为两种情况，若流体通道中未发生氢正仲催化转化过程，则流体单元中的流体与隔板固体域表面和翅片固体域表面换热，若流体通道中发生氢正仲催化反应，流体单元中的流体还与催化剂固体颗粒表面换热，而且氢正仲催化转化热也为流体提供热量；流体单元中的这几部分热量来源满足能量守恒，再结合前面的式(18)、(19)推导得到计算流体温度节点的方程组，若流体沿着+x方向流动，则推导得到方程组的形式为：
[0092][0093]
式中的影响系数b
i,k
，b
w,i,k
,b
w,i+1,k
,b
i,k+1
的计算方法是：
[0094][0095][0096]bi,k+1
＝b
i,k
+b
w,i,k
+b
w,i+1,k
ꢀꢀ
(29)
[0097]
式中括号外的下标为(i,k),则括号内m是mi,代表了编号为i的流体通道中的质量流量，c
p
写作c
p,i,k
，是由温度节点t
i,k
和t
i,k+1
所定义的流体单元的比热容；若在流体通道中发生氢正仲催化转化，则式(29)中的m
e,i,k
需按照式(19)的流体通道中填充催化剂时的计算方法来计算，此时式(28)不封闭，需补充仲氢的组分输运方程：
[0098][0099]
若流体沿着+x方向流动，且流体通道中发生氢正仲催化转化，则根据式(28)、(30)沿着+x方向逐点确定流体温度节点中储存的仲氢浓度值；
[0100]
若流体通道沿着-x方向，计算流体单元温度节点的方程组是：
[0101][0102]
影响系数的计算方法是:
[0103][0104][0105]bi,k
＝b
i,k+1
+b
w,i,k
+b
w,i+1,k
ꢀꢀ
(32)
[0106]
而流体单元中仲氢的组分输运方程是：
[0107][0108]
通过迭代的方法，式(31)、(33)确定流动是-x方向的流体通道中的温度场和仲氢浓度场分布；
[0109]
7)温度场残差和仲氢浓度场残差计算：
[0110]
将步骤5)计算得到的隔板温度场分布与前一个外迭代轮次确定的隔板温度场分布进行比较，计算隔板温度场残差，将步骤6)计算得到的流体温度场分布与前一个外迭代轮次确定的流体温度场分布进行比较，计算流体温度场残差，取隔板温度场残差和流体温度场残差的最大值作为温度场残差；取步骤6)计算得到的仲氢浓度场与前一个外迭代轮次确定的仲氢浓度场分布进行比较，计算仲氢浓度场残差；若温度场残差和仲氢浓度场残差有一个不满足收敛性条件，则继续重复步骤4)～7)，作为新一轮外迭代,直至两种残差都满足收敛性条件，即输出计算结果，计算结果包括了板翅式换热器整体的温度场分布，流体通道中的压力场分布和仲氢浓度分布；温度场包括隔板温度场和流体温度场。
[0111]
所述的步骤5)中对于隔板温度节点，按需要引入逐次亚松弛：
[0112]
t
n0
＝t
n0-1
+α
sur
(t
n0-t
n0-1
),0＜α
sur
≤1
ꢀꢀ
(23)
[0113]
t
n0
是第n0外迭代轮次中的隔板温度场，α
sur
是松弛因子，该参数小于1时为逐次亚松弛迭代。
[0114]
所述的步骤6)中根据式(28)、(30)沿着+x方向逐点确定流体温度节点中储存的仲氢浓度值，对于由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元，已知和t
i,k
，具体步骤如下:
[0115]
a.预估出口仲氢浓度
[0116]
b.计算(计算所需的t
i,k+1
来源于前一个外迭代轮次所确定的流体温度场)，再根据式(26)、(27)计算
[0117]
c.根据式(28)、(30)计算和t
i,k+1
；
[0118]
d.比较新计算的和旧的若残差不满足收敛性条件，则继续从步骤
b开始进行新一轮计算，直至的残差满足收敛性条件；
[0119]
当确定了和t
i,k+1
，继续确定下一个流体单元出口处的和t
i,k+2
，直至计算得到编号为i的流体通道的出口和t
i,n+1
。
[0120]
与现有技术相比，本发明具有以下的有益效果：
[0121]
1.本发明采用的氢正仲催化转化与板翅式换热器流动换热一体化的技术，可以增加换热，缩小设备体积，优化换热流程，减小系统组成部件的复杂程度。
[0122]
2.本发明是耦合氢正仲催化转化的板翅式换热的校核方法，考虑了流体的物性变化和隔板的轴向导热效应，在发生催化转化的通道中引入多孔介质压降源项、氢正仲催化转化热源项和仲氢的组分输运过程，计算速度快且精度高。
[0123]
3.本发明方法可写作一个通用的计算程序，可针对任意排列的耦合氢正仲催化转化的多股流逆流板翅式换热器进行校核计算，克服了传统cfd无法针对几百层翅片通道排布的板翅式换热器进行数值模拟的缺陷，也克服了传统板翅式换热器设计软件如htri、aspen muse和apsen edr等无法在换热器计算中引入氢正仲催化转化的问题，计算的结果包括板翅式换热器整体的压力场、温度场和仲氢浓度场的分布。
附图说明
[0124]
图1为填充氢正仲转化催化剂的多股流逆流板翅式换热器的三维结构(采用了锯齿型翅片)示意图。
[0125]
图2为板翅式换热器的二维简化模型以及局部的温度节点布置示意图。
[0126]
图3为板翅式换热器填充了催化剂的流体通道中的导热示意图。
[0127]
图4为板翅式换热器隔板温度节点分块迭代计算示意图。
[0128]
图5为本发明的流程示意图。
[0129]
图6为本发明实施例计算的其中三层通道中流体沿着流动方向的温度分布示意图。
[0130]
图7为本发明实施例计算的各层流体通道出口温度以及出口仲氢浓度示意图。
[0131]
图8为本发明实施例计算的某一层流体通道中仲氢浓度和平衡氢浓度沿着流动方向的分布示意图。
具体实施方式
[0132]
下面结合附图和实施例对本发明作进一步详细说明。
[0133]
一种耦合氢正仲催化转化过程的低温板翅式换热器的校核方法，包括以下步骤：
[0134]
1)给定参数：
[0135]
选取一个三股流的板翅式换热器，换热器的排布方式是(cbcbcaa*16)cbcbc，一共有117层流体通道，如图1所示，换热器芯体的长l是800mm，宽w是550mm；a流体是氢气，布置了32层流体通道，其中填充了氢正仲转化的催化剂，沿着+x方向流动，质量流量是0.064kg/s，入口的仲氢浓度是45％，入口温度是78.2k，入口压力是2.4mpa，a流体通道中翅片的结构参数为：打孔翅片，翅高hi是9.5mm,翅距si是1.4mm，翅厚ti是0.3mm,孔隙率φi是0.05；b流体是氦气，布置了34层流体通道，沿着+x方向流动，质量流量是1.285kg/s，入口温度是78.5k，入口压力是2.1mpa，b流体通道中翅片的结构参数为：锯齿翅片，翅高hi是6.3mm,翅
距si是1.2mm，翅厚ti是0.2mm，节距li是3mm；c流体是氦气，布置了51层流体通道，沿着-x方向流动，质量流量是1.285kg/s，入口温度是67.5k，入口压力是0.42mpa，c流体通道中翅片的结构参数为：锯齿翅片，翅高hi是6.35mm,翅距si是1.2mm，翅厚ti是0.2mm，节距li是3mm；
[0136]
a流体通道中填充的氢正仲转化催化剂颗粒的参数是：空隙率ε是0.35，空隙率的定义是填料层内颗粒间空隙体积占总体积的百分比，颗粒直径是0.371mm，催化剂颗粒的导热系数是0.58w
·
m-1
·
k-1
；
[0137]
2)模型简化及初始化：
[0138]
为了方便构建模型，假设流体通道中流体沿着z方向均匀分布，即可将三维模型(图1)简化至二维平面(图2)内研究，另外，在计算中只考虑稳态，认为换热器与周围环境绝热，在流体计算中忽略了导热项,模型的构建中考虑了流体域、翅片和隔板固体域，以及催化剂颗粒固体域，在发生氢正仲催化反应的通道中，采用假设：氢气是由正氢和仲氢组成的单相可压缩气体，通道内的催化剂填料均匀，颗粒近似为球体且直径一致；
[0139]
如图2所示，一共有n+1层隔板和n层翅片通道，在隔板中沿着流动方向(图2中x方向)布置n个节点，而在流体通道中沿着流动方向布置n+1个节点，由隔板中的每个隔板温度节点t
w,i,k
定义一个隔板单元，而由流体通道中的每两个流体温度节点t
i,k
和t
i,k+1
定义一个流体单元，流体温度节点还用来储存压力数据p
i,k
和p
i,k+1
，若流体通道中还发生氢的正仲催化转化，则在相应的温度节点处储存仲氢浓度pa
i,k
和pa
i,k+1
，在一个板翅式换热器中一共设置了n
×
(n+1)个隔板单元和n
×
n个流体单元(图2所示)；
[0140]
整个计算过程是迭代的，即每一个外迭代轮次是根据前一个外迭代轮次得到的旧场来计算出一个新场，因此在第一个外迭代轮次的计算之前，需要给定一个初始化的温度场、压力场和仲氢浓度场；每层流体通道中的流体温度节点初始化为步骤1)中给定的入口温度，温度节点中储存的压力值初始化为步骤1)中给定的入口压力，在发生氢正仲催化转化的流体通道中，温度节点中还储存了仲氢浓度值，初始化为步骤1)中给定的入口仲氢浓度；
[0141]
3)流体物性拟合：
[0142]
所需要的流体物性包括比热容、导热系数、动力粘度、密度和普朗特数，物性数据从nist数据库中提取，一般会把流体的密度与温度和压力的关系拟合在样条曲面上，其他类型流体物性则分别与温度拟合在样条曲线上，对于固体的导热系数，也存在与温度的函数关系，可拟合成导热系数-温度样条曲线，最后将所有拟合好的样条曲线和样条曲面的参数加载在内存中，供后面的计算调用；
[0143]
在发生氢正仲催化转化的翅片通道中，正氢和仲氢的物性需分开拟合，在后续计算中，若需要计算某仲氢浓度下氢流体的物性时，可根据仲拟合好的正氢和仲氢的物性插值计算；
[0144]
4)计算板翅式换热器流体通道中的压力场分布：
[0145]
为得到流体通道中的压力场分布，需要先计算得到流体通道中每个流体单元的压力降，根据流体通道中是否发生氢的正仲催化转化，将压力降的计算方式分为两种：
[0146][0147]
式(1)是计算由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元中的压力降，若该层流体通道中未填充催化剂颗粒，即未发生氢正仲催化转化，则采用第一种方法，通过该流体单元中的摩擦因子f
i,k
来计算，另外gi是每层翅片通道中流体的质量流速，δl
x
是流体单元沿着流动方向的长度，di该层流体通道中翅片结构的当量直径，ρ
i,k+1
和ρ
i,k
由流体温度节点t
i,k+1
和t
i,k
以及储存在其中的压力值p
i,k+1
和p
i,k
所确定的密度，是是流体单元的密度，f
i,k
和可根据该流体单元中的定性温度和定性压力来确定：
[0148][0149]
在流体通道中可能会布置平直翅片、打孔翅片或者锯齿翅片，对锯齿翅片采用如下的方法来计算摩擦因子：
[0150][0151]
式中的α,θ,γ分别是：
[0152][0153]
对于打孔翅片采用如下的方法计算摩擦因子：
[0154]
lnf
i,k
＝28.78906-12.31399lnre
i,k
+1.565191(lnre
i,k
)
2-0.06736098(lnre
i,k
)3ꢀꢀ
(5)
[0155]
对于平直翅片采用如下的方法计算摩擦因子：
[0156][0157]
re
i,k
是雷诺数：
[0158][0159]
是由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元的动力黏度，根据其中的定性温度和定性压力来确定；
[0160]
若该层流体通道中填充了氢正仲催化转化的催化剂颗粒，则需要采用式(1)中的第二种方法来计算流体单元中的压力降，式(1)中的是该流体单元中的流体表观速度，即不考虑流体通道中填充的催化剂颗粒所计算得到的流体单元中的平均速度，α是催化剂颗粒的渗透率，c1是催化剂颗粒的惯性阻力系数，两个参数的计算方法分别是：
[0161][0162][0163]
式中d
p
是催化剂颗粒的直径，ε是空隙率；
[0164]
根据式(1)计算得到每个流体单元的压力降，在给定每层流体通道入口压力的情况下，即可沿着流动方向逐点更新储存在每个温度节点中的压力值；将更新得到的压力场分布与旧的压力场分布进行比较，计算压力场残差，若残差不满足收敛性条件，则重复步骤2)，即在旧压力场的基础上继续更新压力场，直至压力场残差满足收敛性条件，进入步骤5)，在后续步骤的计算中采用了此步更新得到的压力场分布作为各个流体单元的定性压力；
[0165]
5)板翅式换热器隔板温度场分布计算：
[0166]
当隔板相邻的流体通道中填充了催化剂颗粒，催化剂颗粒扩展了流体和固体表面的换热面积，此时流体向隔板传导的热量分为三部分，一部分是流体热量通过翅片表面传至翅片固体域，然后热量经翅片固体域传导至隔板固体域中；另一部分是流体热量通过催化剂表面传至催化剂固体域，然后通过催化剂固体域传导至隔板固体域中；第三部分是流体的热量直接通过隔板表面传至隔板固体域中。若隔板相邻的流体通道未填充催化剂，则流体向隔板传导的热量仅包含第一部分和第三部分。
[0167]
对于翅片和催化剂颗粒固体域，为了简化计算，不考虑沿着轴向(+x或者-x方向)的导热，如图3所示，那么由温度节点t
i,k
和t
i,k+1
所定义的流体单元中的翅片和催化剂颗粒固体域的一维导热方程是：
[0168][0169]
式(10)可求得翅片和催化剂固体域的温度t
s,i,k
沿着y方向的温度分布，λ
s,i,k
和λ
p,i,k
分别是翅片和催化剂颗粒固体域的导热系数，是一维导热方程中的换热源项，包括了翅片表面与流体之间的换热，以及催化剂颗粒表面与流体之间的换热：
[0170][0171]
式(11)中λ
s,i,k
和λ
p,i,k
分别由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元中翅片和催化剂颗粒固体域的导热系数，α
p
是催化剂颗粒的比表面积，在假设催化剂颗粒大小均匀且均为正球体的情况下，计算方法是：
[0172]
α
p
＝6(1-ε)/d
p
ꢀꢀ
(12)
[0173]hc,i,k
是流体通道中流体和固体域表面的对流换热系数，当流体通道中未填充催化剂时，h
c,i,k
的计算方法是：
[0174][0175]
式中pr
i,k
是由温度节点t
i,k
和t
i,k+1
所定义的流体单元的普朗特数，λ
fl,i.k
是该流体单元的导热系数，j
i,k
该流体单元中用来评价翅片表面换热能力的评价指标，对锯齿型翅片，采用如下的j因子关联式：
[0176][0177]
对于打孔翅片，采用关联式：
[0178]
lnj
i,k
＝34.57583-15.92678lnre
i,k
+2.137607(lnre
i,k
)
2-0.09544151
×
(lnre
i,k
)3ꢀꢀ
(15)
[0179]
对于平直翅片，采用关联式：
[0180][0181]
若流体通道中填充了催化剂颗粒，则流体单元中固体域表面的对流换热系数h
c,i,k
的计算方法是：
[0182][0183]
式中re
p
是取颗粒直径为特征尺寸，以及流体单元中表观速度所计算得到的雷诺数；
[0184]
根据式(10)、(11)可以推导得到由温度节点t
i,k
和t
i,k+1
所定义的流体单元中固体域沿着y方向的温度分布：
[0185][0186]
式中，sh和ch分别是双曲余弦和双曲正弦，θ
i,k
、θ
i-i,k
和θ
i-(i+1),k
为翅片、翅片顶部和根部的过余温度；m
e,i,k
、θ
i,k
，θ
i-i,k
和θ
i-(i+1),j
表示如下：
[0187][0188]
当得到流体单元中包含的固体域中的温度分布，即可通过傅里叶导热定律求得流体单元包含的固体域与隔板固体域的导热量，实际上就是步骤5)开头讨论的三部分换热量中的前两部分；对于一个隔板单元，除了这三部分来源于相邻流体单元的热量，还有第四部分热量，是该隔板单元与相邻隔板单元的导热量，隔板单元中的这四部分热量满足能量守恒，经推导得到计算隔板温度节点的代数方程组：
[0189][0190][0191]ai-1,k+1
＝a
i-1,k
[0192][0193]ai,k+1
＝a
i,k
[0194][0195][0196][0197]
[0198]aw,i,k
＝a
i-1,k
+a
i-1,k+1
+a
i,k
+a
i,k+1
+a
w,i-1,k
+a
w,i+1,k
+a
w,i,k-1
+a
w,i,k+1
ꢀꢀ
(21)
[0199]
式(21)中t
sp
是板翅式换热器的隔板厚度，λw是隔板的导热系数，式中在括号外标注了下标(i,k),则括号内的参数的下标标注分为三种情况，若是参数δl
x
,w,t
sp
,ε,α
p
等，则不需要有下标，若是参数s,t,h，则下标为i，若是参数λw,λs,hc,me等，则下标为(i,k)；若中括号外的下标是(i-1,k)或者(i,k-1)，则中括号内的参数按照同样的思路标注，总之括号内的参数标注与文中其它公式中出现的该参数的标注方法相同；λ
w,i,k
是隔板温度节点t
w,i,k
和t
w,i,k+1
的界面上的导热系数，若k＝n,则直接根据隔板温度节点t
w,i,n
来确定λ
w,i,n
,λ
w,i,k-1
是隔板温度节点t
w,i,k-1
和t
w,i,k
的界面上的导热系数，若k＝1，则直接根据隔板温度节点t
w,i,1
来确定λ
w,i,1
，另外，式(21)的中m
e,i,k
的计算要根据式(19)分为两种情况；
[0200]
对每一个由隔板温度节点所定义的隔板单元，都要求解形如式(20)所示的代数方程，一共需要求解n
×
(n+1)个方程，为了加快隔板温度节点的计算速度，采用结合了高斯-赛德尔(gauss-seidel)法的线迭代法(line iteration)，把隔板固体域求解分解为平行y轴的多条线，每条线只包含一层温度节点，在同一条线上各节点的值是用代数方程的直接解法来获得的，图4中展示了aa'线上的隔板温度节点所构成的块，同一块内各节点的值是以隐含的方式联系着，但是沿着+x方向从一条线向另一条的推进是用迭代的方式进行：
[0201][0202]
式中t
p
代表了所要求的隔板温度节点，tn,tw,ts,te(north,west,south,east)是相邻的隔板温度节点,n0可理解为外迭代轮次的标识符，b代表了不能归类为相邻隔板温度节点的影响的换热源项，由于扫描方向是沿着+x方向，即从左至右，因此tn,tw,ts采用了已经更新的隔板温度，而te采用了前一个外迭代轮次更新得到的温度值，对于一条线上的每一个隔板温度节点，都要求解形如式(22)的代数方程组，该条线上一共包含了n+1方程，可以采用基于gauss消元法的thomas算法对方程组进行联立求解(参见陶文铨《数值传热学》第二版的章节4.4.1.1)；对于隔板温度节点，有时需要引入逐次亚松弛：
[0203]
t
n0
＝t
n0-1
+α
sur
(t
n0-t
n0-1
),0＜α
sur
≤1
ꢀꢀ
(23)
[0204]
t
n0
是第n0外迭代轮次中的隔板温度场，α
sur
是松弛因子，该参数小于1时为逐次亚松弛迭代；
[0205]
当从左至右更新了所有隔板温度节点，即可进入步骤6)；
[0206]
6)流体通道中温度分布及仲氢浓度分布计算：
[0207]
若编号为i的流体通道中发生氢正仲催化转化过程，则温度节点t
i,k
和t
i,k+1
所定义的流体单元中的氢正仲催化转化热是：
[0208][0209]
式中s
o-p,i,k
是由流体温度节点t
i,k
和t
i,k+1
所定义的流体单元中的氢的正仲转化热源项，m是氢分子的摩尔质量，δh
i,k
是该流体单元中氢正仲催化反应的反应热，一般来说，反应热是温度的函数：
[0210]
[0211]
单位体积内氢正仲转化的速率根据elovich模型计算：
[0212][0213]
式中，pc是氢的临界压力，是由温度节点t
i,k
和t
i,k+1
所定义的流体单元中的仲氢浓度，为elovich计算模型的反应速率常数；a、b1、c和d为实验数据拟合系数，取值分别是1.0924,59.7mol
·
m-3
·
s-1
，-253.9mol
·
m-3
·
s-1
，-11.6mol
·
m-3
·
s-1
。为平衡氢中仲氢含量，计算方法是：
[0214][0215]
式中tc是氢的临界温度；
[0216]
流体通道中温度节点的计算分为两种情况，若流体通道中未发生氢正仲催化转化过程，则流体单元中的流体与隔板固体域表面和翅片固体域表面换热，若流体通道中发生氢正仲催化反应，则除了前述的两部分热量来源，流体单元中的流体还与催化剂固体颗粒表面换热，而且氢正仲催化转化热也为流体提供热量；流体单元中的这几部分热量来源应当满足能量守恒，再结合前面的式(18)、(19)可推导得到计算流体温度节点的方程组，若流体沿着+x方向流动，则可推导得到方程组的形式为：
[0217][0218]
式中的影响系数b
i,k
，b
w,i,k
,b
w,i+1,k
,b
i,k+1
的计算方法是：
[0219][0220][0221]bi,k+1
＝b
i,k
+b
w,i,k
+b
w,i+1,k
ꢀꢀ
(29)
[0222]
式中括号外的下标为(i,k),则括号内m是mi,代表了编号为i的流体通道中的质量流量，c
p
可写作c
p,i,k
，是由温度节点t
i,k
和t
i,k+1
所定义的流体单元的比热容；若在流体通道中发生氢正仲催化转化，则式(29)中的m
e,i,k
需按照式(19)所介绍的流体通道中填充催化剂时的计算方法来计算，此时式(28)不封闭，需补充仲氢的组分输运方程：
[0223][0224]
若流体沿着+x方向流动，且流体通道中发生氢正仲催化转化，则根据式(28)、(30)可以沿着+x方向逐点确定流体温度节点中储存的仲氢浓度值，以由流体温度节点t
i,k
和
t
i,k+1
所定义的流体单元为例，具体步骤如下(已知和t
i,k
):
[0225]
a.预估出口仲氢浓度
[0226]
b.计算(计算所需的t
i,k+1
来源于前一个外迭代轮次所确定的流体温度场)，再根据式(26)、(27)计算
[0227]
c.根据式(28)、(30)计算和t
i,k+1
；
[0228]
d.比较新计算的和旧的若残差不满足收敛性条件，则继续从步骤b开始进行新一轮计算，直至的残差满足收敛性条件；
[0229]
当确定了和t
i,k+1
，可以继续确定下一个流体单元出口处的和t
i,k+2
，直至计算得到编号为i的流体通道的出口和t
i,n+1
；
[0230]
若流体通道沿着-x方向，计算流体单元温度节点的方程组是：
[0231][0232]
影响系数的计算方法是:
[0233][0234][0235]bi,k
＝b
i,k+1
+b
w,i,k
+b
w,i+1,k
ꢀꢀ
(32)
[0236]
而流体单元中仲氢的组分输运方程是：
[0237][0238]
通过前面介绍的一般迭代的方法，式(31)、(33)可确定流动是-x方向的流体通道中的温度场和仲氢浓度场分布；
[0239]
7)温度场残差和仲氢浓度场残差计算：
[0240]
将步骤5)计算得到的隔板温度场分布与前一个外迭代轮次确定的隔板温度场分布进行比较，计算隔板温度场残差，将步骤6)计算得到的流体温度场分布与前一个外迭代轮次确定的流体温度场分布进行比较，计算流体温度场残差，取隔板温度场残差和流体温度场残差的最大值作为温度场残差；取步骤6)计算得到的仲氢浓度场与前一个外迭代轮次确定的仲氢浓度场分布进行比较，计算仲氢浓度场残差；若温度场残差和仲氢浓度场残差有一个不满足收敛性条件，则继续重复步骤4)～7)，作为新一轮外迭代,直至两种残差都满足收敛性条件，即可输出计算结果，计算结果包括了板翅式换热器整体的温度场(隔板温度场和流体温度场)分布，流体通道中的压力场分布和仲氢浓度分布。
[0241]
本发明流程见图5，步骤4)至步骤7)是一个完整的外迭代轮次，其中的步骤4)是压力场的内迭代计算，步骤5)是采用结合高斯-赛德尔法的线迭代法来计算隔板中的温度场，
步骤6)是步进计算流体温度场和仲氢浓度场的分布，在步骤6)的每一个计算步内，需要先预估仲氢浓度，然后引入一般迭代法的思想来计算得到每个流体单元的出口仲氢浓度和出口温度。
[0242]
图6、图7和图8是本实施例的部分计算结果，图6是板翅式换热器其中三层流体通道中的流体温度分布，图7是计算得到的每一层流体通道的出口温度(一共117层)，以及a流体的出口仲氢浓度，图8是板翅式换热器第6层流体通道中的仲氢浓度和平衡氢浓度沿着流体方向的分布。