“双碳”背景下，氢能因其能量密度高、最终产物（水）无污染、资源广泛（工业上通过电解水制氢且地球上水资源丰富）等优点得到了人们的重点关注。质子交换膜燃料电池（PEMFC）是氢能的应用方式之一，具有零污染、转化过程噪音低、电池效率高的优点^［1-2］。由于电池处于变工况的工作环境下，通过供应过量氢气的方法来保证其输出功率的稳定，导致部分氢气未被完全消耗，因此需要将这些氢气回收利用。目前PEMFC常利用机械泵或引射器将燃料电池未消耗完的氢气重新循环到阳极进行二次利用。与机械泵相比，引射器无需消耗燃料电池的功率，噪声小，运维成本低，更适合在PEMFC氢气再循环系统中使用^［3］。

引射器的结构是根据特定工况参数设计，且燃料电池处于变工况的工作环境，因此引射器的工作性能并不稳定，会受到结构尺寸、操作条件的影响^［3–5］。由此，国内外学者通过改变引射器的结构尺寸来提高引射器的性能。尹燕等^［3］对3种不同尺寸的引射器进行建模、划分网格、模拟计算后，指出同一混合室直径下，随着扩散室长度的增大，最优扩散室角度减小，且该角度所对应的回流比增大。贠海涛等^［6］先设计得到引射器的结构尺寸，接着对4种不同尺寸的引射器进行数值模拟，研究结果表明引射器存在最佳的等压混合室收敛角、等容混合室长度。Metin等^［7］通过对3种不同尺寸的引射器进行对比研究，认为优化引射器喷嘴位置可使引射系数提高6%。Feng等^［8］设计了4种不同喉部直径的引射器方案，研究表明引射系数随着喷嘴喉部直径的减小而增大。

目前，PEMFC引射器设计大部分遵照结构参数计算、计算域建模、网格划分、数值模拟验证，然后迭代计算这一流程实现。如果需优化的结构参数有多个并且考虑这些参数之间的相互影响，那么设计方案的数量将会呈几何倍数增加，时间成本将是巨大的。因此，本文以二维PEMFC引射器为研究对象，通过Python编程语言编写计算引射器结构尺寸的程序、调用OpenFOAM软件中的blockMesh网格划分工具和rhoSimpleFoam求解器，来实现自动执行建模、网格划分、数值模拟仿真这一流程，建立一套自动仿真设计工具，从而显著降低设计引射器所需的时间。

1 引射器结构参数确定 1.1 引射器基本工作原理

图1、2分别为 PEMFC氢气循环流程图和引射器示意图。来自高压氢气罐的气体（工作气体）通过压力调节器达到适当压力后在喷嘴中被加速，经过喷嘴出口后会产生一个低压区，在压差的作用下来自燃料电池的未被消耗的氢气（引射气体）到达吸入腔，两者在混合室充分混合后进入扩压室，经减速扩压达到允许压力后通过电池阳极侧入口进入电池，重新参与反应^［9］。

通常用引射系数R_E来评价引射器的性能^［9］，它表示引射气体质量流量与工作气体质量流量之比。

$\qquad {R}_{{\mathrm{E}}}=\dfrac{{m}_{{\mathrm{s}}}}{{m}_{{\mathrm{p}}}} $

(1)

式中：m_s为引射气体质量流量，kg·s⁻¹；m_p为工作气体质量流量，kg·s⁻¹。

1.2 引射器结构尺寸计算流程

本文采用索科洛夫引射器设计方法^［9］，利用能量守恒、动量守恒和质量守恒定理并结合经验公式对PEMFC引射器进行结构设计，流程如下：

（1）由已知参数工作气体压力p_p、工作气体温度t_p、引射气体压力p_h、引射气体温度t_h、引射器背压p_c、气体绝热指数k，先计算出相应的气体动力函数，再计算出引射器所能达到的最大引射系数$ {R}_{{\mathrm{E{\text{，}}max}}} $。

$\qquad {R}_{{{{\mathrm{E}}{\text{，}}{\mathrm{max}}}}}=\dfrac{{K}_{1}{\lambda }_{{\mathrm{ph}}}-{K}_{3}{\lambda }_{{{\mathrm{c}}}_{3}}}{\sqrt{\theta }\left({K}_{4}{\lambda }_{{{\mathrm{c}}}_{3}}-{K}_{2}{\lambda }_{{{\mathrm{h}}}_{2}}\right)} $

(2)

其中

$\qquad {K}_{1}={\varphi }_{1}{\varphi }_{2}{\varphi }_{3} $

(3)

$\qquad {K}_{2}={\varphi }_{2}{\varphi }_{3}{\varphi }_{4} $

(4)

