锂电解槽极间通道内两相流场数值模拟

金属锂广泛应用于新能源汽车、火箭燃料、合金等领域[1-4]，尤其是电池材料行业对金属锂的需求日益增长

在金属锂制备过程中，由于Li+的还原电势远低于H+，而常见的锂化合物在热力学上过于稳定，以至于无法在温和条件下生成锂金属，因此金属锂的制备通常采用熔盐电解法[5]

锂电解槽是电解质、金属锂和氯气构成的液-液-气三相流动体系，其中金属锂以极低的速度上升无法对速度场进行干扰，研究重点为氯气与电解质形成的两相湍流流动

电解质液相流动会影响电解槽内部温度分布，而氯气以气泡的形式运动在流场中影响电解槽中二次反应的强度，两者均会导致电解槽流场不稳定

一旦流场不稳定发生，则槽内循环流动被干扰，氯气无法顺利排出，造成锂的二次氧化，影响锂盐混合物的溶解，最终导致电解槽无法运行

因此，流场是电解槽中重要的物理场之一，直接决定电解槽效率

流场内气泡微观运动捕捉及速度场测定的困难使得实验研究难以进行，而数值模拟方法具有高度精确性及稳定性，已广泛应用于电解过程中复杂流场预测

在数值模拟中，预测气泡运动的方法可以分为三种：流体体积法(VOF)、Euler-Euler法和Euler-Lagrange(或DPM)法

VOF方法主要用来捕获气泡-电解质形成的界面，但它需要足够精细的网格以追踪较小的气泡，计算资源需求较大，不适合对工业电解槽中的小气泡群进行跟踪

目前，该方法主要用于分析预测单气泡或尺寸较大的气泡，研究其对于电解槽内流场分布及电解质与金属界面稳定性的影响[6-10]

在Euler-Euler法中，连续相和分散相均被视为连续相，体积分数被认为是空间的连续函数[11-13]

Euler-Lagrange法将气泡作为离散质点，分析气泡群对流场的影响[11-15]

相对于Euler-Euler法，Euler-Lagrange法所需的建模工作更少，因为其通常不考虑气泡间的相互作用

该方法仅适用于稀悬浮液，即离散相的体积分数在10%~15%之间[16-19]

以欧拉-拉格朗日方法为基础建立的DPM模型属于离散模型，对于连续相采用欧拉坐标系描述，将气泡视为离散相，通过拉格朗日微分方程来求解离散相颗粒的轨道

在工业电解过程中，气泡多尺度运动使得数值模拟具有挑战性

POZZETTI等[20]考虑到体积和流体的细尺度，开发了一种多尺度的DEM-VOF方法，利用VOF将DEM求解与多相流体求解耦合起来，利用双网格方法解决了不同尺度下的多相流问题，粗网格用于在体积尺度上进行CFD-DEM耦合

在流体细尺度上，采用细网格法求解流体相控制方程

DEM法是一种显式求解的数值分析方法，单元之间的力可以通过力和位移求出

核心是颗粒与颗粒之间相互作用的接触模型

其中，所有颗粒的平移、旋转、都是通过求解牛顿运动方程来完成

部分学者[11-13, 15, 21-25]分析了不同电解槽结构及操作参数对气泡层厚度、气液两相流及气泡分布特征的影响

在早期大多数学者采用二维模型研究气体和电磁力对流场的影响[26-28]

随着计算方法的进步，数值模拟已经趋于三维模型全耦合建模[29-31]

SUN等[29]通过CFD仿真研究了镁电解槽中电磁场作用下的液态镁、熔融电解质和氯气三相流行为

SUN等[30]在单阳极实验室电解槽模型中，针对气泡行为和熔池金属MHD提出了多尺度多相流模型：通过VOF模型模拟气泡-电解质-金属三相流，采用DPM法跟踪阳极底部形成的微气泡，结果表明，熔池金属流动可以加速气泡运动并降低气体覆盖率

詹水清等[31]研究了气泡及磁场对流场的影响，发现气体在局部位置对氧化铝的均匀分布有一定影响，但电磁力对全槽氧化铝浓度分布起到了决定性作用

在锂电解槽数值计算研究中，SUN等[32]测试了用于选择和验证合适CFD模型的PIV实验

根据优化结果设计了工业锂电解电池的冷模型，并进行了冷模型实验

OLIAII等[33-35]使用OpenFOAM求解器研究了气举式锂电解槽的传质和耗能状态，通过建立电场-流场-物质传递耦合模型，研究阴阳极距离对电流密度和离子分布的影响；之后又建立带有隔膜的二维轴对称模型，研究了隔膜长度、位置以及孔隙率对电场强度和两相流的影响，尽管隔膜的存在会带来额外的电压降，但有利于阻止电解槽中的二次反应

本文以工业中应用的30 kA锂电解槽为研究对象，建立电解槽内电解质-气泡两相流场，重点关注侧面极间通道和中间极间通道中的y向电解质速度和体积分数分布，同时考虑电流密度分布不均匀的影响，考察极距、阳极半径等设计参数和电解质高度、电流等操作参数对流场分布的影响

1模型建立1.1物理模型的建立锂电解槽物理模型包括石墨阳极，阴极，高、低水平连接板，阴极绑带，电解质，内衬，耐火砖，槽壳，槽盖，槽盖上设有出锂口和排氯口，具体三维结构如图1所示

图1锂电解槽三维模型



Fig. 13D model of lithium electrolytic cell工业锂电解槽一般为多阳极结构，构造复杂

