CN113745591A - 具有过冷却水结冰机理的燃料电池低温冷启动建模方法 - Google Patents

Abstract

本发明提供了一种具有过冷却水结冰机理的燃料电池低温冷启动建模方法，包括以下步骤：S1、分别对膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程以及模型输出电压性能公式进行求解；S2、基于步骤S1中的各个方程建立具有过冷却水结冰机理的燃料电池低温冷启动模型。本发明所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法计算了电池内部的“水‑热‑气‑电”耦合传输机理，过冷却水结冰机理与阴极流道中结冰现象，能够探究单电池结构设计参数以及环境因素的影响规律，并且给出单电池内部的参数分布情况。

Description

具有过冷却水结冰机理的燃料电池低温冷启动建模方法

技术领域

本发明属于质子交换膜燃料电池领域，尤其是涉及具有过冷却水结冰机理的燃料电池低温冷启动建模方法。

背景技术

质子交换膜燃料电池(PEMFC)是氢能产业发展的核心技术之一，具有高功率密度、高能量转换效率、零排放、低噪声等优点，在车用动力源、固定电站以及移动电源的领域具有重要应用前景。

零下温度启动问题是PEMFC实现大规模商业化之前需要解决的技术难题之一。当PEMFC从零下温度启动时，氢氧电化学反应生成的水会结冰，催化层(CL)中的冰会覆盖三相反应界面，导致有效反应面积的减小，从而造成输出性能的下降。气体扩散层(GDL)中的冰会占据气体传输通道，阻碍反应气体从流道往催化层中传输。当催化层中反应界面完全被冰覆盖或者反应气体无法抵达催化层中的反应界面时，低温启动过程被视为失败。由于冰的密度低于水的密度，水结冰时体积会膨胀，因而结冰/融化循环会对电池的微观结构造成破坏，从而严重影响电池的耐久性与寿命。

根据美国能源部(DOE)提出的80kW车载燃料电池动力系统的最终技术目标，系统需要成功实现-30℃的无辅助措施启动以及-40℃的辅助加热启动。此外，从-20℃环境温度启动时，燃料电池系统需要在30秒以内达到额定功率的50％。关于启动/停机的耐久性指标，在DOE设定的特定测试工况下，要求经过5000次循环后输出性能的衰减应小于5％。人们已经针对冷启动问题开展了实验与仿真研究，实验研究聚焦于操作工况与结构设计参数的影响，仿真模型包括传热模型、解析模型、一维、二维和三维模型，但是上述技术指标目前仍然难以实现。由于透明可视化实验装置与红外成像技术观察到了冷启动过程中燃料电池内催化层表面与流道中出现的过冷却水水滴，而且高分辨率中子成像技术同样证实了过冷却水形态的存在，而已有的仿真研究却缺乏过冷却水形态的存在与阴极气体流道中的结冰情况，因此极大降低了模拟结果可靠性与准确性，我们有必要加入上述现象搭建更全面的冷启动仿真模型，以满足燃料电池低温冷启动技术开发需求。

发明内容

有鉴于此，本发明旨在提出具有过冷却水结冰机理的燃料电池低温冷启动建模方法，以弥补已有冷启动模型缺乏过冷却水结冰机理与流道中的结冰现象而导致模型准确度低的短板，解决三维模型仿真资源需求高且计算效率低的问题，能够为低温冷启动技术开发提供仿真数据支持，极大的减小实验成本及研发周期。

为达到上述目的，本发明的技术方案是这样实现的：

具有过冷却水结冰机理的燃料电池低温冷启动建模方法，包括以下步骤：

S1、分别对膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程以及模型输出电压性能公式进行求解；

S2、基于步骤S1中的各个方程以及过冷却水结冰机理建立具有过冷却水结冰机理的燃料电池低温冷启动模型。

进一步的，在步骤S1中的所述能量守恒方程分别与膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程之间数据耦合，所述液压守恒方程分别与冰守恒方程、水蒸气守恒方程之间数据耦合，所述冰守恒方程与水蒸气守恒方程之间数据耦合，所述膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程均控制耦合电池输出电压性能公式；

所述过冷却水结冰机理分别与膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程之间数据耦合，膜态水守恒方程中膜态水源项、液压守恒方程中过冷却水源项、冰守恒方程中冰源项、水蒸气守恒方程中水蒸气源项、能量守恒方程中热源项的计算均基于过冷却水结冰机理进行计算；

所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程、过冷却水结冰机理、输出电压性能公式组成燃料电池低温冷启动模型。

进一步的，在步骤S1中的所述膜态水守恒方程如下：

式中ρ_MEM表示质子交换膜密度，EW表示质子交换膜当量质量，ω表示聚合物体积分数，λ表示膜态水含量，其物理定义为每个质子传导团簇所携带的平均水分子数，t表示时间，

表示有效膜态水扩散系数，S_mw表示膜态水源项；

所述S_mw的计算表达式如下：

式中S_mw表示膜态水源项，S_react表示氢氧电化学反应产物水，S_m-v表示膜态水与水蒸气之间的相变转化，S_m-l表示膜态水与过冷却水之间的相变转化，S_EOD表示电拖拽效应，S_per表示水力渗透作用；

所述有效膜态水扩散系数的参数值基于能量守恒方程的温度参数值进行计算；

