一种用于随机重建气体扩散层的框架，该框架能够增强纤维形态的调谐效果，并实现添加剂的混合融合

《Materials Today Communications》：A framework for stochastic reconstruction of gas diffusion layers with enhanced fiber morphology modulation and additive incorporation

　　杨达南|Himani Garg|Martin Andersson
瑞典隆德大学工程学院能源科学系，邮政信箱118，隆德，SE-221 00

**摘要**
气体扩散层（GDLs）是质子交换膜燃料电池和水电解系统中的关键纤维多孔组件。通过对GDL微结构的随机重建，可以从经验性的试错方法转向基于性能的材料设计。然而，目前的重建方法往往忽略了商用GDL所固有的复杂形态特征，如纤维弯曲度、粘合剂/PTFE的结合以及局部孔隙率分布。此外，在生成过程中精确计算孔隙率仍然具有挑战性，这通常受到简化分析模型（忽略纤维交叉点）和计算密集型体素基方法之间的权衡限制。为了克服这些限制，本研究提出了一种新的、统一的GDL重建框架，该框架包含三个关键创新点：（1）一种通用的纤维生成方案，可以利用多种曲线来构建从直线到高度弯曲的碳纤维的各种几何形状；（2）一种混合堆叠策略，可以减少局部孔隙率波动并精确控制垂直于平面的各向异性；（3）一种双阶段孔隙率计算算法，在保持精确度的同时提高计算效率。此外，通过图像形态技术集成添加剂成分，以模拟真实的炭化粘合剂和疏水涂层分布。最后，流体传输特性的模拟验证了重建GDL与先前实验数据和分析模型的一致性，确立了该框架作为下一代质子交换膜燃料电池和水电解系统靶向设计的强大工具。

**1. 引言**
气体扩散层（GDLs）是先进电化学能量转换设备中的关键功能组件，特别是在质子交换膜燃料电池[1]、水电解器[2]和金属-空气电池[3]中。作为多相传输的主要场所，GDLs同时介导反应气体扩散、产物水去除、电子传输和热耗散。无论其具体的形态分类（即碳纸、碳毡或碳布），这些不同类型的多孔介质都决定了电化学设备内的复杂耦合传输现象。因此，深入阐明其微观结构特性对于最大化整体能量转换效率、避免关键淹没状态和提高体积功率密度至关重要。
为了阐明这些基本的结构-性能关系，研究人员采用了各种方法论途径。基于同步辐射的X射线断层扫描等实验表征方法能够实现操作条件下的高分辨率、原位观察[4]。然而，这些技术经常受到高昂设备成本、耗时过程以及有限参数研究或逆向设计能力的限制。为了补充实验数据，基于严格数学推导的分析模型已被开发出来，以建立直接关联[5]、[6]、[7]。例如，GDL的渗透性已被分析与孔隙率和纤维形状等微观结构特征相关联[8]。尽管这些理论和实验工具很有价值，但数值模拟已成为阐明复杂多物理传输现象（如流体、热、电）在复杂纤维微结构中的最不可或缺的方法[9]、[10]。利用计算流体动力学和有限元分析，这些方法建立了在实验上难以获取的定量关系。然而，任何此类模拟的可靠性从根本上取决于一个关键前提，即获得高保真度、计算上可行的纤维多孔结构的数字表示。

基于几何的重建方法通过逐步将数学生成参数与真实的形态特征相结合而发展起来。早期的努力主要集中在基本的几何输入上，如纤维直径[9]、[11]、体积孔隙率[12]和厚度[13]、[14]，以建立碳骨架的基线物理框架。与简化的均质模型不同，这些高度详细的纤维结构有助于研究GDL孔隙空间内的两相传输[9]、[14]，从而能够明确考虑纤维表面润湿性和孔径分布。随后，意识到实际制造过程会产生高度各向异性的结构，研究人员明确纳入了物理纤维取向[15]、[16]、[17]、[18]。这是通过引入概率密度函数来精确调整纤维垂直于平面的取向角度实现的，反映了碳纸制造的现实世界约束。此外，商用GDLs本身含有二次添加剂，如用于机械/电固化的炭化粘合剂和用于液体水管理的疏水聚合物（例如PTFE）。为了忠实模拟炭化树脂在纤维交叉处的物理积累，基于形态处理方法的复杂粘合剂分布被集成到数字模型中[19]、[20]、[21]、[22]、[24]。这种特定的形态细化有效地弥合了理想化模型与实际几何形状之间的差距，解决了模拟与实验测量属性之间的历史差异。因此，这些结构改进直接使得对耦合热和质量传输的数值预测更加准确。

虽然这些改进显著提高了数字结构的保真度，但它们基本上依赖于理想化的几何假设，留下了两个关键限制未得到解决。首先，许多现有方法将纤维近似为直圆柱[19]。这种假设可能适用于Toray系列的GDLs，但无法再现其他商业产品（如Freudenberg H系列碳纸）所固有的纤维弯曲特性。因此，使用直圆柱模型的模拟无法捕捉实际材料的复杂交织和增强的迂回性，不可避免地导致毛细压力和两相流动路径的预测不准确。其次，大多数现有算法依赖于简化的纤维堆叠策略，通常在整个领域内采用随机分布以满足规定的体积孔隙率[15]、[16]、[17]。这些方法缺乏调节空间孔隙率变化的能力，因此无法重建对驱动液体水远离反应区域和防止灾难性淹没至关重要的工程化局部或垂直于平面的孔隙率梯度。

我们之前的工作[25]试图通过使用杆状周期性表面模型生成弯曲纤维，并实施逐层堆叠策略来施加特定于各层的孔隙率分布，以解决这些差距。据我们所知，这是第一个同时结合纤维弯曲度和对局部孔隙率分布的确定性控制的统一重建模型。尽管如此，这个初步模型仍存在几个限制。最关键的是，它在堆叠过程中没有考虑纤维的垂直于平面的取向，导致相邻层之间的孔隙率出现明显振荡。此外，由于弯曲度受到单一数学公式的限制，可获得的纤维几何形状的多样性仍然有限。

此外，孔隙率评估在指导纤维堆叠决策中起着核心作用，这与孔隙或纤维体积密切相关。然而，大多数现有研究并没有明确报告体积的计算方法，这使得很难评估它们重建程序的可靠性和可重复性。在提供明确方法的研究中，所报告的方法通常可以分为两类，即分析性圆柱体积计算[9]和基于体素的随机计算[26]、[27]。然而，分析方法无法准确确定不规则纤维几何形状或粘合剂结构的体积，并且在迭代堆叠过程中经常导致纤维重叠区域的重复计数。基于体素的方法虽然更为鲁棒，能够处理通用几何形状，但由于需要精细的空间分辨率以实现足够的精度，通常会承受较高的计算成本。

为了应对这些限制，本研究提出了一种新的、全面的重建框架，能够捕捉纤维形态的全部变化性，同时保持对孔隙率特性的精确控制。首先，建立了一种通用的纤维生成方案，可以构建从直线到高度弯曲的各种几何形状，具有可自定义的截面，从而扩展了传统模型之外的几何多样性。这通过与混合堆叠策略的结合，同时调节纤维方向的各向异性并消除人为的孔隙率波动，确保精确的局部结构控制。为了协调计算精度和成本之间的权衡，开发了一种先进的混合孔隙率计算算法，它结合了分析估计和基于体素的校正。此外，该框架通过集成图像处理技术来重建真实的粘合剂和添加剂相，超越了骨架结构。最后，通过质量和动量传输的数值模拟严格验证了生成微结构的保真度，其中预测的渗透率和扩散率与实验数据和经典分析模型进行了对比。通过成功结合多样化的纤维形态并增强对局部孔隙率分布的捕捉，该框架为缩小理想化几何形状与实际GDL物理特性之间的差距提供了可靠的解决方案，为下一代燃料电池的逆向设计和材料优化提供了必要的基础。

**2. 模型开发**
在本研究中，GDL微结构的重建是通过多步骤计算程序实现的，该程序在Python环境中使用现有的开源工具完成。过程从生成和堆叠碳纤维开始，形成多孔纤维骨架框架，然后可选地加入添加剂材料以模拟一些真实的组成特征。

**2.1. 结构关键参数准备和假设**
在开始重建过程之前，需要根据商用GDL规格定义一些关键结构参数。之前的一项研究[9]回顾了一些GDL的典型属性范围，包括纤维直径为7–11 μm、体积孔隙率为0.54-0.9以及厚度为100–400 μm。此外，垂直于平面和沿厚度方向的孔隙率分布可以根据实验扫描图像通过图像处理来确定[28]。所有这些关键几何特征通常来自公开的商业数据表、实验表征研究和高分辨率成像数据库。适当选择GDL类型及其相关参数对于提高随机重建微结构的真实性至关重要。在本研究中，我们主要关注非织造GDL产品，基于以下关键假设：
•每根碳纤维都具有特定的截面形状，其核心几何参数（包括纤维直径、长度和弯曲度）可以定义为预设的固定值，或者作为从用户定义的连续概率分布中抽样的随机变量，以反映材料的真实异质性。
•纤维在垂直于平面和平行于平面的方向上都有明确的取向。
•允许纤维交叉，因为在之前的研究中这非常有限且通常被忽略[16]、[24]。
•添加剂成分（即PTFE和粘合剂）优先分布在纤维连接处，这与实验观察到的积累行为一致[29]。