采用全模型流场计算需要大量资源，因此沿阳极中心线对三维电解槽进行侧面剖分，侧面剖分方式见图2

图2三维电解槽侧面剖分展示



Fig. 23D model of lithium electrolytic cell side division show：(a) Side view; (b) Top view; (c) Left view图2红线表示相交线，蓝色平面为切割平面

所得二维模型如图3所示

图3锂电解槽二维模型Fig.32D model of lithium electrolysis cell图3中用蓝色虚线所强调的边界1~4为气体出口边界，用橙色虚线所强调的边界5~10为气体入口边界

其他未强调的黑色实线为其他边界

1.2参数确定1.2.1电解质及氯气密度、黏度及表面张力的确定电解时，电解槽内发生氧化还原反应，即氯化锂离解为氯离子和锂离子

在电场的作用下，锂离子向阴极运动，并在阴极发生还原反应生成金属锂单质，由于其密度远小于电解质密度，金属锂聚集在电解质表面

氯离子向阳极运动，在阳极发生氧化反应生成氯气，气泡逐渐变大，上升过程中发生聚并，最后破裂逸出

由于气泡上升产生推动力，电解质能够在电解槽中循环流动

流场中流体的密度、黏度和表面张力是数值计算中不可或缺的参数，计算方法如下

流体黏度受温度和压强的共同影响

随着压力的增大，流体黏度逐渐增加，但总体上变化不大；黏度受温度的影响较为明显，且不同流体状态变化不同，气体黏度系数随温度升高而增大，液体黏度系数随温度升高而减小

目前工业生产中采用的方法是420 ℃时电解熔融LiCl(60%)-KCl(40%)混合熔盐，VAN ARTSDALEN等[36]对不同比例LiCl-KCl体系的熔盐密度进行了实验研究，得到了密度与温度的关系，在温度为395~593 ℃时，密度与温度的关系式如下：(1) 当t=429 ℃时，熔盐电解质的密度为1663.55 kg/m3

熔融状态下KCl和LiCl的表面张力分别为95.8 mN/m和123 mN/m，由此可得电解质的表面张力为：(2) mN/m氯气及电解质的动力学参数如表1所示

表1氯气及电解质的动力学参数Table 1Kinetic parameters of chlorine and electrolyteComponentDensity/(kg?m-3)Dynamic viscosity/(Pa?s)Surface tension/(mN?m-1)Chlorine[37]3.22.95×10-521.9Electrolyt16641.2×10-3112.121.2.2主要导体电导率1) 阴极及阴极绑臂电导率阴极及阴极绑臂材料为低碳钢，其电导率符合线性电导率[38]，如式(3)所示：(3) 式中：为普通钢在参考温度下的电阻率(Ω?cm)，普通钢在20 ℃时的电阻率(0.13~0.3)×10-4 Ω?cm，这里取0.3×10-4 Ω?cm；为参考温度(℃)；为电阻率温度系数(℃-1)，这里取=0.005 ℃-1；t为阴极及阴极绑臂的温度(℃)

式(3)可以简化为：(4)}主要导体电导率汇总如表2所示

表2主要导体的电导率Table 2Conductivity of main conductorParameterConductor/(S?m-1)Molten salt electrolyte187Anode1.36×105Cathode jointing plate4.032×1061.2.3气泡直径的确定电解条件下产生气泡的机理是：电活性氯离子转移到电极，转化为分子物质氯气，且在反应器的操作条件下为气态

当其局部浓度过饱和时通过非均相成核形成气泡：电极表面预先存在的气核由于周围过饱和电解液中溶解气体的转移而生成气泡[17]

气泡内部压强由开尔文方程计算：(5) 式中：pa是与气泡相接触的液体压强(N/m)；σ是气泡的表面张力(N/m)；r是气泡半径(m)

电解质中气泡形态各异，计算脱离阳极表面的半径时，假设气泡呈半球形附着在阳极表面[39]

气泡沿着阳极表面向上运动的过程中，受到的表面张力向下分量假设为[40]：(2σcos45°)/r

假设气泡离开阳极时半径为r(mm)，受到的浮力为：(6) N/cm当时，(7) (8) mm2数学模型的建立锂电解槽中，垂直电极间界定了一个具有高气体体积分数的区域，气相可视为连续相

因此，选择欧拉模型进行数值模拟

湍流流动的数值模拟可分为直接数值模拟(DNS)、雷诺平均纳维-斯托克斯(RANS)方程模拟以及大涡模拟(LES)

DNS直接对方程进行求解，会解析所有尺寸大小的涡，但网格需要极为细致

湍流的本质是三维、非定常、多尺度，因此所需要的计算成本高

LES同DNS均模拟非定常的湍流，网格大小则需适中，能够将小尺度流动过滤掉

RANS将非定常的湍流问题转化为定常问题，对速度和压力场取时间平均值，只寻求平均意义下的流动

工业中广泛应用的模型有、、剪切应力输运(SST)湍流模型以及 Spalart-Allmaras 湍流模型

在众多RANS湍流模型中，最常用的是湍流模型

该模型不仅具有较好的稳定性，还能有效节省计算资源，因此常被用于工业应用领域

模型采用壁面函数来描述近壁流动，因此其计算成本低

此电解槽中湍流流动的数值模拟采用标准模型[41]

2.1控制方程由于电解质的流动情况十分复杂，计算过程中无法考虑全部因素，所以必须进行适当简化：1)材料的物理性质保持稳定；2)电解液为不可压缩流体；3)忽略气泡的聚并、破裂等，认为直径是常数；4)电解质-氯气两相之间不存在质量交换；基于以上假设，流体的质量守恒方程为：(9) 式中：下标k为c时表示连续相(continuous phase);下标k为d时表示不连续相(disperse phase)

