技术领域
[0001] 本
发明涉及
电解池模型的技术领域，特别涉及用于质子交换膜电解池特性分析的二维稳态模型建立方法。
背景技术
[0002] 虽然
潮汐能、
风能和
太阳能等
可再生能源是有前途的能源，但是它们受间歇性和地域性影响，并不是十分可靠。为了使
可再生能源技术得到广泛而可靠的应用，迫切需要清洁和可持续的能源技术来解决严重的环境问题并满足人类的需求。
[0003] 氢气是可再生能源存储的有前途的能源载体，过量的可再生能源可用于驱动电解池产生氢气，也可当再生能源不足时，通过
燃料电池将其转化为
电能。此外，氢气是
燃料电池汽车实现低排放和智能运输的理想燃料。
[0004] 质子交换膜电解池(PEMEC)是一种低温电化学电解池，是电解
水制氢气最广泛的一种方法。但是由于
电解质膜需要高水含量来维持膜的高质子传导性，因此除非将系统加压以维持膜的
含水量，否则操作
温度通常低于100℃。然而，在低于100℃的温度下输入PEMEC的
能量是电，并且
热能的贡献非常低。更重要的是，
电极反应迟钝需要使用昂贵的催化剂，例如Pt，这使得PEMEC非常昂贵。随着可替代电解质膜的发展，PEMEC可以在高于100℃的温度下运行，这种高温条件对于制氢是非常理想的。首先，电极动
力学随着温度的升高而增加，从而降低了电极的活化损失并使得能够使用成本更低的催化剂。其次，在130℃下进行水
蒸汽电解时，水分解需要的热量(ΔH)为243kJ mol-1，低于在80℃进行水分解需要的热量(284kJ mol-1)。随后，通过将电解池的温度从80℃(1.18V)升高到130℃(1.16V)，可以稍微降低电解池的可逆
电压。除此之外，高温质子交换膜电解池需要更多的热能输入，则表明更多的工业废热可用于高温质子交换膜电解池的制氢。
[0005] 而在工程实践上，实际的电解池流道通常较长，由于电解池运行时反应气体的消耗，沿流道方向必然存在反应气浓度的下降，从而导致催化层反应速率的变化。因此，对于长流道电解池中的特性分析往往不够精准，影响实际工程上的应用。
发明内容
[0006] 本发明的目的是提供用于质子交换膜电解池特性分析的二维稳态模型建立方法，以解决长流道的电解池特性分析不精准的问题。
[0007] 包括以下步骤：
[0008] S1，建立一个质子交换膜电解池的物理模型，提供电解池各结构的物理参数，以模拟电化学反应；
[0009] S2，根据所述物理模型建立用于计算电化学反应中相关稳态参数的计算模型；
[0010] S3，向所述物理模型内输入边界条件，以启动电化学反应，所述物理模型沿其流道方向被划分为多个网格层，每个网格层内被划分为多个网格
块，不同网格块内的反应物浓度不同，一个网格块内的反应物浓度视为均一；
[0011] S4，所述计算模型将边界条件作为初始值，使用有限元的
迭代方式进行计算，根据不同网格块内的反应物浓度逐个得到所有网格块内与边界条件对应的稳态参数值。
[0012] 上述技术方案中，将物理模型划分为网格状后，增加了沿流道方向的实际反应物浓度的分布考虑，根据每个网格块内的反应度浓度来获得该网格块内的稳态参数值，能够更好地反映电解池运行的局部特性，使得模型物理过程更加完整，提高了准确性和应用价值，有助于电堆及系统级电解池模型的开发，推动电解池实际应用领域的发展，并实现对电解池运行情况进行预测。
[0013] 进一步地，所述S1建立的物理模型的流道内各结构从左至右依次包括
阴极、气体扩散电极、催化层、电解质膜、催化层、气体扩散层、
阳极，流道的上端为气体入口，下端为气体出口；
[0014] 所述物理模型提供的物理参数包括所述催化层的孔隙率、所述气体扩散层的孔隙率、所述阴极与阳极的电极渗透率、所述气体扩散层的渗透率；所述S4中计算模型计算稳态参数值时根据不同网格块所处的
位置分别使用对应的物理参数。
[0015] 进一步地，所述网格块呈矩形，所述网格块的
密度从所述阴极与阳极中间分别向两侧逐渐减小。
[0016] 进一步地，所述S3中输入的边界条件包括操作电压、阳极反应物的
质量分数、阳极气体流速、阴极反应物的质量分数、阴极气体流速、操作压强、操作温度。
[0017] 进一步地，所述S4迭代计算结束后，得到的稳态参数值包括有与操作电压对应的
电流密度值、与操作压强对应的压强值、与阳极反应物的质量分数和阴极反应物的物质分数对应的摩尔分数值、与操作温度对应的温度值。
[0018] 进一步地，所述S2中建立的计算模型包括质量与动量传递模块、电化学反应模块以及热传递模块，所述S4中计算模型对各稳态参数的计算顺序依次为：
[0019] 用于质子交换膜电解池特性分析的二维稳态模型建立方法，所述质量与动量传递模块先逐网格块地迭代出所有网格块的压强值；
[0020] b，所述质量与动量传递模块根据各网格块中的压强值逐网格块地迭代出所有网格块的摩尔分数值，以获得各网格块内的反应物浓度；
[0021] c，所述电化学模块根据各网格块中的反应物浓度沿流道方向逐网格层地迭代出所有截面内的电流密度值；
[0022] d，所述热传递模块根据电流密度值获取各网格层内的热量变化，再逐网格块地迭代出所有网格块的温度值。
[0023] 进一步地，所述S3中不同网格块内的反应物浓度通过所述计算模型计算出的摩尔分数值获得，反应过程中反应物气体从前一个网格块进入到当前网格块的传递过程通过所述质量与动量传递模块中的公式来表达：
[0024]
[0025] 公式中，Ni表示物质运输的通量，R表示通用气体常数，T表示温度值，P表示压强，B0为多孔电极的渗透率，μ表示气体
粘度，yi是组分i的摩尔分数；
[0026] 是组分i的总有效扩散系数，通过下述公式获得：
[0027]
[0028] 其中， 为组分i的克努森扩散系数， 为组分i的分子扩散系数；
[0029]
[0030] 公式中，ci是组分i的摩尔浓度，与组分i的摩尔分数有关，Ri为组分i的质量源项。
附图说明
[0031] 图1为质子交换膜电解池的物理模型结构示意图一；
[0032] 图2为质子交换膜电解池的物理模型结构示意图二；
[0033] 图3为物理模型中网格块的划分示意图一；
[0034] 图4为物理模型中网格块的划分示意图二；
[0035] 图5为电解池阴极侧压强沿流道方向的分布图；
[0036] 图6为电解池阳极侧压强沿流道方向的分布图；
[0037] 图7为电解池阳极侧水的摩尔分数沿流道方向的分布图；
[0038] 图8为电解池阴极侧水的摩尔分数沿流道方向的分布图；
[0039] 图9为电解池内电流密度值沿流道方向的分布图；
[0040] 图10为电解池内温度沿流道方向的分布图。
具体实施方式
[0041] 为了使本发明要解决的技术问题、技术方案及有益效果更为清楚明白，以下结合附图及
实施例，对本发明进行进一步详细说明。应当理解此处描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0042] 一、建立物理模型
[0043] 参照图1和图2，通过COMSOL
软件搭建，电解池流道呈长方形，从左至右依次包括阴极、气体扩散电极、催化层、电解质膜、催化层、气体扩散层、阳极，电解池流道的上端为气体入口，下端为气体出口。
[0044] 电解池长20mm，气体流道高度1mm，气体扩散层厚度0.38mm，催化层厚度0.05mm，电解质膜厚度0.1mm，催化层孔隙率为0.3，气体扩散层孔隙率为0.4，电极渗透率为2.36×10-12m2，气体扩散层渗透率为1.18×10-11m2。物理模型提供各个结构的物理参数，包括催化层的孔隙率、气体扩散层的孔隙率、阴极与阳极的电极渗透率、气体扩散层的渗透率，用于模拟电化学反应。
[0045] 向物理模型中输入边界条件，边界条件包括操作电压、阳极反应物的质量分数、阳极气体流速、阴极反应物的质量分数、阴极气体流速、操作压强、操作温度。输入边界条件后，可以启动电化学反应。
[0046] 二、建立计算模型
[0047] 本实施例公开了一种质子交换膜电解池的计算模型，基于实施例一的物理模型，通过COMSOL软件搭建计算模型，向计算模型内输入边界条件并设定相对容差，计算模型根据输入的边界条件，进行有限元的迭代计算，直到相邻两次计算的结果之间的差值小于等于相对容差时，停止迭代，最终得到与操作电压对应的电流密度值、与操作压强对应的压强值、与阳极气体流速和阴极气体流速对应的气体速度值、与阳极反应物的质量分数和阴极反应物的物质分数对应的摩尔分数值、与操作温度对应的温度值。
[0048] 通过计算模型获得的上述各参数值均为电解池的稳态参数值，计算模型中包含有计算单元和判断单元，首先计算单元把输入的边界条件作为初始值，以开始有限元的迭代计算，每一次迭代均得到与初始值对应的结果值，然后再将结果值作为下一次迭代计算的初始值。
[0049] 判断单元则用于保存预先交互设定的相对容差，并获取每一次迭代的结果值，与前一次迭代的结果值进行比较，直到相邻两次迭代的结果值之间的差值小于等于相对容差时，输出停止指令给计算模块，从而实现停止计算，最后一次计算得到的结果值则为电解池的稳态参数值。
[0050] 计算模型包括三个计算模块：质量与动量传递模块、电化学反应模块以及热传递模块。
[0051] 质量与动量传递模块包括以下公式：
[0052]
[0053] 公式中，Ni表示物质运输的通量，P表示压强值，B0为多孔电极的物理参数渗透率，与位点有关，μ表示气体粘度，yi是组分i的摩尔分数值；
[0054] 是组分i的总有效扩散系数，通过下述公式获得：
[0055]
[0056] 为组分i的克努森扩散系数， 为组分i的分子扩散系数；
[0057]
[0058] 公式中，ci是组分i的摩尔浓度，即反应物浓度，与组分i的摩尔分数值有关，Ri为组分i的质量源项。
[0059]
[0060] ε是计算位置处的孔隙率，τ表示弯曲系数，ρ为气体密度，u为气体速度，P表示压强值，μ为气体粘度，T表示矩阵转置。
[0061] 电化学反应模块包括以下公式：
[0062] V＝E+ηact,an+ηact,ca+ηohmic；……(5)
[0063] V表示交互输入的操作电压，E表示电解池在当前运行条件下的平衡电压，ηact,an表示阳极的活化过电势，ηact,ca表示阴极的活化过电势，ηohmic表示由质子和
电子传导引起的欧姆过电势；
[0064]
[0065] 是标准状态下的平衡电压，R是通用气体常数，T表示电解池的操作温度，F为法拉第常数， 和 分别表示在不同位点上H2、H2O和O2的局部分压，与质量与动量传递模块获得压强值有关；
[0066] 阳极的活化过电势和阴极的活化过电势均通过以下公式得到：
[0067]
[0068] i表示工作电流密度，i0表示交换电流密度，α是电荷传输系数，n是每摩尔电化学反应转移的电子数，γ为指数前因子，Eact表示活化能；
[0069] 欧姆过电势通过欧姆定律获得：
[0070]
[0071] 表示质子导电率，Φs表示质子电势，il为电流密度值， 表示质子导电率，Φl表示质子电势；其中， 随迭代过程中的温度值变化，Φl随交互设定的操作电压变化。
[0072] 热传递模块包括以下公式：
[0073]
[0074] T表示温度值，ρ表示密度；Cp为
流体热容；u为质量与动量传递模块得到的气体速度值，λeff是有效导热系数；Q是热源项，表示电化学反应或过电压损失而消耗或产生的热量。
[0075] λeff＝(1-ε)λs+ελl；……(10)
[0076] λeff是有效导热系数，λs表示固体相导热系数；λl表示液体相导热系数，ε是物理参数孔隙率，与位点有关，不同位点的有效导热系数，得到不同位点处的温度值。
[0077] 以下为电流密度值、压强值、气体速度值、摩尔分数值、温度值的计算原理：
[0078] 参照图3和图4，物理模型沿其流道方向被划分成多个网格层，每个网格层内划分为多个网格块，不同网格块内的反应物浓度不同，一个网格块内的反应物浓度视为均一。网格块呈矩形，网格块的密度从阴极与阳极的中间开始、分别向阴极与阳极两侧逐渐减小。
[0079] a，先从第一个网格块开始、逐网格块地算出所有网格块的压强值与气体速度值，单个网格块内的压强值与气体速度值通过公式(1)至(4)进行有限元的迭代计算获得。由于不同的网格块位于电解池中的不同位置，因此针对不同位置，使用该位置处的物理参数，逐个遍历所有网格块后，得到电解池内压强与气体速度沿流道方向的分布。
[0080] b，反应过程中反应物气体从前一个网格块进入到当前网格块的传递过程通过公式(1)至(3)来表达，因此，根据各网格块中的压强值，可逐网格块地算出所有网格块内的摩尔分数值，从而得到各网格块内的反应物浓度。由于不同的网格块位于电解池中的不同位置，因此针对不同位置，使用该位置处的物理参数，逐个遍历所有网格块后，得到电解池内反应物浓度沿流道方向的分布。
[0081] c，根据各网格块中的反应物浓度沿流道方向逐网格层地迭代出所有截面内的电流密度值，单个网格层内的电流密度值通过公式(5)至(8)进行有限元的迭代计算获得，逐个遍历所有网格层后，得到电解池内电流密度值沿流道的分布。
[0082] d，根据电流密度值得到各网格层内的热量变化，即获得公式(9)中的热源项，然后根据热源项逐网格块地迭代出所有网格块的温度值。由于不同的网格块位于电解池中的不同位置，因此针对不同位置，使用该位置处的物理参数，逐个遍历所有网格块后，得到电解池内温度沿流道方向的分布。
[0083] 三、用于电解池特性分析
[0084] 包括以下步骤：
[0085] S1，建立一中公开的物理模型，提供电解池中各结构的物理参数，包括催化层的孔隙率、气体扩散层的孔隙率、阴极与阳极的电极渗透率、气体扩散层的渗透率。
[0086] S2，根据S1中的物理模型建立二中公开的计算模型。
[0087] S3，向S1的物理模型内输入边界条件：操作电压为1.4V、阳极入口处反应物的质量分数为1、阳极入口处气体流速为0.1m/s、阴极入口处反应物的质量分数为1、阴极入口处气体流速为0.4m/s、操作压强为1atm、操作温度为403.15k。以启动电化学反应。其中，反应物为水，设定相对容差为0.001。
[0088] S4，将S1中输入的边界条件作为初始值，通过实施例二公开的计算模型进行计算，将各网格块的数据收集
整理后，分别获得如图5和图6所示的电解池内压强沿流道方向的分布图；如图7所示的电解池阳极侧反应物摩尔分数沿流道方向的分布图；如图8所示的电解池阴极侧反应物摩尔分数沿流道方向的分布图；如图9所示的电解池内电流密度值沿流道方向的分布图；如图10所示的电解池内温度沿流道方向的分布图。
[0089] 通过该方法，实现了对不同工作状态下的质子交换膜电解池性能和沿流道方向的各稳态参数分布进行求解。由于电解池运行时反应气体的消耗，沿流道方向必然存在反应气浓度的下降，从而导致催化层反应速率的变化，表现为电解池沿流道方向各
节点电流密度不同。在工程实践上，实际的单电解池流道为长流道，流道长度及流道结构对电池的影响往往很显著。而本方法则考虑了沿流道方向的实际反应气体浓度分布，提高了二维稳态模型的准确性和应用价值，特别有助于系统级电解池模型的开发和电解池动态模型研究，推动质子交换膜电解池实际应用领域的发展，并实现对电解池运行情况进行预测。
[0090] 以上仅为本发明的若干个优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
