质子交换膜燃料电池（proton exchange membrane fuel cell， PEMFC）是理想的能量转换装置，可以直接将氢气中的化学能转化为电能，在以新能源汽车为代表的移动领域应用前景广阔^[1-3]. 甲醇蒸汽重整制氢微反应器凭借紧凑、轻量化的结构与高能量密度，可以为PEMFC提供更为安全、高效、低成本的氢源，是有效的移动制氢解决方案^[4-7].

催化剂载体是微反应器中实现甲醇蒸汽重整制氢的关键构件，其结构设计是提升制氢性能的重要课题. 因此，新型微通道、泡沫金属、多孔金属纤维等新型催化剂载体层出不穷^[8-14]. 其中，华南理工大学Tang等^[15]结合切削和固相烧结法开发的多孔金属纤维结构工艺简单、比表面积高，表现出良好的催化反应性能（以下将其简称为纤维载体）. 更进一步，Zhou等^[16]开发了梯度纤维载体，在提升制氢性能的同时，通过引入梯度界面和孔隙率分布2个变量，提高其结构设计的自由度. 但该设计方式依赖于实验，缺乏明确的设计依据.

催化剂载体要求具有较好的传热传质能力^[17-18]，而这与载体中流体的流动分布密切相关. 流动不均会导致反应物停留时间不均，严重影响反应效率；均匀的体积流量分布有利于反应速率和目标产物选择性的提高. 因此均匀的流动分布是催化剂载体结构优化的重要目标^[19-22].

计算流体力学（CFD）仿真是研究微流道中流动分布的重要手段^{[9， 23-26]}，但其计算量大、求解速度慢. 为此，基于等效电阻网络的方法在微流道、微翅片结构的流动分析中得到了发展^{[21， 27-29]}. 该方法借鉴基尔霍夫定律，在建立微流道模型的基础上，将其中的体积流量和流动阻力分别等效为电流和电阻，并根据电压和电流守恒建立微流道的近似压降模型，进而对流速场快速求解. 尽管效果较好，但要将等效电阻网络方法应用于纤维载体仍然存在极大挑战. 原因在于，纤维载体的微流道结构是三维连通的，流道形状不规则、数目众多，长度和方向都具有随机性. 这使得纤维载体的流道形状难以描述、流道阻力难以确定，无法直接应用已有的等效电阻网络模型.

在随机分布纤维结构研究中，研究者们发展了一种统计分析方法，其核心是将纤维结构简化为规则的统计网络^[30-32]. 但该方法面向的是纤维实体结构，而非其中的孔隙结构. 为了发展面向纤维载体流速场快速分析的等效电阻网络模型， 本研究借鉴这一统计网络方法对纤维载体的流道结构进行简化；基于基尔霍夫定律，建立有效的等效电阻网络模型，进行高效求解；通过实例验证方法的有效性.

纤维载体由长度为10~20 mm、当量直径约为100 μm、随机分布的铜纤维段经模压和固相烧结工艺制得（见图1 （a）、（b））^[15]. 对于其内部复杂的孔隙通道（见图1（c）、（d））缺乏统一的描述方法；切削加工工艺使得铜纤维表面形貌丰富（见图1 （d）），但也使得纤维表面对流动阻力的影响难以量化.

纤维载体负载催化剂之后，即可用于组装微反应器. 如图2所示，附有催化剂的纤维载体被嵌于反应腔中；当微反应器工作时，甲醇与水的混合反应物经由蒸发腔充分蒸发，以气态形式流入反应腔，并在纤维载体中充分扩散，同时借助负载的催化剂发生重整制氢反应.

研究表明，反应物在纤维载体中的扩散速度、停留时间以及与催化剂的接触面积都会影响甲醇重整制氢反应效率^[27]. 因此，该区域的流动分布状态对微反应器性能具有重要影响. 鉴于纤维载体直接与进出口腔相连，本研究中纤维载体的整体流道结构包括其本身的流体微通道和进出口腔流道.

纤维随机分布的纤维结构内部纤维搭接点数N_c符合以下统计关系^[30-31]：

(1) ${N_{\rm{c}}} = \frac{1}{{2{\text{π}} S}}\left[ {1 + {{\rm{exp}}\;\;{({{{ - {N_{\rm{f}}}\lambda \omega } \mathord{\left/ {\vphantom {{ - {N_{\rm{f}}}\lambda \omega } S}} \right. } S}})}}} \right]N_{\rm{f}}^2{\lambda ^2}.$

