基于模拟退火算法的真实多孔电极中热-质传递的研究

作者：黄盼练成刘洪来来源：《化工学报》日期：2022-08-30人气：1483

电场作用下电极中电化学反应-热质传递现象是典型的多尺度问题[1-3]，如图1所示。微观尺度(电子-离子迁移、晶格稳定性)、介观尺度(界面热/动力学、热-质传递、电流电压分布)和宏观尺度(散热性能、充放电管理)的传递和反应特性，直接决定系统中的浓度分布和反应速率，最终影响储能转化效率。而系统中发生传递和反应过程的场所主要在多孔电极中，且多孔电极的不规则表面引起的限域效应显著影响着其内的传质和传热过程，所以深入研究介观尺度下多孔电极中的热-质传递现象及其耦合机制，对高性能电化学储能器件的设计有重要意义[3]。但是，与平板电极不同，对多孔电极中热-质传递现象的研究存在两个关键问题：一方面，多孔电极丰富的孔道体系、多样化的表面形态和分布复杂的催化活性位点，难以用传统数学工具去准确描述其结构[4-6]；另一方面，多孔电极中发生的过程涉及在受限空间中的离子迁移、热量产生和传递、双电层形成、电子传递的耦合，极大地影响了离子和热量在孔道中的传递过程，并难以直接用实验进行原位表征[7-8]。因此，多孔电极中的热-质传递的研究一直以来都难以获得突破性进展。

图1

图1 储能和转换过程中的电极尺度问题及与其他尺度问题的联系

Fig.1 The electrode scale issues in energy storage and conversion processes and their relation to other scale issues

目前，研究者们可以采用密度泛函理论 (DFT)[9-10]和分子动力学(MD)模拟[11-12]等方法研究纳米级多孔结构中的演化过程和充电机理。例如纳米多孔电极与离子液体的结合是提高超级电容器能量密度的重要手段。然而，这总是伴随着功率密度的降低，特别是考虑到高黏度和大空间位阻的离子液体。Gan等[11]利用MD模拟发现在阈电位作用下，具有疏离子孔的电极内的离子呈现出一种新的充电机制，即离子吸附。进一步得到了充电时间/电容与电压/耐离子性能之间的定量关系，以评价协同提高超级电容器能量密度和功率密度的临界条件。虽然DFT和MD模拟可以考虑微观粒子间的作用力，但是受限于计算能力的瓶颈，只能模拟离子浓度和温度在含周期性边界条件的高度有序的纳米级多孔结构(如分子筛和金属有机框架)中纳秒级的演化过程，并不能全面描述真实多孔电极中，尤其是常见的无定形碳电极中的离子和热量的传递过程[9]。连续介质模型通过偏微分方程组描述离子浓度和电势的状态演化及其本构关系，可以建立纳米级及以上空间尺度的复杂的多物理场耦合关系[13]。例如，Tao等[14]采用基于Poisson-Nernst-Planck(PNP)方程和Navier-Stokes方程的组合方法，通过控制氧化石墨烯膜的偏置角度，研究了不同层间距分布的氧化石墨烯膜中离子的输运现象。d′Entremont等[15]从第一性原理推导出控制能量方程，并与修正PNP模型相结合，从而得到由于离子扩散、空间效应和混合熵的变化而产生的不可逆焦耳热和可逆热产生速率。但是，目前连续介质模型在多孔电极中的大规模运用仍存在两个难点：准确描述多孔介质的复杂结构和快速求解含有复杂边界条件的偏微分方程组。因此，多种简化的多孔电极表示方法和经验模型被广泛运用。天然形成的和部分人工制造的多孔介质的微结构具有分形特征，可用分形维数表示。Sakaguchi等[16]利用PNP方程研究了分形多孔电极中的充电过程，并发现分形多孔介质充注过程的时间演化遵循幂律，指数与分形维数有关。目前大量研究都集中在单个孔隙中的充放电过程，因此可以通过等效介质近似(effective medium approximation,EMA)方法[17]，用假设的、具有相同导电性的单个孔隙组成的均质网络来替代无序多孔介质的非均质网络中的孔隙，同时保持孔隙网络的连通性。Lian等[18]利用该方法研究了平均孔径、孔径分布和孔连通性对多孔炭电极离子输运特性(包括电导和电导率)的影响。更进一步地，孔隙网络模型(pore network modeling, PNM)采用不同大小的球体和圆柱代表多孔结构中的空腔和吼道[19]。因此，PNM可以考虑多孔电极中的孔径分布和孔道连通性等性质，且计算成本非常低，可以在合理的时间内通过PNP方程计算储能设备的电化学性能和寻找最佳电极结构[20-21]。此外，Lian等[8]报道了一种简单且基于物理的堆叠电极模型来表示多孔电极，成功用PNP方程和等效电路模型解释了超级电容器的缓慢充电动力学，并发现充电过程可以分为两个不同弛豫时间尺度的阶段。虽然上述简化的多孔电极模型被成功运用在解释储能过程中的各种现象，但这些方法只能在一定程度上接近真实的多孔电极的孔径分布，很难表示多样化的表面形态和分布复杂的催化活性位点，因此限制了这些简化模型在多孔电极中的应用。