所述膜态水源项的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算。

进一步的，所述过冷却水结冰机理包括以下步骤：

A1、判断催化层中膜态水含量是否达到饱和，若是则进行下一步，否则直接进入步骤A3；

A2、催化层中膜态水含量转变为过冷却水与水蒸气；

A3、催化层中膜态水与水蒸气之间相互转化；

A4、判断催化层、微孔层、气体扩散层与流道中水蒸气是否达到饱和浓度，若是则进行下一步，否则进入步骤A6；

A5、催化层、微孔层、气体扩散层与流道中水蒸气达液化为过冷却水，并进行下一步；

A6、过冷却水从催化层往流道传输，并同时发生结冰现象。进一步的，在步骤S1中的所述液压守恒方程的求解操作包括以下步骤：

B1、对液压守恒方程进行求解；

B2、当液压守恒方程得到求解后，对过冷却水体积分数公式以及水力渗透公式计算。

进一步的，在步骤B1中的所述液压守恒方程如下：

式中ε表示孔隙率，ρ_lq表示过冷却水密度，s_lq表示过冷却水体积分数，t表示时间，表示过冷却水渗透率，μ_lq表示动力粘度，p_l表示液压，S_lq表示过冷却水源项，u_lq表示过冷却水流速；

所述S_lq的计算表达式如下：

式中S_lq表示过冷却水源项，S_m-l表示膜态水与过冷却水之间的相变转化，S_v-l表示过冷却水与水蒸气之间的相变转化，S_l-i表示过冷却水与冰之间的相变转化，S_per表示水力渗透作用，

表示水的摩尔质量；

所述过冷却水源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中液压参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算；

在步骤B2中的所述过冷却水体积分数公式如下：

p_c＝p_g-p_l (5)；

式中p_g表示气压，p_c表示毛细压力，σ表示表面张力系数，θ表示接触角，ε表示孔隙率，K⁰表示固有渗透率，s_lq表示过冷却水体积分数；

在步骤B2中的所述水力渗透公式如下：

式中S_per表示水力渗透作用，ρ_lq表示过冷却水密度，K_per表示水力渗透系数，p_l，CLc表示阴极催化层的液压，p_l，CLa表示阳极催化层的液压，μ_lq表示动力粘度，δ_CL表示催化层厚度，δ_MEM表示质子交换膜厚度，

表示水的摩尔质量。

进一步的，在步骤S1中的所述冰守恒方程的求解操作包括以下步骤：

C1、对于催化层、微孔层与气体扩散层内的实际结冰点温度进行修正操作；

C2、对冰守恒方程进行求解；

在步骤C1中催化层、微孔层与气体扩散层内的实际结冰点温度修正公式如下：

式中T_FPD表示结冰点温度修正数值，T_N表示正常结冰点温度，σ表示表面张力系数，θ表示接触角，ρ_ice表示冰密度，h_fusn表示融化潜热，r_{CL，MPL，GDL}分别表示催化层、微孔层和气体扩散层内孔径尺寸；

在步骤C2中的所述冰守恒方程如下：

式中s_ice表示冰体积分数，S_ice表示冰源项，ε表示孔隙率，ρ_ice表示冰密度，t表示时间；

所述冰源项的参数值基于液压守恒方程中过冷却水体积分数参数值、能量守恒方程中温度参数值进行计算。

进一步的，在步骤S1中的所述水蒸气守恒方程如下：

式中

分别表示t时刻催化层、微孔层、气体扩散层的水蒸气浓度，

分别表示t-Δt时刻催化层、微孔层、气体扩散层的水蒸气浓度，

分别表示催化层与微孔层之间、微孔层与气体扩散层之间的有效气体扩散系数，δ_CL表示催化层厚度，δ_MPL表示微孔层厚度，δ_GDL表示气体扩散层厚度，S_vp表示水蒸气源项，ε_CL表示催化层孔隙率，ε_MPL表示微孔层孔隙率，Δt表示时间步长；

所述催化层与微孔层之间、微孔层与气体扩散层之间有效气体扩散系数的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、能量守恒方程中温度参数值进行计算；

所述水蒸气源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、能量守恒方程中温度参数值进行计算。

进一步的，在步骤S1中的所述能量守恒方程如下：

式中

表示有效密度比热容，T表示温度，

表示有效导热率，S_T表示热源项，ASR_CL表示催化层的单位面积电阻，ASR表示单位面积电阻，δ_CL表示催化层厚度，δ表示厚度，I表示电流密度，V_act，a，V_act，c分别表示阳极与阴极的活化损失电压，ΔS表示熵变，F表示法拉第常数，S_pc表示相变潜热，t表示时间；

所述有效密度比热容的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算；

所述有效导热率的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算；

所述热源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算。

进一步的，在步骤S1中的所述模型输出电压性能公式如下：

V_out＝V_Nemst-V_act-V_ohmic-V_conc (16)；

V_ohmic＝I·ASR (19)；

式中V_out表示输出电压，V_Nernst表示能斯特电压，V_act表示活化损失电压，V_ohmic表示欧姆损失电压，V_conc表示浓差损失电压，R表示理想气体常数，T表示温度，

