基于多种湍流模型的电解加工温度多场耦合仿真

2022-04-29陈远龙林华陈培譞李回归

华南理工大学学报（自然科学版） 2022年3期

陈远龙林华,2 陈培譞李回归

(1.合肥工业大学机械工程学院,安徽合肥 230009；2.皖西学院机械与车辆工程学院,安徽六安 237000)

电解加工由于其独特的加工优势，被广泛应用于航空、航天、兵器等国防军工及高端装备领域的难加工材料、复杂型面零件的精密加工等，有效支撑和推动了我国航空航天事业的蓬勃发展[1- 4]。电解加工过程是多物理场综合作用的过程，加工间隙内电解液的温度分布不均，对电导率分布有重要影响，使得电场分布不均，最终影响着工件的成型精度。基于单物理场和多物理场耦合的仿真技术在电解加工技术研究中有着越来越广泛的应用，有效提升了加工精度，缩短了研发周期[5- 6]。

多物理场比单物理场的仿真分析更能表达电解加工过程中多场间错综复杂的关系，目前电解加工过程仿真模型已由纯电场或流场的单场仿真向多场耦合仿真的方向发展[7- 9]。Deconinck等[10]基于层流和多离子传输模型研究了电解加工的温度分布和材料去除规律。周小超等[5]基于欧拉双流体气液两相流模型进行了叶片电解加工过程多场耦合研究。陈远龙等[11]基于湍流k-ε模型进行了高频脉冲电化学加工过程多场耦合仿真研究。刘国强等[12]对小孔内扩孔的电解加工过程进行了多场耦合仿真研究。刘强等[13]基于k-ε模型对航空发动机中叶片式扩压器电解加工的流场进行了优化设计。朱彦伟等[14]对采用不同湍流模型对涡轮叶片表面换热计算的影响进行了比较，并指出在不同计算域需采用不同湍流模型和壁面处理函数进行计算。可见，在仿真中使用不同的湍流模型，流场求解结果也不尽相同，尤其是不同流场模型在耦合电场、温度场等物理场时，计算的结果差别更大，很大程度上影响着整个电解加工过程和加工精度分析。

本文针对型面电解加工过程中的间隙温度分布不均匀、难以预测和测量的问题，建立基于两相流场模型并耦合电场、流场和温度场的温度多场耦合仿真模型，采用多种流场模型和壁处理方式进行间隙流速和温度分布的求解，分析其求解精度并与试验结果进行对比，以获得具有较高求解精度的电解加工温度多场耦合模型及其求解方法。

1 几何模型

假设加工进入平衡状态，建立型面电解加工仿真的二维几何模型如图1所示。采用侧流式供液方式，左侧为电解液入口，右侧为电解液出口。工件阳极尺寸为100 mm×90 mm，工具阴极尺寸为100 mm×40 mm，厚度尺寸为12.7 mm，曲率半径为70 mm，加工间隙为0.5 mm，Ω为电解液区域。在加工间隙中设置6个温度观测点1～6，各点间距18 mm，点1和点6距左右端面各5 mm，均高出阴极表面0.05 mm。O1O2为间隙内流速观察线。

图1 电解加工仿真几何模型

2 多场耦合分析

2.1 电场

电场模型如图1所示，Γ1为工件阳极边界，Γ2为工具阴极边界，Γ3=U，Γ4=0，U为加工电压。根据电荷守恒定律和欧姆定律，由电场控制方程，电解液中电位满足[2]:

∇·(κ∇φ)=0

(1)

式中：κ为电解液电导率，φ为电解液电位。在电解液入口和出口满足

∂φ/∂n=0

(2)

式中：n为电位梯度方向上的单位矢量。电解加工中电解液电导率会受到电解液温升和气泡率的影响，其模型可描述为[2]：

κ=κ0[1+ξ(T-T0)](1-β)m

(3)

式中：κ0、T0为电解液初始电导率和温度；ξ为温度系数，取0.02；β为气泡率；m为气泡率影响指数，取2。

2.2 流场

电解加工中，流动的电解液要足以排出间隙中的电解产物与所产生的热量，要能减小电极附近的浓差极化，因此其必须具有一定的速度，呈湍流状态，且流速均匀性要好。为简化计算，假设电解液为理想状态的液体，不含气泡、固体颗粒等，为连续不可压缩黏性流体。由质量守恒定律和动量守恒定律，电解液流体的流动应满足[5]：

(4)