为了对真实的多孔电极(尤其是非均质、各向异性材料)中的热-质传递过程进行深入的理解，需要对孔隙和骨架结构进行高分辨率三维表示。目前可以通过微尺度X射线计算机断层扫描[22-23](computed tomography, CT)或扫描电子显微镜[24](scanning electron microscope, SEM)以非破坏性的方式获取多孔介质的孔隙结构。但是，CT不能分辨微米以下的结构，SEM的缺点是它只提供二维信息。获得多孔电极的三维结构表示的一种流行的替代方法是随机重构。总地来说，随机重构包括处理可用的信息，并生成符合所导出的介质特性的随机结构两个步骤。基本上有两种重构方法。第一种依赖于基于流程的建模。这种方法试图模拟自然材料在其原始环境中的形成过程[25]，或创造工程材料的过程[26]。然而，现有的基于过程的方法只考虑了某些方面，难以准确描述材料结构。一种试图解决这些问题的方法是相场方法[27]。第二种类型的重构方法采用了一种更加数据驱动的方法，并使用了统计数据。这些统计数据主要来自通过图像分割识别出孔隙相和骨架相的各类图片，使用两点统计信息，如自相关或概率函数[28]、线性路径函数[29]和聚类函数[26]。最常用的两点统计量重构方法是模拟退火(simulated annealing, SA)算法[30-32]。与可能非常复杂，并且可能需要大量超参数调优的基于多点统计[33]的方法相反，模拟退火算法的简单性和准确性在实际运用中很有吸引力。Wu等[34]采用SA算法重构商业锂离子电池阴极LiCoO2的三维微观结构，包括活性材料相、孔隙相和添加剂相。重构阴极的表征提供了重要的结构和输运性质，包括两点相关函数、体积比表面积、曲折度和单个相的几何连通性。Habte等[35-36]研究了正极材料的微观结构形态对锂离子电池性能的影响，首先通过模拟退火算法生成了球形电极的结构，然后计算出结构参数并代入阻抗谱经验公式中，最后得到的阻抗值与实验结果符合较好。目前研究大多都停留在多孔电极的重构、孔结构参数[34,37]和等效传递系数[38-39]的计算上，真实多孔结构的形貌和孔隙结构对离子传递和电极导热的影响尚未被研究过。

本文提出一套研究介观尺度下真实多孔电极中的热-质传递的研究框架：首先采用一种基于改进的状态更新的随机重建方法和动态退火系数相结合的模拟退火算法，将图像分割后的二维SEM图重构为真实介观尺度的三维多孔电极，重构生成的多孔电极的结构和真实多孔电极截面上的结构在统计学意义上是一致的，并在此基础上建立真实多孔电极中的离子传递和电极导热模型。

1 研究思路

本文的研究思路如图2所示。通过改进的模拟退火算法和有限元计算，研究了真实多孔电板中的离子传递和电极导热现象。

图2

图2 基于模拟退火算法重构的真实多孔电极中的离子传递和电极导热

Fig.2 The ion transport and electrode heat conduction in real porous electrode based on simulated annealing algorithm

2 真实多孔电极的三维重构及其结构参数计算

如图3所示，多孔电极三维重构过程如下：首先利用SEM扫描多孔电极，得到参考图像；然后对其进行预处理，即通过图像分割识别孔隙相和骨架相，得到二值化图像序列，并将二值化图像序列作为参考模型；然后生成与参考模型的孔隙率相同的随机初始点云，利用改进的模拟退火算法进行随机重构，将两点相关函数作为优化目标，采用一种改进的体素交换方法更新状态，得到重构的多孔电极的点云；最后将点云实体化为立方体构成的空隙相和骨架相的结构，并计算了其结构参数。

图3

图3 多孔电极三维重构流程图

Fig.3 The flow chart of three-dimensional reconstruction of porous media

2.1 图像预处理

目前研究和应用较多的电极材料主要集中在多孔炭材料、过渡金属氧化物和导电高分子等方面。其中多孔炭材料具有比表面积大、导电性好、结构稳定、资源丰富、价格低廉等优势，既可直接作为电极，又可与其他材料复合，起到传输电子、舒缓体积膨胀以及优化界面反应等作用，被广泛用为电极材料[40]。多孔炭材料的SEM图来自于Kunanusont等[24]的工作，如图3(a)所示。他们采用超临界二氧化碳干燥法，在10.0~20.0 MPa和40℃条件下制备了乙炔炭黑和聚偏氟乙烯基多孔电极，并研究了压力对电极形态、孔隙率和电化学性能的影响。

通过SEM扫描得到的灰度图像序列中，孔隙相和骨架相的灰度值差异明显。因此，首先通过阈值分割算法将灰度图像序列转化为只有孔隙相(黑色区域，灰度值为0)和骨架相(白色区域，灰度值为255)的二值化图像序列。采用SEM图中所有像素的灰度值的平均值(本文为144)作为阈值。当某处的灰度值高于阈值时，则将此处的灰度值设为255；当灰度值小于等于平均值时，则设灰度值为0。处理后的图像中含有斑点噪声和大面积的孔洞，采用形态学运算填充大面积孔洞，同时去除图像中的杂质和亮斑等噪声干扰[41]，从而得到预处理之后的二值化图像序列，如图3(b)所示。