表示参考反应气体浓度，α_a，α_c分别表示阳极与阴极电荷传递系数，n_a，n_c分别表示阳极与阴极传递电子数，j_a，j_c分别表示阳极与阴极反应速率，p_a，p_c分别表示阳极与阴极气体压强，ST_a，ST_c分别表示阳极与阴极化学计量比，δ_CLa，δ_CLc分别表示阳极催化层与阴极催化层的厚度，I表示电流密度，ASR表示单位面积电阻，I_lim，a，I_lim，c分别表示阳极与阴极极限电流密度；

所述输出电压的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气浓度守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算。

相对于现有技术，本发明所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法具有以下优势：

(1)本发明所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，计算了电池内部的“水-热-气-电”耦合传输机理，过冷却水结冰机理与阴极流道中结冰现象，能够探究单电池结构设计参数(如多孔介质层憎水性、质子交换膜厚度、流道肋宽比)以及环境因素(如冷启动温度、换热系数、结冰速率)的影响规律，并且给出单电池内部的参数分布情况(如水含量、冰含量、水蒸气浓度、温度)。

(2)本发明所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法弥补了已有冷启动模型缺乏过冷却水结冰机理与流道中的结冰现象而导致模型准确度低的短板，同时解决了三维模型仿真资源需求高且计算效率低的问题，能够为低温冷启动技术开发提供仿真数据支持，极大的减小实验成本及研发周期。

附图说明

构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为本发明实施例所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法相变转化关系示意图；

图2为本发明实施例所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法其失败启动工况下模型仿真结果与实验数据验证图；

图3为本发明实施例所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法其成功启动工况下模型仿真结果与实验数据验证图；

图4为本发明实施例所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法其催化层中冰体积分数随时间的变化情况图；

图5为本发明实施例所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法其催化层中水体积分数随时间的变化情况图；

图6为本发明实施例所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法其电池平均温度随时间的变化情况图；

图7为本发明实施例所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法各个守恒方程之间的数据耦合原理图。

附图标记说明：

1、质子交换膜；2、催化层；3、微孔层。

具体实施方式

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明中的具体含义。

下面将参考附图并结合实施例来详细说明本发明。

如图1至图7所示，具有过冷却水结冰机理的燃料电池低温冷启动建模方法，包括以下步骤：

S1、分别对膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程以及模型输出电压性能公式进行求解；

S2、基于步骤S1中的各个方程以及过冷却水结冰机理建立具有过冷却水结冰机理的燃料电池低温冷启动模型。本建模方法充分计算电池内部的“水-热-气-电”耦合传输机理，过冷却水结冰机理与阴极流道中结冰现象，能够探究单电池结构设计参数(如多孔介质层憎水性、质子交换膜1厚度、流道肋宽比)以及环境因素如温度、换热系数、结冰速率的影响规律，并且详细的给出单电池内部的参数分布情况(如水含量、冰含量、水蒸气浓度、温度)。

在步骤S1中的所述能量守恒方程分别与膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程之间数据耦合，所述液压守恒方程分别与冰守恒方程、水蒸气守恒方程之间数据耦合，所述冰守恒方程与水蒸气守恒方程之间数据耦合，所述膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程均控制耦合电池输出电压性能公式；

所述过冷却水结冰机理分别与膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程之间数据耦合，膜态水守恒方程中膜态水源项、液压守恒方程中过冷却水源项、冰守恒方程中冰源项、水蒸气守恒方程中水蒸气源项、能量守恒方程中热源项的计算均基于过冷却水结冰机理进行计算；

所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程、过冷却水结冰机理、输出电压性能公式组成燃料电池低温冷启动模型。

在步骤S1中的所述膜态水守恒方程(在质子交换膜1与催化层2中求解)如下：

式中ρ_MEM表示质子交换膜1密度，EW表示质子交换膜1当量质量，ω表示聚合物体积分数，λ表示膜态水含量，其物理定义为每个质子传导团簇所携带的平均水分子数，t表示时间，

表示有效膜态水扩散系数，S_mw表示膜态水源项。

所述S_mw的计算表达式如下：

式中S_mw表示膜态水源项，S_react表示氢氧电化学反应产物水，S_m-v表示膜态水与水蒸气之间的相变转化，S_m-l表示膜态水与过冷却水之间的相变转化，S_EOD表示电拖拽效应，S_per表示水力渗透作用；

所述有效膜态水扩散系数的参数值基于能量守恒方程的温度参数值进行计算；

所述膜态水源项的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算。

所述过冷却水结冰机理包括以下步骤：

A1、判断催化层中膜态水含量是否达到饱和，若是则进行下一步，否则直接进入步骤A3；

A2、催化层中膜态水含量转变为过冷却水与水蒸气；

A3、催化层中膜态水与水蒸气之间相互转化；

A4、判断催化层、微孔层、气体扩散层与流道中水蒸气是否达到饱和浓度，若是则进行下一步，否则进入步骤A6；

A5、催化层、微孔层、气体扩散层与流道中水蒸气达液化为过冷却水，并进行下一步；