**2.2. 纤维生成**
现有的GDL随机重建主要使用直圆柱形纤维几何形状[9]，这适用于基于碳纸的基底，如Toray TGP-H-060/090/120系列和某些SGL SIGRACET?产品（例如31BA）。底层算法已经很成熟，通常涉及在两点之间随机放置中心轴，然后生成具有指定直径的圆柱表面。虽然这种方法适用于主要含有直线纤维的材料建模，但它无法捕捉非织造GDL（如Freudenberg H系列（例如H23C5）和某些SGL基于纤维的基底（例如SIGRACET? 10CA）[30]的制造过程中固有的纤维弯曲特性。很少有研究关注这一特征。在我们之前的研究中，使用杆状周期性表面模型生成了弯曲纤维[25]。然而，该模型产生的弯曲特性是规则性的，缺乏随机变化性，无法生成具有自然随机弯曲的纤维。此外，生成具有极端弯曲的纤维带来了额外的挑战，限制了模型在表示完整纤维几何形状范围方面的适用性。所有这些限制都强调了需要更多几何上多样化的纤维生成技术来准确表示更广泛的商用GDL微结构。

本研究介绍了一种改进的方法，用于生成直线和弯曲纤维，以解决上述问题。首先生成一个直的或弯曲的中心轴，然后围绕该轴构建具有定义半径的圆柱表面。对于每个单独生成的纤维，可以通过从概率密度函数（例如高斯分布或均匀分布）中抽样来动态分配该半径。这种随机分配允许框架在同一区域内明确生成具有高度可变和异质纤维直径的微结构。可以使用已建立的几何公式生成曲线形中心轴，包括插值（例如三次样条）、逼近（例如贝塞尔曲线、B样条和非均匀有理B样条）以及随机过程（如随机行走和几何螺旋）。这些选择大大增加了模仿真实弯曲纤维特征的可能性。在这些方法中，一个关键因素是定义控制点的数量，这些控制点不仅区分了各自的核心算法，更重要的是，它们直接控制了所得曲线的关键属性，包括其分割曲率、自由度，以及最终其对预期路径的几何忠实度。本研究采用贝塞尔曲线来生成直纤维和曲纤维，如图1所示，该框架可以很容易地扩展到未来的工作中，以适应其他或混合曲线公式。贝塞尔曲线的构建需要指定控制点，包括端点（A0和A1）以及一个或多个中间点（B1, B2, …, BM）。端点主要决定了纤维的整体方向和长度，而中间点则影响其曲率和局部几何特征。在模拟域内直接指定两个端点会限制对端点相对位置和间距的控制能力。为了提高可控性和灵活性，引入了三个参数来定义曲线端点：一个位于域内的基点B0、一个方向向量v，以及一个预期的纤维长度L。

基点B0是在一个固定域内确定的，该域的尺寸为[具体尺寸未提供]。首先随机初始化一个点P0，然后使用二维Sobol序列更新其x-y坐标，以实现均匀的空间分布，并避免纯随机采样可能导致的纤维聚集，如图A1所示。方向向量v定义为：
v = (θ0, θ1)，其中θ0表示平面内方向角，θ1表示垂直于平面的角度。θ0和θ1的值根据各自的范围内的特定规则来选择。在本研究中，θ0是随机选择的，而θ1取决于一个概率分布函数，这将在下一节中详细讨论。可选地，为了避免生成小于预设阈值的过短纤维，首先计算由v和方向向量确定的直线与域边界的交点P1和P2。如果交点之间的距离小于预设阈值，则重新生成方向向量v以满足这一限制。

直接控制纤维的曲率长度需要迭代调整控制点，直到达到目标长度，这在计算上效率低下且难以收敛。考虑到实际纤维长度与其两个端点之间的距离密切相关，我们改为调节端点之间的距离d，从而间接控制纤维长度。这种方法提供了一种计算效率高且几何上一致的生成不同曲率纤维的方式。应当注意的是，一旦引入曲率，实际曲率长度总是会超过d；然而，这种偏差在整个样本集合中统计上是一致的。在实际应用中，也很难在代表性的样本区域内精确控制纤维长度。尽管如此，所提出的方法通过统计控制端点距离，有效地捕捉了短纤维、中等长度或长纤维的主要长度特征。目标端点距离d可以固定为一个常数值，或者从指定的分布函数中采样，以适应不同的结构特性。最终端点将通过以下步骤确定。

一方面，如果A0和A1被指定为纤维中心轴的方向，则相应的纤维会因为域的大小而部分截断。指定这些端点可以保持纤维轴在整体材料中的连续性，并确保纤维方向和连通性的统计特性能够代表整体材料。另一方面，如果使用A1、B1、C1和D1，则根据以下方程计算两个新点P2和P3：
P2(x, y) = A1 + v * (x - A0) + D1 * (y - A1)
此外，检查P2和P3是否位于固定域内。如果任何坐标超出了域的限制，则将其夹紧到相应的边界值，即：
P2(x, y) = min(Ax, Ay) if x ≤ Ax and y ≤ Ay
P2(x, y) = max(Bx, By) if x ≥ Ax and y ≤ By
这确保了所有生成的纤维端点都位于计算域内。然后，将得到的受限点P2和P3分别指定为纤维的中心点A0和A1。

沿纤维中心轴的中间控制点则是通过考虑曲率控制来确定的。对于直纤维，两个确定的端点足以定义直的中心轴。但对于曲纤维，至少需要一个中间控制点来定义贝塞尔曲线的形状。控制点的数量、位置以及最终的曲率可以通过多种方式设置。在本研究中，我们采用以下方法来生成中间控制点：
(4)
在方程(4)中，前两项定义了连接端点的弦线上的点序列P3, P4, …, Pn，其中n是控制点的数量。k和l是权重系数，满足条件k + l = n，并且k * l ≤ L（L是纤维的总长度）。k和l直接影响纤维的曲率频率。过多的控制点和过小的k可能会导致不现实的纤维结果。中间控制点的数量n由总纤维长度与参考长度L的比例决定，参考长度L定义为L的上限值。因此，k的增加会导致离散化密度的降低。所以，k的值应该相应地调整。为了使这些控制点沿弦线均匀分布，定义了权重系数ω(k, l) = n / L，这些系数将段长度按比例划分。最后，方程(4)中的最后一项引入了一个垂直于参考轴的横向偏移量，从而控制纤维的整体曲率。偏移量的大小定义为：
Δθ = k * (ln(L / n)) / L
其中Δθ表示曲率偏移因子，可以定义为常数或遵循特定概率密度函数的随机变量。Δθ的值为0对应于完全直的纤维。使用一个随机符号运算符σ来使偏移点在弦线上振荡，以确保无偏的双侧分布。引入参数β来限制控制点从连接两个端点的弦线的横向位移，β是一个缩放因子。具体来说，对于给定的k，结果曲率与β成正比。

单位向量e0垂直于纤维向量v0，并定义为：
e0 = (cos(θ0), sin(θ0))
当k = 0时，所有曲线点都位于一个平面内；而非零的β引入了一个垂直于平面的旋转，从而产生三维的曲率轨迹。

接着，使用一组控制点P3, P4, …, Pn来构建复合贝塞尔曲线。由于n阶贝塞尔曲线需要n个控制点，直接使用高阶曲线可能导致数值不稳定和曲线的不必要振荡。为了减轻这一问题，我们结合使用二次和三次贝塞尔段。具体来说，连接到端点A0和A1的终端段由二次贝塞尔曲线定义，而中间区间由三次贝塞尔曲线覆盖[31]。为了确保连接处的平滑性和连续性，为每个内部点P3, P4, …, Pn生成两个辅助支持点P3', P4'，如图2所示：
(7)
(8)
其中γ用于控制局部曲率，在本研究中γ设置为0.5。

因此，v0和v1分别用于两个二次贝塞尔曲线：
v0 = (cos(θ0), sin(θ0))
v1 = (cos(θ1), sin(θ1))
v2, v3, …, vn分别用于生成三次贝塞尔曲线：
v2 = (cos(θ1 + δ), sin(θ1 + δ))
v3 = (cos(θ2 + δ), sin(θ2 + δ))
vn = (cos(θn + δ), sin(θn + δ))