2.2 图像结构描述符计算

多孔电极的结构与其结构描述符密切相关，在对多孔电极重构时，需要考虑其结构描述符等信息，使重构模型与参考模型更加符合[42]。采用特征函数 $I_{r}^{j}$ 定义不同的结构描述符。在多相系统中，假设第 $j$ 相所占的区域为 $V_{j}$ ，其在系统中的浓度为 $φ_{j}$ ，则第 $j$ 相的特征函数为

$I_{r}^{j} = \{\begin{matrix} 1 \begin{matrix} r \in V_{j} \end{matrix} \\ 0 \begin{matrix} r \notin V_{j} \end{matrix} \end{matrix}$ (1)

式中， $r$ 是系统中任意一点。因为本文中的多孔电极仅由孔隙相和骨架相两相构成，所以令 $j = 1$ 表示骨架相， $j = 2$ 为孔隙相，则多孔电极的孔隙率为 $φ = 1 - φ_{1} = φ_{2} = \bar{I^{2} (r)}$ 。

本文采用的第一个结构描述符是两点概率函数 $S_{2}^{j} (r_{1}, r_{2})$ ，简写为 $S_{2}$ ，它给出了位置 $r_{1}$ 和 $r_{2}$ 的两个像素属于同一相位的概率[30]。 $S_{2}^{j} (r_{1}, r_{2})$ 定义为

$S_{2}^{j} (r_{1}, r_{2}) = 〈I_{r_{1}}^{j} I_{r_{2}}^{j}〉$ (2)

在平衡态系统中，粒子之间由于相互作用而存在一定的相关性。两点概率函数作为一种常用的结构描述符，可以很好地描述系统的微观结构特征。第二个描述符是两点线性路径函数 $L_{2}^{j} (r_{1}, r_{2})$ ，简写为 $L_{2}$ 。这个函数计算从 $r_{1}$ 到 $r_{2}$ 的线段完全位于第 $j$ 相的概率[29]。值得注意的是， $S_{2}$ 和 $L_{2}$ 仅在执行重构时使用。而为了评估重构质量，本文使用两点聚类函数 $C_{2}^{j} (r_{1}, r_{2})$ ，简写为 $C_{2}$ ，它计算在位置 $r_{1}$ 和 $r_{2}$ 上的两个像素属于同一簇的概率。簇被定义为属于同一相并相互连接的像素的集合。

2.3 模拟退火算法

模拟退火算法是一种通用的概率优化算法，起源于金属的退火过程，用于在一个很大的搜寻空间中寻找出最优解，可以求解复杂的非线性优化问题。模拟退火算法其实是一种贪心算法，每次都选择一个当前最优解，因此只能搜索到局部的最优值。但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出局部的最优解，达到全局的最优解。

2.3.1 模拟退火算法基本框架

图4是由Yeong等[31]首次提出的模拟退火算法运用在多孔介质重构中的基本框架。首先获得参考模型 $w_{0}$ [图3(b)]的能量函数 $E (w_{0})$ ,设置初始温度为 $T = T_{0}$ 。然后随机产生与 $w_{0}$ 的孔隙率相同的初始结构 $w$ ，并计算其能量函数 $E (w)$ 。通过像素交换的方式产生新的结构 $w'$ ，计算其能量函数 $E (w')$ 以及能量增量 $Δ E = E (w') - E (w)$ 。若 $Δ E < 0$ ，则接受新结构，即 $w = w'$ , $E (w) = E (w')$ ；否则以概率 $P = e^{- Δ E / T}$ 接受新结构。在当前温度下新结构被多次拒绝后，通过退火过程降低温度 $T$ ，然后重新搜索最优结构，直到满足算法的终止条件。

图4

图4 模拟退火算法流程图

Fig.4 The flow chart of simulated annealing algorithm

退出模拟退火算法的条件如下：①系统的能量函数小于阈值 $E_{t h}$ ；②相邻两步能量下降的绝对值小于 $Δ E_{t h}$ ；③连续 $N_{c o n}$ 次结构无更新，即连续拒绝随机产生的新结构 $N_{c o n}$ 次；④达到最大的迭代次数 $N_{i t e r}$ 。

2.3.2 能量函数

原则上，任何结构描述符或描述符的组合都可以用来计算能量函数 $E$ 。一般来说，定义 $E$ 时考虑的描述符越多，计算量就越大，且过多描述符的加入并不能显著提高重建质量。因此，本文仅使用 $S_{2}$ 和 $L_{2}$ 的组合 $S_{2} - L_{2}$ 来计算 $E$ [31-32]。设能量函数 $E$ 为重构模型与参考模型在各个方向上 $S_{2} - L_{2}$ 的差值平方和的范数，计算方法为

$E = \sqrt[]{\sum_{d} {\{\sum_{r} {[(S_{2, d}^{j} (r) - L_{2, d}^{j} (r)) - (S_{2 r e f, d}^{j} (r) - L_{2 r e f, d}^{j} (r))]}^{2}\}}^{2}}$ (3)

式中，d指模型维度。对于二维图像，d包含4个维度，即沿 $x$ 轴和 $y$ 轴，以及两个主对角线方向 $d x y$ 和 $d y x$ ；对于三维图像，d包含9个维度，除了二维的4个方向，还有沿 $z$ 轴，以及 $d x z$ 、 $d z x$ 、 $d y z$ 和 $d z y$ 方向。

2.3.3 基于DPN值的像素交换

模拟退火算法的最重要的过程是在保持多孔电极孔隙率 $ϕ$ 不变的前提下，通过一定规则产生新结构。将孔隙相和骨架相的像素进行位置交换，等效于产生新结构的过程。

利用传统的随机像素交换方法产生新结构时，图像中所有像素的交换概率相同，这种方法容易破坏已经重构好的局部结构，从而增加了寻找最优结构的时间。同时，图像中孤立的像素点可以被视为噪点，使用传统方法也难以去除。本文引入DPN(different phase neighbors)[43]的概念来改进传统的随机像素交换方法。像素的DPN值指在该像素的邻域内与该像素不同相的像素个数。像素的DPN值反映了该像素在图像中的孤立程度， $D P N = 0$ 表明该像素邻域内所有的像素都与其同相。对于二维图像，选取8邻域；对于三维图像序列，选取26邻域。图5列举了二维图像中几种不同的 $3 \times 3$ 大小的结构，并指出了中心像素的DPN值。

图5

图5 中心像素在不同结构中的DPN值

Fig.5 DPN values of center pixel in different structures

2.3.4 动态退火系数

退火过程就是体系温度下降的过程。温度越低，体系跳出当前解的概率就越低。温度降低太快，容易陷入局部最优；温度降低太慢，算法运行效率太低。因此，为了平衡精度和速度，需要选取合适的退火方法。常用的退火方法为 $T_{k} = λ T_{k - 1}$ ，其中， $k$ 为迭代的次数， $λ$ 为退火系数。传统模拟退火算法的退火系数是常数，无法考虑系统在不同重构阶段中能量函数的变化，容易陷入非最优结构或者增加寻优时间。本文采用基于动态退火系数的退火过程，考虑模拟退火过程中系统能量的变化情况，退火系数的计算方法[44]为

$λ = \{\begin{array}{l} λ_{m i n} E_{m i n} / E_{a v} < λ_{m i n} \\ λ_{m a x} E_{m i n} / E_{a v} > λ_{m a x} \\ E_{m i n} / E_{a v} o t h e r w i s e \end{array}$ (4)

式中， $λ_{m i n}$ 和 $λ_{m a x}$ 分别是退火系数的最小值和最大值，用于对退火系数限幅； $E_{m i n}$ 和 $E_{a v}$ 是迭代过程中能量函数的最小值和平均值。

2.3.5 参数设置、点云实体化和运行环境

根据前文模拟退火算法的退出条件，设置具体参数如表1所示。

表1 模拟退火算法退出迭代的参数设置

Table 1 Parameter setting of simulated annealing algorithm to exit the iteration

Eth	ΔEth	Ncon	Niter
1×10-6	1×10-8	500	1×106

新窗口打开| 下载CSV

通过模拟退火算法产生的点云无法通过有限元方法模拟其内的离子传递和电极导热过程，需要被实体化为由面构成的三维结构。本文采用Blender自带的API(application program interface)，将点云实体化为由小立方组成的多孔电极，最后提取表面结构。图3(c)、(d)分别展示了 $200 \times 200 \times 200$ 的多孔电极的骨架相和孔隙相，其中一个小立方体的边长为 $8 . 2^{} n m$ 。

本文的计算平台为4 x Intel(R) Xeon(R) Platinum 9242 CPU at 2.30 GHz和384G内存的服务器，模拟退火算法和实体化算法通过Python 3.8编程实现。

2.4 多孔电极的结构参数

孔隙率 $φ$ 是保证模拟退火算法随机构造的多孔电极结构符合实际的前提之一。 $φ$ 的计算方法是孔隙相的立方体数量 $N_{p}$ 与多孔结构所有立方体数量 $N_{t}$ 的比值，即 $φ = N_{p} / N_{t}$ 。因为本文构建的多孔结构为实际结构的一小部分，所以构建的多孔结构的六个方向的面为截面而非表面。因此，本文比表面积 $S_{V}$ 的计算方法为除去截面面积的骨架相表面积与骨架相的体积之比。

多孔电极是分形物质，它的孔隙结构复杂度与其孔隙分形维数等特征参数密切相关。目前常用的计算多孔结构分形维数的方法有盒计数法[45]、豪斯道夫方法[46]和随机游走法[47]。盒计数法由于其精度更高、运算速度更快而被广泛用于计算分形维数。对于多孔电极的孔隙体积 $V_{p}$ ，有幂律关系 $V_{p} (ε) \propto ε^{d}$ ，其中， $V_{p} (ε)$ 是在对应长度尺度 $ε$ 下测得的孔隙体积， $d$ 为多孔介质孔隙相的分形维数。对于三维多孔电极, 用尺度为 $ε$ 的立方体测量孔隙体积，则孔隙分形维数计算方法为 $d = l n N_{p} (ε) / l n ε$ ，其中， $N_{p} (ε)$ 为在长度尺度 $ε$ 下识别出孔隙的立方体数量。

在实际的多孔电极中，孔道的形态一般是弯曲的，离子的传输路径总长度为 $l$ ，一般远大于电极的厚度 $H$ 。曲折因子 $τ$ 表示结构对扩散过程的影响程度，定义为 $τ = l / H$ 。本文采用皮尔逊随机行走模拟[48]来计算曲折因子。

所有的结构参数均采用Python 3.8编程计算。表2展示了尺寸分别为 $50 \times 50 \times 50$ 、 $100 \times 100 \times 100$ 和 $200 \times 200 \times 200$ 个立方体的多孔电极的结构参数。结果表明多孔电极的尺寸对结构参数影响不大，验证了重构的多孔电极结构是稳定且有效的。

表2 通过模拟退火算法重构的多孔电极的结构参数

Table 2 Structural parameters of porous electrode reconstructed by simulated annealing algorithm

多孔电极尺寸/ (个立方体)	边长/ $μ m$	结构参数
多孔电极尺寸/ (个立方体)	边长/ $μ m$	孔隙率	比表面积/ $m^{- 1}$	分形维数	曲折因子
$50 \times 50 \times 50$	0.41	0.545	1.493×107	3.28	2.36
$100 \times 100 \times 100$	0.82	0.545	1.223×107	3.19	2.01
$200 \times 200 \times 200$	1.64	0.545	1.306×107	3.11	2.13

新窗口打开| 下载CSV

3 三维多孔电极中的传质与传热原理

为了研究多孔电极中的传质和传热现象，参考实际超级电容器的工作原理[8]，构建了如图6所示的两个模型。受限于计算机的计算资源，离子传递模型采用两个 $100 \times 100 \times 100$ 大小的立方体多孔电极(边长约为 $0 . 82^{} μ m$ )的孔隙相与 $200 \times 100 \times 100$ 大小的长方体体相相接。值得注意的是，两个孔隙相互为镜像关系，即右边的多孔电极是左边多孔电极经过x、y、z三个方向的镜像处理后生成的。而电极导热模型采用 $100 \times 100 \times 100$ 个格子大小的立方体多孔电极的骨架相。

图6

图6 真实多孔电极中的离子传递和电极导热模型

Fig.6 The model of ion transfer and electrode heat conduction in real porous electrodes

3.1 充电过程

刚开始时整个体系的离子密度分布是均匀的，电势为零。当左右两个多孔电极施加相反的电势后，体相和孔隙相中的阴阳离子受到电场力的作用，分别朝着相反的方向定向移动到多孔电极表面，从而形成电流并最终形成双电层(electrical double layer, EDL)储存能量。双电层的厚度通过德拜长度来表示。

$λ_{D} = \sqrt[]{(ε_{0} ε_{r} k_{B} T) / (2 e^{2} ρ_{b u l k})}$ (5)

式中， $ε_{0}$ 是真空介电常数，F·m-1; $ε_{r}$ 是相对介电常数； $k_{B}$ 是Boltzmann常数，J·K-1；T是温度，K;e是基本电荷，C; $ρ_{b u l k}$ 是体相离子数密度,m-3。由于在这个过程中离子电迁移导致电解液中的离子浓度不均匀，体系存在浓度梯度，所以充电过程中离子的运动还受到扩散的影响。

本文考虑到生成的双电层不会发生重叠以及求解的复杂度两个因素，于是采用Poisson-Nernst-Planck(PNP)方程来描述上述复杂孔隙结构中的离子充电动力学。电势 $ϕ (V)$ 和离子的数密度 $ρ ({m^{-}}^{3})$ 是PNP方程要求解的变量，其控制方程形式为：

$\frac{\partial ρ_{\pm}}{\partial t} + \nabla \cdot J_{\pm} = 0$ (6) $J_{\pm} = - D_{\pm} (\nabla ρ_{\pm} + \frac{z_{\pm} e ρ_{\pm}}{k_{B} T} \nabla ϕ)$ (7) $ε_{0} ε_{r} \nabla \cdot (- \nabla ϕ) = e q$ (8)

式中， $±$ 分别代表阴阳离子；t是时间,s; $\nabla = \frac{\partial}{\partial x} i + \frac{\partial}{\partial y} j + \frac{\partial}{\partial z} k$ ，是哈密顿算子;J是离子通量，m-3·s-1;D是离子扩散系数，m2·s-1; $z$ 是离子的价态; $q = z_{+} ρ_{+} + z_{-} ρ_{-}$ ，是电荷密度。式(7)左边的项为体系的离子通量，它由右边第一项的扩散通量和第二项的电场中电迁移通量共同组成。

设由两个多孔电极的孔隙相和体相组成的离子运动空间为 $Ω$ 。刚开始时整个空间的离子数密度 $ρ$ 是均匀分布的，电势 $ϕ$ 为零，即

I.C. $ρ_{\pm} (Ω, t < 0) = N_{A} c_{0}$

$ϕ (Ω, t < 0) = 0$ $\partial_{x} ϕ (Ω, t < 0) = 0$

式中， $N_{A}$ 是阿伏伽德罗常数。

因为多孔炭电极是良好的导体，所以可以假设整个多孔炭电极是一个等势体，即多孔电极和电解液的界面处的电势相等。图6中红色和蓝色区域代表电解液和多孔电极的界面以及电解液和集电极的界面，分别带正电势 $V_{m}$ 和负电势 $- V_{m}$ ，其中 $V_{m} = k_{B} T / e$ 为热电势。由于重构出的多孔电极并不是周期性结构，不能采用周期性边界条件，所以本文采用不可渗透且电势为零的边界条件，即灰色区域。因此。设带正电的边界为 $S_{+}$ ，带负电的边界为 $S_{-}$ ，其他边界为 $S_{0}$ ，则

B.C. $ϕ (S_{+}, t) = + V_{m}$

$ϕ (S_{-}, t) = - V_{m}$ $\partial_{x} ϕ (S_{+}, t) = 0$ $\partial_{x} ϕ (S_{-}, t) = 0$ $J_{i} (S_{0}, t) = 0$

本文的研究重点是在 $T = 298 . 15^{} K$ 时NaCl水溶液中的电解质。钠离子( $+$ )和氯离子( $-$ )的价态分别为 $z_{+} = 1$ 和 $z_{-} = 1$ 。本文使用以下参数设置[49-50]： $D_{\pm} = 1.6 \times 10^{- 9}^{} m^{2} \cdot s^{- 1}$ 、 $ε_{r} = 71$ 、 $c_{0} = 10^{} m o l \cdot m^{- 3}$ 。因此， $V_{m} = 0.02526 V$ ， $λ_{D} = 2 . 87^{} n m$ ， $τ_{R C} = 2.94 \times 10^{- 6}^{} s$ 。模拟时间为( $10^{- 3} ~ 10^{- 1}) τ_{R C}$ ，其中 $τ_{R C} = (H + L) λ_{D} / D$ 为表征平板电极充电完成的弛豫时间[10]， $H = 0 . 82^{} μ m$ 和 $L = 0 . 82^{} μ m$ 分别代表电极的厚度和体相的半长。

本文基于COMSOL Multiphysics 5.4软件，采用有限元法求解上述含有复杂边界条件的PNP方程，从而得到离子在多孔电极中的迁移行为。计算平台为4 x Intel(R) Xeon(R) Platinum 9242 CPU at 2.30 GHz和384G内存的服务器。

3.2 导热过程

在实际的充放电过程中，由于电流的生成和离子的重排等因素[15]，整个体系会产生热量。体系温度主要受到热生成、电解液内部的热传导和热对流，以及多孔电极内的热传导影响。而温度又会影响离子的输运参数和电极的导电性能，从而影响整个充电过程。充电、产热和传热三种现象相互耦合，极大地增加了计算量。此外，本文研究体系的电极导热的弛豫时间远小于充电的弛豫时间(或产热的弛豫时间)，且模拟时间较短，可以将充电现象和电极导热现象解耦，从而分开研究。因此，为了简化电极导热模型，只考虑电极中的热传导。本文采用傅里叶定律描述热传导过程。

$\frac{\partial T}{\partial t} = \frac{k}{ρ_{f} c_{p}} \nabla^{2} T = a \nabla^{2} T$ (9)

式中，k是热导率，W·m-1·K-1; $ρ_{f}$ 是骨架质量密度，kg·m-3; $c_{p}$ 是骨架相的比定压热容，J·kg-1·K-1;a是热扩散系数，m2·s-1。

设骨架相空间为 $γ$ ，骨架相的初始温度为 $T_{0}$ ，即

I.C. $T (γ, t < 0) = T_{0}$ (10)

热量通过电解液和电极的界面从电解液传入电极中，设电解液和电极的界面为 $Λ_{T}$ ，并假设 $Λ_{T}$ 上产生不随时间和位置发生变化的温升 $Δ T_{0}$ 。同2.1节，因为不能设置周期性边界条件，所以假设其他面 $Λ_{0}$ 是热绝缘边界。因此边界条件为

B.C. $T (Λ_{T}, t) = T_{0} + Δ T_{0}$ (11) $- n \cdot [- k \nabla T (Λ_{0}, t)] = 0$ (12)

式中， $n$ 是法线方向的单位向量。

Kunanusont等[24]只测量了多孔炭电极的电导率。因此，根据文献[51]，设置多孔炭电极的 $k = 0 . 81^{} W \cdot m^{- 1} \cdot K^{- 1}$ 。假设 $ρ_{f} = 1100^{} k g \cdot m^{- 3}$ ， $c_{p} = 1100 J \cdot k g^{- 1} \cdot K^{- 1}$ ，则 $a = 6.69 \times 10^{- 7}^{} m^{2} \cdot s^{- 1}$ ， $τ_{a} = 1.01 \times 10^{- 6}^{} s$ ， $Δ T_{0} = 0 . 1^{} K$ 。模拟时间为 $(10^{- 3} ~ 10^{- 1})^{} τ_{a}$ ，其中 $τ_{a} = H^{2} / a$ 为表征厚度为 $H$ 的平板电极在表面加温后整个电极内部的温度达到均一时的时间[51]。

电极导热过程同样基于COMSOL Multiphysics 5.4软件，采用有限元法求解。计算平台为4 x Intel(R) Xeon(R) Platinum 9242 CPU at 2.30 GHz和384G内存的服务器。

4 三维多孔电极中的传质与传热原理

本节首先验证了改进后的模拟退火算法的高效性，然后模拟了如图6所示的多孔电极模型中的离子传递和电极导热现象。

4.1 改进的模拟退火算法结果

为了验证改进后的模拟退火算法重构的准确性，图7展示了重构结构w的能量函数E(w)随迭代次数的变化曲线。结果表明相较于随机的像素交换策略和静态的退火系数的经典模拟退火算法，采用DPN值的像素交换策略和采用动态的退火系数的改进模拟退火算法在相同的迭代次数下的重构精度更高(即E(w)更小)。但由于每次迭代需要计算每个像素点的DPN值，所以花费的总时间也更长。但是，由于改进的模拟退火算法在每次迭代时E(w)下降的值更大(即斜率更大)，所以单位时间内E(w)的减小量可能也更大（即重构效率更高）。

图7

图7 重构结构的能量函数随迭代次数的变化

Fig.7 The energy function of the reconstructed structure varies with the number of iterations

4.2 离子传递

阴阳离子在电势的驱动下，分别朝着正极和负极移动。 $10^{- 1}^{} τ_{R C}$ 时整个超级电容器的电势分布 $ϕ / V_{m}$ 如图8(a)所示，正极部分的孔隙相为高电势，负极部分的孔隙相为低电势，体相的电势则整体线性变化。图8(b)、(c)分别展示了正极部分的孔隙相与体相的截面处( $x = - L$ )以及此截面前一个小立方体边长距离的截面处( $x = - L + 0.01 H$ )的 $ϕ / V_{m}$ 。如边界条件设置一样，多孔电极的骨架相和体相的接触面的电势 $ϕ$ 始终为 $V_{m}$ ，即图8(b)中红色区域 $ϕ / V_{m} = 1$ 。此外，电势在yz平面并不是均匀的，骨架相附近的电势较高，且电势随着x的变化而存在较大的变化。

图8

图8 $10^{- 1} τ_{R C}$ 时不同位置的电势分布 $ϕ / V_{m}$

Fig.8 The potential distribution $ϕ / V_{m}$ at different position at $10^{- 1} τ_{R C}$

图9展示了阴阳离子浓度分布的归一化平均值 $(c_{+} + c_{-}) / (2 c_{0})$ 。因为模拟时间只有 $10^{- 1} τ_{R C}$ ，所以阴阳离子主要朝与体相接触的骨架相截面迁移，较少地往孔中迁移。此外，结果还表明，阴阳离子更倾向于从截面的边缘往中心迁移，如图9(b)所示，边缘的 $(c_{+} + c_{-}) / (2 c_{0})$ 大于中心的。而迁移的离子来自于体相对应的位置，如图9(c)所示，边缘的 $(c_{+} + c_{-}) / (2 c_{0})$ 小于中心的。

图9

图9 $10^{- 1} τ_{R C}$ 时不同位置的浓度分布 $(c_{+} + c_{-}) / (2 c_{0})$

Fig.9 The concentration distribution $(c_{+} + c_{-}) / (2 c_{0})$ at different position at $10^{- 1} τ_{R C}$

4.3 电极导热

图10为多孔电极中的归一化温度分布 $[(T - T_{0})] / Δ T_{0}$ 随时间变化的等值面图，包括在 $10^{- 3} τ_{a}$ 、 $10^{- 2} τ_{a}$ 和 $10^{- 1} τ_{a}$ 时刻下。图中黄色区域，即 $[(T - T_{0})] / Δ T_{0} = 0$ ，代表低温;白色区域，即 $(T - T_{0}) / Δ T_{0} = 1$ ，代表高温。当 $(T - T_{0}) / Δ T_{0} = 0$ 时，表明热量还没传导到这个位置；当 $(T - T_{0}) / Δ T_{0} = 1$ 时，说明此时温度已经上升到能达到的最高温度 $T_{0} + Δ T_{0}$ 。结果表明：在 $10^{- 1} τ_{a}$ 时，多孔电极的整个骨架相几乎都是白色，表明从导热边界传递给骨架相的热量已经完全传递到各个角落，骨架相的温度已经均匀分布。此时的时间 $10^{- 1} τ_{a}$ 远小于平板电极的温度弛豫时间 $τ_{a}$ ，原因是热量不仅仅是从骨架相与体相接触的界面传导入骨架相的，还包括孔道内部的壁面，从而减小了实际的导热距离，最终使在多孔电极中的热弛豫时间远小于平板的弛豫时间 $τ_{a}$ 。

图10

图10 不同时刻下的温度分布 $(T - T_{0}) / Δ T_{0}$ 的等值面图

Fig.10 Theisosurface diagram of temperature distribution $(T - T_{0}) / Δ T_{0}$ at different time

5 结论

（1）一种改进的模拟退火算法重构出的多孔电极结构是稳定且有效的。

（2）多孔电极内的电势分布不均匀，体相中的电势沿x方向整体线性变化。

（3）当充电时间为 $10^{- 1} τ_{R C}$ 时，阴阳离子主要吸附在多孔电极骨架相与体相的接触面上，且阴阳离子从截面边缘往中心迁移。

（4）由于电解液深入多孔电极内部，导致实际的导热距离远小于多孔电极厚度，于是多孔电极中的热弛豫时间远小于平板的弛豫时间 $τ_{a}$ 。

采用更加真实的物理模型和求解方法用于真实多孔电极中的离子传递和电极导热现象是未来的重点研究方向。PNP方程并没有考虑离子的体积排阻效应，而是假定所有离子是点粒子。离子体积在多孔电极内，尤其是高浓度和高电压时，会极大地影响离子传递现象。为此，Kilic等[52]推导出包含体积排阻效应的MPNP方程。进一步地，为了更全面地弥补连续介质模型在介观尺度下对微观信息的缺失，还可以采用多尺度多物理场经典密度泛函理论代替MPNP方程[7]。目前通过有限元方法求解含有复杂边界条件的偏微分方程组需要消耗大量计算资源，难以深入研究多孔电极结构对热-质传递现象的影响。因此，亟需快速准确求解含有复杂边界条件的偏微分方程组的方法，比如格子Boltzmann方法[53](lattice Boltzmann method，LBM)和神经网络(neural network，NN)[54]。

符号说明

a	热扩散系数，m2·s-1
$C2j$	第j相的两点聚类函数
cp	骨架相的比定压热容，J·kg-1·K-1
c+,c-,c0	分别为阳离子、阴离子的浓度和初始浓度，mol·m-3
D	离子扩散系数，m2·s-1
d	分形维数
E,Eth,Emin,Eav	分别为能量函数、退出迭代的能量阈值、迭代过程中能量函数的最小值和平均值
ΔE,ΔEth	分别为能量增量和退出迭代的能量增量阈值
e	基本电荷， $C$
H	电极厚度， $m$
$Irj$	多孔介质第j相中 $r$ 点上的特征函数
J	离子通量，m3·s-1
k	热导率，W·m-1·K-1
kB	Boltzmann常数，J·K-1
L	体相半长，m
$L2j$	第j相的两点概率函数
l	孔隙相的传输路径总长度， $m$
NA	阿伏伽德罗常数
Ncon,Niter	分别为连续拒绝随机产生的新结构的次数的上限值和模拟退火算法可以迭代的最大次数
Np	孔隙相立方体的数量
Nt	多孔电极所有立方体的数量
n	法线方向的单位向量
q	电荷数密度，m-3
$S2j$	第j相的两点线性路径函数
S+,S-,S0	分别为带正电、带负电和不带电的边界
T,T0	分别为局部温度和体系初始温度， $K$
Tk-1,Tk	分别为第k-1、k次迭代时的温度， $K$
ΔT0	温升， $K$
t	时间， $s$
V	多孔介质某个相所占的区域
Vm	热电势， $V$
Vp	孔隙相体积， $m3$
w0,w,w′	分别为参考模型的二维结构、模拟退火算法的初始三维结构和下一步生成的三维结构
z	离子的价态
γ	骨架相空间
ε	长度尺度， $m$
εr	相对介电常数，F·m-1
ε0	真空介电常数，F·m-1
ΛT, Λ0	分别为温升边界和绝热边界
λ,λmin,λmax	分别为退火系数、退火系数的最小值和最大值
λD	德拜长度， $m$
ρbulk	体相离子数密度，m-3
ρf	骨架质量密度，kg·m-3
τ	曲折因子
τa	单板电极中热传导的弛豫时间， $s$
τRC	单板电极中离子传递的弛豫时间， $s$
ϕ	电势， $V$
φj,φ	分别为多孔介质第j相的浓度和多孔介质的孔隙率
Ω	离子运动空间
上角标
j	多孔介质某个相的序号
下角标
av	平均值
con	新结构连续被拒绝
f	骨架相
iter	迭代
k	迭代次数
min	最小值
r	相对
r	图像中任意一点的位置
th	阈值
-	带负电
+	带正电
±	阴阳离子

关键字：优秀论文

上一篇：气固流化系统多尺度跨流域EMMS建模
下一篇：复杂物联网联盟链系统通信机制研究

栏目分类

热门排行

推荐信息

期刊知识