A6、过冷却水从催化层往流道传输，并同时发生结冰现象；所述过冷却水结冰机理如下，当催化层中膜态水含量达到饱和时，其转变为过冷却水与水蒸气。当水蒸气达到饱和浓度后，其液化为过冷却水。过冷却水从多孔介质中往流道中传输，与此同时，过冷却水发生结冰现象。具体的，本发明对于冷启动中的相变关系进行了相应完善并且计算了阴极气体流道中的结冰情况，如图1所示，当催化层2中膜态水含量达到饱和时，其转变为过冷却水形态与水蒸气形态。当水蒸气达到饱和水蒸气浓度后，其液化为过冷却水。过冷却水在液压梯度的作用下从多孔介质中往流道中传输，与此同时，过冷却水以一定的结冰速率发生结冰现象。

在步骤S1中的所述液压守恒方程的求解操作包括以下步骤：

B1、对液压守恒方程进行求解；在实际建模时，过冷却水体积分数基于多孔介质中液压连续进行求解。

B2、当液压守恒方程得到求解后，对过冷却水体积分数公式以及水力渗透公式计算。

在步骤B1中的所述液压守恒方程如下：

式中ε表示孔隙率，ρ_lq表示过冷却水密度，s_lq表示过冷却水体积分数，t表示时间，表示过冷却水渗透率，μ_lq表示动力粘度，p_l表示液压，S_lq表示过冷却水源项，u_lq表示过冷却水流速；在实际建模时，过冷却水体积分数基于多孔介质中液压连续进行求解。

所述S_lq的计算表达式如下：

式中S_lq表示过冷却水源项，S_m-l表示膜态水与过冷却水之间的相变转化，S_v-l表示过冷却水与水蒸气之间的相变转化，S_l-i表示过冷却水与冰之间的相变转化，S_per表示水力渗透作用，

表示水的摩尔质量；

所述过冷却水源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中液压参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算；

在步骤B2中的所述过冷却水体积分数公式如下：

p_c＝p_g-p_l (5)；

式中p_g表示气压，p_c表示毛细压力，σ表示表面张力系数，θ表示接触角，ε表示孔隙率，K⁰表示固有渗透率，s_lq表示过冷却水体积分数；

在步骤B2中的所述水力渗透公式如下：

当液压的数值得到求解后，则水力渗透作用计算如下：

式中S_per表示水力渗透作用，ρ_lq表示过冷却水密度，K_per表示水力渗透系数，p_l，CLc表示阴极催化层的液压，p_l，CLa表示阳极催化层的液压，μ_lq表示动力粘度，δ_CL表示催化层2厚度，δ_MEM表示质子交换膜1厚度，

表示水的摩尔质量。

在步骤S1中的所述冰守恒方程的求解操作包括以下步骤：

C1、对于催化层2、微孔层3与气体扩散层内的实际结冰点温度进行修正操作；

C2、对冰守恒方程进行求解；

由于多孔介质中结冰点温度比正常结冰点温度略低，因此，为了准确表征相变转换过程，本发明对于催化层2、微孔层3与气体扩散层内的实际结冰点温度进行了修正，在步骤C1中催化层2、微孔层3与气体扩散层内的实际结冰点温度修正公式如下：

式中T_FPD表示结冰点温度修正数值，T_N表示正常结冰点温度(273.15K)，σ表示表面张力系数，θ表示接触角，ρ_ice表示冰密度，h_fusn表示融化潜热，r_{CL，MPL，GDL}分别表示催化层2、微孔层3和气体扩散层内孔径尺寸；

在实际建模时，冰体积分数在多孔介质与流道中求解，由于冰是固态，故不涉及到流动现象，在步骤C2中的所述冰守恒方程如下：

式中s_ice表示冰体积分数，S_ice表示冰源项，ε表示孔隙率，ρ_ice表示冰密度，t表示时间；

所述冰源项的参数值基于液压守恒方程中过冷却水体积分数参数值、能量守恒方程中温度参数值进行计算。

在步骤S1中的所述水蒸气守恒方程如下：

本发明中，多孔介质与流道中水蒸气的浓度在每一个时间步进行更新，即采用显式格式更新算法，为提高计算效率，层内不再细分网格，水蒸气浓度在上述各层的中心处求解，离散化处理后计算表达式如下：

式中

分别表示t时刻催化层、微孔层、气体扩散层的水蒸气浓度，

分别表示t-Δt时刻催化层、微孔层、气体扩散层的水蒸气浓度，

分别表示催化层与微孔层之间、微孔层与气体扩散层之间的有效气体扩散系数，δ_CL表示催化层厚度，δ_MPL表示微孔层厚度，δ_GDL表示气体扩散层厚度，S_vp表示水蒸气源项，ε_CL表示催化层孔隙率，ε_MPL表示微孔层孔隙率，Δt表示时间步长；

所述催化层与微孔层之间、微孔层与气体扩散层之间有效气体扩散系数的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、能量守恒方程中温度参数值进行计算；

所述水蒸气源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、能量守恒方程中温度参数值进行计算。

在步骤S1中的所述能量守恒方程如下：

能量守恒方程在整个燃料电池中求解：

式中

表示有效密度比热容，T表示温度，

表示有效导热率，S_T表示热源项，ASR_CL表示催化层的单位面积电阻，ASR表示单位面积电阻，δ_CL表示催化层厚度，δ表示厚度，I表示电流密度，V_act，a，V_act，c分别表示阳极与阴极的活化损失电压，ΔS表示熵变，F表示法拉第常数，S_pc表示相变潜热，t表示时间；

