3.3　曲率的计算

曲率的计算分两步进行：

首先重构SOFC电极的几何结构，由于孔在电极中是随机分布的，颗粒随机堆积法是一种有效的微结构重建手段。由于真实电极中的孔是复杂的无规则结构，采用不同形状的颗粒堆积得到的孔隙结构差别又很大，所以我们很难预先确定采用哪种形状的颗粒来堆积能最好的预测曲率。但从公开的文献报道中，通过比较，我们至少能确定3D球堆的模型在对孔隙曲率的预测上是不准确的。图3-3即为3D球堆模型预测的曲率与实验值的比较，很明显可以看出，模型的结果低于实验值。除了球堆积，也有学者采用3D立方体堆积方法重构电极微观结构，并用以研究SOFC其他方面的问题，这提供了一个很好的思路。于是本章将试着采用这种结构，其有效性将在后面的结果中得到证明。另外，用立方体堆积将极大地简化模型，提高计算效率，这样就可以研究更大的模型。

有了几何模型后，第二步是通过求解扩散方程获得流量并进一步计算出曲率。需要指出的是，这里将要用到的方法是一套相对成熟的方案，并已经在3D球堆积模型和重构模型中广为采用。

3.3.1　3D立方体堆积

如前文所述，本次研究确定采用3D立方体来堆积重构电极结构。除了颗粒形状的选择，还需确定堆积方式。由于电极中孔隙形貌和分布的不规则性，我们可以设想从电极中取出代表体元，将其均匀地分成立方体后，如果立方体足够小，每个立方体既可以近似地被视作孔隙或者固体，而无论是孔隙还是固体的立方体将是随机分布在这个代表体元中的。遵循相反的思路，不妨先用立方体堆积成代表体元，然后根据孔隙率随机地从总数中取出相应数量立方体，赋予它们孔隙的本征特性，如本征扩散率D0，而剩下的则赋予为固体的本征特性，对流体绝缘。举例来说，对于一个x，y，z方向均为10个立方体的代表体元来说，总共有1000个立方体。如果孔隙率为0.3，则将有300个立方体被随机指定为孔隙，剩下的700个则为固体，如图3-4所示。通过这种方式，就很容易但又不失准确性地构建出电极的微观结构。这种方法的一个附加优点是立方体的结构很容易用结构网格进行离散，而且代表体元尺寸可以通过堆积立方体的数量轻松调整，而这些都是3D球堆积或者实验重构法所不具备的。

3.3.2　扩散模拟

重构电极结构后，第二步是在这个结构中模拟扩散过程。尽管电极中的扩散过程十分复杂，包括了分子扩散、Knudsen扩散和黏滞流三种机制。但是考虑到曲率是一个几何参数，所以必须是独立于气体种类和扩散机理的[9]。因此为了简化计算，本文采用最简单的扩散方程即FM来描述电极微结构中的扩散，这套方案广泛见于文献[29～31]。FM的控制方程见式（3-1）。图3-5中展示的是模型的边界设置，为了简化起见，上下面的浓度差Δc被预设为1mol·m-3，其中上边界为ctop=1mol·m-3，下边界为cbottom=0，其他面则为绝缘的零流量边界。有限元软件COMSOL Multiphysics将被用来求解方程（3-1）。求解域（代表体元）也为立方体形状，其边长为L，L是小立方体的边长l的整数倍。

通过求解方程（3-1）获得有效扩散摩尔流量Jeff后，有效扩散系数Deff就可以通过下式获得：

　　（3-15）

式中，A是求解域的横截面积。

将方程（3-15）和方程（3-8）结合，得到曲率的计算式：

　　（3-16）

从上文的论述可知，由于孔隙在电极中是随机分布的，对于一个给定的孔隙率，如果代表体元选取合适，分别计算若干个孔隙率相同的随机模型，由于随机模型本身的因素，所得曲率值也会略微不同，但必然在一个固定值周围波动。所以，为了将每次模型计算的随机差异考虑进来，下面对于每一个孔隙率下曲率的计算将取五个随机模型的平均值，以减少单次计算的不确定性。

3.3.3　模型验证与计算结果分析

（1）网格独立性验证　由于采用了有限元手段来求解扩散方程，为了保证模型计算的准确性，网格独立性验证尤为重要。由于本模型采用的是正六面体网格，所以只需改变每个小立方体边上的节点数就可以很轻松地调整网格数量。通过对不同网格数量的系列化计算，网格独立性得到了仔细验证。图3-6即为曲率值与网格数量的依赖关系，其中每个数据点都添加了1%的误差棒，可以看出当单个小立方体网格数分别为3375和4096时，相对误差已经小于1%。同时考虑到计算时间，本模型将采用3375的网格数量。

（2）代表体元尺寸　由于不可能重构整个电极的微观结构，而只是取出代表体元作为求解域，则取出来的区域必须足够大才能具有代表性，使得任意的两个模型的计算值在正常的误差范围内。更重要的是，由于真实的电极结构是各向同性的，所以构建的模型还要保证从x、y、z三个方向计算的曲率值也要在正常误差范围内。也即取出的代表体元要保证足够大并各向同性。图3-7中（a）和（b）分别展示了对代表体元大小和x、y、z三个方向性质的验证，其中横坐标是代表体元（求解域）边长L与小立方边长l的比值。从图3-7（a）中我们可以看出当L/l=15时，模型就足够大了，因为L/l=15和L/l=17之间的曲率相差不到1%；同样从图3-7（b）中可以看出三个方向的性质在L/l<15时有明显的各向异性，L/l=15时三个方向的计算结果基本一致，说明在L/l=15时已足够保证各向同性。这个结论十分重要，也和很多公开报道文献相一致，前面提到，有学者采用3D球堆模型来研究其他电极参数，也有相同的发现[12，32]，而采用实验重构法的也有类似结论[22]。这也从侧面证明了本模型采用的堆积方式有一定合理性。这个结论为采用FIB-SEM和X射线等实验手段的重构体积提出了要求，重构的模型要达到孔隙特征尺寸的15倍。

（3）孔隙的影响　在公开报道的文献中，有很多关于计算多孔介质曲率的关系式，如表3-1所示。从这些关系式中可以得出两个结论，首先，曲率主要由孔隙率决定；其次，对于不同的多孔介质，曲率与孔隙率之间的关系是不同的。所以下面我们将研究SOFC电极中曲率与孔隙率之间的关系。

　　实际电极中可能出现的孔隙率范围是0.2～0.5。因此，本文中将采用这个孔隙率范围。图3-8中展示了本模型预测的曲率和孔隙率的依赖关系以及和实验值的对比，其中的实验值所针对的电极是通过传统粉末烧结所得，这就意味着孔隙在这种电极中是随机分布的，用3D立方体堆积来构建其电极微观结构是合理的。从图3-8中可以看出本计算的曲率与实验值十分吻合，说明本模型能够用来预测SOFC电极曲率。从图3-8中还可以看出，对于0.2～0.50的孔隙率范围，曲率典型值的为1.3～4，这比很多模型所采用的值要低。同时可以看出曲率值与孔隙率是紧密相关的，而不应该像很多文献中那样将其视作一个可以任意调整的经验值；而且曲率和孔隙率的关系不是简单的线性关系，当曲率较小时，随着孔隙率增大，曲率下降得快，而孔隙率较大时，下降趋势更平缓。