动量守恒方程为：(10) (11) (12) 式中：是应力张量；是湍流黏度；是动力黏度；Fb是体积力，通常包括阻力、升力和虚拟质量力

与阻力相比，升力和虚拟质量力的影响微不足道[42]，由式（13）、（14）、（15）给出：(13) (14) (15) 此外，标准方法通过两个附加变量的解来表示流体湍流特性[43]，附加变量分别是湍动能k和动能耗散率ε

湍动能k方程为：(16) 动能耗散率ε方程为：(17) 式中：(18) (19) (20) 上述方程中出现的常数如表3所示

表中Ck、Cε为模型常数，这两个值具有很高的相关性，可通过实验数据进行调整，以获得实验和仿真结果较好的一致性

若Ck等于零表示忽略气泡引起的湍流

表3标准k-ε方法中的常数Table 3Constants of standard k-ε approachCε1Cε2CμσkσεCkCε1.441.920.0911.31.001.462.2边界条件二维流场的边界条件如表4所示，表中的边界序号对应图3中的边界编号

其中：ud表示气体速度；uc表示液体速度

表4流场边界条件Table 4Boundary conditions for velocity fieldBoundary No.Boundary conditions1-4ud= ud, uc?n=05-10Othersud=0, uc=0 (wall)1) 边界1~4为气体出口，设定液相滑移边界条件，气体速度不变，ud=ud；液体速度法向为零，

2) 边界5~10为气体进口，阳极表面氯气气泡生成量与阳极表面电流密度成正比，氯气生成量Q根据法拉第定律计算：(21) 式中：M是氯气的相对摩尔质量，M=7.1×10-2 kg/mol；z是电荷数，z=2；F是法拉第常数，F=96485 C/mol；j为阳极电流表面密度(A/m2)

3) 其他边界为壁面，设定无滑移边界条件，气体速度ud=0；液相速度uc=0

2.3网格独立性检验数值模拟的实质是已知初值和边值后求解偏微分方程，由于计算机的位宽限制及采用有穷表示无穷，因此数值模拟会分别产生舍入误差和离散误差

划分网格是将连续问题离散化，最终得到精确解，理论上网格越密，离散误差越小；总误差是离散误差和舍入误差的总和[44]，但由于舍入误差持续增加，且网格不能无限加密，因此，需要在离散误差和舍入误差之间找到一个平衡点，使得总误差最小，即实现网格的独立性

为了解析气泡浓度的陡峭程度，入口附近的网格必须非常细化，所以添加了边界层网格进行求解

由于氯气的存在形成了比较复杂的流动模式，因此反应器内部的网格也应该细化

整体采用结构化网格对二维流场模型求解域进行划分，由于电解质中阳极表面附近流场分布梯度大，因此在此处添加边界层网格，边界层数为8层，第一层厚度为1.3 mm，厚度调节因子为1.2

最终网格最大控制尺寸为5 mm，网格数量为59968

具体划分结果如图4所示

图4二维流场区域网格



Fig. 42D flow field area grid网格最大尺寸的长度从20 mm减小至5 mm，相应的网格数量从14608增至94159

将80 s后的电解质最大速度作为判定依据，详细数据见表5

表5电解质最大速度与网格数的关系Table 5Relationship between maximum electrolyte velocity and number of gridsNo.Max size/mmMin size/mmMesh quantityMaximum electrolyte velocity/(m?s-1)1204146080.5472104241490.5343102404350.522454599680.511552941590.510当网格数量增加至59968时，电解质最大速度逐渐趋于稳定；网格数量从59968继续增加至94159时，电解质最大速度变化很小，即当网格总数达到59968时，速度精度并不会随着网格总数的增长有显著提升，说明此时网格数量已满足计算要求

选用COMSOL Multiphysics中气泡流接口进行求解，此物理场接口不考虑单个气泡的界面，而是追踪气相的局部体积分数，最终求得液相速度

整体采用瞬态求解器分步计算，具体到因变量时采用GMRES求解器，再次开始前的迭代次数为50；对于多重网格采用平滑聚集AMG求解器，V循环方式，迭代次数为2

总求解时间设为80 s，步长设为0.1 s，时间步进采用向后差分公式

3计算结果与分析3.1阳极电流分布图5所示为电解液中的电流密度分布，可以发现绝大多数电流都集中在阴阳极之间，即电流从阳极工作面流经极间通道直接流至阴极工作面；只有一小部分电流从阳极流出，经由电解质，分别从顶部、侧面和底部回到阴极

图5中数字1~3为模型中阳极编号

图5电解质内电流密度分布



Fig. 5Current density in electrolyte根据法拉第定律可知，气泡的产生与阳极电流密度相关，故提取三维电热场耦合模型阳极的边线电流，各阳极电流密度如图7所示

图5中数字1~3分别对应图7中Anode 1~Anode 3

以Anode 1为例说明，见图6

图6x轴坐标选取示意图



Fig. 6Sketch map of selectioncoordinate x图7阳极电流密度



Fig. 7Current density in anode: (a) AB side; (b) CD side为表示图7中阳极表面AB侧和CD侧随阳极高度变化的电流密度，在此，将x轴如图6处理，图7中x坐标与图6中黄色箭头标示的距离一致

由图7(a)和(b)可见，不同阳极电流密度变化趋势一致，可以忽略由于阳极位置造成的电流密度差异