式中：ρ为电解液密度；p为电解液压力；v为电解液流速；μ为电解液动力黏度；μT为湍流黏性系数。

2.2.1 单相流场模型

湍流的计算方法主要分为3类：雷诺时均模拟、尺度解析模拟和直接数值模拟，前两类可以看成是非直接模拟。常用的湍流模型有SA模型、k-ε模型、k-ω模型、SST模型、雷诺应力模型和大涡模拟等[14]。

(1)Spalart-Allmaras模型

湍流Spalart-Allmaras(简称SA模型)模型是单方程模型，只需求解湍流黏性的输运方程。该模型通常用于模拟中等复杂的内流和外流以及压力梯度下的边界层流流动。

(2)标准k-ε模型

湍流标准k-ε模型(简称k-ε模型)是两个方程模型，要解两个变量，即速度和长度尺度，其在现代工程流场计算中的应用最为广泛，但由于其假设流动为完全湍流，故只适用于完全湍流的流动过程模型。模型中湍动能输运方程是通过精确的方程推导得到的，而耗散率方程则是通过物理推理及数学上模拟相似原型方程得到。k-ε模型液相湍动能及耗散率方程为[11]：

(5)

式中：k为湍动能；ε为湍流耗散率；σk、σε、C1ε和C2ε为模型常数；Pk为平均速度梯度产生的湍动能。在标准k-ε模型中，模型的常数取值为

C1ε=1.44，C2ε=1.92，σk=1.0，σε=1.3。

(3)标准k-ω模型

湍流标准k-ω模型(简称k-ω模型)基于Wilcoxk-ω模型，并考虑低雷诺数、可压缩性和剪切流传播对其修改得到[13]。模型预测了自由剪切流传播速率，可以应用于墙壁束缚流动和自由剪切流动。k-ω方程克服了k-ε方程在近壁流动模拟的缺陷，对于近壁流动或存在逆压梯度流动的湍流计算时具有较大优势，同时还考虑到了湍流剪应力的影响修改了湍流黏性公式。在近壁流动模型中k-ω模型精度较高，其适用于近壁低雷诺数区域的模拟，但比k-ε模型收敛困难，计算结果对初始条件较敏感。

(4)SSTk-ω模型

湍流SSTk-ω模型[15](简称SST模型)综合了k-ω模型在近壁处边界层内部的模拟和k-ε模型在外部区域高雷诺数流动模拟的优点，比标准k-ω模型应用更为广泛。

2.2.2 气液两相流场模型

分散气泡流模型的理论依据是欧拉-欧拉模型，主要求解3个方程[16- 17]，其中流体动量方程为

Φlρlg+F

(6)

混合物连续性方程为

(7)

气相输运方程为

(8)

式中：下标为l和g表示液相和气相；Φ为体积分数；mgl表示液相到气相的传质速率。

气泡流速度为

(9)

阴极产生氢气的质量通量为

(10)

式中,NH2为氢气质量通量，M为氢气的摩尔质量，i为局部电流密度，η为电流效率，F为法拉第常数。

2.3 温度场

工具阴极和工件阳极的电阻很小，可忽略其加工中所产生的焦耳热，电解加工过程中产生的热量主要由电解液焦耳热和电解液/电极边界处的电化学反应热组成。电解液温度分布可由流体对流扩散方程描述[10- 11]：

(11)

Qbulk=Qj+Qr

(12)

式中,cp为电解液常压比热容,kt为电解液热导率,T为电解液的温度,Qbulk为电解加工产生的热量,Qj为焦耳热,Qr为电化学反应热。

电解液中产生的焦耳热为

Qj=κ(∇φ)2

(13)

工具阴极和工具阳极与电解液的边界电化学反应热为

Qr=Ukik

(14)

式中：Uk为过电位，ik为双电层的电流密度。

假设工具阴极和阳极的边界Γ4和Γ3与空气自然对流，其余边界均为热绝缘。

2.4 多物理场耦合分析

在电解加工中，阳极的溶解速率直接取决于加工间隙内各点法向电流密度。当施加的外电压一定时，间隙电场分布受到电解液电导率和几何结构的影响，电解液的电导率受温度和气泡的共同影响，而温度分布又受到两相流场和电场等的作用。因此，电解加工过程间隙的温度分布是多物理场耦合的结果[12]。