根据获得的贝塞尔曲线，可以通过改变参数t从0到1来生成一系列样本点。在本研究中，采用了一种均匀采样策略，总共生成H个样本点。默认的采样点数量设置为：
N = ceil(L / n)
其中N是向上取整的函数。包括N是为了自动调整点数，以减少计算成本。因此，纤维中心轴C由一系列离散的三维点定义。为了生成具有精确曲率特性的纤维表面网格，沿离散的中心轴构建了一个平行传输框架。这种方法避免了传统Frenet-Serret框架在曲率接近零时出现的奇异性，确保了局部正交坐标系的连续性和稳定性。

每个中心轴点的切向量如下确定。对于内部点，切向量是相邻线段的平均值：
tk = (v0 * (x - x0) + v1 * (y - y0)) / 2
对于端点，切向量由相邻线段定义：
tk = v0 * (x - x0) + v1 * (y - y0)

为了消除尺度依赖性，将tk规范化以获得最终的单位切向量：
t?k = (tk / |v0||v1|)^(1/2)

局部右手法线向量n?k和副法线向量n?k通过平行传输构建，以沿中心轴保持平滑变化的正交框架。为了保持法向量场的连续性，我们使用旋转最小化（或平行传输）方法来更新法向量。选择一个初始法向量n?0垂直于v0。对于每个点P3, P4, …, Pn，计算从n?0到n?k的最小旋转角度θ，以及旋转轴和旋转角度θ：
θ = arccos((v0 * (x3 - x0) - v1 * (y3 - y0)) / (v1 * (x4 - x0) - v2 * (y4 - y0))
然后构建表示旋转的的四元数q：
q = exp(i * θ)

接着，使用四元数q获得旋转后的法向量n?k：
n?k = q * e0

通过上述信息，可以生成表面网格。例如，为了生成圆柱形纤维，在每个中心轴点，使用局部框架定义一个半径为r的圆形截面：
r?k = r * cos(θ0), sin(θ0)
在区间[θ0, θn]内将r?k离散化为N个点。使用三角形补片连接所有截面的顶点。具体来说，对于相邻的截面P3, P4, …, PN，对于每对相邻顶点(P3, P4)和(P4', P5)，生成两种类型的三角形：
(23)

通过沿中心轴重复此过程，生成了一个平滑的高质量管状三角化表面网格。通过修改最后一步IV-c，可以用其他类型的表面形状（如三角形、矩形等）来调整纤维形状。对于椭圆形截面，指定两个不同的半径r1和r2，以及角坐标α1, α2。对于矩形截面，将α1设置为4，并为角坐标α1, α2分配值。通过调整α1和α2，可以灵活控制矩形的宽高比。同样，对于三角形截面，也有类似的公式。这种通用公式使得在统一的几何框架内生成具有多样化截面轮廓的纤维成为可能。

2.3. 纤维堆叠

在纤维生成过程之后，一个重要的后续步骤是纤维的空间堆叠，以构建整体结构。这个过程需要考虑几个关键因素，包括目标整体孔隙率、局部孔隙率分布和纤维方向特性。在文献中，广泛采用了两种通用的堆叠策略，例如随机堆叠和逐层堆叠。这两种方法之间的主要区别在于是否明确控制了局部孔隙率的贯穿平面变化。虽然逐层方法允许改善沿贯穿平面方向的结构异质性的调节，但仍存在一些限制。具体来说，大多数现有研究仅考虑了逐层的孔隙率，其中每层的厚度大约等于纤维直径。然而，当在更小的切片间隔下检查孔隙率时，由于单个纤维的圆柱形状，相邻层之间经常会出现明显的波动。此外，在逐层堆叠中，几乎不考虑贯穿平面的方向。为了解决这些问题，本文首次提出了一种混合堆叠策略，以实现对局部孔隙度的更好地控制，同时保留纤维介质的随机特性。所提出的混合堆叠策略包括以下步骤。为了清晰地展示这些调控机制，图3展示了混合堆叠策略的完整动态反馈控制循环。控制逻辑由三个核心自适应机制管理：（1）孔隙度偏差加权排序，优先在孔隙度缺陷最大的层中沉积纤维；（2）在定义的全局偏差阈值（0.1）触发双阶段孔隙度计算切换；（3）定向规则化协议，在连续生成失败或总层数的百分比达到某个值时，严格强制水平纤维放置。首先，根据参考结构确定整体孔隙度和平面切片孔隙度分布，以确保平均切片孔隙度与规定的整体孔隙度相匹配。纤维域的总厚度定义为，每层的厚度计算为，其中表示层数，通常由纤维厚度决定。因此，通过平均相应厚度区间内的切片级孔隙度值来获得每层的孔隙度。

下载：下载高分辨率图片（631KB）
下载：下载全尺寸图片

图3. 混合堆叠策略的动态反馈控制循环流程图。

整体孔隙度和层孔隙度均初始化为1，在计算域内每一根纤维生成后都会进行更新。然后应用孔隙度偏差加权调整方案来细化z方向上的采样范围。在每次迭代中，评估并排名所有层中目标孔隙度值与当前孔隙度值之间的偏差。偏差最大的层被选为后续纤维生成的优先区域。一旦某层达到其指定的孔隙度，该层就被“冻结”，意味着不允许在该层内添加额外的纤维端点或交叉点。

为了抑制切片级孔隙度的高频振荡，在纤维堆叠过程中引入了两种调控机制。首先，允许一定程度的平面内纤维倾斜，以增强局部结构连续性和物理真实性。其次，采用混合堆叠规则。具体来说，纤维最初随机堆叠，直到指定比例的层达到目标孔隙度。在此之后，新生成的纤维优先沉积在孔隙度低于期望阈值的层中，同时逐渐限制纤维的平面内取向角度。当整体孔隙度和层孔隙度在预定义的容差范围内收敛时，生成过程终止。

为了实现上述细节，有一些关键部分需要详细说明，包括孔隙度计算和纤维取向控制，这些在之前的研究[32]、[33]中描述得不够充分。

2.3.1 孔隙度计算
孔隙度表示在添加固体纤维后指定域内剩余的空隙空间。因此，可以通过以下公式计算孔隙度：
（25）
其中表示纤维体积，表示指定域的体积。因此，核心方法是使用圆柱管的体积计算公式，即截面积与纤维长度的乘积，这在上述纤维生成过程中可以轻松计算。然而，这种方法难以考虑纤维的允许交叉区域。此外，这种方法也不适用于局部层状区域以及随后不规则附加结构的体积计算。另外，在之前的研究[16]中也采用了一种基于体素的方法。纤维生成域被离散成具有指定空间分辨率的体素。一旦在域内生成纤维，就会以某种方式识别出被该纤维占据的体素，从而能够通过考虑目标体素来估计整体和局部孔隙度。当后续纤维与已生成的纤维相交时，重叠的体素只计算一次，从而消除了交叉区域的影响。此外，这些标记的体素可以很容易地划分为层，以确定各层的孔隙度分布。

然而，这种基于体素的方法面临两个主要挑战。首先是计算成本与准确度之间的权衡。更细的体素分辨率可以提高孔隙度估计的准确性，但也会导致计算负载呈立方级增加。当微结构的特征长度（例如纤维直径）显著小于整个域的大小时，需要更高的分辨率。第二个挑战是识别被纤维相占据的体素。对于圆柱形纤维，可以通过计算每个纤维采样点与其相邻体素之间的欧几里得距离来实现。满足最小距离标准（等于纤维半径）的体素被标记为属于纤维相。在整个域中，被识别为纤维相的体素总数为。然而，这种体素标记策略对于非圆柱形或不规则纤维几何形状来说是不够的。

为了克服这些限制，本研究首次开发了一种混合孔隙度计算方法，以实现计算准确性和效率之间的平衡。该方法分为两个阶段。在纤维堆叠的早期阶段，基于表面三角剖分体积计算纤维体积，即计算由连接三角形表示的三维表面所包围的空间体积[34]。然而，在堆叠过程的早期阶段，忽略了纤维交叉体积对孔隙度的影响，因为总固体体积只是简单地通过累加所有单个纤维或其相应剪裁部分的体积来获得。结果，由于交叉区域被重复计算，孔隙度被低估了。虽然通过网格合并操作生成单一封闭表面是一个潜在的解决方案，但其实际应用受到了过度计算开销和在处理数百根纤维时可能出现灾难性故障的风险的限制。因此，在第二阶段使用基于体素的计算方法来解决纤维重叠问题，该方法涉及对纤维占据的体素进行高级评估。为了连接这两个重建阶段，引入了一个与实际整体孔隙度相关的判断指数，并始终存储生成纤维的多边形数据。判断指数（表示为阈值）定义为实际整体孔隙度与目标孔隙度之间的差异，在本研究中设定为0.1。默认值0.1被确定为在计算效率和几何准确度之间严格平衡的最佳经验阈值。具体来说，如果指数大于0.1，则会在堆叠过程的早期触发计算密集型的基于体素的计算，这会显著且不必要地增加总模拟时间。相反，如果指数过小，则会过度延迟这一触发，有可能导致算法发散；这是因为从第一阶段（表面三角剖分）继承下来的未校正的固体体积累积成在指定域内无法数学修复的显著孔隙度偏差。

每个阶段是如何进行的？具体来说，在纤维表面生成过程中，每个纤维表面由一个三角剖分的网格表示，其两端由三角形元素封闭，形成一个多面体结构。然后应用散度定理仅使用表面三角剖分来计算这个多面体几何体的体积。这种方法适用于任意或不规则形状的纤维，前提是表面是封闭的，这为各种情况提供了更高的灵活性。因此，一旦纤维在域内生成，就可以直接计算整体纤维体积，并且整体孔隙度可以表示为：
（26）
其中表示纤维的总数，表示第根纤维的体积。表示组成单个纤维的三角剖分网格面数，而、和分别对应于质心位置向量、外法向量以及三角元素的表面积。

为了实现层属性计算，在纤维生成过程中需要同步执行两个必要步骤。首先，每根生成的纤维根据预定义的层厚度和取向进一步划分为一系列指定层。在剪裁操作之后，通过额外的三角剖分表面封闭纤维的开放截面，以确保每个分割层形成一个不透水的多面体结构。因此，可以使用与整体孔隙度相同的公式计算每层内的孔隙度，相关参数更新以反映特定层的几何形状。其次，每根生成的纤维都被体素化并合并成相同大小的体素。需要注意的是，如果未激活体素计数，此操作不会导致大量的计算。此外，这样可以避免与一次性体素化相关的额外存储和计算成本。因此，这里的局部目标是一种高效且稳健的体素化方法。

与其使用体素质心与纤维中心轴之间的距离来标记体素，不如使用VTK框架[35]中实现的基于模板的方法来进行复杂多边形几何体的体素化。具体来说，使用vtkPolyDataToImageStencil算法将多边形表面转换为定义预定义体素网格内外区域的空间掩模。然后将得到的模板应用于空的vtkImageData体积，通过vtkImageStencilToImage生成二值体素场，其中每个体素值表示它位于几何体的内部还是外部。这种方法确保了几何保真性和稳健的内部-外部分类，为后续分析提供了准确的体积表示。然而，为了实现这种方法，必须适当指定体素分辨率，以准确捕捉纤维表面的主要特征，通常优选较高的体素分辨率。

为了减轻高分辨率体素网格相关的计算成本，在孔隙度计算的第二阶段，采用了一种基于域分解和共享内存并行性的并行体素计数策略，以加速大型三维数据集上的体积计算。具体来说，离散的三维体素空间表示为。引入了一个指示函数，以区分域内的纤维相和孔隙相：
（27）
利用共享内存并行性，沿平面方向（例如y轴）将域分解为M个连续且不重叠的子域，即。通过汇总每个子域内独立评估的局部体素计数，即对中的所有体素求和，获得整体固体纤维体积。假设每个 fiber 层有 层体素，通过汇总每个子域内独立评估的局部体素计数，获得整体固体体积：
（28）
因此，层孔隙度和整体孔隙度定义为：
（29）

在实践中，这种策略是在Python环境中使用joblib库和基于线程的并行性实现的。所有并行工作者直接操作表示体素域的共享NumPy数组，从而避免了内存重复并最小化了通信开销。这种共享内存并行化框架使得大规模三维数据集的工作负载分配和性能可扩展。一旦激活基于体素的方法，就会利用不断更新的体素表示，并与随后生成的纤维合并。通过提出的混合方法，重建过程显著加速，同时保持了高精度，这将在结果和讨论部分详细讨论。

2.3.2 纤维平移和旋转
纤维堆叠过程的另一个关键方面是控制纤维在平面内和平面内的取向。纤维的取向强烈影响纤维结构的固有各向异性特性。实际上，堆叠的纤维完全平行于一个平面（即本研究中考虑的平面内方向）是不现实的。因此，一定程度的平面内取向是不可避免的。然而，在之前的逐层方法中，这种效应经常被忽略，因为忽略它可以简化平面内孔隙度分布的控制。因此，本研究首次提出了一种两阶段纤维堆叠过程，包括初始随机堆叠和随后的逐层堆叠阶段。需要注意的是，这两个阶段与之前的孔隙度计算阶段是独立的。在第一个纤维堆叠阶段，激活了随机堆叠策略，并应用了Stoyan等人[36]提出的概率密度函数，以确保纤维的透射平面方向角度落在指定的合理范围内，如下所示（30）。这里的参数用于调整纤维的各向异性。当参数为1时，所有纤维都与透射平面方向对齐；而当参数大于1时，更多的纤维会与平面方向对齐。该角度是z轴与纤维中心轴矢量之间的夹角，具体取决于纤维中心轴的端点，其范围在[...]之间，不同参数下的概率密度函数分布图可见于图4。下载：下载高分辨率图像（182KB）下载：下载全尺寸图像图4。在不同各向异性参数（= 1, 3, 5, 7, 9, 11）下，纤维透射平面方向的概率分布函数。此外，在早期步骤中已经确定了纤维在平面方向上的 orientation。纤维的方向可以视为相对于x轴或y轴的范围[0,...]。由于纤维的x和y坐标已经给出，为了随机堆叠它们，我们还需要调节起始点的z坐标，而终点的z坐标将相应地根据给定的参数进行调整（31）。为了在孔隙率约束下加速纤维堆叠过程，提出了一种基于孔隙率差异加权的策略来确定起始位置。在每生成一层纤维后，会更新实际的层孔隙率分布，并计算每层实际孔隙率与目标孔隙率之间的差异；然后对所有层进行排序。选择孔隙率偏差最大的层，并在其z方向上的上下限范围内随机分配纤维。另外，一旦某层达到目标孔隙率，该层将被“冻结”，即任何可能导致新纤维穿过该层的随机选定的方向角都会被拒绝。这个角度选择过程会重复进行，直到找到一个可行的方向。随着更多层被冻结，拒绝候选角度的概率会增加。当达到指定比例的层（）达到其目标孔隙率值时，就过渡到逐层堆叠阶段。在本研究中，这个比例为0.65。这个阈值确保了初始65%的结构保持其固有的随机各向异性，而剩余的35%则通过严格的透射平面方向控制来快速消除剩余的孔隙率偏差，而不会产生过多的计算拒绝。在这个后期阶段，由公式（30）控制的透射平面方向逐渐限制在更小的角度范围内，以确保随后生成的纤维优先填充剩余的未充分发展的层。最后，当全局和层的孔隙率都在规定的容忍范围内收敛时，堆叠过程终止。重建的几何结构随后被导出为立体光刻文件，然后是最终的体素数据，这些数据可以直接用于后续的图像处理，例如增材处理。

2.4 增材处理
通过上述两个主要步骤，可以生成纤维骨架。然而，在实际应用中，这类材料通常会经过二次处理或化学处理以增强特定的性能特性，如机械强度或化学功能。这些后处理不可避免地会在多个尺度上引入或放大不规则的结构特征。一个显著的例子是应用化学试剂或涂层，这通常会导致纤维表面上的非均匀吸附、沉积或包裹。这种异质性导致表面形态和化学成分的局部变化，通常表现为不均匀的颗粒聚集或不连续的薄膜层。一个典型的例子是在质子交换膜燃料电池中对纤维基材（GDL）进行后处理，其中分别添加聚合物粘合剂和聚四氟乙烯（PTFE）以提高机械强度和疏水性。粘合剂，例如碳化树脂，通过形成碳纤维之间的刚性固定连接来提高压缩强度和结构完整性，通常通过高温固化或石墨化来稳定。PTFE处理是一个关键的疏水化过程，包括将GDL基底浸渍在水性PTFE乳液中，然后在327°C以上进行烧结。这个热处理步骤促使亚微米级的PTFE颗粒聚集和融合，在碳纤维上形成低表面能的包覆层，接触角通常超过120度。已经确认粘合剂结构主要聚集在纤维交叉且紧密堆积的区域[22]。为了重建这种结构，常见的方法是基于图像形态学处理[37]，例如闭合或开放。在本研究中，采用了与之前基于体素数据的方法相结合的方式。在引入粘合剂结构时，相应的也会改变孔隙率。没有使用三维形态学处理方法，而是采用了二维方法，沿着x、y和z方向依次处理体素切片，使用基于圆盘形结构元素进行闭合处理。在二值化的体素数据集中，1表示固体纤维相，0表示孔隙空间。圆盘元素的半径最初设置为1个体素。在对每个切片进行形态学闭合操作后，会重新计算局部和全局孔隙率。这种程序提供了控制重建结构孔隙率的灵活性（或域中粘合剂的量）。在所有三个方向上的切片都迭代完成后，如果未达到目标孔隙率，则将结构元素的半径增加一个体素，并重复上述过程。更多细节可以在我们之前的研究[38]中找到。与粘合剂结构不同，目前用于PTFE重建的方法表现出很大的多样性，没有标准化的方法，这可能是由于目前关于GDL中PTFE结构部分的研究非常有限。此外，只有少数研究描述了他们的PTFE重建方法，这对后续的复制带来了挑战。例如，Daino等人[22]通过对现有的粘合剂和纤维应用形态学闭合操作来生成PTFE结构。然而，这种方法可能导致PTFE无法牢固地附着在纤维上，从而影响模型的物理真实性。Molaeimanesh等人[36]采用了一种不同的方法，在粘合剂和纤维的表面随机生成特定直径的PTFE颗粒，这些颗粒随机堆叠直到达到目标PTFE含量。受到参考文献[36]中方法的启发，本研究提出并实现了一种基于局部几何竞争和动态形态学修复的并行随机生长算法，以准确重建GDL中PTFE的非均匀聚集分布。这种方法的核心逻辑摒弃了传统的均匀厚度假设，转而采用受毛细力启发的随机选择机制。首先，算法通过利用每个空洞体素的局部固体邻居计数作为评估指标来建立评分系统。在这个系统中，纤维交叉点和狭窄的裂缝被赋予更高的成核概率，因为它们的邻居占用率较高（即高曲率特性），从而模仿了乳液在浸渍过程中倾向于在微孔喉部优先积聚的物理倾向。为了进一步增强PTFE分布的不均匀特性，引入了对候选体素采样的严格限制，确保生长仅发生在具有极端几何优势的位置。此外，为了模拟干燥过程中由表面张力驱动的聚集和固化效应，将动态半径形态学闭合机制集成到迭代循环中。该机制在给定的元素半径变化范围内自主演化结构元素半径，有效地修复PTFE簇内的内部随机缺陷，同时避免不必要的全局均匀扩散。最后，为了应对大规模体素模型带来的计算压力，该算法采用了基于幽灵单元重叠的并行化卷积框架，显著提高了动态孔隙率评估和邻居检索的效率。需要注意的是，这种方法需要足够小的体素分辨率以避免厚涂层。

Niu等人[19]发现添加剂的分布受到干燥方法的影响。空气干燥可能导致添加剂的分布比真空干燥方法更加不均匀。为了证明添加剂分布的可控性，在该框架中实现了两种不同的透射平面生成策略，即均匀分布和非均匀分布。均匀分布是通过在整个计算域上应用不受约束的形态学闭合操作来实现的，使得固体添加剂在纤维交叉处均匀积累。相反，非均匀分布旨在模仿制造干燥过程中添加剂在毛细力驱动下向GDL边界侧的物理迁移。这是通过将形态学操作与空间概率密度函数耦合在计算上实现的。具体来说，在迭代生成过程中引入了蒙特卡洛过滤机制，新形成的添加剂体素只有在满足预定义的双峰轮廓定义的局部概率约束时才被保留。

3 结果与讨论
本研究的目的是提供一个更灵活的框架来重建GDL，以匹配其物理结构。为了与燃料电池应用中GDL的特征纤维尺度相匹配，采用了1微米作为基本长度单位。

3.1 与现有纤维堆叠方法的比较
在本节中，将混合纤维堆叠策略与之前的随机堆叠和逐层堆叠策略在控制能力和准确性方面进行了比较。此外，还评估了考虑和不考虑纤维透射平面方向的随机堆叠。图5显示了四种方法之间的比较。每种堆叠策略重复五次，以使用提出的孔隙率计算方法生成五个独立样本。这里是用于比较的默认重建设置。纤维中心线使用贝塞尔曲线插值，而纤维横截面被建模为一个直径为10微米的圆。纤维生成域被定义为...。体素长度分辨率设置为每个体素...。目标纤维端点距离从一个截断的高斯分布中采样，该分布被截断到区间...。为了更好地比较结果，曲率偏移因子被设置为常数值0.2，同时...= 150。纤维各向异性使用参数...进行控制。在纤维堆叠过程中，实现目标层状孔隙率和全局孔隙率的容忍度分别设定为0.004和0.002。对于涉及逐层堆叠的情况，透射平面维度被划分为22层。图6显示了逐层堆叠和混合堆叠策略中使用的透射平面孔隙率分布，这与之前研究[39]、[40]中显示的GDL孔隙率分布特征一致。本节没有考虑增材处理，这将在后面展示。

下载：下载高分辨率图像（681KB）下载：下载全尺寸图像
图5. 四种纤维堆叠算法的算法比较流程图，明确对比了它们在透射平面空间约束（Z轴）和平面外方向（）方面的不同生成逻辑。

下载：下载高分辨率图像（109KB）下载：下载全尺寸图像
图6. 默认的透射平面层状孔隙率分布。每个数据点对应于每一层的中间位置。
如图7(a–b)所示，两种随机堆叠策略在五个独立样本中的透射平面孔隙率上表现出显著的变异性。这种明显的随机波动表明纤维的随机放置导致了较高的结构不确定性，表明结果的代表性较差。从统计角度来看，这种变异性意味着单个微观结构实现不足以捕捉材料的整体属性，需要一个显著更大的代表性基本体积或大量样本的集合平均值才能达到收敛。然而，引入纤维的透射平面方向作为一个关键的稳定因素。如图7(b)所示，加入这些有方向的纤维有效地平滑了沿透射平面方向的局部孔隙率梯度。这种平滑效果可能是由于倾斜的纤维 bridged 层与层之间的间隙，减少了纯水平随机分布中常见的大局部空洞的形成。

下载：下载高分辨率图像（680KB）下载：下载全尺寸图像
图7.四种不同纤维堆叠策略的透射面孔隙率分布：(a) 无透射面纤维方向的随机堆叠；(b) 有透射面纤维方向的随机堆叠；(c) 层层堆叠；以及 (d) 混合堆叠。相比之下，图7(c)中展示的层层堆叠策略表现出不同的权衡。它在所有五个样本的局部孔隙率分布上实现了显著的一致性，大大降低了样本间的方差。然而，在目标孔隙率分布周围引入了显著且持续的振荡，最大振荡幅度达到了0.36。这些周期性波动是离散层叠过程的固有现象，其中孔隙率在界面处达到峰值，在纤维层中心处达到最低值。虽然这些振荡是一致的，但它们可能导致对局部传输现象的预测不准确。最终，图7(d)中展示的混合堆叠策略被认为是最佳方法。通过结合层叠的结构规则性和定向纤维的平滑效果，它实现了高度一致的分布，同时大大减弱了有害的振荡。这种双重优势展示了更高级的微观结构代表性，为后续的物理建模和性能预测提供了更可靠的基础。总体而言，无透射面纤维方向的随机堆叠具有较差的代表性。有透射面纤维方向的随机堆叠虽然平滑了梯度，但仍然具有高度的随机性。层层堆叠提供了样本一致性，但引入了人为的物理振荡。相比之下，混合方法独特地解决了这些冲突，同时实现了一致性和物理真实性。图8展示了通过四种不同纤维堆叠策略生成的初始样本的横截面视图。这些可视化结果沿着两个主要正交平面捕获：纵向中平面和横向中心平面，提供了对内部微观结构的双重视角分析。关于横截面，视觉证据清楚地表明，规定的纤维方向是塑造内部结构的关键因素。它不仅决定了空间孔隙率分布，还从根本上改变了局部孔隙形态。例如，某些堆叠序列导致高度拉长的间隙空隙，而其他序列则产生更加各向同性的、分隔的孔隙。这种形态的变化至关重要，因为它直接影响材料的渗透性和局部应力集中区。同时，横截面分析显示样本间的孔隙率梯度存在显著差异。这些梯度并非偶然；它们源自每种堆叠策略特定的纤维分配逻辑。一些策略在整个深度上促进均匀的纤维密度，而其他策略则由于其层层沉积特性而引入了有意的或系统性的纤维浓度变化。图8中的这些见解强调了如何利用堆叠策略的选择来调整复合材料的功能梯度。

图8. 在四种不同纤维堆叠策略下，z = 115和y = 250平面内的横截面图片：(a) 无透射面纤维方向的随机堆叠；(b) 有透射面纤维方向的随机堆叠；(c) 层层堆叠；以及 (d) 混合堆叠。

3.2. 透射面孔隙率分布
最近，梯度孔隙率GDL结构在燃料电池研究领域受到了广泛关注[41][42]。这种兴趣源于一个基本认识，即传输属性（包括气体渗透性、有效扩散率和液水毛细压）与局部孔隙结构密切相关。通过调整孔隙率分布，可以优化质量传输和欧姆电阻之间的平衡。在这项研究中，我们使用了一种为精细控制局部形态而开发的随机重建算法，生成了三种不同的GDL变体来评估我们建模方法的稳健性。这些变体包括均匀分布、线性梯度分布和非单调的中低（U形）分布，如图9所示。此外，模拟的透射面孔隙率分布与GDL的相对厚度进行了对比。目标（虚线）和实际（实线）孔隙率之间的严格比较表明，重建的几何形状与预期分布高度匹配。即使对于更复杂的“中低”分布，算法也能很好地捕捉到从高孔隙率区域到密集核心的过渡。观察到的局部振荡幅度很小，最大为0.08。这种振荡源于纤维放置的离散性和空间重叠。与前面提到的两种纤维堆叠策略相比，整体上对目标曲线的遵循情况证明了重建过程的高保真度和可靠性。这种能力确保了生成的微观结构可以作为未来孔隙尺度传输模拟的可靠基础。

图9. 具有三种不同透射面孔隙率的重建结构。

3.3. 孔隙率计算成本和准确性
在本节中，我们将比较三种方法的计算效率和数值准确性，即纯表面三角剖分体积法、纯体素基方法和提出的混合方法。为了建立一个可靠的基准，所有方法都在相同的设置下对一个代表性的GDL微观结构进行了五次随机重建测试。所有计算都在配备有Intel Core i7-14700 CPU（2.10 GHz，8个性能核心和12个高效核心）和32 GB RAM的台式计算机上进行，而并行计算仅限于使用16个逻辑线程的8个性能核心。为了比较三种方法，从表面三角剖分计算的初始孔隙率使用体素基方法进行校正，以获得重建结构的最终全局孔隙率。分辨率为0.5/体素。

选择了三个目标全局孔隙率：0.65、0.75和0.85，分别对应于均匀的层孔隙率分布。在每种孔隙率平均五个重建结果中，混合方法与纯体素基方法的精度相匹配（两者相对于目标的绝对误差均低于0.0007），同时处理时间分别减少了47.69%（5分钟18秒对比10分钟8秒）、45.1%（3分钟45秒对比6分钟50秒）和44%（2分钟10秒对比3分钟52秒）。相比之下，基于表面三角剖分的方法虽然更快（4分钟28秒、2分钟50秒和1分钟52秒），但其结果（0.651、0.752和0.852）与体素化真实值（0.709、0.766和0.862）有显著偏差，这突显了其在GDL建模中的不足准确性。这种差异表明，提出的混合方法有效地平衡了计算效率和几何准确性。通过降低目标全局孔隙率并提高体素分辨率，三种方法在孔隙率准确性和效率方面的差异变得更大。在不需要精确孔隙率计算的情况下，仅使用第一阶段就足够快速重建。为了直观展示算法的优越性并提供明确的事实依据，图10总结了三种评估方法的定量性能。如图所示，混合方法成功地保持了计算量大的体素方法的严格几何准确性，同时显著减少了所有测试孔隙率范围内的总模拟时间。

图10. 三种孔隙率计算算法在不同目标孔隙率下的性能比较。(a) 计算时间，突出显示混合方法相比纯体素基方法提高了45%的效率。(b) 几何准确性，表明混合方法保持了严格的精度，而表面三角剖分方法则表现出显著偏差。

3.4. 参数化结构重建
为了展示所提出的随机重建框架的多功能性和稳健性，我们通过改变关键几何参数系统地生成了一系列GDL微观结构。这些参数代表GDL制造中最关键的自由度，并显著影响最终的传输性能。通过隔离这些因素，该模型提供了一个经济高效的虚拟实验室，用于研究GDL性能对其内部形态的敏感性。这些参数变化的综合结果汇总在图6中，展示了所提算法处理多种纤维几何形状的能力，包括横截面形状、直径、曲率和长度。尽管大多数GDL的随机重建假设纤维是圆柱形的以简化问题，但真实的碳纤维由于制造过程通常具有非圆形横截面。这种简化可能会忽略对传输和机械性能的关键微观结构影响。在这里，我们提出了一个通用的重建框架，该框架能够纳入不同的纤维横截面，从而系统地研究横截面形态如何影响GDL以及其他纤维材料的多物理性能[43]。在图11(a)中，生成了三种具有不同纤维横截面形状的纤维结构，即圆形、矩形和三角形，所有三种结构的中心到顶点的距离（半径）均为10。这种能力提供了研究边缘锐度和横截面形状如何影响局部热和质量传输的可能性。

图11. 随机重建的GDL微观结构，从上到下排列不同的形态参数，即纤维横截面形状（圆形、矩形和三角形）、纤维直径（=5、7和10）、纤维曲率（从截断的高斯分布中采样，均值分别为0、0.6和1.2，以及相同的变体0.2）和纤维长度（从截断的高斯分布中采样预期的终点距离，均值分别为50、350和580，以及相同的变体20）。此外，纤维直径在GDL[44]以及其他纤维材料（如静电纺丝制成的材料[45]）中已被广泛研究。图11(b)显示了随着纤维直径从5变化到10的重建结果。可以看出，在恒定孔隙率下，较小的直径会导致更高的比表面积和更复杂的、曲折的孔隙路径，这对于表征GDL内的气体扩散和毛细压至关重要。需要注意的是，在这些特定参数情况下使用统一直径是为了严格隔离纤维厚度的影响，但所提出的框架本身支持在单个生成的结构中分配不同直径。通过在初始化步骤中从实验测量的直径分布中采样，模型可以轻松适应高度异质的 porous 媒体，在这些介质中纤维直径表现出显著的自然变化。

正如我们之前讨论的，纤维曲率是GDL中普遍存在但量化不足的形态特征，它很少被纳入随机重建模型中。我们的框架能够明确包含曲率作为随机变量，从而系统地研究其对传输、机械和电化学性能的影响。图11(c)展示了三种具有不同纤维曲率的纤维结构，由不同的（=0.2、0.6和1）和相同的（=150）生成。在GDL制造过程中，纤维长度主要通过切割过程在纤维堆叠之前进行控制[46]。此外，实验观察到在某个压缩应力后纤维会断裂[47]，这可能导致纤维长度分布的变化。除了GDL外，纤维长度在其他增强纤维聚合物复合材料中也是一个重要因素[48]。因此，我们使用开发的方法生成了具有不同纤维长度特性的纤维结构，如图11(d)所示。通过从三个截断的高斯分布中采样预期的纤维终点距离，生成了三种纤维结构，分别如图12所示。结果，中心轴的距离可以用来间接控制弯曲的纤维长度。具体来说，对于直纤维，终点距离分布就是纤维长度分布。

图12. 纤维终点距离的高斯分布函数，三个平均值分别为50、350和580，共同的标准偏差值为20。

3.5. 添加剂处理
在GDL制造过程中，通常会使用粘合剂处理来增强机械强度并调节传输特性。为了方便对这些添加剂的影响进行稳健的研究，我们开发的框架包含了一个专门的粘合剂相集成模块。图13展示了从具有基线体积孔隙率0.85的相同纤维骨架重建的形态，这些骨架系统地改变了粘合剂体积分数()。这种受控的比较突显了所开发算法在精确调节添加剂负荷的同时保持初级纤维网络结构一致性的能力。随着粘合剂含量的增加，可以观察到孔隙几何形状和相互连接性的渐进性变化，这通过粘合剂在纤维接头处的积累得到了证明。下载：下载高分辨率图像（352KB）下载：下载全尺寸图像 图13. 使用相同的GDL纤维骨架重建的结构，其体积孔隙率为0.85，但添加剂比例系统性地变化，分别为：(a) 0，(b) 0.05，以及 (c) 0.1。图14(a)和(b)分别展示了纤维-粘合剂二元系统和纤维-粘合剂-PTFE三元系统的3D重建形态。尽管两种结构都保持了一致的纤维网络拓扑结构，但PTFE的引入带来了明显的形态变化，尤其是在黄色虚线椭圆标出的区域。具体来说，PTFE相表现出双重分布特性：它要么作为薄层均匀包裹现有的纤维-粘合剂骨架，要么选择性地填充间隙微孔。这种选择性沉积有效改变了纤维基体的局部表面粗糙度和孔隙喉部几何形状。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像 图14. 使用相同的纤维和粘合剂体积分数但不同PTFE体积分数重建的GDL结构：(a) 0，(b) 0.05。(c) 包括PTFE相的重建结构截面视图。为了更细致地检查内部相的排列情况，图14(c)提供了三元系统的正交截面视图。在这些切片中，每个组成相都被赋予了不同的颜色，以方便观察其空间关联性和界面边界。值得注意的是，PTFE主要聚集在高曲率区域和狭窄的孔隙喉部。这种优先积累并非随机现象，而是与润湿性添加剂在浸渍或涂层过程中的物理行为一致，其中毛细力将材料驱动到受限的几何形状中。这种对相分布的真实再现展示了该框架模拟复杂多相相互作用的能力，为后续的运输性质模拟提供了可靠的几何基础。图15显示了添加剂处理前后重建GDL的平面穿透孔隙率分布。为了与参考文献[19]中报告的特征一致，比较了两种不同的添加剂分布：一种是相对均匀的分布，另一种是异质的分布。与参考文献[19]中的研究类似，仅基于GDL纤维骨架重建PTFE，以清楚地展示对添加剂分布的精确控制。如图所示，在均匀分布的情况下，孔隙率在整个厚度上相对均匀地降低，表明添加剂自然填充了孔隙喉部；相比之下，异质策略产生了特征性的U形平面穿透孔隙率分布。这种分布成功地限制了核心区域的过量添加剂积累，同时促进了顶部和底部边界附近的局部富集。这种分布完美地再现了在蒸发和干燥过程中由毛细力驱动的粘合剂分离的物理现实，同时严格保持了规定的目标体积孔隙率。下载：下载高分辨率图像（416KB）下载：下载全尺寸图像 图15. 添加剂处理前后平面穿透分布的比较。(a) 裸露GDL纤维骨架（黑色点划线）、添加均匀添加剂的GDL（蓝色实线）以及添加异质添加剂的GDL（绿色虚线）的平面穿透孔隙率分布，并根据实验参考数据[11]（红色圆圈）进行了验证。(b) 均匀和异质分布策略下沿相对厚度的相应添加剂比例，说明了粘合剂向边界表面的物理U形分离。关于图中颜色参考的解释，请参阅本文的网络版本。3.6. 重建GDL运输性质的评估 随机重建算法的准确性主要是通过将微观结构特征与物理GDL样本进行比较来验证的。图16(a)展示了重建GDL的孔径分布与商业Freudenberg H2315 GDL [49]的实验结果进行对比。两种基材具有相同的宏观属性，尺寸为500×500×223，体积孔隙率为...，平均纤维直径为...。原始数据源[49]中没有局部平面穿透孔隙率曲线，因此采用了均匀分布来匹配提供的体积孔隙率。这种方法与参考文献[15]中的实验观察结果一致，这些观察结果表明商业H2315 GDL样本具有几乎均匀的平面穿透孔隙率分布。考虑到Freudenberg H2315 GDL中纤维曲率的极端随机性，在重建中包括了变化范围内的这一特征。观察到重建的PSD捕捉到了特征性的单峰分布，其中一个主导的孔径峰值约为8-10微米，其中重建GDL显示出8.5%更高的概率密度。它们的总体相对偏差小于5%。这一一致性证实了纳入纤维曲率参数()使得算法能够再现现实中碳纤维复杂的随机缠绕。尾部区域的轻微偏差()可以归因于重建中使用的均匀局部孔隙率假设，这在参考的物理GDL中是未知的。尽管如此，统计一致性验证了生成微观结构对后续运输模拟的代表性质。下载：下载高分辨率图像（806KB）下载：下载全尺寸图像 图16. 重建GDL的形态和运输性质的评估。(a) 重建微观结构和Freudenberg H2315 GDL [49]实验数据之间的孔径分布(PSD)比较。插图显示了随机纤维网络。(b) 平面渗透性()随孔隙率的演变，参照Bruggeman [50]和Tomadakis-Sotirchos [51]模型进行了基准测试。(c) 与LBM模拟[23]、实验测量[52],[53],[54]以及分析相关性[55],[56]进行的内在平面渗透性()验证。(d) 与各种经验模型[50],[51],[57],[58],[59],[60]和液压曲折度近似()相比的标准化平面有效扩散率()。误差条表示来自五个独立随机重建的标准偏差。图(b-d)中展示的模拟结果是使用相同重建参数生成的五个独立GDL样本的集合平均值。为了进一步验证所提出的重建框架，通过使用OpenFOAM软件进行计算流体动力学模拟来评估宏观运输性质，特别是渗透率和扩散率。所使用的GDL结构是根据...和...生成的，结果与既定文献数据进行了基准测试。模拟在0.55到0.95的广泛孔隙率范围内进行。由于模拟尺度的适用性限制，省略了较低体积孔隙率的案例。在这些孔隙率范围内，特征孔径大小与网格分辨率限制相当，使得基于连续体的模拟方法可能不适合准确解析控制物理过程。此外，在低孔隙率水平下，增强的孔壁相互作用和明显的瓶颈效应显著主导了运输阻力，需要更复杂的子网格尺度建模。为了确保统计独立性并最小化重建过程中的随机不确定性，图16(b)-(d)中呈现的每个数据点代表使用不同随机种子生成的五个独立重建的集合平均平面结果。误差条（可见时）表示标准偏差，表明在指定约束范围内的结构变异性。这些误差条的狭窄范围进一步证明了所提出算法在再现具有一致运输特性的代表性基本体积方面的稳定性。附录B提供了数值模型和参数提取方法的详细描述。平面内结果和相关模拟结果已附加在支持信息中。图16(b)中的曲折度分析进一步阐明了控制流动阻力的潜在结构复杂性。从体积平均速度场得出的液压曲折度()用于比较。这种差异突出了流体包裹体受到动态限制的物理现实。模拟结果与两种经典分析模型进行了比较，即Bruggeman相关性() [50]和Tomadakis和Sotirchos提出的随机纤维模型[51]。如图所示，我们的重建GDL与这些相关性存在偏差，尽管它们都显示出随着孔隙率增加而减少的趋势。这一总体趋势是由关键孔隙喉部的逐渐扩展和固相堵塞的减少所驱动的，有效地减轻了粘性阻力。与Bruggeman模型的偏离是预期的，因为该模型最初是为各向同性的球形颗粒堆积提出的。GDL的纤维性质导致了与球体相比显著不同的障碍形状和连通性，从而导致了更高的各向异性阻力。关于与Tomadakis和Sotirchos模型的比较，观察到的差异可以归因于我们重建的具体形态特征。虽然Tomadakis模型假设了具有特定重叠概率的无限长、直立的圆柱形纤维，但我们的重建考虑了纤维曲率和各向异性。这些有限长度和曲率的纤维特性可能创造了与理论预测的无限长直纤维不同的局部连接路径，特别是在高孔隙率范围内，其中“通道效应”变得更加明显。从机制上讲，这种通道效应允许流体优先通过阻力最小的路径，从而显著改变了宏观液压曲折度。随着...的收敛，这些值趋近于 unity，证实了在高度多孔介质中结构障碍的消失。为了评估动量传输能力，在0.55到0.95的孔隙率范围内评估了内在渗透性()。如图16(c)所示，模拟的渗透率随着孔隙率的增加而呈指数级增长。与孔隙空间的简单线性扩展不同，这种流体传导的指数级飞跃是由于孤立空隙随着纤维密度的降低逐渐合并成连续的、低阻力的宏观通道。这一趋势反映了液压通道的快速扩展和粘性阻力的减少，这与Tomadakis和Robertson [55]提出的粘性阻力模型很好地吻合。我们的重建准确性与Gao等人[23]的Lattice Boltzmann模拟数据以及Becker等人[53]、Kirsch等人[52]和Tamayol等人[54]对商用Toray TGP-H GDL的实验测量结果的高度一致性进一步得到了证实。关于我们的结果与Shou等人[56]的分析模型之间的差异，这些偏差可以归因于我们随机模型中复杂的各向异性纤维排列，这种模型在理想化的分析推导中经常被简化。理想化模型通常简化了纤维排列，而我们的随机重建捕捉了由随机纤维交叉和不同方向引起的微观尺度流动分离和局部形式的阻力。这些复杂的几何特征作为额外的动量吸收体，导致了与高度理想化预测的观察偏差。图16(d)展示了作为孔隙率函数的归一化有效扩散率()。为了评估经典近似关系的有效性，同时绘制了这一指标进行比较。有趣的是，观察到了一种交叉行为，即在低孔隙率时略微高估了扩散率，而在孔隙率增加时低估了它。这种差异可以归因于液压传输机制和扩散传输机制之间的根本差异。实际的质量传输不仅受到路径长度延长（曲折度）的阻碍，还受到孔隙截面面积变化的“瓶颈效应”的显著影响，这一现象通过解决浓度扩散方程得到了严格捕捉，但在液压曲折度中只是部分反映。与压力驱动的液压流动不同，后者在通过狭窄部分时动态加速，而浓度驱动的分子扩散依赖于在这些微观结构陷阱中的随机游走，这不可避免地导致了偏差。尽管存在这些细微差异，模拟结果仍然与现有的经验相关性[50],[51],[57],[58]和之前的数值研究[59]一致。特别是，与Shou等人[60]提出的分析模型的强一致性证实了重建GDL非常适合预测PEMFC应用中的气体物种传输。值得注意的是，这些传输性质在机械压缩下可能会进一步退化，正如Rashapov等人[61]所证明的。图17展示了在重建GDL内随着液水饱和度()的变化而变化的两相相对渗透率()。通过使用OpenPNM [62]中的孔隙网络建模(PNM)获得了和之间的关系。为了比较[63]中的结果，我们使用了提供的体积孔隙率（0.79）和纤维直径（7.5）。由于数据有限，假设穿透平面的孔隙分布是均匀的。为了全面评估微观结构的随机性对流体传输特性的影响，结果是在五个独立生成的随机样本上进行了集合平均。在图中，散点（非润湿水相的蓝色实心三角形和润湿空气相的红色空心圆圈）代表这些样本的统计平均值，而相应的误差条表示由局部形态波动引起的标准偏差。随着液态水饱和度的增加，水相的相对渗透性表现出典型的非线性增长趋势，而空气相的相对渗透性由于孔隙通道被水逐渐占据而显著减弱。从多相流动的角度来看，因为水在疏水性气体扩散层（GDL）环境中是是非润湿相，毛细力决定了水优先侵入最大、阻力最小的孔隙空间。这种选择性侵入迅速阻塞了空气传输的主要宏观通道，导致空气渗透性急剧下降。同时，水渗透性的非线性激增从根本上受到渗透力学的影响；孤立的液体团簇必须首先聚集并合并成连续的、贯穿系统的网络，然后才能发生大量液体流动。相对较窄的误差带表明，在指定的孔隙率和特征参数下，重建算法生成的微观拓扑具有优异的均匀性和统计代表性。

下载：下载高分辨率图像（263KB）
下载：下载全尺寸图像

图17. 重建的GDLs的相对渗透性（）与液态水饱和度（）的比较（体积孔隙率为0.79，纤维直径为8）。带有误差条的符号代表五种随机实现的孔隙网络建模（PNM）预测的统计平均值和标准偏差（水为蓝色三角形，空气为红色空心圆圈）。这些结果通过LBM模拟（水为青色星号，空气为绿色正方形）和经典的经验幂律模型（用点划线表示）进行了验证。（关于此图例中颜色的解释，请参阅本文的网页版本。）

PNM模拟结果与先前的LBM结果[63]和经典的经验幂律模型[ ]进行了比较。与理想的经验模型不同，后者假设在绝对零度或完全饱和之前流动路径是连续的，模拟结果揭示了复杂纤维多孔介质特有的非线性传输特性。润湿相（空气）的相对渗透性在液态水饱和度达到1.0之前就接近于零。这种早期截断现象在物理上代表了“断开”机制，即空气相的传输路径被侵入的液体在水中的复杂多孔骨架内阻挡，进一步影响了电化学反应。更显著的是，非润湿相（水）的相对渗透性在低至中等饱和度时显著超过了经典的立方理论模型[ ]，这与LBM结果相似。这种正向偏差在物理上表明在GDL内部发生了优先流动通道现象，其中重叠的纤维和狭窄的通道显著增加了传输阻力和毛细屏障。由于GDL具有高度疏水性，侵入的液态水遵循最小毛细阻力的路径，选择性地占据了最大的连通孔隙（即毛细指形）。因此，建立了一条连续的、高导性的通道，所需流体体积最小（低），从而导致相对渗透性显著上升。这表明，与理想的各向同性多孔介质相比，GDL的异质纤维拓扑促进了非润湿流体的更有效渗透。PNM预测与LBM模拟在成功捕捉这些关键机制（即断开、连接限制和优先流动通道）方面高度一致，不仅验证了重建GDL形态的可靠性，还为后续的水管理优化提供了坚实的理论基础。

3.7. 扩展到其他纤维材料的可行性
尽管所提出的重建框架主要是针对GDLs进行验证的，但其底层数学算法几乎是材料无关的。这种灵活性为将该框架扩展到更广泛的纤维多孔材料提供了可行的技术基础，例如静电纺丝微/纳米纤维垫、玻璃棉和生物支架。虽然GDLs和静电纺丝垫在基本拓扑上具有相似性，都是高度多孔的、相互连接的纤维网络，但它们在形态上存在显著差异。GDLs通常由刚性的、微米级的碳纤维组成，这些碳纤维相对较短且略微弯曲。相比之下，静电纺丝垫的特点是超细的、连续的、高度柔性的纳米纤维，由于制造过程中的混沌摆动不稳定性而形成极其曲折的轨迹[64]。

这种扩展性的核心在于通用纤维生成方案的形态适应性。通过克服传统的直圆柱限制，我们的框架可以明确控制纤维的曲率、截面和长度分布。因此，它可以无缝模拟静电纺丝喷射和不规则非织造织物的高度曲折、连续蜿蜒的特性。此外，所提出的结构组装和评估策略对于极端结构情况下非常有利。通过结合混合堆叠策略和双阶段孔隙率计算，该框架可以轻松处理超高孔隙率的气凝胶或异常密集的过滤垫，有效降低传统基于体素的方法的严重计算成本和重复计算错误。将此模型应用于特定替代材料主要需要尺寸缩放（例如，按比例调整纤维直径值并将计算单位从微米改为纳米）和物理验证，考虑尺度依赖的传输现象（例如，克努森扩散）。然而，仅从结构生成的角度来看，建立的框架为多样化的纤维材料的数字孪生重建提供了高度可靠的蓝图。

4. 结论
本研究旨在开发一个用于重建多孔纤维结构的综合框架，特别是能量转换系统中的气体扩散层。该框架包括三个步骤，即纤维生成、堆叠和添加结构的整合。提出了一种综合的纤维生成方案，以提高重建框架的通用性和适应性，涵盖了从中心轴的数学定义到纤维表面构建的整个工作流程。通过采用各种已建立的纤维轴几何公式（例如，贝塞尔曲线、三次样条、B样条），这种方法能够精确且独立地调节关键形态参数，包括纤维曲率、厚度、方向和截面拓扑（例如，圆形、三角形或矩形）。此外，在纤维轴生成过程中，通过使用Sobol采样策略提高了纤维分布的空间均匀性。通过引入纤维端点间距离的概率密度函数，有效控制了纤维长度。

为了克服现有重建方法的局限性，引入了一种新颖的混合堆叠策略，将随机沉积与基于孔隙率加权的逐层填充相结合。这种方法解决了调节穿透平面孔隙梯度与保持真实纤维各向异性之间的权衡，有效消除了传统分层方法常出现的人为孔隙率振荡。此外，通过开发双阶段混合孔隙率计算方法，显著提高了重建过程的计算效率。通过结合分析几何估计与并行化的基于体素的校正，所提出的算法将计算时间大幅减少，与纯基于体素的方法相比大约减少了45%，同时保持了严格的数值精度。除了纤维骨架之外，该框架还通过整合先进的随机生长算法扩展到了多相重建，用于添加处理。结果模型真实地捕捉了粘合剂在纤维连接处的优先积累以及在高曲率区域中疏水性PTFE剂的异质分布。通过针对实验基准和经典分析模型的全面验证，确认了重建微结构的真实性，特别是在孔隙尺寸分布、渗透性和有效扩散性方面。最终，这项工作提供了一个强大且灵活的计算工具，用于建立定量结构-性能关系，为高性能燃料电池和水电解器的下一代电极材料的逆向设计和优化铺平了道路。

与现有的GDL重建方法相比，当前的框架包含了一组更全面的物理结构特征。尽管取得了这些进展，但仍存在某些局限性。首先，添加处理步骤目前受到高计算负载的限制，这是由于需要 fine 体素分辨率所致。为了解决这个问题并进一步改进该方法，未来的工作应侧重于整合通过不同商业GDL系列获得的真实结构特征信息。此外，纤维重叠仍然是随机重建过程的固有局限性。在真实材料中，固体纤维会相互偏转或压缩；因此，模拟的重叠特征人为改变了纤维之间的连接拓扑。虽然这在评估固相传输（例如，电导和热传导）和机械性能时是一个关键因素，但其对本文分析的流体传输特性的影响预计是有限的。正如先前研究[24]所示，重叠和非重叠的纤维结构具有相似的平均孔隙大小和整体连通性。在GDLs中，气体扩散主要受液体流动路径的影响，这些路径主要由大的疏水孔隙控制。添加剂的加入通过填充纤维交叉点的小孔隙进一步减少了重叠伪影。虽然通过随机游走算法或迭代距离检查生成严格非重叠结构在理论上是可能的，但在低孔隙率下计算成本很高。严格执行严格的不重叠约束也会严重限制所提出框架的核心优势，例如调节纤维形态和局部孔隙分布的灵活性。因此，对于流体传输特性的预测，当前框架在计算效率、形态灵活性和物理准确性之间达到了实用的必要平衡。尽管如此，为了提高纤维结构重建的现实性和普遍性，我们承认仍有改进的空间。

cRediT作者贡献声明
Danan Yang：撰写——原始草稿、可视化、验证、软件、方法论、形式分析、数据管理、概念化。
Himani Garg：撰写——审阅与编辑、监督。
Martin Andersson：撰写——审阅与编辑、监督、资源获取、资金获取。