每个阳极表面电流密度不均匀，随着阳极浸入电解质的深度增加，电流密度逐渐下降

电解质中阴阳极正对区域为工作区域，电流密度在工作区域的范围为7000~6000 A/m2，在工作面顶部达最大

但是，事实上阳极顶端A、D所在水平位置比底端B、C水平位置气泡聚集更多，前者电极间电阻更大

在流场仿真中，并未考虑气泡对电解液电阻的影响，所以出现了A、D处电流密度却比B、C处更大

由于在阳极边缘处存在尖端放电效应，因此阳极边缘的电流密度值远大于工作区域的电流密度值

3.2液相速度分布和体积分数分布在气体逸出电极上，实际电流密度始终大于计算电流密度，实际电流流过的面积小于几何电极表面积，这是因为部分电极表面被黏附的气泡占据，并且具有电化学惰性[42]

根据法拉第定律，阳极表面电流密度影响氯气的生成速率

为研究电流密度对两相流的影响，决定建立不同模型：分别将阳极电流密度设为均匀边界和不均匀边界

均匀电流密度值由式(22)计算得到，不均匀电流密度值根据图7得到

(22) (A/m2)图8所示分别为均匀模型和不均匀模型中，锂电解槽二维模型的液相速度分布和体积分数分布云图

如图8(a)所示，不同通道间的不同位置用缩写表示，比如C1L表示侧面极间通道底部位置

分别研究侧面极间通道(C1)和中间极间通道(C2)的底部(L)、中部(M)和顶部(H)的电解质速度和体积分数

与C1和C2有关的坐标见表6

表6中，(x1, Δx, x2)表示x轴坐标起点是x1，x轴坐标终点是x2，跨度为Δx

图8电解质速度和体积分数的分布



Fig. 8Distributions of electrolyte velocity((a), (b)) and volume fraction((c), (d)): (a), (c) Uniform model; (b), (d) Non-uniform model表6参数坐标Table 6Parameter coordinatesChannelCoordinate of x(Abscissa range)/mmCoordinate of y(Height)/mmC1H(325,40,365)700C1M(325,40,365)350C1L(325,40,365)0C2H(810,40,850)700C2M(810,40,850)350C3L(810,40,850)0对比图8(a)和(b)发现，气泡上升时，阳极附近电解质随之向上运动，两模型中的电解质速度分布云图相似

在侧面极间通道C1中，电解质从通道流出到达自由液面后，转向180°后沿着侧壁向下运动，在阴极与侧壁组成的通道间形成循环

中间极间通道(C2)中电解质流动与C1完全不同，电解质从C2H流出到达自由液面后，沿阴极壁循环进入同一通道中

且电解质在C2中的流动速度比C1中慢，这是因为C2中同时存在向上和向下的运动

对比图8(c)和(d)可知，两图中体积分数分布云图相似，但均匀电流密度模型的体积分数最大值大于非均匀电流密度模型，这是两者电流密度值的差异所造成的

由于两通道间电解质速度不同，两图中通道C2的体积分数明显大于通道C1

为了清楚地反映不同通道间不同高度电解质速度分布的区别，C1、C2中L、M和H高度的液相速度曲线如图9所示

图9通道C1和C2中均匀模型和非均匀模型的电解质速度



Fig. 9Electrolyte velocities of two channels in uniform and non-uniform models: (a) C1; (b) C2对比图9(a)和(b)可以发现，任意高度处均匀模型的液相速度均大于非均匀模型的液相速度

两个通道内的速度模型完全不同：C1中，由于L处体积分数很小，产生的曳力不足以使电解质流动，通道间的循环作用是L处液相流动的主要原因；液相在y方向的最大速度出现在通道中间，但速度值不超过0.35 m/s；而在C1M和C1H中，液相在y方向的最大速度出现在阳极附近，H位置的最大速度值接近0.5 m/s，M位置的最大速度值接近0.45 m/s，H和M处的最大速度均大于L处

C2的速度曲线分布完全不同于C1，C2中存在负y方向电解质速度，且C2M位置液相在y方向的最大速度达到0.15 m/s，远高于C2H的最大速度值0.05 m/s

这是因为在C2H上方液相形成一个循环，电解质在H处的流动受到影响；C2L中，液相速度几乎为0，一方面是由于气体含量低，曳力作用不足，另一方面是由于没有相邻通道间的电解质作为补充，所以电解质速度很低

图10所示为不同模型两通道中不同位置的体积分数，均匀模型的体积分数高于不均匀模型，但两模型中体积分数分布曲线形状几乎一致

两通道中，随着通道高度的降低，即从H到L处，体积分数随之降低，到达L处时，体积分数几乎为0

在同一高度处，越靠近阳极，体积分数值越大，到达阳极附近时，两通道中的体积分数达到最大值，为12%

在M和H高度处，C2的体积分数均大于C1

电解质的速度决定气泡的停留时间，C2中由于电解质存在向下的运动，速度较低，氯气停留时间较C1长，因此均匀模型和非均匀模型的体积分数之间的差异更加明显

图10通道C1和C2处均匀模型和非均匀模型的体积分数



Fig. 10Volume fractions of two channels in uniform and non-uniform models: (a) C1; (b) C2根据上述分析可知，阳极表面电流密度影响电解槽内电解质速度和体积分数的分布，后续研究将通过电场得到的电流密度作为流场初始边界条件进行计算

3.3工艺参数对两相流场的影响如今工业生产中锂电解基本采用30 kA电流进行操作，因此下面以电流30 kA为例，分析极距、阳极半径及电解质液面高度对两相流场的影响