式中：S为纤维结构的面积，λ为纤维的平均长度，ω为纤维的当量直径，N_f为纤维的根数.

基于式（1），可以将纤维随机分布的纤维结构简化为阶数为 $N_{\rm{c}}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}$的规则统计网络^[32]，如图3 （c）所示. 其中，网格节点之间是纤维段.

统计网络模型面向的是纤维实体结构，以其作为流道边界，即可得到纤维载体的内部流体微通道结构. 由于纤维载体在厚度方向是40 mm×70 mm的长方形，引入划分阶数n，并按纤维载体的宽长比，将其结构简化为如图4所示的4n×7n的规则纤维网络. 对于同一个纤维载体，其内部纤维搭接点数目应该是相同的，可以得到

(2) $ {N_{\rm{c}}} = {\rm{ }}\left( {4n + 1} \right)^2 . $

须注意的是，由于纤维载体的厚度（2 mm）远小于其长度和宽度，同时，其中的流动一般呈低速层流状态^[16]，以单层纤维结构的流道划分结果作为纤维载体的整体流道结构. 该层纤维结构处在纤维载体厚度方向中心面上. 对于该层纤维结构，本研究不加推导地给出式（1）中的纤维根数为

(3) $ {N_{\rm{f}}} = \frac{{4S\left( {1 - E} \right)}}{{{\text{π}} \omega \lambda }}. $

式中：E为纤维载体的孔隙率，可以按文献[15]计算.

显然，在如图4所示的流道结构中，每个网格单元都代表由纤维段包围形成的孔隙. 由于网格单元满足均匀性标准^[33]，将反应物在其中的流动阻力等价为孔隙中纤维段产生的平均阻力.

进出口腔流道的划分方法已经较为成熟^[27]，本研究亦采用相同的方法. 由于进出口腔具有中心对称性，以下仅以出口腔为例对该方法进行简要说明.

1）出口中心O在出口边界角∠MGN的角平分线GH上（见图5 （a））. 以点O为圆心、OH为半径作辅助圆，与出口腔的边界线交于点E、F. 此时，反应物从EH、FH这2条线扩散到点O受到的阻力近似相等，区域GEHF内的流速分布对纤维载体流速场不会造成直接影响，可以忽略其内部压降.

2）纤维载体的流道边界与出口腔边界之间存在一系列交点，假设其数量为m+1；加上点H，总计有m+2个点. 用A₀，A₁，···，A_m+1对这些点进行标记（见图5 （a））. 其中，点H的标记是A_b，下标m和纤维载体划分阶数n的关系为

(4) $ m = {\rm{ }}6n. $

3）对于点H左侧的点A_k，过点A_k作EH的平行线段，使其与出口腔边界相交于点B_k−1；对于点H右侧的点A_k，过点A_k作FH的平行线段，使其相交于点B_k，从而将出口腔划分为m+1条流道（见图5 （b））. 令点E为点B_b−1，点F为点B_b，则点H左侧第k条流道为 ${{{A}}_k}{{{A}}_{k{\rm{ + }}1}}{{{B}}_k}{{{B}}_{k{\rm{ - }}1}}$，右侧为 ${{{A}}_k}{{{A}}_{k{\rm{ + }}1}}{{{B}}_{k{\rm{ + }}1}}{{{B}}_k}$.

4）由于出口腔的流道阻力远小于纤维载体流道，可以将上述流道近似为矩形（见图5 （c））. 其中，流道k的长 ${l_{k\left( 1 \right)}}$和宽 ${l_{k\left( 2 \right)}}$分别如下：