$ \qquad {K}_{3}=1 + {\varphi }_{3}\dfrac{{p}_{{\mathrm{c}}}}{{p}_{{\mathrm{p}}}}\frac{{\varPi }_{{\mathrm{c}}3}-\dfrac{{p}_{{\mathrm{h}}}}{{p}_{{\mathrm{c}}}}}{k{\varPi }^{*}{\lambda }_{{\mathrm{c}}3}{q}_{{\mathrm{ph}}}} $

(5)

$ \qquad {K}_{4}=1 + {\varphi }_{3}\dfrac{{p}_{{\mathrm{c}}}}{{p}_{{\mathrm{h}}}}\dfrac{{\varPi }_{{\mathrm{c}}3}-{\varPi }_{{\mathrm{c}}2}}{k{\varPi }^{*}{\lambda }_{{\mathrm{c}}3}{q}_{{\mathrm{h}}2}} $

(6)

式中：$ {\lambda }_{{\mathrm{ph}}} $为工作气体在喷嘴处的折算等熵速度；$ {\lambda }_{{{\mathrm{c}}}_{3}} $为混合气体在混合室出口截面上的折算等熵速度；$ {\lambda }_{{{\mathrm{h}}}_{2}} $为引射气体在混合室入口处的折算等熵速度；$ \theta $为引射气体温度与工作气体温度之比；$ {K}_{1} $、$ {K}_{2} $、$ {K}_{3} $、$ {K}_{4} $为系数；$ {\varphi }_{1} $为喷嘴的速度系数，取0.95；$ {\varphi }_{2} $为混合室的速度系数，取0.975；$ {\varphi }_{3} $为扩压室的速度系数，取0.9；$ {\varphi }_{4} $为引射气体进入混合室之前的速度系数，取0.925；$ {\varPi }_{{\mathrm{c}}2} $为混合流体在混合室入口截面上的相对压力；$ {\varPi }_{{\mathrm{c}}3} $为混合气体在混合室出口截面上的相对压力；$ {\varPi }^{*} $为临界压比；p_c为引射器出口压力，Pa；q_ph为工作流体在喷嘴出口截面上的折算质量速度；q_h2为引射流体在混合室入口截面上的折算质量速度。

（2）根据所确定的最大引射系数$ {R}_{{\mathrm{E{\text{，}}max}}} $及其对应的气体动力函数$ {\lambda }_{{{\mathrm{c}}}_{3}} $、$ {q}_{{{\mathrm{h}}}_{2}} $和经验公式得出引射器的结构数据。

$\qquad {f}_{{p}^{*}}=\frac{{G}_{p} {a}_{{p}^{*}}}{k{\varPi }^{*}{p}_{p}} $

(7)

$\qquad {f}_{3}=\frac{{f}_{{p}^{*}}{p}_{{{p}}}\left(1 + {R}_{{\mathrm{E{\text{，}}max}}}\sqrt{\theta }\right)}{{p}_{{\mathrm{c}}} {q}_{{\mathrm{c}}3}}$

(8)

$ \qquad {l}_{{\mathrm{h}}}=\left(6 \sim 10\right){d}_{3} $

(9)

式中：$ {f}_{{p}^{*}} $为喷嘴喉部面积，$ {\mathrm{m}\mathrm{m}}^{2} $；$ {f}_{3} $为混合室出口面积，$ {\mathrm{m}\mathrm{m}}^{2} $；$ {l}_{{\mathrm{h}}} $为混合室长度，一般取6～10倍的$ {d}_{3} $，$ {d}_{3} $为混合室出口直径，$ \mathrm{m}\mathrm{m} $；$ {a}_{{p}^{*}} $为工作气体临界速度，$ {\mathrm{m}}\cdot {\mathrm{s}}^{-1} $；$ {q}_{{\mathrm{c}}3} $为混合气体在混合室出口截面上的折算质量速度；q_h3为引射流体在混合室出口截面上的折算质量速度；G_p为工作流体入口处的质量流量，kg·s⁻¹。

2 数值模拟仿真

基于1.2节中引射器结构尺寸计算流程计算得到引射器的结构参数，采用计算流动力学（CFD）仿真软件OpenFOAM中的blockMesh工具、rhoSimpleFoam求解器对引射器流场进行计算域建模、网格划分、数值模拟。blockMesh是CFD仿真软件OpenFOAM中的网格生成工具之一，它依赖于字典文件 blockMeshDict。该文件主要由5个列表组成，分别为vertices、edges、blocks、boundary和mergePatchPairs。

2.1 计算域建模

vertices列表用于定义几何区域的顶点坐标，是构建复杂几何形状和划分网格的基础数之一。对于二维引射器，可通过vertices描述二维模型实现建模功能。图3(a)为vertices列表局部信息示例。

2.2 网格划分

列表blocks的作用是通过对列表vertices中的点进行有序定义，将计算域分解成多个三维的六面体块，并指定块方向上的网格数量、变化率以及网格的疏密程度。