在多场耦合仿真中，通过流体传热的对流扩散方程中的对流项来耦合流场和温度场的影响；通过热源项(电解液焦耳热和电化学反应热)来耦合电场的影响。在模型中可通过设置阴极/电解液边界的气体通量公式实现电场和气相流场的耦合；通过设置电解液电导率公式来实现电场、温度场和两相流场的耦合。温度多场耦合仿真流程如图2所示。

图2 温度多场耦合仿真流程图

2.5 仿真参数

基于多物理场有限元仿真软件COMSOL Multiphysics进行求解。将几何模型导入仿真软件中，添加多相流场、电场和流体传热等模块，设置边界条件并进行网格划分，单元尺寸选择极细化。工具阴极、工件阳极和电解液的物理特性设置如表1所示。

表1 工具阴极、工件阳极和电解液的物理特性

3 仿真分析

近壁区流速的求解精度影响着电流密度分布的求解精度，流场分布最终影响着温度场分布。近壁区流场的求解通常有两种方法：一种是使用壁函数将壁面上的物理量与湍流核心区内的相应物理量联系起来，能大大提高计算速度和收敛性，但由于简化了边界层内流动的计算，未模拟缓冲区中的流动，因此该方法对边界层流场计算精度不高，但在工程中是很实用，应用广泛；另一种是使用低雷诺数壁处理模型对近壁区黏性底层和过渡层进行求解，精度较高，但计算成本也相对较高。

3.1 不同湍流模型的间隙流动分析

设置入口流速为1 m/s，出口压力为0 MPa(相对于1个大气压)，分别在湍流SA、k-ε、k-ω、SST和低雷诺数k-ε等5个模型中计算间隙内流速分布。在仿真中，求解近壁区流动时，k-ε模型使用壁函数，另4个模型使用自动壁处理和低雷诺数壁处理，求得O1O2的流速分布如图3所示。

从图3(a)和图3(b)中可以看出，5种流场模型采用不同壁处理方式时仿真得到的流速分布有区别，尤其是在近壁区。在近壁区的流速值，k-ε模型的仿真值最高，而SA、SST和低雷诺数k-ε模型的仿真值则很接近。在湍流中心区的流速值，k-ε模型和k-ω模型比另3个模型要低。从图3(a)可知，采用自动壁处理的近壁区流速仿真值差别较

(a)自动壁处理

(b)低雷诺数壁处理

大，其中k-ω模型为4 m/s，SA、SST和低雷诺数k-ε模型为1.2 m/s，k-ε模型为11 m/s。从图3(b)中可知，除k-ε模型以外，其它4个采用低雷诺数壁求解的速度在边界层内均是从0 m/s开始，一直到湍流中心区均有求解，其中SA、SST和低雷诺数k-ε模型计算结果相近。

造成图3(a)和图3(b)中的模型求解结果的差异是因为选择自动壁处理时，仿真软件会根据网格精细程度自动选择壁函数或低雷诺数壁处理，当网格足够精细时自动调用求解精度高的低雷诺数壁处理。k-ε模型是高雷诺数模型，只对充分发展的湍流有效，对流动情况变化大的近壁区黏性底层和缓冲过渡层不进行精确求解，只使用壁函数进行简化计算。

SA、SST、k-ω和低雷诺数k-ε等4个模型是低雷诺数模型，均可对边界黏性底层和缓冲过渡层进行求解，在近壁区的计算精度相对较高[18]。这是因为：SA模型增加了一个额外的无衰减运动学涡流粘度变量，它可求解实体壁之内的整个流场；ω模型是近壁湍流模型，求解的是动能耗散的具体速率，对近壁区边界黏性底层低雷诺数流动模拟较好；SST模型在近壁处采用k-ω模型对近壁区边界层流速较低的黏性底层和缓冲过渡层进行计算；低雷诺数k-ε模型是对标准k-ε模型适应于求解低雷诺数流动而进行修正的模型。在本模型中划分网络时尺寸选择的是极细化，因而使用自动壁和低雷诺数壁处理的仿真结果相差很小。

对于电解加工过程的流体传热研究而言，提高近壁区的流场求解精度无疑将提高温度场分布求解结果的可信度。SA模型常用于空气动力学流动分析，对墙壁束缚流动模拟效果较好，而对电解加工中流动尺度变换比较大的情况并不太适合；k-ω模型非线性程度大，对初值敏感且难以收敛。SST模型结合了k-ε模型和k-ω模型的优点，求解精度高，求解成本不高，其与低雷诺数k-ε模型一样是电解加工仿真中值得研究的模型。

3.2 间隙温度演变分析