考虑到工业生产中为了提高电流效率而增大电流，为此，探究了34 kA、38 kA电流对两相流场的影响，同时，并与30 kA电流共同绘图进行比较

3.3.1极距通道C1中不同高度处的体积分数和电解质速度随极距的变化如图11所示

图11(a)、(b)和(c)分别显示了通道C1中H、M和L高度的体积分数随极距的变化规律

随着极距的缩小，体积分数随之增大：C1H中，体积分数从11.51%增加至11.67%；C1M中，体积分数从11.71%增加至11.83%；C1L中，体积分数从0.011%增加至0.013%；由于不同极距给定的电流密度相同，产生的气体总量保持一定，所以随着极距的减小，即空间缩小，通道内的体积分数增加

图11(a′)、(b′)和(c′) 分别显示了通道C1中H、M和L高度的电解质y方向速度随极距的变化情况

随着极距的减小，电解质y方向速度也在增加：C1H中，速度从0.403 m/s增加至0.408 m/s；C1M 中，速度从0.312m/s增加至0.339 m/s；C1L中，速度从0.0256 m/s增加至0.0362 m/s；原因是两者动量交换增加，而阳极气体抽空引起的液体能耗也随之减少，引起速度增加

图11通道C1中体积分数和电解质速度随极距的变化



Fig. 11Changes of volume fraction and electrolyte velocity in channel C1 with ACD: (a), (a′) C1H; (b), (b′) C1M; (c), (c′) C1L通道C2中不同高度处的体积分数和电解质速度随极距(ACD)的变化如图12所示

图12(a)、(b)和(c)分别显示了通道C2中H、M和L高度的体积分数随极距的变化规律

同样地，通道C2中体积分数随极距的减小而增加：C2H中，体积分数最大值从10.78%增加至11.46%；C2M中，体积分数最大值从11.05%增加至11.19%；C2L中，体积分数最大值从0.12%增至0.125%；与通道C1中体积分数的规律呈现一致

对于通道C2中的电解质y向速度，图12(a′)、(b′)和(c′) 分别显示了通道C2中H、M和L高度的电解质y方向速度随极距的变化情况

在C2H和C2M中靠近阴极的位置出现y轴负方向的速度，说明电解质会在同一通道中循环流动

不同高度处电解质y方向最大速度都随着极距的减小而减小，C2H中，最大速度从0.136 m/s降低至0.0791 m/s；C2M中，最大速度从0.255 m/s降低至0.178 m/s；C2L中，最大速度从0.00448 m/s降低至0.00397 m/s

极距的减小导致气泡在电解质内部移动困难，从而使得气泡本身的运动速度减小，而气泡运动产生的浮力是电解质运动的主要外力，因此电解质的运动速度随之减小

图12通道C2中体积分数和电解质速度随极距的变化



Fig. 12Changes of volume fraction and electrolyte velocity in channel C2 with ACD: (a), (a′) C2H; (b), (b′) C2M; (c), (c′) C2L3.3.2阳极半径通道C1中不同高度处的体积分数和电解质速度随阳极半径的变化如图13所示

由图13(a)、(b)和(c)可知，随着阳极半径的增加，不同高度的体积分数均有所下降：C1H中，体积分数最大值由11.84%下降至11.17%；C1M中，体积分数最大值由12.02%下降至11.53%；C1L中，体积分数最大值由0.024%下降至0.011%；体积分数与电流密度大小成正比，电流密度在总电流一定的情况下随着阳极半径的增加而下降，所以体积分数也随之减少

在靠近阳极位置，H高度处体积分数降幅最大，这是因为H处离自由液面最近，气泡能够比在阴阳极间通道内空间更大的位置上浮

由图13(a′)、(b′)和(c′)可知，由于体积分数的减少，气液两相之间的动量交换减少，所以电解质y方向速度随着阳极半径的增加而减小：C1H中，最大速度由0.461 m/s下降至0.425 m/s；C1M中，最大速度由0.396 m/s下降至0.367 m/s；C1L中，电解质y方向最大速度由0.306 m/s下降至0.283 m/s；且不同高度处速度变化较为明显的地方发生在通道的中心位置

图13通道C1中体积分数和电解质速度随阳极半径的变化



Fig. 13Changes of volume fraction and electrolyte velocity in channel C1 with r_anode: (a), (a′) C1H; (b), (b′) C1M; (c), (c′) C1L图14所示为通道C2中体积分数和电解质y方向速度随阳极半径的变化，可以发现，通道C2与C1呈现相同的变化规律，随着阳极半径从150 mm增加至190 mm时，不同高度处的体积分数和电解质y方向速度都随之减少：H处的体积分数最大值从11.61%降低至11.05%，电解质y方向最大速度从0.0562 m/s降低至0.0538 m/s；M处的体积分数最大值从11.34%降低至10.36%，电解质y方向最大速度从0.143 m/s降低至0.140 m/s；L处的体积分数最大值0.18%降低至0.033%，y轴负方向电解质最大速度从0.00268 m/s降低至0.00223 m/s

在C2H和C2M中，电解质速度减小幅度很小，不同阳极半径下的速度变化几乎重合

一方面是因为动量交换减少，导致电解质速度下降；另一方面是因为能量耗散减少，导致电解质速度上升；两者共同作用使得电解质y方向速度变化不是很明显

图14通道C2中体积分数和电解质速度随阳极半径的变化



Fig. 14Changes of volume fraction and electrolyte velocity in channel C2 with r_anode: (a), (a′) C2H; (b), (b′) C2M; (c), (c′) C2L3.3.3电解质液面高度通道C1中H、M和L位置的体积分数和电解质速度随电解质液面高度的变化如图15所示