(5) ${l_{k\left( 1 \right)}} = \left\{ \begin{array}{l} {L_{_{{{{A}}_{k{\rm{ + }}1}}{{{B}}_k}}}},\;k \leqslant b - 1 ; \\ {L_{_{{{{A}}_k}{{{B}}_k}}}},\;\;\;\;k \geqslant b .\\ \end{array} \right.$

(6) ${l_{k\left( 2 \right)}} = \left\{ \begin{array}{l} {L_{_{{{{A}}_k}{{{A}}_{k{\rm{ + }}1}}}}}\sin \;\left( { - {{\tan }^{ - 1}}\;{{{k}}_{{{{A}}_{k + 1}}{{{B}}_k}}}} \right),\; k \leqslant b - 1 ;\\ {L_{_{{{{A}}_k}{{{A}}_{k{\rm{ + }}1}}}}}\sin \;\left( {{{\tan }^{ - 1}}\;{{{k}}_{{{{A}}_k}{{{B}}_k}}}} \right),\;\;\;\;\;\;\;k \geqslant b .\\ \end{array} \right.$

式中： $L$为线段的长度， ${{k}}$为线段的斜率，其计算方式可以参考文献[27]；下标A、B如图5所示.

最终，纤维载体整体流道网络模型如图6所示.

纤维载体等效电阻网络建模的任务是在其流道划分的基础上，将流道单元等效为纯电阻电路，进而有效地模拟流体在流道中受到的流动阻力.

根据纤维载体流体微通道网格单元的特点（见图7 （a）），假设每个网格单元4条边界的中点是流体在相邻网格单元间流动的端口，即相邻网格单元通过4个端口进行等效电流传输. 端口从下至上按逆时针方向标记为A、B、C、D，如图7 （b）所示.

流体在任意2个端口间流动都会受到纤维及其表面的阻力作用，这一作用可以等效为端口之间的电阻. 这样，确定一个网格单元的电阻电路如图7 （c）所示. 其中，4个等效电阻阻值相同，被标记为1、2、3、4. 采用这种电阻电路，可以通过端口之间的组合来模拟流体在孔隙中不同的流动方向. 为了分析方便，进一步将纤维载体流体微通道的等效电阻整理为如图7 （d）所示的形式.

须说明的是，每个网格单元4个等效电阻的阻值等于AC或BD端的电阻. 因此，流体沿AC或BD端口流动受到的阻力等价为该单元电路中单个等效电阻的阻值. 另外，由于部分边缘网格单元与反应腔壁面相接，部分端口的等效电流为0.

与纤维载体流体微通道类似，进出口腔流道中划分的矩形单元被视作流体流动通道. 不同的是，由于在进出口腔的边界不存在流体流动，假设每个矩形单元3条边界的中点是流体在相邻流道间流动的端口，并根据端点在H的左侧还是右侧（见图8 （a）），分别按顺时针或逆时针方向将矩形单元中的端口标记为A、B、C （见图8 （b））. 其中，端口A是进出口腔矩形流道与纤维载体流体微通道连接的通道；BC端口间的电流等效为相邻矩形流道之间的流动，AC端口间的电流等效为矩形流道和纤维载体流体微通道之间的流动；由于反应物几乎不会沿远离进/出口的方向流动，AB端口间的等效电流可以忽略. 另外，由于直接与纤维载体流体微通道相连，进出口腔最左侧和最右侧的2条矩形流道都缺少1条与之相邻的矩形流道，只有AC端口之间存在等效电流.

根据上述分析，确定单个矩形流道的纯电阻电路如图8 （b）所示. 其中，平行于矩形单元长边的电阻标记为1，垂直于矩形单元长边的电阻标记为2. 最终，出口腔的等效电阻网络可以整理成如图8 （c）所示的形式.可以看出，与前文对应，进出口腔最左侧和最右侧的矩形流道的等效电阻电路中只有电阻1，没有电阻2.

基于上述工作，得到纤维载体的整体等效电阻网络如图9所示. 图中，R为电阻，上标p、i、o分别表示纤维载体、进口腔和出口腔，各单元的电阻标记分别由图7 （c）和图8 （c）确定，比如， $R_{k(2)}^{\rm{o}}$为出口腔第k个矩形单元中的第2个电阻， $R_{i, j(3)}^{\rm{p}}$为纤维载体（i， j）个网格单元中的第3个电阻. 矩形单元索引k和网格单元索引j之间的关系为

(7) $k = \left\{ \begin{array}{l} {{\left( {26n - 1{\rm{ + }}{{\left( { - 1} \right)}^n}} \right)} \mathord{\left/ {\vphantom {{\left( {26n - 1{\rm{ + }}{{\left( { - 1} \right)}^n}} \right)} 4}} \right. } 4} - j,\;\forall {\varOmega ^{\rm{i}}} ; \\ j - {{\left( {2n - 1{\rm{ + }}{{\left( { - 1} \right)}^n}} \right)} \mathord{\left/ {\vphantom {{\left( {2n - 1{\rm{ + }}{{\left( { - 1} \right)}^n}} \right)} 4}} \right. } 4},\;\;\;\forall {\varOmega ^{\rm{o}}} . \\ \end{array} \right.$

式中： ${\varOmega ^{\rm{i}}}$、 ${\varOmega ^{\rm{o}}}$分别表示进口腔和出口腔区域. j的取值范围如下：

(8) $j \in \left[ {{{\left( {2n - 1{\rm{ + }}{{\left( { - 1} \right)}^n}} \right)} \mathord{\left/ {\vphantom {{\left( {2n - 1{\rm{ + }}{{\left( { - 1} \right)}^n}} \right)} 4}} \right. } 4},\;{{\left( {26n - 1{\rm{ + }}{{\left( { - 1} \right)}^n}} \right)} \mathord{\left/ {\vphantom {{\left( {26n - 1{\rm{ + }}{{\left( { - 1} \right)}^n}} \right)} 4}} \right. } 4}} \right].$

对于如图9所示的等效电阻网络，设定向左和向上为正. 在此基础上，借鉴基尔霍夫定律，将反应物的体积流量、流动阻力和压力降分别等效为电路中的电流、电阻和电压降，即可建立反应物流动的质量和压力守恒方程，进而形成完整的纤维载体等效电阻网络模型.

对于进口腔，总输入流量为

(9) $q_{\rm{t}}^{\rm{i}}{\rm{ = }}q_{b({\rm{1}})}^{\rm{i}} - q_{b(2)}^{\rm{i}}{\rm{ + }}q_{b - 1({\rm{1}})}^{\rm{i}}{\rm{ + }}q_{b{\rm{ - 1}}(2)}^{\rm{i}}.$

式中： $q_{i(j)}^{\rm{i}}$ （i = b−1，j = 1， 2）为流经电阻 $R_{i（j）}^{\rm{i}}$（见图9）的流量.

在进口腔内部，各节点的质量守恒方程为

(10) $q_{k(2)}^{\rm{i}}{\rm{ = }}\left\{ \begin{array}{l} q_{k - {\rm{1}}(1)}^{\rm{i}}, \;k = 1 ;\\ q_{k - {\rm{1}}(1)}^{\rm{i}}{\rm{ + }}q_{k - 1(2)}^{\rm{i}},\; k \in \left[ {2,\; b - 1} \right] ;\\ q_{k{\rm{ + }}1({\rm{1}})}^{\rm{i}} - q_{k{\rm{ + 1}}({\rm{2}})}^{\rm{i}},\; k \in \left[ {b,\; m - 2} \right] ;\\ - q_{k + {\rm{1}}(1)}^{\rm{i}},\; k = m - 1 . \\ \end{array} \right.$

式中：m可以由式（4）计算得出，j由式（8）约束.

$q_{k\left( 1 \right)}^{\rm{i}}$表达式为

(11) $\left. { \begin{array}{l} q_{k(1)}^{\rm{i}}{\rm{ = }}q_{1, j（(4)}^{\rm{p}} - q_{1, j(1)}^{\rm{p}};\; k \ne b ,\\ q_{k - 1(1)}^{\rm{i}}{\rm{ + }}q_{k(1)}^{\rm{i}}{\rm{ = }}q_{1, j（(4)}^{\rm{p}} - q_{1, j(1)}^{\rm{p}};\; k = b .\\ \end{array}} \right\}$

出口腔的质量守恒方程与入口腔的类似，表达式如下：

(12) $q_{\rm{t}}^{\rm{o}}{\rm{ = }}q_{b(1)}^{\rm{o}} - q_{b(2)}^{\rm{o}}{\rm{ + }}q_{b - 1(1)}^{\rm{o}} + q_{b - 1(2)}^{\rm{o}},$

(13) $q_{k(2)}^{\rm{o}}{\rm{ = }}\left\{ \begin{array}{l} q_{k - 1(1)}^{\rm{o}},\; k = 1 ; \\ q_{k - 1(1)}^{\rm{o}}{\rm{ + }}q_{k - 1(2)}^{\rm{o}},\; k \in \left[ {2,\; b - 1} \right] ; \\ q_{k{\rm{ + }}1(1)}^{\rm{o}} - q_{k + 1(2)}^{\rm{o}},\; k \in \left[ {b,\; m - 2} \right] ; \\ - q_{k + 1(1)}^{\rm{o}}, \;k = m - 1 . \\ \end{array} \right.$

(14) $\left. { \begin{array}{l} q_{k(1)}^{\rm{o}}{\rm{ = }}q_{4n, j(2)}^{\rm{p}} - q_{4n, j(3)}^{\rm{p}}; \;k \ne b , \\ q_{k - 1(1)}^{\rm{o}}{\rm{ + }}q_{k(1)}^{\rm{o}}{\rm{ = }}q_{4n, j(2)}^{\rm{p}} - q_{4n, j(3)}^{\rm{p}};\; k = b . \\ \end{array}} \right\}$

对于纤维载体，其最左和最右侧、最上和最下侧电阻电路各节点的质量守恒方程为

(15) $\left. { \begin{array}{l} q_{i, 1(2)}^{\rm{p}}{\rm{ = }} - q_{i, 1(1)}^{\rm{p}} ,\\ q_{i, 7n（(4)}^{\rm{p}}{\rm{ = }} - q_{i, 7n(3)}^{\rm{p}} , \\ \end{array}} \right\}$

(16) $ \begin{array}{l} q_{1, j(1)}^{\rm{p}}{\rm{ = }}q_{1, j（(4)}^{\rm{p}} ,\; q_{4n, j(2)}^{\rm{p}}{\rm{ = }}q_{4n, j(3)}^{\rm{p}};\; j \ne k. \\ \end{array}$

式中：k由式（7）确定.

除此之外，纤维载体内部纵向、横向电阻电路各节点的质量守恒方程分别如下：

(17) $q_{i, j(2)}^{\rm{p}} + q_{i + 1, j(1)}^{\rm{p}}{\rm{ = }}q_{i, j(3)}^{\rm{p}} + q_{i + 1, j（(4)}^{\rm{p}}, $

(18) $q_{i, j(3)}^{\rm{p}} + q_{i, j + 1(1)}^{\rm{p}}{\rm{ = }}q_{i, j（(4)}^{\rm{p}} + q_{i, j + 1(2)}^{\rm{p}}.$

与质量守恒方程建立方式类似，纤维载体进、出口腔和孔隙区域的压力平衡方程如下：

(19) $\left. { \begin{array}{l} {{U}}_{k\left( 1 \right)}^{\rm{i}} + {{U}}_{1, j + 1(1)}^{\rm{p}}{{ = U}}_{k + 1}^{\rm{i}} + {{U}}_{1, j（(4)}^{\rm{p}};\; k \leqslant b - 2 ,\\ {{U}}_{k\left( 1 \right)}^{\rm{i}} = {{U}}_{k + 1\left( 1 \right)}^{\rm{i}};\; k = b - 1 , \\ {{U}}_k^{\rm{i}} + {{U}}_{1, j + 1（(4)}^{\rm{p}}{{ = U}}_{k + 1\left( 1 \right)}^{\rm{i}} + {{U}}_{1, j（(4)}^{\rm{p}};\; k \geqslant b ,\\ \end{array} } \right\} $

(20) $\left. { \begin{array}{l} {{U}}_{k\left( 1 \right)}^{\rm{o}} + {{U}}_{4n, j(3)}^{\rm{p}}{{ = U}}_{k + 1}^{\rm{o}} + {{U}}_{4n, j + 1(2)}^{\rm{p}};\; k \leqslant b - 2 , \\ {{U}}_{k\left( 1 \right)}^{\rm{o}} = {{U}}_{k + 1\left( 1 \right)}^{\rm{o}};\; k = b - 1 ,\\ {{U}}_k^{\rm{o}} + {{U}}_{4n, j(3)}^{\rm{p}}{{ = U}}_{k + 1\left( 1 \right)}^{\rm{o}} + {{U}}_{4n, j + 1(2)}^{\rm{p}}; \;k \geqslant b , \\ \end{array}} \right\} $