选择低雷诺数k-ε湍流模型，采用低雷诺数壁处理，进行温度多场耦合仿真，设置入口流速为1 m/s，出口压力为0 MPa，加工电压为20 V。在0.1 s时，模型的温度场分布如图4所示。间隙内各测量点在0～0.1 s时间段内的温升变化如图5所示，为更清晰展示温度随时间的变化规律，采用对数横坐标。

图4 温度分布

从图4可以看出，热量在电解液和边界处产生，大部分热量由电解液带离加工间隙，在靠近出口处温升最大，部分热量在电极内进行热传导。

由图5可知，在同一时间点上，加工间隙内沿程1～6点的温度均呈上升态势。随着加工时间的增加，各点温度急剧上升后又急速下降，经过一段时间后逐渐趋于稳定。在加工刚开始时，电流密度很大，瞬间温升很快，在流场的作用下，温度沿程累积(图5中点6较为明显)，各点在0.02 s后温升值逐渐稳定。若记进入稳态温升的±(2%～5%)范围内为进入准稳态，在当前加工参数下，点6在0.05 s左右进入准稳态，其余各点的时间更短。因此，后文可取0.05 s时的物理场分布进行分析。

图5 点1～6的温升变化

3.3 不同湍流模型下间隙流动传热分析

基于不同湍流模型(k-ε、SST和低雷诺数k-ε模型)，以及耦合气泡的低雷诺数k-ε模型对间隙温度进行多场耦合仿真。设置入口流速为1 m/s，出口压力为0 MPa，加工电压为20 V。在0.05 s时，点1～6的温升值如图6所示。

图6 不同流场模型的温度分布

由图6可以看出，不同流场模型对间隙温度分布求解结果不尽相同，沿程各点温度值呈上升趋势。SST和低雷诺数k-ε模型求解值近似相等，且高于k-ε模型求解值，这是因为前两个模型的计算流速值近似相等且低于后者(如图3(b)所示)，电解液流速大能在相同时间内带走更多的热量，因而温升相对较小。根据式(3)可知，电解液中气泡的存在会造成电解液电导率、电流密度的下降，因此耦合气泡的低雷诺数k-ε模型求解的沿程各点温升会有所减小。

3.4 不同加工电压和入口流速时的温度分布

选择耦合气泡的低雷诺数k-ε模型进行电解加工温度多场耦合仿真。图7为加工电压分别为16、20和24 V时的沿程温度分布。图8为入口流速分别为0.6、1.0和1.4 m/s时的沿程温度分布。

图7 不同加工电压时沿程温度分布

图8 不同入口流速时沿程温度分布

从图7和图8中可以看出，在电解加工中，加工电压越大，单位时间内产生的热量越多，在电解液的带动下，温度沿程累积，在出口处最大温升分别为16、26和40 K左右。增加入口电解液流速或压力，可使得加工间隙内电解液的流速增加，带走热量的能力增强，因而沿程温升变小。

4 试验验证

为更进一步验证前文温度模型求解的精度以及不同模型求解的差异性，与相关试验结果进行对比分析。试验数据来源见文献[19]。仿真中，调整几何模型与文献所述的试验装置一致，流场模型选择低雷诺数k-ε模型(低雷诺数壁处理)、k-ε模型(壁函数壁处理)，以及耦合气泡的低雷诺数k-ε模型等3种模型进行温度多场耦合仿真。仿真参数的设置与试验值一致：间隙入口质量流速为11.3×10-5m3/s，出口背压为0.17 MPa，加工电流密度为16 A/cm2，电解液质量分数为10%的NaCl溶液，初始温度为293 K，电导率为13 S/m，阴极进给速度为0.35 mm/min，加工时间为3 min。1～6点仿真温升数据与试验数据对比如图9所示。

从图9可以看出，仿真值与试验值的变化趋势一致，其中低雷诺数k-ε模型仿真值比k-ε模型低，这是因为前者在测量点处的流速相对较高。还可以看出，耦合气泡的低雷诺数k-ε模型的仿真值整体上与试验值更为接近，在出口处的温度近似相等，但仿真和试验值的差值在流程的前半程比后半程大，其主要原因是试验装置中测温传感器的安装是伸出阴极表面的，在间隙入口处氢气泡在测温头附近更容易累积，使得测温头附近气液混合流体的电阻率更大，以致测得的温度会比实际值高，试验误差较大，而在后半程由于气泡的充分累积，试验误差变小。