英语原文共 11 页
用多相多重弛豫时间格子玻尔兹曼模型研究质子交换膜燃料电池中气体扩散层水-气输运的过程
摘要
为了研究质子交换膜（PEM）燃料电池系统中气体扩散层（GDL）中的水-气输运过程，我们提出了一种新的计算方法。本文提出了多相、多重弛豫时间格子玻尔兹曼模型。该模型基于平均场扩散界面。理论上可以处理密度比大、粘度大的多相流。使用标准的反弹边界条件以及在GDL中详细的液气输送的防滑和湿润边界墙的近似平均方案，模拟以及精确的边界条件都很难实现。与大多数基于Bhtanagar–Gross–Krook碰撞算子的格子玻尔兹曼方法不同，该模型显示了一个与粘度无关的速度场，这在模拟多种粘度共存的多相流动中非常重要。我们通过模拟湿壁上的静态液滴来验证我们的模型，并与理论预测相比较。然后，我们模拟了质子交换膜燃料电池GDL中的水-气流动，并研究了在不同的条件下与饱和相关的输运性能。结果与以往的数值和理论研究结果在定性上是一致的。
介绍
质子交换膜（PEM）燃料电池被认为是寻找新能源的最佳方法之一。在PEM燃料电池中，两种反应气体氢和氧在约50mu;m厚的膜上结合，覆盖有催化层。膜被约200mu;m厚的气体扩散层（GDL）包围。在这种设置中，GDL起着重要作用，因为它具有几个特定的功能，例如提供反应气体的连续传输和阳极与阴极之间的电子传导性。然而，由于PEM燃料电池通常在较低温度下用加湿的反应气体操作，因此可能引起所谓的“水淹”现象，特别是在阴极侧，水蒸汽凝结并阻塞GDL的孔隙。因此降低了燃料电池的性能。因此，研究GDL中的液气输送并研究输送特性是至关重要的，以便找到具有最佳性能的扩散介质。
在文献中，为了模拟GDL中的液-气两相流动，通常对两种具有特性的相都应用达西定律，尤其是两相效应。其性质是毛细管压力，即两相之间的压差，以及相对渗透率，它调节每相达西定律的渗透率。这些性质可以是液体饱和度的函数，即液体体积与孔隙体积之比。由于GDL通常是一个复杂的随机微观结构，因此很难精确描述其几何结构，从而计算其输运性质。这些困难可以在图1中说明，图1显示了作为典型示例的Toray TGP-H-60的显微照片，并分别通过虚拟材料设计和耗散粒子动力学（DPD）可视化重建的三维（3D）几何。显然，我们我们可以清楚地看到介质结构的强各向异性。
传统上，为了阐明GDL中的两相传输现象，已经开发了各种数值模型。最简单的方法是单相模型，其中气体和液体被认为是单流体混合物，因此共享相同的速度场。而且，完全忽略了界面张力效应。在这种情况下，通过求解单个方程而不区分水蒸气和液态水，可以获得总水量。一旦获得总水浓度场，可以允许水浓度超过饱和水平，基本上假设气相中的超饱和[3-5]。在GDL中，对液态水输送更为严格的方法是将经典的两相模型重新表述为一个单一的方程，并且燃料电池成功运行所必需的界面张力效应和GDL润湿性都是充分的[6-10]。该模型可以有效地生成大部分GDL特征，但忽略了局部结构细节。另一种方法是体积平均或孔隙水平模型[11-16]。这些模型假设局部界面平衡，即孔隙水平的电，化学和热平衡。仔细定义了局部界面平衡的有效性条件。严格地说，所有上述模型都是宏观模型，尽管在某些工作中可能存在理论上的不一致性，并且本质上两个流体相之间的任何界面边界都是介观的。因此，期望开发一种介观模型来研究双流体流动。
近年来，作为一种介观模型，格子玻尔兹曼（LB）方法已成为一种很有前途的多相流模拟工具[17-23]。LB法起源于晶格气体元胞自动机（LGCA）[24,25]。由于LB方法的动力学来源，它具有一些与宏观模型不同的特点。这些特征包括编程简单、内在并行性和复杂边界和多种流体的直接解析，后两个特征在模拟类GDL结构中的多相流方面非常有吸引力。在这些情况下，流体 - 固体界面通过锯齿形方法近似，以便可以直接应用标准反弹边界条件，这通过模拟非粒子 - 固体相互作用来逆转流体粒子与实心边界碰撞的动量。然而，用于多相流建模的LB方法大多是基于Bhatnagar–Gross–Krook（BGK）碰撞模型[26,27]，它们经常遇到数值不稳定性和粘性相关速度场[23,28]等问题。由于两相流中通常同时存在不同的粘度，因此上述限制对数值模拟的精度至关重要。为了克服这些缺点，开发了多弛豫时间（MRT）lb模型[22,23,28]。已经证明，在模拟多相流时，MRTLB模型比BGK模型具有更好的稳定性和准确性[22,23]。然而，大多数MRT多相模型是基于单组分势模型[18]，该模型通过调节相互作用势来模拟两相流，仅对密度比较小（le;100）的两相流有效。最近提出了一个双分布函数MRTLB模型，它克服了以往BGK和MRT模型的大部分局限性，但需要界面张力的隐式处理[29]。
本文建立了一个大密度比、不同粘度的多相流的MRTLB模型，并将其应用于PEM燃料电池GDL中的水气输送模拟。该模型基于扩散界面理论[30,31]。为了处理较大的密度比和不同的粘度，该模型采用了两种分布函数，一种用于模拟不同相位总密度的速度场，另一种用于跟踪界面，包括不同相位密度差的影响。与以往的MRT模型不同，在MRT模型中，界面张力和润湿边界条件必须根据固体壁上的静态液滴提前计算，目前的模型明确地结合了这两个因素。为了捕捉边界效应，除了模拟非滑动边界条件的标准反弹条件外，还提出了一种基于泰勒级数展开的近似平均方法来模拟湿化边界。
本文的其余部分安排如下。 在第2节中，介绍了一种多相MRTLB模型。 我们首先简要介绍扩散界面理论，然后介绍MRTLB方程，然后是非滑移和润湿边界条件的数值实现。 第3节包括两部分。 在第一部分中，通过模拟润湿壁上的静态液滴来验证该模型，并将结果与理论解决方案进行比较。 第二部分，对GDL试验模型中的液气流进行了数值模拟，研究了输运性质。 我们在第4节得出结论。
多相多重弛豫时间格子玻尔兹曼模型
2.1扩散界面理论
我们的模型基于扩散界面理论[29,30]。在这里，我们考虑具有两种流体或相的流动，例如气体和液体，其质量密度分别为rho;g和rho;l，粘度分别为eta;g和eta;l。 整个流动系统的特征在于质量密度rho;=rho;l rho;g，粘度eta;=eta;l eta;g，局部有序参数phi;表示两相的密度差。 然后可以通过以下自由能函数来描述系统的热力学行为
其中V是控制体积，S是润湿壁的表面积。 右侧的第一个积分是标准的Ginzburg-Landau表达式，由[31]定义的体积自由能密度psi;
其中参数k和A确定界面张力sigma;=4radic;2kA/ 3，界面宽度xi;=radic;2k/ A [32]。 利用体积自由能的这种形式，系统将分别贡献两个平衡状态-1和 1，相应的低密度和高密度场。
式（1）中的第二项说明了由于表面能psi;s，固体-流体边界处的特定润湿相互作用。要求f的函数最小化表面热力学-液相和气相可与表面发生不同的相互作用。
（3）
在实心墙上评估，其中n是指向流体的表面法线[33]。 该条件导致在没有流动的情况下在平坦壁处的静态接触角theta;，使得润湿电位gamma;具有相互作用电位
其中beta;= arccos（sin2theta;）。 给定方程（1）的自由能函数，两相流体的动态行为由Navier-Stokes方程控制，
和Cahn-Hilliard方程，
u是宏观流体速度，Fb是体力，M是流动性。 压力张量P的梯度和化学势mu;分别写为[30-32]
cs是后来给出的声速。 系统的体积压力可以计算为
2.2 MRTLB method
在n离散速度的MRT格子玻尔兹曼方法的框架中[28,29]，方程（5）和（6）可通过以下两个方程求解：
其中i = 1,2，...，N，f和g是晶格位置r和时间t处的密度和阶次参数分布函数的相应矢量。 Lambda;alpha;（alpha;= f，g）是由12给出的对角弛豫矩阵
Q是一个N*N矩阵，它将分布函数f和g分别线性转换为速度矩mf和mg，
主体符号f，g，m和G表示尺寸列向量并且具有
这里T代表转置算子，G代表强制体，其中的构成要素是
式中， mu;▽ϕ来自压力张量（见等式（7）），wi是依赖于离散速度的权重系数
模型。对于三维19离散速度（d3q19）模型cs=1/radic;3（delta;r=delta;t=delta;=1），离散速度和加权系数为：
表1给出了and;alpha;中的19个平衡力矩及其相应的松弛系数。
需要注意的是，在meqgi的平衡中，所有非线性速度项都可以省略，因为在验证Cahn–Hilliard方程时不需要四阶各向同性晶格张量（等式（6））[32]。本文选择的弛豫参数通过对MRT模型[34]的分析，来减小粘度依赖速度场。d3q19模型的变换矩阵Q在[34]中给出。
局部密度rho;，阶次参数phi;和动量j分别被计算为
应该指出，在我们的工作中定义的密度rho;可以被认为是标称密度，实际局部密度应该根据阶数参数获得rho;（phi;）=rho;g （phi; - phi;g）（rho;l - rho;g） ）/（phi;1-phi;g），其中phi;1和phi;g是阶次参数的上限和下限，可以通过麦克斯韦等面积规则确定。我们在工作中使用标称密度的原因是它在整个流场中的变化很小，因此满足了LBM的不可压缩极限。 另一方面，密度差反映在有序参数phi;及其相关的化学势mu;中，其充当力并且被添加在流场中。 该力的大小仅由表面张力和界面厚度控制。 因此，在该模型中，可以容易地处理大密度差异情况并且稳定性不是关键。除此之外，可以注意到Q中的所有元素都是整数，矩阵Q的倒数可以通过Q-1 = QT /（Q·QT）来计算，因为Q是正交矩阵，因此，程序可以有效手动编码以减少耗时的矩阵运算，从而提高运行程序的计算效率。 实施本方法的主要步骤可归纳如下：
bull;第一步：根据流体密度和速度计算力矩{meqfi}和{meqgi}（初始设置中mfi=meqfi,）
bull;第一步:在矩空间中计算{meqfi}和{meqgi}
bull;第三步：用公式计算碰撞后密度分布函数（10）和（11）。
bull;第四步：平分分布函数{fi}和{gi}。
表1
平衡矩及其相应的松弛系数f和g（rho;0是密度的初始值，tau;f=3eta;/（rho;delta;t） 0.5和tau;g= M /Gamma; 0.5和Gamma;是一个自由参数控制的迁移率）
2.3衍生物的边界条件和数值评估
与传统的Navier-Stokes解算器相比，LB方法的一个优点是在复杂几何体（如图1所示的GDL结构）上实现边界条件。理论上，由于实验测量或数学描述精确边界位置的困难，这种几何图形只能用锯齿形网格（见图2）来近似。在这种情况下，LB方法的标准反弹边界条件提供了一种根据分布函数fi建立防滑条件模型的有效而简单的方法。如图2所示，来自实体部分的未知分布函数（黑色虚线箭头）可以简单地设置为其对应方向相反的分布函数（蓝色实线箭头）的已知值，实际的防滑边界可以在最后一个流体节点以外的半个网格间距处实现[23]
图2：复杂几何边界条件实现的二维草图（流体点：蓝色实心圆；固体点：黑色实心圆；最近和下一个最近的流体侧点：绿色短划线椭圆圆内）
由于湿润边界条件（式（4））是一个平衡条件，因此将平衡{meqgi}[35]施加它是合适的。从方程（4），为了设置平衡边界条件，我们通常需要计算边界的法向量，以便获得边界节点上的阶数参数phi;及其二阶导数nabla;2phi;。然而，正如我们之前强调的，对于像GDL这样的几何图形来说是很困难的。为了避免这种困难，在这项工作中，我们使用近似策略来评估它们。在图2中，边界点P（红色实心圆）被一些流体点（蓝色实心圆）和一些实心点（黑色实心圆）包围。为了获得点P上的phi;和nabla;2phi;的值，在点P的最近和下一个最近的相邻流体侧点（见图2）上进行了phi;的泰勒级数展开，对于每个展开点，我们得到了phi;p=phi;f，lminus;delta;part;r_ O（delta;2），可进一步近似为phi;p=phi;f，lminus;delta;part;nphi;。因此，平均而言，点P上的有序参数ϕ及其二阶导数nabla;2phi;可近似为：
图3。在theta;=60°和120°时，扩散和收缩液滴形状的时间变化。等值面表示阶数参数a=0。
其中delta;是单位网格空间，phi;f，l是点B的最近和最近相邻流体侧点上的有序参数; l = 1,2，...，N，其中N代表点P的最近和下一个最近的流体侧点的总数; gamma;可以从方程式（4）直接计算公式20和21可以通过格子模板直接使用，并且可以通过离散速度的流动方向判断和计数流体侧点，无论其是否指向流体场。
在LB框架中，流体场中phi;的一阶和二阶导数也可以通过使用晶格模板常规地计算，从而减少离散误差。 在这项工作中，使用以下模板方案：
其中j = 1,2,3表示尺寸。
结果与讨论
3.1模型验证
为了验证当前模型，模拟保持在墙壁上的静态液滴，并测试接触角。 考虑将半径为R的半球形液滴置于壁上，并且在施加于壁上的湿润电势gamma;的作用下使其在静止时达到平衡状态。
图4.计算的接触角与理论预测的比较（方程（4））（三角形，理论预测;圆形，模拟结果）
图5.三维气体扩散层模型，其边缘平行于从左到右的主流。 切片显示平行于主流方向的中平面处的速度矢量（delta;p= 1.7times;10-3，theta;= 120°，eta;l/eta;g= 100）。
图6.计算的GDL绝对渗透率与流体动力粘度的函数的比较
计算域分为51times;
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