所述有效密度比热容的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算；

所述有效导热率的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算；

所述热源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算。

在步骤S1中的所述模型输出电压性能公式如下：

模型输出电压通过塔菲尔公式进行计算：

V_our＝V_Nernst-V_act-V_ohmic-V_conc (16)；

V_ohmic＝I·ASR (19)；

式中V_out表示输出电压，V_Nernst表示能斯特电压，V_act表示活化损失电压，V_ohmic表示欧姆损失电压，V_conc表示浓差损失电压，R表示理想气体常数，T表示温度，

表示参考反应气体浓度，α_a，α_c分别表示阳极与阴极电荷传递系数，n_a，n_c分别表示阳极与阴极传递电子数，j_a，j_c分别表示阳极与阴极反应速率，p_a，p_c分别表示阳极与阴极气体压强，ST_a，ST_c分别表示阳极与阴极化学计量比，δ_CLa，δ_CLc分别表示阳极催化层与阴极催化层的厚度，I表示电流密度，ASR表示单位面积电阻，I_lim，a，I_lim，c分别表示阳极与阴极极限电流密度；

所述输出电压的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气浓度守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算。

实施例1

以下通过具体计算实施例1对发明的方法以及模型建立的具体步骤进行详细的说明。

单电池结构参数：

有效反应面积：25cm²。极板厚度：2mm；流道厚度：1mm；气体扩散层厚度：0.2mm；微孔层3厚度：0.02mm；催化层2厚度：0.01mm；膜厚度：0.03mm；流道长度：0.15m；流道与肋板宽的比值：1。

质子交换膜1密度：1980kg m^-3；催化层2密度：1000kg m^-3；微孔层3密度：1000kg m^-3；气体扩散层密度：1000kg m^-3；极板密度：1000kg m^-3。质子交换膜1比热容：833J kg^-1 K^-1；催化层2比热容：3300J kg^-1 K^-1；微孔层3比热容：568J kg^-1 K^-1；气体扩散层比热容：2000J kg^-1 K^-1；极板比热容：1580J kg^-1 K^-1。质子交换膜1电导率：300S m^-1；催化层2电导率：300S m^-1；微孔层3电导率：300S m^-1；气体扩散层电导率：300S m^-1；极板电导率：20000Sm^-1。质子交换膜1导热系数：0.95W m^-1 K^-1；催化层2导热系数：1.0W m^-1 K^-1；微孔层3导热系数：1.0W m^-1 K^-1；气体扩散层导热系数：1.0W m^-1 K^-1；极板导热系数：20W m^-1 K^-1；

催化层2孔隙率：0.3；微孔层3孔隙率：0.4；气体扩散层孔隙率：0.7。催化层2接触角：100°；微孔层3接触角：100°；气体扩散层接触角：120°。催化层2平均孔径尺寸：1.2×10^{- 8}m；微孔层3平均孔径尺寸：1.0×10^-7m；气体扩散层平均孔径尺寸：3.9×10^-5m。催化层2固有过冷却水渗透率：6.2×10^-13m²；微孔层3固有过冷却水渗透率：1.0×10^-12m²；气体扩散层固有过冷却水渗透率：6.2×10^-12m²。

质子交换膜1当量质量：1.1kg mol^-1。催化层2中聚合物体积分数：0.4。单电池与冷却液之间的换热系数：20W m^-2K^-1。

操作工况参数：

运行电流密度：0.1A cm^-2。阴阳极入口气体压强：2.0atm；1.2atm。阴阳极入口相对湿度：0；0。阴阳极化学计量比：2.0；1.5。阴阳极入口气体温度：-20℃；-20℃。环境温度：-20℃。初始膜态水含量：6.4。初始冰体积分数：0。初始过冷却水体积分数：0。时间步长：0.05s。

以下选取0.95s到1.0s的时间段进行本发明实施说明。

(1)膜态水守恒方程求解

式中ρ_MEM表示质子交换膜1密度，EW表示质子交换膜1当量质量，ω表示聚合物体积分数，λ表示膜态水含量，t表示时间，

表示有效膜态水扩散系数，S_mw表示膜态水源项。

在0.95s时刻，阳极催化层、质子交换膜1、阴极催化层中平均膜态水含量分别为5.606，6.404，7.289，平均温度分别为253.923K，253.948K，253.967K，从而通过下式计算阳极催化层、质子交换膜1、阴极催化层中膜态水有效扩散系数分别为7.717×10^-12m² s^-1，2.907×10^-11m² s^-1，7.056×10^-12m² s^-1。

膜态水源项的计算表达式如下：

式中S_react表示氢氧电化学反应产物水，S_m-v表示膜态水与水蒸气之间的相变转化，S_m-1表示膜态水与过冷却水之间的相变转化，S_EOD表示电拖拽效应，S_per表示水力渗透作用，从而计算得出阳极催化层、质子交换膜1、阴极催化层中膜态水源项平均值为-0.660kmolm^-3 s^-1，0kmol m^-3 s^-1，0.725kmol m^-3 s^-1。

经过计算得出，在1.0s时刻，阳极催化层、质子交换膜1、阴极催化层中平均膜态水含量分别为5.571，6.404，7.328。

(2)液压守恒方程求解