由图15(a)、(b)和(c)可知，不同高度处体积分数不随电解质液面高度变化而发生改变

由图15(a′)、(b′)和(c′)可知，电解质y方向速度随电解质液面高度的增加而增加：当电解质液面总高度由1140 mm增加至1240 mm时，H处的电解质y方向最大速度由0.457 m/s增加至0.466 m/s；M处的电解质y方向最大速度由0.393 m/s增加至0.398 m/s；L处的电解质y方向最大速度由0.301 m/s增加至0.310 m/s

电解质液面高度增加，气泡从极间通道流出具有更大的空间上升至自由液面，产生更大的曳力来加速电解质的流动速度

图15通道C1中体积分数和电解质速度随电解质液面高度的变化



Fig. 15Changes of volume fraction and electrolyte velocity in C1 with h_electrolyte: (a), (a′) C1H; (b), (b′) C1M; (c), (c′) C1L图16(a)、(b)和(c)显示了通道C2中不同高度处体积分数随电解质液面高度的变化规律，与C1一致，随着电解质液面高度增加，体积分数不发生改变

由图16(a′)、(b′)和(c′)可知，通道C2中H、M和L处电解质y方向速度不随电解质液面高度的增加而发生改变，且在H和M处同样存在电解质负y方向速度，说明通道C2中总是存在电解质的回流

电解质速度不变是由于电解质回流导致的湍流能量耗散引起的速度下降和相互抵消，因此改变电解质高度对于通道C2的体积分数和电解质速度几乎没有影响

气泡具有更大的空间，所产生的曳力会引起电解质速度增加

图16通道C2中体积分数和电解质速度随电解质高度的变化



Fig. 16Changes of volume fraction and electrolyte velocity in C2 with h_electrolyte (a), (a′) C2H; (b), (b′) C2M; (c), (c′) C2L3.3.4电流提高锂金属产量最直接有效的方法是增大锂电解槽的电流，在不改变其他条件的情况下，阳极电流密度随着电流增大而增大

分别提取30~38kA下锂电解槽阳极电流密度作为两相流场的初始边界条件，不同电流下阳极电流密度分布曲线如图17所示

图17(a)和(b)所示分别为阳极两侧的电流密度，从图中可以发现，电流的增大主要影响阳极工作区域的电流密度，相比之下，对阳极末端电流密度的影响很小

图17不同电流下阳极表面的电流密度



Fig. 17Current density along anode with different currents通道C1中H、M和L位置的体积分数和液相速度随电流的变化如图18所示

图18(a)、(b)和(c) 分别显示了通道C1中H、M和L处体积分数的分布

随着电流的增大，电流密度增大，根据法拉第定律，气体总量与电流密度成正比，这意味着阳极周围的气泡会逐渐积聚，体积分数因此增加

当电流从30 kA增大至38 kA时，H位置的体积分数最大值由11.84%增加至12.09%，增幅为2.11%；M位置的体积分数最大值由12.02%增加至12.33%，增幅为2.58%；L处的体积分数最大值由0.011%增加0.013%，增幅为18.18%

图18(a′)、(b′)和(c′)分别显示了通道C1中H、M和L处电解质y方向速度的分布

随着电流增大，液相y方向最大速度也在增大：H处的电解质y方向最大速度由0.462 m/s增加至0.503 m/s，增幅为8.80%；M处的电解质y方向最大速度由0.396 m/s增加至0.433 m/s，增幅为9.31%；L处的电解质y方向最大速度由0.306 m/s增加至0.338 m/s，增幅为10.67%

由此可见，对于通道C1而言，电流的增加对L处的体积分数及电解质y方向速度影响最大

这是因为具有较高体积分数的气泡析出可以增强气泡驱动的电解质流动

图18通道C1中体积分数和电解质速度随电流的变化



Fig. 18Changes of volume fraction and electrolyte velocity in channel C1 with current: (a), (a′) C1H; (b), (b′) C1M; (c), (c′) C1L由图19(a)、(b)和(c)可知，通道C2不同高度处的体积分数随电流的增大而增加，且C2中的体积分数随电流的变化远大于C1：H处的体积分数最大值由11.61%增加至11.94%，增幅为2.84%；M处的体积分数最大值由11.34%增加至11.98%，增幅为5.64%；L处的体积分数最大值由0.12%增加至0.16%，增幅为33.33%

由图19(a′)、(b′)和(c′)可知，通道C2中的电解质y方向速度变化远小于C1，且在H处靠近阳极的位置出现了相反趋势，即随着电流的增大，电解质y方向速度反而下降

一方面是因为此处聚集的气泡较多，气泡移动困难，气泡速度下降；另一方面是因为电解质回流会导致耗散掉更多的湍流能量，从而导致较低的速度幅值

图19通道C2中体积分数和电解质速度随电流的变化



Fig. 19Changes of volume fraction and electrolyte velocity in channel C2 with current: (a), (a′) C2H; (b), (b′) C2M; (c), (c′) C2L4结论1) 建立电流密度均匀模型和不均匀模型，对比发现不均匀模型的体积分数和电解质y向速度均小于均匀模型，在后续研究中全部采用不均匀电流密度模型

2) 随着极距降低，两通道中体积分数均增加，电解质y向速度在通道C1中增加，在C2中降低

3) 随着阳极半径增加，两通道内体积分数和电解质y向速度均减少

4) 随着电解质液面高度增加，两通道中体积分数不改变，电解质y向速度在C1中增加，在C2中不发生变化