(21) $\left. { \begin{array}{l} {{U}}_{i, j(1)}^{\rm{p}}{{ + U}}_{i, j（(4)}^{\rm{p}}{{ = U}}_{i, j(2)}^{\rm{p}}{{ + U}}_{i, j(3)}^{\rm{p}} , \\ {{U}}_{i, j(3)}^{\rm{p}}{{ + U}}_{i, j + 1(2)}^{\rm{p}}{{ = U}}_{i + 1, j（(4)}^{\rm{p}}{{ + U}}_{i + 1, j + 1(1)}^{\rm{p}} . \\ \end{array} } \right\}$

若令X指代图9中某个矩形单元或网格单元的等效电路（如第k个矩形单元或第（i， j）个网格单元），则式（19）~（21）中的 ${{{U}}_{{X} \left( x \right)}} = {q_{{X} \left( x \right)}}{R_{{X} \left( x \right)}}$，表示单元X中第x个电阻上的电压， ${{{U}}_{{X}}} = \displaystyle\sum\nolimits_x {{q_{{{X}}\left( x \right)}}{R_{{{X}}\left( x \right)}}}$，表示单元X中所有电阻上的电压之和.

纤维载体中的流动阻力与压力降有关. 由于流体在其中的流动为层流^[34]，可以利用达西定律来描述纤维载体的整体压力降^{[33， 35]}：