过冷却水体积分数基于多孔介质中液压连续进行求解：

式中ε表示孔隙率，ρ_lq表示过冷却水密度，s_lq表示过冷却水体积分数，K_lq表示过冷却水渗透率，μ_lq表示动力粘度，p_l表示液压，S_lq表示源项，u_lq表示过冷却水流速，S_v-l表示过冷却水与水蒸气之间的相变转化，S_l-i表示过冷却水与冰之间的相变转化，

表示水的摩尔质量。

在0.95s时刻，阴极催化层、微孔层3、气体扩散层与流道中的液压数值分别为202783.919Pa，202650Pa，202650Pa，202650Pa，液态水源项分别为3.84kg m^-3 s^-1，0kg m^-3s^-1，0kg m^-3 s^-1，0kg m^-3 s^-1，从而计算得出在1.0s时刻，阴极催化层、微孔层3、气体扩散层与流道中的液压数值分别为202791.703Pa，202650.356Pa，202650Pa，202650Pa。

当液压连续方程得到求解后，水体积分数根据毛细压力计算得出：

p_c＝p_g-p_l

式中p_g表示气压，p_c表示毛细压力。

阴极气体压强的数值为202650Pa，通过计算能够得出，在1.0s时刻，阴极催化层、微孔层3、气体扩散层与流道中的毛细压力分别为141.703Pa，0.356Pa，0Pa，0Pa，过冷却水体积分数通过如下表达式计算：

；

式中σ表示表面张力系数，θ表示接触角，K⁰表示固有渗透率。

通过计算能够得出，在1.0s时刻，阴极催化层、微孔层3、气体扩散层与流道中的过冷却水体积分数分别为0.011，2.888×10^-5，0，0。

(3)冰守恒方程求解

由于多孔介质中结冰点温度比正常结冰点温度略低，因此，为了准确表征相变转换过程，本发明对于催化层2、微孔层3与气体扩散层内的实际结冰点温度进行了修正：

；

式中T_FPD表示结冰点温度修正数值，T_N表示正常结冰点温度(273.15K)，ρ_ice表示冰密度，h_fusn表示融化潜热，r_{CL，MPL，GDL}分别表示催化层2、微孔层3和气体扩散层内平均孔径尺寸，通过计算得出，催化层2、微孔层3与气体扩散层内的实际结冰点温度分别为272.240K，273.041K，273.150K。

冰体积分数在多孔介质与流道中求解，由于冰是固态，故不涉及到流动现象，控制方程如下：

；

式中s_ice表示冰体积分数，S_ice表示冰源项。

在0.95s时刻，阴极催化层、微孔层3、气体扩散层与流道中冰体积分数分别为3.959×10^-3，2.271×10^-6，0，0，冰源项分别为2.993kg m^-3 s^-1，0.008kg m^-3 s^-1，0kg m^-3 s^-1，从而计算得出，在0.95s时刻，阴极催化层、微孔层3、气体扩散层与流道中冰体积分数分别为4.501×10^-3，3.418×10^-6，0，0。

(4)水蒸气守恒方程求解

多孔介质与流道中水蒸气的浓度在每一个时间步进行更新，即采用显式格式更新算法，为提高计算效率，水蒸气浓度在上述各层的中心处求解，离散化处理后计算表达式如下：

；

；

式中

分别表示t时刻催化层2、微孔层3、气体扩散层及流道中的水蒸气浓度，上标t-Δt则代表前一时刻的水蒸气浓度，

分别表示催化层2与微孔层3之间、微孔层3与气体扩散层之间的有效气体扩散系数，δ_MPL表示微孔层3厚度，δ_GDL表示气体扩散层厚度，δ_CH表示流道厚度，S_vp表示水蒸气源项，r_CH表示气体扩散层和流道界面处的有效传质面积系数，Δt表示时间步长。

在0.95s时刻，阴极催化层、微孔层3、气体扩散层的水蒸气浓度分别为6.321×10^{- 2}mol m^-3，5.245×10^-2mol m^-3，2.351×10^-2mol m^-3，1.923×10^-2mol m^-3，水蒸气的源项分别为0.454kmol m^-3 s^-1，0kmol m^-3 s^-1，0kmol m^-3 s^-1，0kmol m^-3 s^-1，催化层2与微孔层3之间、微孔层3与气体扩散层之间的有效气体扩散系数分别为1.901×10^-6m² s^-1，5.186×10^{- 6}m² s^-1，从而能够计算得出，在1.0s时刻，阴极催化层、微孔层3、气体扩散层的水蒸气浓度分别为6.315×10^-2mol m^-3，5.240×10^-2mol m^-3，2.349×10^-2mol m^-3，1.928×10^-2mol m^-3。

(5)能量守恒方程求解

能量守恒方程在整个燃料电池中求解：

；

；

式中

表示有效密度比热容，T表示温度，

表示有效导热率，S_T表示热源项，ASR表示单位面积电阻，I表示电流密度，V_act，a，V_act，c分别表示阳极与阴极的活化损失电压，ΔS表示熵变，F表示法拉第常数，S_pc表示相变潜热。