值得注意的是，在OpenFOAM软件中即使进行二维数值模拟，网格也必须划分成三维，只需在定义边界条件时将无关方向的两个面略加更改即可。通过修改列表blocks可实现对网格划分信息的修改。图3(b)为blocks列表局部信息示例。

2.3 定义边界条件

列表boundary主要是对列表blocks中所形成的六面体块中各个面进行边界条件的定义。通过修改列表boundary中的信息实现定义边界条件的功能。图3（c）为boundary列表局部信息示例。

对vertices、blocks、boundary列表进行修改后，通过blockMesh命令会生成一个长1 m、宽0.1 m、高1 m的长方体。该长方体与长对应的四条边会被均匀地划分为若干个网格，与宽对应的四条边上会被划分为1个网格，与高对应的四条边会被均匀地划分为若干个网格。

2.4 求解器

OpenFOAM软件中有许多自带的求解器，例如：icoFoam求解器，可用于模拟层流流动；rhoSimpleFoam求解器，主要用于求解可压缩流动以及求解亚音速或超音速流动；interFoam求解器，可用于求解多相流的界面问题。用户还可以根据自身的需求编写相应功能的求解器。针对引射器中的湍流、超音速、流体的压缩性等流动问题，本文采用rhoSimpleFoam求解器进行求解。图4为基于Python的引射器自动设计流程。

2.5 数值仿真验证

为验证数值仿真程序的可靠性，采用经典二维圆柱绕流案例进行验证，计算雷诺数Re=200，圆柱直径D=1m，计算域为长20D，宽15D的矩形，并对部分区域进行加密，圆柱周围添加边界。图5为Re=200的圆柱绕流数值计算网格，并将计算结果与其他文献结果进行对比。

2.6 计算结果

表1为二维圆柱绕流数值验证的结果。将数值模拟得到的平均阻力系数$ {\overline{C}}_{{\mathrm{D}}} $、斯特劳哈尔数St与文献[10–12]中的数据进行对比发现，数值模拟结果与文献中的结论基本一致。由于模型尺寸、网格质量、求解器的差异会引起一定的误差，但误差在允许范围之内，由此可见本程序具有较高的可靠性。

当流体流经圆柱体时，在特定条件下会出现不稳定的边界层分离现象，使流体从圆柱体两侧剥离，并在圆柱体下游两侧产生两道非对称排列、交替的旋涡。图6为圆柱绕流瞬时涡量场。

3 引射器自动设计实例

结合前文讨论，以给定的一组引射器工况参数作为初始条件，然后利用本文开发的程序，按照图4所示流程，自动完成结构参数计算、计算域建模、网格划分、数值模拟等步骤。

3.1 计算引射器结构参数

表2为引射器设计工况点参数，其中：p为压力；T为温度。输入设计工况点参数后，程序自动计算得到的引射器主要结构数据如表3所示。

3.2 计算域建模

根据获得的引射器结构数据，程序自动将其转换成点的坐标，并写入vertices列表中，生成引射器计算域模型。图7为二维引射器计算域模型。

3.3 网格划分

通过Python软件向blocks列表中写入计算域所需分成的六面体块信息，并指定网格数量和网格变化率以达到加密网格的目的；接着调用blockMesh命令生成网格。图8为二维引射器喷嘴处局部网格。

3.4 数值模拟计算

通过Python软件调用rhoSimpleFoam求解器进行模拟计算，得到引射器内部流场结构。图9、10分别为引射器内部速度云图、引射器扩压室流线图。

图9中工作气体经喷嘴加速进入混合室后出现偏离混合室中轴线的现象。这是因为引射气体从下方入口端被吸入混合室时，会在混合室入口处（图9中虚线框出部分）向上挤压工作气体导致工作气体偏离混合室的中轴线。

图10中引射器扩压室出现明显的非对称性流动，这可能是由以下几个因素造成：①非对称现象可能与引射器的几何结构有关^［13］。②由于引射器出口压力过大，不利于引射过程，会在扩压室中出现非对称现象^［14］。③混合气体从超音速到亚音速的减速过程会导致边界层内的气流分离，超音速射流上的高剪切应力会导致其波动 ^[15］。

4 结论

本文针对质子交换膜燃料电池引射器设计时所面临的耗时问题，利用Python软件开发了一套PEMFC引射器自动仿真设计工具，该工具可自动计算引射器结构参数，并进行计算域建模、网格划分以及数值模拟。

针对经典二维圆柱绕流问题对求解器进行了数值验证，结果与其他文献结果基本一致，验证了该求解器的可靠性。给出了利用该工具进行引射器设计的一般流程和应用实例，展示了该工具的应用效果。此外，未来可进一步提高网格质量，提高数值模拟精度，并将优化过程集成到Python程序中，以进一步提升该设计工具的性能。