(22) $\frac{{\Delta p}}{L}{\rm{ = }} - \frac{\mu }{K}{{V}}. $

式中：Δp为沿长度L方向的压降；μ为流体的动力黏度；K为渗透因子；V为流体流过该段长度的表观速度，负号表示压降方向和速度方向相反.

对比基尔霍夫定律中电压、电流与电阻之间的关系，式（22）可以改写为

(23) $\Delta p{\rm{ = }}\frac{{\mu L}}{{K{{A}}_{\rm{t}}^{\rm{p}}}}q_{\rm{t}}^{\rm{p}} = R_{\rm{t}}^{\rm{p}}q_{\rm{t}}^{\rm{p}}. $

式中： $q_{\rm{t}}^{\rm{p}}$为流体流过纤维载体孔隙区域的总流量， $A_{\rm{t}}^{\rm{p}}$为垂直于流动方向的横截面积； $R_{\rm{t}}^{\rm{p}}$为流体流过纤维载体时受到的总流动阻力.

类似地，可以根据纤维载体及流体参数确定流体微通道等效电阻网络模型中第（i， j）个网格单元中第x个电阻（见图9）的阻值，即流体的流动阻力：

(24) $R_{i, j(x)}^{\rm{p}} = \frac{{\mu l}}{{KA_{i, j}^{\rm{p}}}};\;x = 1,2,3,4.$

