本发明属于质子交换膜燃料电池领域，尤其是涉及一种采用显式格式更新算法的燃料电池瞬态建模方法。
质子交换膜燃料电池（pemfc）具有高功率密度、高能量转换效率、零排放等优点，被认为是未来交通运输行业的清洁能源之一。
典型的燃料电池结构包括质子交换膜、催化层、微孔层、气体扩散层、流道与极板。质子交换膜的作用是传输电化学反应的离子，即将阳极生成的质子传递到阴极，而电子无法直接穿过质子交换膜，此外，质子交换膜起到隔离阴阳极反应气体的作用。催化层提供电化学反应发生的场所，氢气在阳极催化层中失去电子发生氧化反应，氧气在阴极催化层中得到电子发生还原反应，并且与质子结合生成产物水。催化层与质子交换膜统称为膜电极结构。微孔层能够改善电池内部的水传输情况。气体扩散层提供结构支撑作用，并且传输反应气体与产物水。流道提供反应气体由电池入口到电池出口的传输路径。极板提供结构支撑作用与收集电流作用。
燃料电池运行过程中，氢氧电化学反应生成水的同时，伴随着热量的产生，若产热无法及时排出，则会造成电池运行温度过高，对膜电极结构造成不可逆破坏，导致电池耐久性的降低。由于质子传导率受到膜电极湿润程度的影响，若反应气体加湿不足，则会造成膜干现象，增大欧姆电压损失并导致输出性能的下降，但是电池内部液态水无法及时排出电池则会占据多孔介质中气体传输通道并覆盖催化层中有效反应位点，严重降低输出性能。因此，燃料电池内部的传热传质过程是“水-气-热-电”强耦合的，其中水热管理问题不仅影响了电池的输出性能，而且显著影响了耐久性与使用寿命。由于燃料电池的膜电极与多孔介质层厚度为微米级别，实验手段无法在线观测到电池内部的水热传输现象，仿真模型为探究电池内部的耦合传热传质机理提供了重要研究方法，而商业化软件中（如simulink，amesim）燃料电池模块为零维模型，仅计算输出电压，无内部的传热传质过程，导致模型准确性欠缺。此外，三维模型计算资源需求高且稳态模型居多，难以满足长时间瞬态仿真分析的需求。因此，建立计算效率高、可靠性强、准确度好的仿真模型对于燃料电池水热管理策略的开发至关重要。
技术实现要素：
有鉴于此，本发明旨在提出一种采用显式格式更新算法的燃料电池瞬态建模方法，通过显式格式更新算法在每层结构的中心处进行物理量求解，层内不再细分网格，从而极大提高模型计算效率,以弥补商业化软件中零维燃料电池模型准确性欠缺的不足，也解决了三维模型难以实现长时间瞬态仿真分析的短板，能够为水热管理策略开发提供仿真数据支持，极大的减小实验成本及研发周期。
为达到上述目的，本发明的技术方案是这样实现的：
一种采用显式格式更新算法的燃料电池瞬态建模方法，包括以下步骤：
s1、建立采用显式格式的膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程以及能量守恒方程，建立电池输出电压性能公式；
s2、分别对步骤s1中的所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程和能量守恒方程之间相互耦合，并分别对膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程进行求解，得到守恒方程更新算法求解结果；
s3、将守恒方程更新算法求解结果输入到电池输出电压性能公式内，得到电池输出电压性能公式求解结果；
s4、基于步骤s3中的守恒方程更新算法求解结果以及电池输出电压性能公式求解结果建立显式格式更新算法的燃料电池瞬态模型。
进一步的，在步骤s3中的所述能量守恒方程分别与膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程之间数据耦合，所述液压守恒方程分别与冰守恒方程、气体守恒方程之间数据耦合，所述冰守恒方程与气体守恒方程之间数据耦合，所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程均控制耦合电池输出电压性能公式，所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程、电池输出电压性能公式组成燃料电池瞬态模型。
进一步的，所述膜态水守恒方程如下：
（1）；
（2）；
（3）；式中，，分别表示时刻阳极催化层、质子交换膜、阴极催化层的膜态水含量，，，分别对应时刻阳极催化层、质子交换膜、阴极催化层的膜态水含量，，分别表示阳极催化层与质子交换膜之间、阴极催化层与质子交换膜之间的有效膜态水传输系数，，，分别表示阳极催化层、质子交换膜、阴极催化层的厚度，表示质子交换膜密度，表示质子交换膜当量质量，表示时间步长大小，，，分别表示阳极催化层、质子交换膜、阴极催化层聚合物体积分数，表示膜态水源项；所述的计算表达式如下：（4）；式中表示膜态水与水蒸气之间的相变源项，表示电化学反应产物水源项，表示膜态水与液态水之间的相变源项，表示电拖拽效应源项；在时刻下，所述有效膜态水传输系数的参数值基于能量守恒方程的催化层温度参数值、能量守恒方程的质子交换膜温度参数值进行计算。
进一步的，在步骤s2中的所述液压守恒方程的求解操作包括以下步骤：
a1、对液压守恒方程进行求解；
a2、当液压守恒方程得到求解后，对多孔介质层液态水体积分数公式以及流道液态水公式计算。
进一步的，在步骤a1中的所述液压守恒方程如下：
（5）；
（6）；
（7）；式中表示液压数值，，，分别表示时刻表示催化层、微孔层、气体扩散层的液压数值，，，分别表示时刻催化层、微孔层、气体扩散层的液压数值，表示液态水体积分数，分别表示微孔层与气体扩散层之间、催化层与微孔层之间、气体扩散层有效液态水传输系数，，，分别表示催化层、微孔层、气体扩散层的厚度，表示时间步长大小，表示液态水源项，表示液态水密度，，，分别表示催化层、微孔层、气体扩散层的孔隙率，表示气体扩散层与流道交界面的传质面积系数；在时刻下，所述有效液态水传输系数的参数值基于液态水体积分数参数值、冰守恒方程的冰体积分数参数值、能量守恒方程的催化层温度参数值、能量守恒方程的微孔层温度参数值、能量守恒方程的气体扩散层温度参数值进行计算；
所述的计算表达式如下：（8）；式中表示液态水与水蒸气之间的相变源项，表示液态水与冰之间的相变源项，表示水力渗透源项，表示流道入口处与出口处液态水流量的差值，表示水的摩尔质量，表示膜态水与液态水之间的相变源项；
步骤a1求解液压守恒方程后，在步骤a2中的所述多孔介质层液态水体积分数公式如下：
（9）；（10）；式中表示液压数值，，分别表示气压，毛细压力数值，表示时刻的毛细压力数值，，分别表示时刻，时刻的液态水体积分数，表示表面张力系数，表示接触角，表示孔隙率，表示固有液态水渗透率；
所述流道液态水公式如下：
（11）；
（12）；式中分别表示时刻气体扩散层、流道的液压数值，，分别表示时刻，时刻流道液态水体积分数，，分别表示气体扩散层、流道厚度，表示表面张力系数，表示接触角，表示固有液态水渗透率，表示孔隙率，表示液态水源项，表示液态水密度，表示气体扩散层有效液态水传输系数，表示时间步长大小，表示时刻流道的毛细压力数值。
进一步的，所述冰守恒方程如下：（13）；式中，分别表示时刻，时刻的冰体积分数，表示冰密度，表示冰源项，表示孔隙率，表示时间步长大小，表示时刻的冰源项；所述如下：（14）；式中表示液态水与冰之间的相变源项，表示冰源项，在时刻下，所述冰源项的参数值基于液态水体积分数的参数值、气体守恒方程的气体浓度参数值、能量守恒方程的微孔层温度参数值、能量守恒方程的气体扩散层温度参数值、能量守恒方程的流道温度参数值进行计算。
进一步的，所述气体守恒方程如下：（15）
（16）；
（17）；
（18）；式中i表示氢气、氧气、氮气与水蒸气，，，，分别表示时刻催化层、微孔层、气体扩散层、流道中的气体浓度，，，，表示时刻的气体浓度，，，分别表示催化层与微孔层之间、微孔层与气体扩散层之间、气体扩散层与流道之间的有效气体扩散系数，表示气体源项，，，，分别表示催化层、微孔层、气体扩散层、流道的厚度，，，分别表示催化层、微孔层、气体扩散层的孔隙率，表示气体扩散层与流道交界面的传质面积系数，表示时间步长大小；在时刻下，所述有效气体扩散系数的参数值基于液态水体积分数的参数值、冰体积分数的参数值、能量守恒方程的催化层温度参数值、能量守恒方程的微孔层温度参数值、能量守恒方程的气体扩散层温度参数值、能量守恒方程的流道温度参数值进行计算。
进一步的，所述能量守恒方程如下：
（20）；
（21）；
（22）；
（23）；
（24）；
（25）；式中，，，，，分别表示时刻极板、流道、气体扩散层、微孔层、催化层、质子交换膜的温度，，，，，，分别表示时刻极板、流道、气体扩散层、微孔层、催化层、质子交换膜的温度，，，，，分别表示极板与流道之间、流道与气体扩散层之间、气体扩散层与微孔层之间、微孔层与催化层之间、催化层与质子交换膜之间的有效导热系数，，分别表示阴极催化层与质子交换膜之间、质子交换膜与阳极催化层之间的有效导热系数，表示极板与环境的换热系数，表示比热容，表示密度，表示热源项，表示流道的厚度，，，分别表示催化层、微孔层、气体扩散层、极板的厚度；所述的计算表达式如下：（26）；式中表示热源项，表示电流密度，，，，，，分别表示催化层、质子交换膜、气体扩散层、微孔层、流道、极板的单位面积电阻，表示气体对流道的热效应，，分别表示阳极与阴极活化损失，表示相变潜热，表示温度，表示熵变，表示法拉第常数，，，，，分别表示催化层、微孔层、气体扩散层、极板、流道、质子交换膜的厚度；在时刻下，所述有效导热系数的参数值基于液态水体积分数的参数值、冰体积分数的参数值进行计算；在时刻下，所述的参数值基于膜态水含量的参数值、液态水体积分数的参数值、冰体积分数的参数值、气体浓度的参数值以及的温度参数值、的温度参数值、的温度参数值、的温度参数值、的温度参数值、的温度参数值进行计算；在时刻下，所述，的参数值基于液态水体积分数的参数值、冰体积分数的参数值、气体浓度的参数值以及的温度参数值、的温度参数值、的温度参数值、的温度参数值、的温度参数值、的温度参数值进行计算；在时刻下，所述的参数值基于膜态水源项的参数值、液态水源项的参数值、冰源项的参数值、气体源项的参数值进行计算。
进一步的，所述电池输出电压性能公式如下：
（27）；（29）；
（30）；（31）；式中表示输出电压，，，，分别表示能斯特电压、活化损失电压、欧姆损失电压、浓差损失电压，表示理想气体常数，表示阳极催化层中氢气浓度，表示阴极催化层中氧气浓度，表示阳极参考氢气浓度，表示阴极参考氧气浓度，，分别表示阳极与阴极电荷传递系数，，分别表示阳极与阴极传递电子数，，分别表示阳极与阴极反应速率，，分别表示阳极与阴极极限电流密度，表示电流密度，表示单位面积电阻，表示法拉第常数，表示温度，，分别表示阳极催化层、阴极催化层的厚度；在时刻下，，的参数值基于液态水体积分数的参数值、冰体积分数的参数值以及温度参数值、温度参数值、温度参数值、温度参数值、温度参数值、温度参数值进行计算；在时刻下，，的参数值基于液态水体积分数的参数值、冰体积分数的参数值、气体浓度的参数值以及温度参数值、温度参数值、温度参数值、温度参数值、温度参数值、温度参数值进行计算。
相对于现有技术，本发明所述的一种采用显式格式更新算法的燃料电池瞬态建模方法具有以下优势：
（1）本发明所述的一种采用显式格式更新算法的燃料电池瞬态建模方法，充分计算电池内部的“水-气-热-电”耦合传输过程以准确反映内部能量流动规律，通过显式格式更新算法在每层结构的中心处进行物理量求解，层内不再细分网格，从而极大提高模型计算效率。本发明能够深入探究燃料电池操作工况与结构设计参数的影响，如电流密度、运行温度、气体加湿程度、多孔介质层憎水性、聚合物/碳载体比值等，同时本发明能够详细给出燃料电池内部各物理量的瞬态变化情况，如膜态水含量、水体积分数、气体浓度、温度等。
（2）本发明所述的一种采用显式格式更新算法的燃料电池瞬态建模方法既弥补了商业化软件中零维燃料电池模型准确性欠缺的不足，也解决了三维模型难以满足长时间瞬态仿真分析的短板，能够为水热管理策略开发提供仿真数据支持，极大减小实验成本及研发周期。
附图说明
构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为本发明实施例所述的一种采用显式格式更新算法的燃料电池瞬态建模方法其模型仿真结果与实验数据验证图；
图2为本发明实施例所述的一种采用显式格式更新算法的燃料电池瞬态建模方法其燃料电池输出电压随时间的变化情况示意图；
图3为本发明实施例所述的一种采用显式格式更新算法的燃料电池瞬态建模方法其质子交换膜中膜态水含量随时间的变化情况示意图；
图4为本发明实施例所述的一种采用显式格式更新算法的燃料电池瞬态建模方法其阴极催化层中氮气浓度随时间的变化情况示意图；
图5为本发明实施例所述的一种采用显式格式更新算法的燃料电池瞬态建模方法其阴极催化层温度随时间的变化情况示意图；
图6为本发明实施例所述的一种采用显式格式更新算法的燃料电池瞬态建模方法各个守恒方程之间的数据耦合原理图。
具体实施方式
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明中的具体含义。
下面将参考附图并结合实施例来详细说明本发明。
名词解释：显式格式：为显式差分格式，显式差分格式属于数值格式中的一种。
如图1至图6所示，一种采用显式格式更新算法的燃料电池瞬态建模方法，包括以下步骤：s1、建立采用显式格式的膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程以及能量守恒方程，建立电池输出电压性能公式；
s2、分别对步骤s1中的所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程和能量守恒方程之间相互耦合，并分别对膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程进行求解，得到守恒方程更新算法求解结果；
s3、将守恒方程更新算法求解结果输入到电池输出电压性能公式内，得到电池输出电压性能公式求解结果；
s4、基于步骤s3中的守恒方程更新算法求解结果以及电池输出电压性能公式求解结果建立显式格式更新算法的燃料电池瞬态模型。本建模方法充分计算电池内部的“水-气-热-电”耦合传输过程，通过显式格式更新算法在每层结构的中心处进行物理量求解，能够探究燃料电池操作工况与结构设计参数的影响，同时能够给出电池内部各物理量的瞬态变化情况；且本文中的更新算法，就是在不同的时刻，利用五个守恒方程进行迭代更新计算的表达式/算法。
在步骤s3中的所述能量守恒方程分别与膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程之间数据耦合，所述液压守恒方程分别与冰守恒方程、气体守恒方程之间数据耦合，所述冰守恒方程与气体守恒方程之间数据耦合，所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程均控制耦合电池输出电压性能公式，所述膜态水守恒方程、液压守恒方程、冰守恒方程、气体守恒方程、能量守恒方程、电池输出电压性能公式组成燃料电池瞬态模型；其中“水”为膜态水守恒方程、液压守恒方程、冰守恒方程，“气”为气体守恒方程，“热”为能量守恒方程，“电”为输出性能计算。
所述膜态水守恒方程如下：（1）；
（2）；
（3）；式中，，分别表示时刻阳极催化层、质子交换膜、阴极催化层的膜态水含量，，，分别对应时刻阳极催化层、质子交换膜、阴极催化层的膜态水含量，，分别表示阳极催化层与质子交换膜之间、阴极催化层与质子交换膜之间的有效膜态水传输系数，，，分别表示阳极催化层、质子交换膜、阴极催化层的厚度，表示质子交换膜密度，表示质子交换膜当量质量，表示时间步长大小，，，分别表示阳极催化层、质子交换膜、阴极催化层聚合物体积分数，表示膜态水源项；
所述的计算表达式如下：（4）；式中表示膜态水与水蒸气之间的相变源项，表示电化学反应产物水源项，表示膜态水与液态水之间的相变源项，表示电拖拽效应源项；在时刻下，所述有效膜态水传输系数的参数值基于能量守恒方程的催化层温度参数值、能量守恒方程的质子交换膜温度参数值进行计算。即有效膜态水传输系数的数值受到阳极催化层温度、质子交换膜温度以及阴极催化层温度的影响。
在步骤s2中的所述液压守恒方程的求解操作包括以下步骤：
a1、对液压守恒方程进行求解；a2、当液压守恒方程得到求解后，对多孔介质层液态水体积分数公式以及流道液态水公式计算。
在步骤a1中的所述液压守恒方程如下：
（5）；
（6）；
（7）；式中表示液压数值，，，分别表示时刻表示催化层、微孔层、气体扩散层的液压数值，，，分别表示时刻催化层、微孔层、气体扩散层的液压数值，表示液态水体积分数，分别表示微孔层与气体扩散层之间、催化层与微孔层之间、气体扩散层有效液态水传输系数，，，分别表示催化层、微孔层、气体扩散层的厚度，表示时间步长大小，表示液态水源项，表示液态水密度，，，分别表示催化层、微孔层、气体扩散层的孔隙率，表示气体扩散层与流道交界面的传质面积系数；在时刻下，所述有效液态水传输系数的参数值基于液态水体积分数参数值、冰守恒方程的冰体积分数参数值、能量守恒方程的催化层温度参数值、能量守恒方程的微孔层温度参数值、能量守恒方程的气体扩散层温度参数值进行计算，其数值受到液态水体积分数、冰体积分数以及温度的影响；
所述的计算表达式如下：（8）；式中表示液态水与水蒸气之间的相变源项，表示液态水与冰之间的相变源项，表示水力渗透源项，表示流道入口处与出口处液态水流量的差值，表示水的摩尔质量，表示膜态水与液态水之间的相变源项；液压守恒方程由于催化层、微孔层、气体扩散层之间的憎水性差异会导致界面处液态水体积分数突变现象，本发明中液态水体积分数基于多孔介质中液压连续进行求解。
步骤a1求解液压守恒方程后，在步骤a2中的所述多孔介质层液态水体积分数公式如下：（9）；（10）；式中表示液压数值，，分别表示气压，毛细压力数值，表示时刻的毛细压力数值，，分别表示时刻，时刻的液态水体积分数，表示表面张力系数，表示接触角，表示孔隙率，表示固有液态水渗透率；所述流道液态水公式如下：
（11）；
（12）；式中分别表示时刻气体扩散层、流道的液压数值，，分别表示时刻，时刻流道液态水体积分数，，分别表示气体扩散层、流道厚度，表示表面张力系数，表示接触角，表示固有液态水渗透率，表示孔隙率，表示液态水源项，表示液态水密度，表示气体扩散层有效液态水传输系数，表示时间步长大小，表示时刻流道的毛细压力数值。为了准确的反映电池内部的水传输过程，本发明不仅求解了多孔介质内的液态水，而且求解了流道中的液态水。
所述冰守恒方程（即催化层、微孔层、气体扩散层与流道中的冰体积分数计算表达式）如下：（13）；式中，分别表示时刻，时刻的冰体积分数，表示冰密度，表示冰源项，表示孔隙率，表示时间步长大小，表示时刻的冰源项；所述如下：（14）；式中表示液态水与冰之间的相变源项，表示冰源项，在时刻下，所述冰源项的参数值基于液态水体积分数的参数值、气体守恒方程的气体浓度参数值、能量守恒方程的微孔层温度参数值、能量守恒方程的气体扩散层温度参数值、能量守恒方程的流道温度参数值进行计算；其数值受到液态水体积分数、气体浓度以及温度的影响；其数值受到液态水体积分数、气体浓度以及温度的影响。
所述气体守恒方程如下：由于扩散作用是气体在多孔介质中沿着垂直于膜方向上传输的主导因素，气体组分（氢气、氧气、氮气、水蒸气）的计算表达式如下：（15）
（16）；
（17）；
（18）；式中i表示氢气、氧气、氮气与水蒸气，，，，分别表示时刻催化层、微孔层、气体扩散层、流道中的气体浓度，，，，表示时刻的气体浓度，，，分别表示催化层与微孔层之间、微孔层与气体扩散层之间、气体扩散层与流道之间的有效气体扩散系数，表示气体源项，，，，分别表示催化层、微孔层、气体扩散层、流道的厚度，，，分别表示催化层、微孔层、气体扩散层的孔隙率，表示气体扩散层与流道交界面的传质面积系数，表示时间步长大小，其中氢气源项与氧气源项受到电化学反应影响，水蒸气源项受到液态水体积分数、冰体积分数、饱和水蒸气压强以及温度影响；在时刻下，所述有效气体扩散系数的参数值基于液态水体积分数的参数值、冰体积分数的参数值、能量守恒方程的催化层温度参数值、能量守恒方程的微孔层温度参数值、能量守恒方程的气体扩散层温度参数值、能量守恒方程的流道温度参数值进行计算；其数值受到液态水体积分数、冰体积分数以及温度的影响。
其中水蒸气源项计算表达式如下：（19）；式中表示流道入口处与出口处水蒸气流量的差值。
所述能量守恒方程如下：由于实际流入电池的反应气体温度可能与电池运行温度存在差异，因此本发明计算了反应气体对电池的热效应。能量守恒方程在燃料电池每层结构中进行求解，计算表达式如下：（20）；
（21）；
（22）；
（23）；
（24）；
（25）；式中，，，，，分别表示时刻极板、流道、气体扩散层、微孔层、催化层、质子交换膜的温度，，，，，，分别表示时刻极板、流道、气体扩散层、微孔层、催化层、质子交换膜的温度，，，，，分别表示极板与流道之间、流道与气体扩散层之间、气体扩散层与微孔层之间、微孔层与催化层之间、催化层与质子交换膜之间的有效导热系数，，分别表示阴极催化层与质子交换膜之间、质子交换膜与阳极催化层之间的有效导热系数，表示极板与环境的换热系数，表示比热容，表示密度，表示热源项，表示流道的厚度，，，分别表示催化层、微孔层、气体扩散层、极板的厚度；
所述的计算表达式如下：
（26）；式中表示热源项，表示电流密度，，，，，，分别表示催化层、质子交换膜、气体扩散层、微孔层、流道、极板的单位面积电阻，表示气体对流道的热效应，，分别表示阳极与阴极活化损失，表示相变潜热，表示温度，表示熵变，表示法拉第常数，，，，，分别表示催化层、微孔层、气体扩散层、极板、流道、质子交换膜的厚度；在时刻下，所述有效导热系数的参数值基于液态水体积分数的参数值、冰体积分数的参数值进行计算；在时刻下，所述的参数值基于膜态水含量的参数值、液态水体积分数的参数值、冰体积分数的参数值、气体浓度的参数值以及的温度参数值、的温度参数值、的温度参数值、的温度参数值、的温度参数值、的温度参数值进行计算，其数值受到膜态水含量、液态水体积分数、冰体积分数、气体浓度以及温度的影响；在时刻下，所述，的参数值基于液态水体积分数的参数值、冰体积分数的参数值、气体浓度的参数值以及的温度参数值、的温度参数值、的温度参数值、的温度参数值、的温度参数值、的温度参数值进行计算；其数值受到液态水体积分数、冰体积分数、气体浓度以及温度的影响；在时刻下，所述的参数值基于膜态水源项的参数值、液态水源项的参数值、冰源项的参数值、气体源项的参数值进行计算；其数值受到膜态水源项、液态水源项、冰源项以及气体源项的影响。
所述电池输出电压性能公式如下：（27）；
（28）；
（29）；
（30）；（31）；式中表示输出电压，，，，分别表示能斯特电压、活化损失电压、欧姆损失电压、浓差损失电压，表示吉布斯自由能，表示参考温度，表示熵变，表示理想气体常数，表示阳极催化层中氢气浓度，表示阴极催化层中氧气浓度，表示阳极参考氢气浓度，表示阴极参考氧气浓度，，分别表示阳极与阴极电荷传递系数，，分别表示阳极与阴极传递电子数，，分别表示阳极与阴极反应速率，，分别表示阳极与阴极极限电流密度，表示电流密度，表示单位面积电阻，表示法拉第常数，表示温度，，分别表示阳极催化层、阴极催化层的厚度；在时刻下，，的参数值基于液态水体积分数的参数值、冰体积分数的参数值以及温度参数值、温度参数值、温度参数值、温度参数值、温度参数值、温度参数值进行计算；其数值受到液态水体积分数、冰体积分数以及温度的影响；在时刻下，，的参数值基于液态水体积分数参数值、冰体积分数参数值、气体浓度参数值以及温度参数值、温度参数值、温度参数值、温度参数值、温度参数值、温度参数值进行计算，其数值受到液态水体积分数、冰体积分数、气体浓度以及温度的影响。
实施例1
以下通过具体计算实施例1对发明的方法以及模型建立的具体步骤进行详细的说明。
有效反应面积：120cm2。极板厚度：2mm；流道厚度：1mm；气体扩散层厚度：0.19mm；微孔层厚度：0.02mm；催化层厚度：0.01mm；膜厚度：0.050mm；流道长度：0.2m；流道与肋板宽的比值：1。质子交换膜密度：1980kgm-3；催化层密度：1000kgm-3；微孔层密度：1000kgm-3；气体扩散层密度：1000kgm-3；极板密度：1000kgm-3。质子交换膜比热容：833jkg-1k-1；催化层比热容：3300jkg-1k-1；微孔层比热容：568jkg-1k-1；气体扩散层比热容：2000jkg-1k-1；极板比热容：1580jkg-1k-1。质子交换膜电导率：5000sm-1；催化层电导率：5000sm-1；微孔层电导率：5000sm-1；气体扩散层电导率：5000sm-1；极板电导率：20000sm-1。质子交换膜导热系数：0.95wm-1k-1；催化层导热系数：1.0wm-1k-1；微孔层导热系数：1.0wm-1k-1；气体扩散层导热系数：1.0wm-1k-1；极板导热系数：20wm-1k-1；催化层孔隙率：0.4；微孔层孔隙率：0.4；气体扩散层孔隙率：0.7。催化层接触角：100º；微孔层接触角：120º；气体扩散层接触角：120º。催化层固有过冷却水渗透率：6.2×10-13m2;微孔层固有过冷却水渗透率：6.2×10-13m2;气体扩散层固有过冷却水渗透率：6.2×10-12m2。质子交换膜当量质量：1.1kgmol-1。催化层中聚合物体积分数：0.4。单电池与冷却液之间的换热系数：20wm-2k-1。运行电流密度：0.2acm-2。阴阳极入口气体压强：1.0atm；1.0atm。阴阳极入口相对湿度：1；1。阴阳极入口气体温度：25℃；25℃。环境温度：25℃。初始膜态水含量：6.0。以下选取0.9s到1.0s的时间段对本实施说明。
（1）膜态水守恒方程求解,；
；
；式中，，分别表示时刻阳极催化层、质子交换膜、阴极催化层的膜态水含量，，，分别对应时刻的膜态水含量，表示阳极催化层与膜之间的有效膜态水传输系数，表示阴极催化层与膜之间的有效膜态水传输系数，，，分别表示阳极催化层、质子交换膜、阴极催化层的厚度，表示质子交换膜密度，表示质子交换膜当量质量，表示时间步长大小，表示聚合物体积分数。在0.9s时刻，阳极催化层、质子交换膜、阴极催化层中膜态水含量，，分别为4.849，5.992，8.237，温度，，分别为298.805k，298.827k，298.846k，从而通过下式计算阳极催化层与质子交换膜之间、质子交换膜与阴极催化层之间有效膜态水扩散系数，分别为1.012×10-10m2s-1，9.207×10-11m2s-1;；
；；；膜态水源项通过如下表达式进行计算：；式中表示膜态水源项，表示膜态水与水蒸气之间的相变源项，表示电化学反应产物水源项，表示膜态水与液态水之间的相变源项，表示电拖拽效应源项，计算表达式如下：；；；；以阴极催化层为例，在0.9s时刻，电化学反应产物水源项为1.036kmolm-3s-1，膜态水与水蒸气之间的相变源项为0.121kmolm-3s-1，膜态水与液态水之间的相变源项为0.121kmolm-3s-1，电拖拽效应源项为1.853kmolm-3s-1，计算得出阴极催化层、阳极催化层的膜态水源项分别为2.646kmolm-3s-1、-1.180kmolm-3s-1。
经过计算得出，在1.0s时刻，阳极催化层、质子交换膜、阴极催化层中膜态水含量，，分别为4.782，5.991，8.427。
（2）液压守恒方程求解；液态水体积分数基于多孔介质中液压连续进行求解：
；
；
；式中，，分别表示时刻表示催化层、微孔层、气体扩散层的液压数值，，，分别表示时刻的液压数值，表示液态水体积分数，表示微孔层与气体扩散层之间的有效液态水传输系数，表示催化层与微孔层之间的有效液态水传输系数，表示气体扩散层的有效液态水传输系数，，，分别表示催化层、微孔层、气体扩散层的厚度，表示液态水源项，表示液态水密度，，，分别表示催化层、微孔层、气体扩散层的孔隙率，表示气体扩散层与流道交界面的耦合传质面积系数，其数值根据肋宽比进行计算，如果气体流道宽度和肋宽度相同，则等于0.5。
；；；；
以阴极为例，在0.9s时刻，催化层、微孔层、气体扩散层与流道中的液压数值，，分别为101386.692pa，101328.100pa，101325.004pa，101325.000pa，催化层、微孔层、气体扩散层的温度，，分别为298.846k，298.827k，298.582k，经过上述表达式计算得出，催化层与微孔层之间的有效液态水传输系数为3.892×10-13m2s-1，微孔层与气体扩散层之间的有效液态水传输系数为3.759×10-12m2s-1，气体扩散层的有效液态水传输系数为6.961×10-12m2s-1。
液态水源项通过如下表达式进行计算:；
式中表示液态水与水蒸气之间的相变源项，表示液态水与冰之间的相变源项，表示水力渗透源项，表示流道入口处与出口处液态水流量的差值，表示水摩尔质量。
；；以阴极催化层为例，在0.9s时刻，水蒸气与液态水之间的相变源项为-0.0076kgm-3s-1，液态水与冰之间相变源项为0，水力渗透源项为0.034kgm-3s-1，经过计算得出阴极催化层中液态水源项为2.183kgm-3s-1。同理可以得出微孔层、气体扩散层、流道中液态水源项分别为0kgm-3s-1，0kgm-3s-1，0kgm-3s-1，0.0504kgm-3s-1，0kgm-3s-1，0kgm-3s-1，0kgm-3s-1；经过计算得出，在1.0s时刻，阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道中液压数值分别为101325.000pa，101325.002pa，101325.513pa，101328.446pa，101380.575pa，101327.401pa，101325.003pa，101325.000pa。
当液压连续方程得到求解后，水体积分数根据毛细压力计算得出：
；；
；
；式中，分别表示气压，毛细压力数值，表示时刻的毛细压力数值，，分别表示时刻，时刻的液态水体积分数，表示表面张力系数，表示接触角，表示固有液态水渗透率，表示时刻流道的液压数值，，分别表示时刻，时刻流道的液态水体积分数，表示流道的厚度。
经过计算得出，在1.0s时刻，阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道中液态水体积分数分别为0，0，1.443×10-5，2.988×10-4，5.067×10-3，7.597×10-5，0，0。
（3）冰守恒方程求解,催化层、微孔层、气体扩散层与流道中的冰体积分数计算表达式如下：；在0.9s时刻，催化层、微孔层、气体扩散层与流道中的冰体积分数分别为0，0，0，0，对应的冰源项分别为0，0，0，0。经过计算得出，在1.0s时刻，催化层、微孔层、气体扩散层与流道中的冰体积分数分别为0，0，0，0。
（4）气体守恒方程.气体组分（i表示氢气、氧气、氮气、水蒸气）的计算表达式如下：；
；
；
；式中，，，分别表示时刻催化层、微孔层、气体扩散层、流道中的气体浓度，，，，表示时刻的气体浓度，，，分别表示催化层与微孔层之间、微孔层与气体扩散层之间、气体扩散层与流道之间的有效气体扩散系数，表示气体源项。
以水蒸气浓度为例，在0.9s时刻，阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道中水蒸气浓度分别为0.904molm-3,0.905molm-3,0.905molm-3,0.905molm-3,1.112molm-3,1.109molm-3,1.099molm-3,1.094molm-3。水蒸气的有效扩散系数计算表达式如下：；；；经过计算得出，在0.9s时刻，阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道中水蒸气有效扩散系数分别为6.849×10-5m2s-1，1.482×10-5m2s-1，6.407×10-6m2s-1，3.370×10-6m2s-1，3.349×10-6m2s-1，6.408×10-6m2s-1，1.482×10-5m2s-1，6.847×10-5m2s-1。
水蒸气源项计算表达式如下：；经计算得出，在0.9s时刻，阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道中水蒸气源项分别为-0.056molm-3,-3.690×10-6molm-3,-0.002molm-3,2.798molm-3,121.280molm-3,-0.006molm-3,-2.862×10-5molm-3,-2.425molm-3；经过计算得出，在1.0s时刻，阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道中水蒸气浓度分别为0.897molm-3,0.897molm-3,0.898molm-3,0.898molm-3,1.123molm-3,1.119molm-3,1.109molm-3,1.104molm-3；同理，经过计算得出，在1.0s时刻，阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层中氢气浓度分别为38.915molm-3,38.902molm-3,38.878molm-3,38.862molm-3,阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道中氧气浓度分别为4.965molm-3,4.984molm-3,5.030molm-3,5.055molm-3。阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道中氮气浓度分别为1.832×10-3molm-3,1.837×10-3molm-3,1.848×10-3molm-3,1.853×10-3molm-3,33.222molm-3,33.222molm-3,33.222molm-3,33.222molm-3。
（5）能量守恒方程求解；能量守恒方程在燃料电池每层结构中进行求解，计算表达式如下：；
；
；
；
；
；式中，，，，，分别表示时刻极板、流道、气体扩散层、微孔层、催化层、质子交换膜的温度，，，，，，分别表示时刻的温度，，，，，分别表示极板与流道之间、流道与气体扩散层之间、气体扩散层与微孔层之间、微孔层与催化层之间、催化层与膜之间的有效导热系数，，分别表示阴极催化层与膜之间、膜与阳极催化层之间的有效导热系数，表示极板与环境的换热系数，表示比热容，表示密度，表示热源项；在0.9s时刻，阳极极板、阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、质子交换膜、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道、阴极极板的温度分别为298.287k，298.347k，298.580k，298.786k，298.805k，298.827k，298.846k，298.827k，298.582k，298.313k，298.262k。以阴极为例，经计算得出，催化层与膜之间、微孔层与催化层之间、气体扩散层与微孔层之间、流道与气体扩散层之间、极板与流道之间的有效导热系数分别为0.864wm-1k-1，0.603wm-1k-1，0.331wm-1k-1，1.697wm-1k-1，15.013wm-1k-1。
热源项计算表达式如下：
；式中表示电流密度，表示单位面积电阻，、分别表示阳极与阴极活化损失，表示相变潜热，表示温度，表示熵变，表示法拉第常数，表示气体对流道的热效应。经计算得出，在0.9s时刻，阳极极板、阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、质子交换膜、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道、阴极极板的热源项分别为-9.834×103wm-3，329.854wm-3，4.869×103wm-3，1.627×103wm-3，2.221×107wm-3，1.543×106wm-3，1.348×108wm-3，1.492×103wm-3，4.868×103wm-3，-6.069×104wm-3，-9.828×103wm-3。
因此，在1.0s时刻，经过计算得出，阳极极板、阳极流道、阳极气体扩散层、阳极微孔层、阳极催化层、质子交换膜、阴极催化层、阴极微孔层、阴极气体扩散层、阴极流道、阴极极板的温度分别为298.306k，298.365k，298.598k，298.803k，298.823k，298.844k，298.864k，298.845k，298.599k，298.328k，298.277k。
（6）输出性能计算；模型输出电压通过塔菲尔公式进行计算：
；；
；；；式中表示输出电压，表示能斯特电压，表示活化损失电压，表示欧姆损失电压，表示浓差损失电压，表示吉布斯自由能，表示参考温度，表示理想气体常数，表示阳极催化层中氢气浓度，表示阴极催化层中氧气浓度，表示电荷传递系数，表示传递电子数，，分别表示阳极与阴极反应速率，，表示参考反应气体浓度，表示极限电流密度，经过计算得出，在1.0s时刻，活化损失电压为0.548v，欧姆损失电压为0.071v，浓差损失电压为3.621×10-3v，输出电压为0.592v。
图1是模型仿真结果与实验数据验证图，仿真结果与实验数据吻合良好，模型有效性得到验证。
图2至图5为本发明实施效果，图2是燃料电池输出电压随时间的变化情况，图3是质子交换膜中膜态水含量随时间的变化情况，图4是阴极催化层中氮气浓度随时间的变化情况，图5是阴极催化层温度随时间的变化情况。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。