在0.95s时刻，阳极极板、阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、质子交换膜1、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道、阴极极板的温度分别为253.342K，253.411K，253.728K，254.012K，254.036K，254.061K，254.081K，254.062K，253.734K，253.378K，253.317K，热源项分别为50W m^-3，99.989W m^-3，2.029×10⁴W m^-3，7.172×10³W m^-3，6.760×10⁵W m^-3，7.566×10⁵W m^-3，1.037×10⁸W m^-3，9.987×10³W m^-3，2.029×10⁴W m^-3，99.989W m^-3，50W m^-3，从而能够计算得出，在1.0s时刻，阳极极板、阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、质子交换膜1、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道、阴极极板的温度分别为253.348K，253.416K，253.734K，254.018K，254.042K，254.067K，254.086K，254.068K，253.750K，253.383K，253.322K。

(6)输出性能计算

模型输出电压通过塔菲尔公式进行计算：

V_out＝V_Nemst-V_act-V_ohmic-V_conc；

；

；

V_ohmic＝I·ASR；

；

式中V_out表示输出电压，V_Nernst表示能斯特电压，V_act表示活化损失电压，V_ohmic表示欧姆损失电压，V_conc表示浓差损失电压，ΔG表示吉布斯自由能，T_ref表示参考温度，R表示理想气体常数，

表示阳极氢气浓度，

表示阴极氧气浓度，

表示参考反应气体浓度，α_a，α_c分别表示阳极与阴极电荷传递系数，n_a，n_c分别表示阳极与阴极传递电子数，j_a，j_c分别表示阳极与阴极反应速率，p_a，p_c分别表示阳极与阴极气体压强，ST_a，ST_c分别表示阳极与阴极化学计量比，I_lim，a，I_lim，c分别表示阳极与阴极极限电流密度。

经过计算得出，在1.0s时刻，活化损失电压为0.806V，欧姆损失电压为0.006V，浓差损失电压为3.928×10^-4V，输出电压为0.398V。

图1为本发明实施例所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法相变转化关系示意图；其中，微孔层3内的相变转化关系也是扩散层以及流道内的相变转化关系。

图2是失败启动工况下模型仿真结果与实验数据验证图，图3是成功启动工况下模型仿真结果与实验数据验证图，仿真结果与实验数据吻合良好，模型有效性得到验证。

图4至图6均为本发明实施效果，图4是催化层2中冰体积分数随时间的变化情况，图5是催化层2中水体积分数随时间的变化情况，图6是电池平均温度随时间的变化情况。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：包括以下步骤：

S1、分别对膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程以及模型输出电压性能公式进行求解；

S2、基于步骤S1中的各个方程以及过冷却水结冰机理建立具有过冷却水结冰机理的燃料电池低温冷启动模型。

2.根据权利要求1所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：在步骤S1中的所述能量守恒方程分别与膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程之间数据耦合，所述液压守恒方程分别与冰守恒方程、水蒸气守恒方程之间数据耦合，所述冰守恒方程与水蒸气守恒方程之间数据耦合，所述膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程均控制耦合电池输出电压性能公式；

所述过冷却水结冰机理分别与膜态水守恒方程、液压守恒方程、冰守恒方程、水蒸气守恒方程、能量守恒方程之间数据耦合，膜态水守恒方程中膜态水源项、液压守恒方程中过冷却水源项、冰守恒方程中冰源项、水蒸气守恒方程中水蒸气源项、能量守恒方程中热源项的计算均基于过冷却水结冰机理进行计算；

所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程、过冷却水结冰机理、输出电压性能公式组成燃料电池低温冷启动模型。

3.根据权利要求2所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：在步骤S1中的所述膜态水守恒方程如下：

式中ρ_MEM表示质子交换膜密度，EW表示质子交换膜当量质量，ω表示聚合物体积分数，λ表示膜态水含量，其物理定义为每个质子传导团簇所携带的平均水分子数，t表示时间，

表示有效膜态水扩散系数，S_mw表示膜态水源项；

所述S_mw的计算表达式如下：

式中S_mw表示膜态水源项，S_react表示氢氧电化学反应产物水，S_m-v表示膜态水与水蒸气之间的相变转化，S_m-l表示膜态水与过冷却水之间的相变转化，S_EOD表示电拖拽效应，S_per表示水力渗透作用；

所述有效膜态水扩散系数的参数值基于能量守恒方程的温度参数值进行计算；

所述膜态水源项的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算。

4.根据权利要求2所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：所述过冷却水结冰机理包括以下步骤：

A1、判断催化层中膜态水含量是否达到饱和，若是则进行下一步，否则直接进入步骤A3；

A2、催化层中膜态水含量转变为过冷却水与水蒸气；

A3、催化层中膜态水与水蒸气之间相互转化；

A4、判断催化层、微孔层、气体扩散层与流道中水蒸气是否达到饱和浓度，若是则进行下一步，否则进入步骤A6；

A5、催化层、微孔层、气体扩散层与流道中水蒸气达液化为过冷却水，并进行下一步；

A6、过冷却水从催化层往流道传输，并同时发生结冰现象。

5.根据权利要求3所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：在步骤S1中的所述液压守恒方程的求解操作包括以下步骤：

B1、对液压守恒方程进行求解；

B2、当液压守恒方程得到求解后，对过冷却水体积分数公式以及水力渗透公式计算。

6.根据权利要求5所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：在步骤B1中的所述液压守恒方程如下：