式中：l为第（i， j）个网格单元的长度， $A_{i, j}^{\rm{p}}$为流体流动方向上的网格单元的横截面积.

对于进出口腔的矩形单元（见图5），由腔体内壁面的剪切应力引起的摩擦损失导致的压力降占主导作用，可以采用哈根-泊肃叶方程进行计算^[27]：

(25) $\Delta p{\rm{ = }} - \frac{{32\mu L\xi }}{{{D^2}}}{{V}} = \frac{{32\mu L\xi }}{{{D^2}A_{\rm{t}}^{{\rm{io}}}}}q_{\rm{t}}^{{\rm{io}}} = R_{\rm{t}}^{{\rm{io}}}q_{\rm{t}}^{{\rm{io}}}.$

式中：ξ、D分别为非圆系数和横截面当量直径，可以通过矩形单元在垂直流动方向的截面深度和宽度表示，具体方法参见文献[36]； $A_{\rm{t}}^{{\rm{io}}} $为进口腔或出口腔面积； $q_{\rm{t}}^{{\rm{io}}} $为流体流过进口腔或出口腔的总流量； $R_{\rm{t}}^{{\rm{io}}} $为流体流过进口腔或出口腔时所受到的总流动阻力.

由于进出口腔的形状不是规则的四边形，式（25）中沿流体流动方向的长度L须按如图8所示的流向计算. 当流体流向下一个矩形单元时，对于第k个矩形单元，分别取ξ、D_k、A_k作为其在该流向的非圆系数、当量直径和横截面积，则流动阻力为

(26) $R_{k(2)}^{{\rm{io}}} = \frac{{32\mu {l_{k\left( 2 \right)}}\xi }}{{D_k^2A_k^{{\rm{io}}}}}.$

式中： ${l_{k\left( 2 \right)}}$由式（6）计算.

当流体从进口腔的矩形单元流入纤维载体孔隙区域，或反过来流向出口腔的矩形单元时，流体流入方向与流出方向垂直，可以沿流动方向对单位长度上的流动阻力进行积分得到其流动阻力：