5) 随着电流增大，两通道内体积分数和电解质y向速度均增加

参考文献

YANG S X, ZHANG F, DING H P, et al. Lithium metal extraction from seawater[J

Joule, 2018, 2(9): 1648-1651.SCHIEMANN M, BERGTHORSON J, FISCHER P, et al. A review on lithium combustion[J

Applied Energy, 2016, 162: 948-965.XU W Q, BIRBILIS N, SHA G, et al. A high specific strength and corrosion resistant magnesium alloy[J

Nature Materials, 2015, 14(12): 1229-1235.WANG L, RAWAL A, QUADIR M Z, et al. Nanoconfined lithium aluminium hydride (LiAlH4) and hydrogen reversibility[J

International Journal of Hydrogen Energy, 2017, 42(20): 14144-14153.ZHANG X, HAN A, YANG Y A. Review on the production of high-purity lithium metal[J

Journal of Materials Chemistry A, 2020, 8(43): 22455-22466.ZHAO Zhi-bin, WANG Zhao-wen, GAO Bing-liang, et al. Anodic bubble behavior and voltage drop in a laboratory transparent aluminum electrolytic cell[J

Metallurgical and Materials Transactions B: Process Metallurgy and Materials Processing Science, 2016, 47(3): 1962-1975.WANG Yu-feng, ZHANG Li-feng, ZUO Xiang-jun. Simulation of the fluid flow-related phenomena in the electrolyte of an aluminum electrolysis cell[J

Metallurgical and Materials Transactions B: Process Metallurgy and Materials Processing Science, 2011, 42(5): 1051-1064.ZHANG Kai-yu, FENG Yu-qing, WITT PJ, et al. SCHWARZ M P. A numerical assessment of bubble-induced electric resistance in aluminium electrolytic cells[J

Journal of Applied Electrochemistry, 2014, 44(10): 1081-1092.ZHANG Kai-yu, FENG Yu-qing, SCHWARZ P, et al. Computational fluid dynamics (CFD) modeling of bubble dynamics in the aluminum smelting process[J

Industrial & Engineering Chemistry Research, 2013, 52(33): 11378-11390.郭学益, 王双, 王亲猛, 等. 氧气底吹熔池熔炼过程气泡生长行为仿真研究[J].中国有色金属学报, 2018, 28(6): 1204-1215.GUO Xue-yi, WANG Shuang, WANG Qin-meng, et al. Simulation research of bubble growth behavior in oxygen bottom blowing smelting process[J

The Chinese Journal of Nonferrous Metals, 2018, 28(6): 1204-1215.ZHAN Shui-qing, WANG Zhen-tao, YANG Jian-hong, et al. 3D numerical simulations of gas-liquid two-phase flows in aluminum electrolysis cells with the coupled model of computational fluid dynamics-population balance model[J

Industrial & Engineering Chemistry Research, 2017, 56(30):8649-8462.ZHAN Shui-qing, HUANG Yu-jie, WANG Zhen-tao, et al. CFD simulations of gas-liquid multiscale flow characteristics in an aluminum electrolysis cell with population balance model: Effect of anode slot configuration[J

JOM, 2019, 71(1):23-33.ZHAN Shui-qing, YANG Jian-hong, WANG Zhen-tao, et al. CFD simulation of effect of interphase forces and turbulence models on gas-liquid two-phase flows in non-industrial aluminum electrolysis cells[J

JOM, 2017, 69(9):1589-1599.HREIZ R, ABDELOUAHED L, FUNFSCHILLING D, et al. Electrogenerated bubbles induced convection in narrow vertical cells: A review[J

Chemical Engineering Research & Design, 2015, 100: 268-281.黄义鹏, 王兆文, 杨酉坚, 等. 铝电解槽阳极形状对气泡排出的影响[J

东北大学学报(自然科学版), 2019, 40(2): 228-233.HUANG Yi-peng, WANG Zhao-wen, YANG You-jian, et al. Effect of anode geometry on bubble behavior in aluminum reduction cell[J

Journal of Northeastern University(Natural Science), 2019, 40(2): 228-233.LITRICO G, VIEIRA C B, ASKARI E, et al. Strongly coupled model for the prediction of the performances of an electrochemical reactor[J

Chemical Engineering Science, 2017, 170: 767-776.HREIZ R, ABDELOUAHED L, FUENFSCHILLING D, et al. Electrogenerated bubbles induced convection in narrow vertical cells: Ply measurements and Euler-Lagrange CFD simulation[J

Chemical Engineering Science, 2015, 134:138-152.WANG Qiang, SUN Mei-jia, LI Bao-kuan, et al. Numerical modeling of effect of slot on bubble motion in aluminum electrolytic process[J

Transactions of Nonferrous Metals Society of China, 2018, 28(8): 1670-1678.LI Lin-min, LI Bao-kuan. Investigation of bubble-slag layer behaviors with hybrid Eulerian-Lagrangian modeling and large eddy simulation[J

JOM, 2016, 68(8): 2160-2169.POZZETTI G, PETERS B. A multiscale DEM-VOF method for the simulation of three-phase flows[J

International Journal of Multiphase Flow, 2018, 99: 186-204.SCHILLINGS J, DOCHE O, DESEURE J. Modeling of electrochemically generated bubbly flow under buoyancy-driven and forced convection[J

International Journal of Heat and Mass Transfer, 2015, 85: 292-299.HAAS T, HILGENDORF S, VOGEL H, et al. A comparison between two cell designs for electrochemical neodymium reduction using numerical simulation[J

Metallurgical and Materials Transactions B: Process Metallurgy and Materials Processing Science, 2017, 48(4): 2187-2194.ZHAN Shui-qing, WANG Jun-feng, WANG Zhen-tao, et al. Computational fluid dynamics-population balance model simulation of effects of cell design and operating parameters on gas-liquid two-phase flows and bubble distribution characteristics in aluminum electrolysis cells[J

JOM, 2018, 70(2): 229-236.HOU Wen-yuan, LI He-song, FENG Yuan, et al. Effects of the application of a perforated anode in an aluminum electrolysis cell on the gas-liquid two-phase flow and bubble distribution characteristics[J

Industrial & Engineering Chemistry Research, 2020, 59(32): 14522-14530.王富强, 刘 伟, 王兆文. 铝电解槽中局部阴极电流增大对电解质-铝液两相流场的影响[J

金属学报, 2020, 56(7): 1047-1056.WANG Fu-qiang, LIU Wei, WANG Zhao-wen. Effect of local cathode current increasing on bath-metal two-phase flow field in aluminum reduction cells[J

Act Metallurgica Sinica, 2020, 56(7): 1047-1056.DOHEIM M A, EL-KERSH A M, ALI M M. Computational modeling of flow in aluminum reduction cells due to gas bubbles and electromagnetic forces[J

Metallurgical and Materials Transactions B: Process Metallurgy and Materials Processing Science, 2007, 38(1): 113-119.XU Y J, ZHANG H L, LI J, et al. A nonlinear shallow-water model combined with gas bubble effect for meltflows and interface instability in aluminum reduction cells[J

JOM, 2013, 65(11): 1459-1466.ZHANG H L, YANG S, ZHANG H H, et al. Numerical simulation of alumina-mixing process with a multicomponent flow model coupled with electromagnetic forces in aluminum reduction cells[J

JOM, 2014, 66(7): 1210-1217.SUN Z, LI P, LU G M, et al. Effect of electromagnetic field on three-phase flow behavior[J

Industrial & Engineering Chemistry Research, 2010, 49(21): 10798-10803.SUN M J, LI B K, LI L M. Multiscale simulation of bubble behavior in aluminum reduction cell using a combined discrete bubble model volume of fluid magnetohydrodynamical method[J

Industrial & Engineering Chemistry Research, 2019, 58(8): 3407-3419.詹水清, 李 茂, 周孑民, 等. 铝电解槽熔体内氧化铝浓度分布的数值模拟[J

中国有色金属学报, 2014, 24(10): 2658-2667.ZHAN Shui-qing, LI Mao, ZHOU Jie-min, et al. Numerical simulation of alumina concentration distribution in melts of aluminum reduction cells[J

The Chinese Journal of Nonferrous Metals, 2014, 24(10): 2658-2667.SUN Ze, NI Hai-ou, CHEN Hang, et al. Designing and optimizing a stirring system for a cold model of a lithium electrolysis cell based on CFD simulations and optical experiments[J

RSC Advances, 2015, 5(103): 84503-84516.OLIAII E, DESILETS M, LANTAGNE G. Effect of the design parameters on mass transfer and energy consumption inside a lithium electrolysis cell[J

Journal of Applied Electrochemistry, 2018, 48(6): 725-737.OLIAII E, LITRICO G, DESILETS M, et al. Mass transport and energy consumption inside a lithium electrolysis cell[J

Electrochimica Acta, 2018, 290: 390-403.OLIAII E, DESILETS M, LANTAGNE G. Numerical analysis of the effect of structural and operational parameters on electric and concentration fields of a lithium electrolysis cell[J

Journal of Applied Electrochemistry, 2017, 47(6): 711-726.ARTSDALEN E R, YAFFE I S. Electrical conductance and density of molten salt systems: KCl-LiCl, KCl-NaCl and KCl-KI[J

The Journal of Physical Chemistry, 1955, 59(2):118-127.有色冶金炉设计手册编委会. 有色冶金炉设计手册[M

北京: 冶金工业出版社, 2000: 1035-1035.Editorial Committee of nonferrous Metallurgy Furnace Design Manual. Non-ferrous metallurgical furnace design manual[M

Peking: Metallurgical Industry Press, 2000: 1035-1035.张永健. 镁电解生产工艺学[M

长沙: 中南大学出版社, 2006: 542-543.ZHANG Yong-jian. Electrolytic metallurgy of magnesium[M

Changsha: Central South University Press, 2006: 542-543.SPICKA P, DIAS MM, LOPES JCB. Gas-liquid flow in a 2D column: Comparison between experimental data and CFD modelling[J

Chemical Engineering Science, 2001, 56(21): 6367-6383.VPGT H, BALZER RJ. The bubble coverage of gas-evolving electrodes in stagnant electrolytes[J

Electrochimica Acta, 2005, 50(10): 2073-2079.SORENSEN L S. An introduction to computational fluid dynamics: The finite volume method[J

Pearson Schweiz Ag, 2007, 20(5): 400.VOGT H. The actual current density of gas-evolving electrodes-Notes on the bubble coverage[J

Electrochimica Acta, 2012, 78: 183-187.DIAZ, ME, IRANZO A, et al. Numerical simulation of the gas-liquid flow in a laboratory scale bubble column: Influence of bubble size distribution and non-drag forces[J

Chemical Engineering Journal, 2008, 139(2): 363-379.LEE M, PARK G, PARK C, et al. Improvement of grid independence test for computational fluid dynamics model of building based on grid resolution[J

Advances in Civil Engineering, 2020, 2020: 8827936.