式中ε表示孔隙率，ρ_lq表示过冷却水密度，s_lq表示过冷却水体积分数，t表示时间，表示过冷却水渗透率，μ_lq表示动力粘度，p_l表示液压，S_lq表示过冷却水源项，u_lq表示过冷却水流速；

所述S_lq的计算表达式如下：

式中S_lq表示过冷却水源项，S_m-l表示膜态水与过冷却水之间的相变转化，S_v-l表示过冷却水与水蒸气之间的相变转化，S_l-i表示过冷却水与冰之间的相变转化，S_per表示水力渗透作用，

表示水的摩尔质量；

所述过冷却水源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中液压参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算；

在步骤B2中的所述过冷却水体积分数公式如下：

p_c＝p_g-p_l (5)；

式中p_g表示气压，p_c表示毛细压力，σ表示表面张力系数，θ表示接触角，ε表示孔隙率，K⁰表示固有渗透率，s_lq表示过冷却水体积分数；

在步骤B2中的所述水力渗透公式如下：

式中S_per表示水力渗透作用，ρ_lq表示过冷却水密度，K_per表示水力渗透系数，p_l，CLc表示阴极催化层的液压，p_l，CLa表示阳极催化层的液压，μ_lq表示动力粘度，δ_CL表示催化层厚度，δ_MEM表示质子交换膜厚度

表示水的摩尔质量。

7.根据权利要求6所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：在步骤S1中的所述冰守恒方程的求解操作包括以下步骤：

C1、对于催化层、微孔层与气体扩散层内的实际结冰点温度进行修正操作；

C2、对冰守恒方程进行求解；

在步骤C1中催化层、微孔层与气体扩散层内的实际结冰点温度修正公式如下：

式中T_FPD表示结冰点温度修正数值，T_N表示正常结冰点温度，σ表示表面张力系数，θ表示接触角，ρ_ice表示冰密度，h_fusn表示融化潜热，r_{CL，MPL，GDL}分别表示催化层、微孔层和气体扩散层内孔径尺寸；

在步骤C2中的所述冰守恒方程如下：

式中s_ice表示冰体积分数，S_ice表示冰源项，ε表示孔隙率，ρ_ice表示冰密度，t表示时间；

所述冰源项的参数值基于液压守恒方程中过冷却水体积分数参数值、能量守恒方程中温度参数值进行计算。

8.根据权利要求7所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：在步骤S1中的所述水蒸气守恒方程如下：

式中

分别表示t时刻催化层、微孔层、气体扩散层的水蒸气浓度，

分别表示t-Δt时刻催化层、微孔层、气体扩散层的水蒸气浓度，

分别表示催化层与微孔层之间、微孔层与气体扩散层之间的有效气体扩散系数，δ_CL表示催化层厚度，δ_MPL表示微孔层厚度，δ_GDL表示气体扩散层厚度，S_vp表示水蒸气源项，ε_CL表示催化层孔隙率，ε_MPL表示微孔层孔隙率，Δt表示时间步长；

所述催化层与微孔层之间、微孔层与气体扩散层之间有效气体扩散系数的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、能量守恒方程中温度参数值进行计算；

所述水蒸气源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、能量守恒方程中温度参数值进行计算。

9.根据权利要求8所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：在步骤S1中的所述能量守恒方程如下：

式中

表示有效密度比热容，T表示温度，

表示有效导热率，S_T表示热源项，ASR_CL表示催化层的单位面积电阻，ASR表示单位面积电阻，δ_CL表示催化层厚度，δ表示厚度，I表示电流密度，V_act，a，V_act，c分别表示阳极与阴极的活化损失电压，ΔS表示熵变，F表示法拉第常数，S_pc表示相变潜热，t表示时间；

所述有效密度比热容的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算；

所述有效导热率的参数值基于液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算；

所述热源项的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气守恒方程中水蒸气浓度参数值进行计算。

10.根据权利要求9所述的具有过冷却水结冰机理的燃料电池低温冷启动建模方法，其特征在于：在步骤S1中的所述模型输出电压性能公式如下：

V_out＝V_Nernst-V_act-V_ohmic-V_conc (16)；

V_ohmic＝I·ASR (19)；

式中V_out表示输出电压，V_Nernst表示能斯特电压，V_act表示活化损失电压，V_ohmic表示欧姆损失电压，V_conc表示浓差损失电压，R表示理想气体常数，T表示温度，

表示参考反应气体浓度，α_a，α_c分别表示阳极与阴极电荷传递系数，n_a，n_c分别表示阳极与阴极传递电子数，j_a，j_c分别表示阳极与阴极反应速率，p_a，p_c分别表示阳极与阴极气体压强，ST_a，ST_c分别表示阳极与阴极化学计量比，δ_CLa，δ_CLC分别表示阳极催化层与阴极催化层的厚度，I表示电流密度，ASR表示单位面积电阻，I_lim，a，I_lim，c分别表示阳极与阴极极限电流密度；

所述输出电压的参数值基于膜态水守恒方程中膜态水含量参数值、液压守恒方程中过冷却水体积分数参数值、冰守恒方程中冰体积分数参数值、水蒸气浓度守恒方程中水蒸气浓度参数值、能量守恒方程中温度参数值进行计算。

Images