(27) $R_{k(1)}^{{\rm{io}}} = \int\limits_0^{{l_{k\left( 1 \right)}}} {\frac{{32\mu {\xi _k}}}{{D_k^2A_k^{{\rm{io}}}{l_{k\left( 1 \right)}}}}} l{\rm{d}}l = \frac{{16\mu {\xi _k}{l_{k\left( 1 \right)}}}}{{D_k^2A_k^{{\rm{io}}}}}.$

式中： ${l_{k\left( 1 \right)}}$由式（5）计算.

由式（22）~（27）可知，如果确定纤维载体的尺寸参数和渗透因子K，即可求得等效电阻网络中的流动阻力. 在此基础上，综合式（9）~（21）即可求解出每个等效电阻中的体积流量，以此除以单元横截面积则可以得到对应的流体速度.

为了验证纤维载体等效电阻网络模型的有效性和求解效率，与前期CFD方法^[34]进行对比. 以孔隙率E=70%、80%、90%的纤维载体为例进行介绍. 其网格划分阶数n按式（1）~（3）计算分别为13、9、5，渗透因子K=1.19×10⁻¹⁰、2.17×10⁻¹⁰、5.92×10⁻¹⁰ m⁻².

由于不同孔隙率纤维载体的分析结果具有相似性，如图10所示仅展示了孔隙率为80%的纤维载体在入口速度为0.096 m/s时2种方法获得的流速场分布. 图中，v为流速. 可以看出，纤维载体的流速场在宽度及长度方向上均近似呈对称分布，且流体流速在与进出口腔边界相交的交点附近有较大起伏，在流速场上出现了4个峰值与4个谷值；纤维载体中心区域的流速变化则较为平缓；对于这些分布特征，2种方法获得的流速场无论是在形式还是数值方面都基本吻合.

如表1所示为皮尔森相关系数r分析结果. 可以看出，2种方法获得的3种孔隙率纤维载体的流速场具有平均约为98%的强相关性，且纤维载体流速场的平均误差η=1.434×10⁻⁴~2.499×10⁻⁴ m/s，平均绝对百分比误差MAPE=0.741%~1.280%，均方误差MSE=4.188×10⁻⁴~4.884×10⁻⁴ m/s. 因此，使用等效电阻网络模型模拟纤维载体流速场的分布是可行的.

如表2所示，进一步对比2种方法获得的纤维载体流速场信息. 表中，v_avg、v_max、v_min分别为平均流速、流速峰值、流速谷值，t为计算时间. 可以看出，2种流速场的平均流速基本一致，因此它们在整体流速上基本没有偏移；2种流速场对应的峰值与谷值 $v_{\text{min}}$之间都有一定的差异，主要原因是不同流速之间的流体存在着摩擦力^[27]. 在将纤维载体孔隙简化为网格单元时，并未考虑单元之间因不同流速流体之间摩擦力造成的压降损失. 这使得等效电阻网络模型的流速分布数据在靠近纤维载体与进出口腔边界交点处的起伏稍大. 尽管如此，对比2种流速场的均方差可知其整体数值及变化趋势仍然具有较高的拟合度.

另外，须注意的是，2种方法所用的时间t相差巨大. 由于CFD方法须迭代求解非线性Navier-Stokes偏微分方程，而等效电阻网络模型仅须求解式（9）~（27）组成的线性方程组，因此后者的求解效率平均提升约2.9×10⁴倍. 从这个角度看，即使CFD方法理论上更为准确，但由于等效电阻网络模型求解效率极高，仍可以用其替代CFD方法进行流速场的前期分析与筛选，后期则可以基于筛选的流速场采用CFD方法对纤维载体进行更为细致的分析.

（1）通过拓展面向多孔纤维结构固相材料的统计网络模型，可以将纤维载体的复杂流体微通道结构简化为规则的流道网络.

（2）基于等效电阻网络模型与采用商用CFD方法获得的纤维载体流速场，无论是在形式还是在数值方面都较吻合，具有约98%的相关性.

（3）由于只须求解线性方程组，而不是非线性Navier-Stokes偏微分方程，基于等效电阻网络模型的纤维载体流速场分析效率约为CFD方法的2.9×10⁴倍.

本研究的研究结果为多孔纤维型重整制氢微反应器流速场的分析和结构优化设计提供了支撑技术. 后续将考虑不同流速流体之间摩擦力造成的压降损失，进一步提升方法的准确性.