工程科学学报，第38卷，第8期：1091-1097,2016年8月 Chinese Journal of Engineering,Vol.38,No.8:1091-1097,August 2016 D0l:10.13374/j.issn2095-9389.2016.08.007:http://journals..ustb.edu.cn 基于数值模拟的气液相界面积计算方法 宋彦坡12)四，刘志高”，陶焰明”，彭小奇12”，陈卓” 1)中南大学能源科学与工程学院，长沙410083：2)湖南第一师范学院信息科学与工程学院，长沙410205 ☒通信作者，E-mail:songyanpo@csu.edu.cm 摘要相界面积对气液两相流中传热、传质、物理化学反应等动力学过程影响重大.为获取这一参数，提出一种根据两相 流数值模拟结果计算相界面积的方法.此方法借鉴分段线性重构界面的思想，在各网格单元内以平面近似真实相界曲面，根 据目标流体的体积分数及其梯度向量将网格内相界面形貌归为五类，进而采用不同的方法分别计算各类相界面的面积.在 铜转炉熔池内两相流数值模拟结果分析中的应用效果表明：该方法能有效提取两相流体系中任意区域的相界面积，从而为体 系动力学特征的定量分析提供依据.利用相界面积数据，进一步计算了氧气利用率并识别出熔池内“高效反应区”，计算和识 别结果与工程实际吻合，证实了该方法的准确性 关键词两相流：界面：面积：数值模拟：计算流体力学 分类号TF801.2:0359.1 Computation method of gas-liquid interfacial area based on numerical simulation results SONG Yan-po),LIU Zhi-gao),TAO Yan-ming,PENG Xiao-i),CHEN Zhuo) 1)School of Energy Science and Engineering,Central South University,Changsha 410083,China 2)School of Information Science and Engineering,Hunan First Normal University,Changsha 410205,China Corresponding author,E-mail:songyanpo@csu.edu.cn ABSTRACT The gas-liquid interfacial area has significant influence on kinetic processes,such as heat transfer,mass transfer and physiochemical reactions in two-phase flow.A new method for getting this parameter is introduced to compute the interfacial area with the numerical simulation results of two-phase flow.Referring to the idea of a piecewise linear interface calculation (PLIC)method,it uses a plane in each cell to approximate the real curved interface between two phases,these planes are then categorized into five types according to the volume fractions of the target fluid and their gradient vectors in each cell,and the interfacial areas are respectively calculated by different equations for different plane types.This method is applied in analyzing the numerical simulation results of a copper converter.It is shown that this method can extract the interfacial area of any spatial region effectively in a two-phase flow system and can be used in analyzing the features of kinetic processes quantitatively in a two-phase dispersion system.Moreover,by using the computed interfacial area,the oxygen utilization ratio is calculated and the "highly efficient reaction zone"in the bath is recognized. The results agree with practical data and experience,indicating the accuracy of the proposed method in some extent. KEY WORDS two-phase flow;interfaces;area:numerical simulation:computational fluid dynamics 气液两相流体系是治金及化工领域中的常见体 行，因此较为准确地获取相界面积对研究这类体系中 系，相间的传热、传质以及化学反应主要通过界面进 两相间的动力学特征具有重要意义.目前，气液相界 收稿日期：201509-11 基金项目：国家自然科学基金资助项目(61105080,61273159,61134006)：湖南省科技重大专项资助项目(2013F1009):湖南省自然科学基金 资助项目(114057)：安徽省科技厅重大科技专项(1301021018)
工程科学学报，第 38 卷，第 8 期: 1091--1097，2016 年 8 月 Chinese Journal of Engineering，Vol． 38，No． 8: 1091--1097，August 2016 DOI: 10． 13374 /j． issn2095--9389． 2016． 08． 007; http: / /journals． ustb． edu． cn 基于数值模拟的气液相界面积计算方法 宋彦坡1，2) ，刘志高1) ，陶焰明1) ，彭小奇1，2) ，陈 卓1) 1) 中南大学能源科学与工程学院，长沙 410083; 2) 湖南第一师范学院信息科学与工程学院，长沙 410205 通信作者，E-mail: songyanpo@ csu． edu． cn 摘 要 相界面积对气液两相流中传热、传质、物理化学反应等动力学过程影响重大． 为获取这一参数，提出一种根据两相 流数值模拟结果计算相界面积的方法． 此方法借鉴分段线性重构界面的思想，在各网格单元内以平面近似真实相界曲面，根 据目标流体的体积分数及其梯度向量将网格内相界面形貌归为五类，进而采用不同的方法分别计算各类相界面的面积． 在 铜转炉熔池内两相流数值模拟结果分析中的应用效果表明: 该方法能有效提取两相流体系中任意区域的相界面积，从而为体 系动力学特征的定量分析提供依据． 利用相界面积数据，进一步计算了氧气利用率并识别出熔池内“高效反应区”，计算和识 别结果与工程实际吻合，证实了该方法的准确性． 关键词 两相流; 界面; 面积; 数值模拟; 计算流体力学 分类号 TF801. 2; O359. 1 Computation method of gas--liquid interfacial area based on numerical simulation results SONG Yan-po1，2) ，LIU Zhi-gao1) ，TAO Yan-ming1) ，PENG Xiao-qi1，2) ，CHEN Zhuo1) 1) School of Energy Science and Engineering，Central South University，Changsha 410083，China 2) School of Information Science and Engineering，Hunan First Normal University，Changsha 410205，China Corresponding author，E-mail: songyanpo@ csu． edu． cn ABSTＲACT The gas--liquid interfacial area has significant influence on kinetic processes，such as heat transfer，mass transfer and physiochemical reactions in two-phase flow． A new method for getting this parameter is introduced to compute the interfacial area with the numerical simulation results of two-phase flow． Ｒeferring to the idea of a piecewise linear interface calculation ( PLIC) method，it uses a plane in each cell to approximate the real curved interface between two phases，these planes are then categorized into five types according to the volume fractions of the target fluid and their gradient vectors in each cell，and the interfacial areas are respectively calculated by different equations for different plane types． This method is applied in analyzing the numerical simulation results of a copper converter． It is shown that this method can extract the interfacial area of any spatial region effectively in a two-phase flow system and can be used in analyzing the features of kinetic processes quantitatively in a two-phase dispersion system． Moreover，by using the computed interfacial area，the oxygen utilization ratio is calculated and the“highly efficient reaction zone”in the bath is recognized． The results agree with practical data and experience，indicating the accuracy of the proposed method in some extent． KEY WOＲDS two-phase flow; interfaces; area; numerical simulation; computational fluid dynamics 收稿日期: 2015--09--11 基金项目: 国家自然科学基金资助项目( 61105080，61273159，61134006) ; 湖南省科技重大专项资助项目( 2013FJ1009) ; 湖南省自然科学基金 资助项目( 11JJ4057) ; 安徽省科技厅重大科技专项( 1301021018) 气液两相流体系是冶金及化工领域中的常见体 系，相间的传热、传质以及化学反应主要通过界面进 行，因此较为准确地获取相界面积对研究这类体系中 两相间的动力学特征具有重要意义． 目前，气液相界
·1092· 工程科学学报，第38卷，第8期 面积的实验测定方法主要有光学仪器直接测量法和化 中是以直线代替曲线)重构界面.在二维网格内，PLIC 学间接测量法两类四.光学仪器直接测量法-网采用 的界面重构原理如图1所示 高速摄影仪等电子光学仪器直接测量相界面的分布， 进一步应用计算机图像处理技术获得相界面积：但该 方法仅适用于连续相可透光的两相流体系，应用面较 窄，且难以获得清晰的流场图片和精确的定量结果 化学法通过向气相与液相中分别加入一定量的溶 质与溶剂物质，如C02与K,C03/KHC0,然后测量一 段时间内溶质与溶剂的浓度，计算出溶质的吸收速率， 最后采用Danckwerts标绘方法可获得相界面积：然而 化学法测定的值与选用的溶质一溶剂物系有关，且实 验测试手段的准确性较差. B H 另一方面，利用现有的数值模拟方法，例如流体体 图1PLIC界面重构示意图 积函数(volume of fluid,.VOF)方法。，可以较为准确 Fig.1 Interface reconstruction scheme of the PLIC method 地模拟气液两相流及其界面特征.VOF方法对流场中 图1中，虚线为网格内的真实相界面，直线EF为 的每个网格定义了流体体积函数（即目标流体与网格 重构的相界面.重构界面的法向量n=(n,n,)与该网 的体积比)，利用该函数值描述多相流各相及界面的 格的体积分数梯度方向一致 分布特征.分段线性化（Piecewise Linear Interface 由重构界面的法向量n=(n,n,),界面的直线方 Calculation,PLIC)则是VOF方法常用的界面重构方 程可表示为 案，采用“以直代曲”的方法近似重构相界面，因具有 nx +n,y=a. (2) 精度高且计算量小的特点而被广泛应用-0.运用 式中，α为待定常数，可由C值（即图1中阴影部分与 PLIC-VOF方法，张勤等四模拟浮法玻璃生产过程中 网格的面积之比)计算得到，从而唯一确定该网格内 玻璃液与气体间的相界面：刘国林等四模拟宽板坯连 的相界面形状及位置. 铸结晶器内的液面波动；Maurya等模拟斜平板上气 PLIC对三维网格内的界面重构方法与其对二维 液降膜流动过程中相界面的演化历程：杨俊英等模 网格的处理方法类似，用平面代替曲面重构相界面，重 拟树脂传递模塑工艺充模过程中流动前沿相界面的位 构界面的单位法向量由下式计算： 置和形状，这些研究表明现有数值模拟方法可以较为 VC 准确地追踪和重构两相流体系的气液相界面.然而， n=Cl (3) 在数值模拟的基础上进一步计算获取相界面积的方法 1.2气液相界面积的计算方法 尚未见文献报道. 在0<C<1的长方体网格中，相界面的法向量可 鉴于相界面积对多相流研究的重要性以及现有获 表示为 取手段的不足，本文借鉴PLIC方法重构界面的思想， no=(n,n'n). (4) 研究提出了一种根据两相流数值模拟结果计算两相界 为了简化三维网格中重构相界面的形貌分类，依 面面积的方法，以LUENT软件为平台开发相应的应 据重构相界面的法向量及网格单元的边长（分别表示 用程序，并将其用于分析铜锍吹炼PS转炉内气液分布 为6，、8，和8)对坐标系进行翻转与对称操作，并保证 特征及部分动力学特征 新坐标系中重构相界面的“新”法向量（表示为八1= 1气液相界面积的计算方法 (abs(n,),abs(n,),abs(n),其中a、b和c互不相同 且均属于集合{x,y,z})满足 1.1PLIC-VOP的界面重构基本原理 abs(n6)≤abs(nδ)≤abs(nδ）. (5) VOF方法定义每个网格中目标流体与网格的体 在sna-as（ab-sln c(简化表示为 积比为目标流体的体积分数C,满足如下控制方程： n。 n。 aCu+v(uC）=0. (1) a'-b'-c)的新坐标系中，网格中的相界面形貌可以 分为五类，直观形状如图2所示，具体判别条件见 式中，t为时间，s;u为流体速度，m·s.显然，仅0< 表1. C。<1的网格内存在相界面.为了确定界面在网格内 表1各式中，为网格体积，C为体积分数较小流体 的形状及位置，PLIC方法认为界面的法方向应该与C: 的体积分数（即C=min(C,C,),下标g与1分别表示气 值的梯度方向一致，并以平面代替曲面（在二维网格 相和液相)，n,=abs(n.),n2=abs(n,),n3=abs(n)
工程科学学报，第 38 卷，第 8 期 面积的实验测定方法主要有光学仪器直接测量法和化 学间接测量法两类［1］． 光学仪器直接测量法［1--2］采用 高速摄影仪等电子光学仪器直接测量相界面的分布， 进一步应用计算机图像处理技术获得相界面积; 但该 方法仅适用于连续相可透光的两相流体系，应用面较 窄，且难以获得清晰的流场图片和精确的定量结果． 化学法［3--4］通过向气相与液相中分别加入一定量的溶 质与溶剂物质，如 CO2 与 K2 CO3 /KHCO3，然后测量一 段时间内溶质与溶剂的浓度，计算出溶质的吸收速率， 最后采用 Danckwerts 标绘方法［5］获得相界面积; 然而 化学法测定的值与选用的溶质--溶剂物系有关，且实 验测试手段的准确性较差． 另一方面，利用现有的数值模拟方法，例如流体体 积函数( volume of fluid，VOF) 方法［5--8］，可以较为准确 地模拟气液两相流及其界面特征． VOF 方法对流场中 的每个网格定义了流体体积函数( 即目标流体与网格 的体积比) ，利用该函数值描述多相流各相及界面的 分布 特 征． 分 段 线 性 化 ( Piecewise Linear Interface Calculation，PLIC) 则是 VOF 方法常用的界面重构方 案，采用“以直代曲”的方法近似重构相界面，因具有 精度高且计算量小的特点而被广泛应用［9--10］． 运用 PLIC--VOF 方法，张勤等［11］模拟浮法玻璃生产过程中 玻璃液与气体间的相界面; 刘国林等［12］模拟宽板坯连 铸结晶器内的液面波动; Maurya 等［13］模拟斜平板上气 液降膜流动过程中相界面的演化历程; 杨俊英等［14］模 拟树脂传递模塑工艺充模过程中流动前沿相界面的位 置和形状，这些研究表明现有数值模拟方法可以较为 准确地追踪和重构两相流体系的气液相界面． 然而， 在数值模拟的基础上进一步计算获取相界面积的方法 尚未见文献报道． 鉴于相界面积对多相流研究的重要性以及现有获 取手段的不足，本文借鉴 PLIC 方法重构界面的思想， 研究提出了一种根据两相流数值模拟结果计算两相界 面面积的方法，以 FLUENT 软件为平台开发相应的应 用程序，并将其用于分析铜锍吹炼 PS 转炉内气液分布 特征及部分动力学特征． 1 气液相界面积的计算方法 1. 1 PLIC--VOF 的界面重构基本原理 VOF 方法定义每个网格中目标流体与网格的体 积比为目标流体的体积分数 Ctf，满足如下控制方程: Ctf t + Δ ·( uCtf ) = 0． ( 1) 式中，t 为时间，s; u 为流体速度，m·s － 1 ． 显然，仅 0 ＜ Ctf ＜ 1 的网格内存在相界面． 为了确定界面在网格内 的形状及位置，PLIC 方法认为界面的法方向应该与 Ctf 值的梯度方向一致，并以平面代替曲面( 在二维网格 中是以直线代替曲线) 重构界面． 在二维网格内，PLIC 的界面重构原理如图 1 所示． 图 1 PLIC 界面重构示意图 Fig． 1 Interface reconstruction scheme of the PLIC method 图 1 中，虚线为网格内的真实相界面，直线 EF 为 重构的相界面． 重构界面的法向量 n = ( nx，ny ) 与该网 格的体积分数梯度方向一致． 由重构界面的法向量 n = ( nx，ny ) ，界面的直线方 程可表示为 nx x + ny y = α． ( 2) 式中，α 为待定常数，可由 Ctf值( 即图 1 中阴影部分与 网格的面积之比) 计算得到，从而唯一确定该网格内 的相界面形状及位置． PLIC 对三维网格内的界面重构方法与其对二维 网格的处理方法类似，用平面代替曲面重构相界面，重 构界面的单位法向量由下式计算［6］: n = Δ Ctf | Δ Ctf | ． ( 3) 1. 2 气液相界面积的计算方法 在 0 ＜ C ＜ 1 的长方体网格中，相界面的法向量可 表示为 n0 = ( nx，ny，nz ) ． ( 4) 为了简化三维网格中重构相界面的形貌分类，依 据重构相界面的法向量及网格单元的边长( 分别表示 为 δx、δy 和 δz ) 对坐标系进行翻转与对称操作，并保证 新坐标系中重构相界面的“新”法向量( 表示为 n1 = ( abs( na ) ，abs( nb ) ，abs( nc ) ) ，其中 a、b 和 c 互不相同 且均属于集合{ x，y，z} ) 满足 abs( na δa ) ≤abs( nb δb ) ≤abs( ncδc) ． ( 5) 在abs( na ) na a － abs( nb ) nb b － abs( nc) nc c( 简化表示为 a' － b' － c') 的新坐标系中，网格中的相界面形貌可以 分为 五 类，直 观 形 状 如 图 2 所 示，具 体 判 别 条 件 见 表 1． 表 1 各式中，v 为网格体积，C 为体积分数较小流体 的体积分数( 即 C = min( Cg，Cl ) ，下标 g 与 l 分别表示气 相和液相) ，n1 = abs( na ) ，n2 = abs( nb ) ，n3 = abs( nc) ． · 2901 ·
宋彦坡等：基于数值模拟的气液相界面积计算方法 ·1093· 类型 类型2 类型3 类型4 类型5 图2五类基本界面网格 Fig.2 Five generic types of interface configuration 表1五类基本界面网格的判别条件 Table 1 Criteria for the five classified interfaces in cells 类型 判决条件 c (n6)3 C()(n6) 6nn2ng 6nnan3t n6.+d,≥n,d且ad)产-d-n6) a且a8-Aa8-a8-Csa6+-aS-a-a成+am 5 4成+4<8且dtn6-a6户-e5)<Cs色i-a6-ad-a80+5-n6.-nm6 6ninan3u 在新坐标系a'-b-c中，重构相界面所在的平 根据“网格被相平面所切块（图2中阴影面以下部 面方程可表示为 分)的体积等于体积分数C与网格体积：之积”这一条 na'+nb'+ngc'=d. (6) 件，可构造五类相平面方程中待定常数d的计算式，而 式中，d为待定常数. 由表1中各不等式可确定d的求解范围，详见表2. 表2五类相平面方程常数d的计算式 Table2 Equations for computing the constantd 类型 关于待定常数d的方程 求解范围 d=6nn2n3C d≤n16。 2 P-(d-n16.)3=6n1n2m3Cm n16。<d≤nd P-(d-n16)3-(d-n2d)3=6n1n2nCm n2d<d≤n3d且d≤md。+nm2i d-(d-m18)3-(d-n6)3-(d-n36)3=6n12n3Cm n36<d≤n16.+n28 dP-(d-n18)3-(d-n28)3+(d-n1d。-nd)3=6n1n2nCcm n16。+n8<d≤n36
宋彦坡等: 基于数值模拟的气液相界面积计算方法 图 2 五类基本界面网格 Fig． 2 Five generic types of interface configuration 表 1 五类基本界面网格的判别条件 Table 1 Criteria for the five classified interfaces in cells 类型 判决条件 1 C≤ ( n1 δa ) 3 6n1 n2 n3 v 2 ( n1 δa ) 3 6n1 n2 n3 v ＜ C≤ ( n2 δb ) 3 － ( n2 δb － n1 δa ) 3 6n1 n2 n3 v ( 3 n1 δa + n2 δb≥n3 δc 且( n2 δb ) 3 － ( n2 δb － n1 δa ) 3 6n1 n2 n3 v ＜ C≤ ( n3 δc) 3 － ( n3 δc － n1 δa ) 3 － ( n3 δc － n2 δb ) 3 6n1 n2 n3 v ) ( or n1 δa + n2 δb ＜ n3 δc 且( n2 δb ) 3 － ( n2 δb － n1 δa ) 3 6n1 n2 n3 v ＜ C≤ ( n1 δa + n2 δb ) 3 － ( n1 δa ) 3 － ( n2 δb ) 3 6n1 n2 n3 v ) 4 n1 δa + n2 δb ＞ n3 δc 且( n3 δc) 3 － ( n3 δc － n1 δa ) 3 － ( n3 δc － n2 δb ) 3 6n1 n2 n3 v ＜ C≤ ( n1 δa + n2 δb ) 3 － ( n1 δa ) 3 － ( n2 δb ) 3 － ( n1 δa + n2 δb － n3 δc) 3 6n1 n2 n3 v 5 n1 δa + n2 δb ＜ n3 δc 且( n1 δa + n2 δb ) 3 － ( n1 δa ) 3 － ( n2 δb ) 3 6n1 n2 n3 v ＜ C≤ ( n3 δc) 3 － ( n3 δc － n1 δa ) 3 － ( n3 δc － n2 δb ) 3 + ( n3 δc － n1 δa － n2 δb ) 3 6n1 n2 n3 v 在新坐标系 a' － b' － c'中，重构相界面所在的平 面方程可表示为 n1 a' + n2 b' + n3 c' = d． ( 6) 式中，d 为待定常数． 根据“网格被相平面所切块( 图 2 中阴影面以下部 分) 的体积等于体积分数 C 与网格体积 v 之积”这一条 件，可构造五类相平面方程中待定常数 d 的计算式，而 由表 1 中各不等式可确定 d 的求解范围，详见表 2． 表 2 五类相平面方程常数 d 的计算式 Table 2 Equations for computing the constant d 类型 关于待定常数 d 的方程 求解范围 1 d3 = 6n1 n2 n3C d≤n1 δa 2 d3 － ( d － n1 δa ) 3 = 6n1 n2 n3Cv n1 δa ＜ d≤n2 δb 3 d3 － ( d － n1 δa ) 3 － ( d － n2 δb ) 3 = 6n1 n2 n3Cv n2 δb ＜ d≤n3 δc 且 d≤n1 δa + n2 δb 4 d3 － ( d － n1 δa ) 3 － ( d － n2 δb ) 3 － ( d － n3 δc) 3 = 6n1 n2 n3Cv n3 δc ＜ d≤n1 δa + n2 δb 5 d3 － ( d － n1 δa ) 3 － ( d － n2 δb ) 3 + ( d － n1 δa － n2 δb ) 3 = 6n1 n2 n3Cv n1 δa + n2 δb ＜ d≤n3 δc · 3901 ·
·1094· 工程科学学报，第38卷，第8期 可以证明，表2中计算式在相应的求解范围内有 开始 唯一解，且可由相关数值方法求取.在此基础上，由 式(6)确定的相平面方程可计算得到网格内相界面 i-1,M=0 (图2中阴影面)的面积，见表3. 表3五类相界面积的计算式 读取第号网格信总： 、òd、C和n Table 3 Calculation formulas for the five interfacial area 类型 相界面积计算公式 山表1的判别条件归类界面形貌 厅+店+店. 2n12y 山表2骑定界平面方科 3 厅+防+心-(d-nd,)鬥 e+】 2n12y 出表3计算网格相界面面积4， 3 +n店+[f-(d-n8,)2-(d-mn)门 2nin2n3 计算总相界面积A,=A,+A 居++[d-（d-n,6)2-(d-d,)2- 2n12n3 -分析区域内的网格总数m? (d-nd)2] 是Y 屏+时+md-(d-ò.)2-(d-m)2+ 返回分析区域内总相界面积4，并结束 2nin2n3 (d-n16。-m2d)2] 图3计算机实现方法流程图 Fig.3 Flow chart of the computing method 运用以上方法，理论上可以较为准确地计算直六 泡的生成过程，计算区域应采用足够细的网格进行离 面体网格单元内的界面面积：但对于斜六面体等其他 散，由于受到网格尺寸以及计算机性能的限制，难以对 类型的网格单元，上述方法并不适用.对于采用结构 整个转炉进行数值模拟，因此本文仅对单个喷嘴区域 化网格划分的仿真，正六面体网格往往占总网格数 进行研究，其几何结构如图4所示 的绝对多数，斜六面体等其他类型网格单元占网格 总数的比例较低并且其中多数扭曲度较小.因此，在 采用结构化网格的仿真中，将斜六面体作为同体积、 同长宽高比的直六面体网格处理不会产生很大的计 D 算误差. 1.3计算机实现方法 本文提出的两相界面积提取方法可以利用LU- ENT软件提供的用户编程接口一用户自定义函数 (user-defined function,UDF)编程实现，也可以将数值 图4PS转炉单喷嘴模型的几何结构 仿真结果中网格单元的几何尺寸、中心坐标、体积、目 Fig.4 Geometry of the single tuyere model of a PS converter 标流体体积分数及其梯度向量等数据导出后利用其他 编程软件编程实现.具体的计算机实现过程如图3 对计算区域采用分块划分六面体结构化网格，并 所示. 对喷嘴附近网格进行加密处理.通过网格独立性分 析，选用1cm×1cm×lcm尺寸网格划分计算区域.网 2算法验证及应用实例分析 格总数约130万，熔池内大部分区域为正六面体网格 2.1铜锍PS转炉造铜期吹炼过程数值模拟 (约占网格总数的80%)，壁面附近有少量斜六面体网 为了验证上节提出的相界面积计算方法并说明其 格（约占网格总数的20%，其中中间角失真大于0.1 用途及意义，本文将之用于分析国内某治炼厂一台PS 仅占约4%)，因此可以采用本文提出的方法计算其界 铜锍吹炼转炉的数值模拟结果.该转炉有59个0°喷 面积 嘴，在转炉一侧沿水平方向等距布置，炉膛截面直径 固体壁面看作无滑移边界条件，壁面附近用标准 D=3.1m,喷嘴间距L=0.152m,喷嘴内半径r。= 壁面函数进行处理：两端圆面为对称边界条件；喷嘴入 22.9mm在造铜期初期，炉内白铜锍高度H=1.3m, 口为质量入口边界条件，鼓入空气由于马赫数Ma> 喷嘴浸入熔体深度S=0.735m.为准确再现熔池内气 0.3,将其看作可压缩气体，质量流量为0.23kgs,湍
工程科学学报，第 38 卷，第 8 期 可以证明，表 2 中计算式在相应的求解范围内有 唯一解，且可由相关数值方法求取． 在此基础上，由 式( 6) 确定的相平面方程可计算得到网格内相界面 ( 图 2 中阴影面) 的面积，见表 3． 表 3 五类相界面积的计算式 Table 3 Calculation formulas for the five interfacial area 类型 相界面积计算公式 1 n2 1 + n2 槡 2 + n2 3 2n1 n2 n3 ·d2 2 n2 1 + n2 槡 2 + n2 3 2n1 n2 n3 ·［d2 － ( d － n1 δa ) 2 ］ 3 n2 1 + n2 槡 2 + n2 3 2n1 n2 n3 ·［d2 － ( d － n1 δa ) 2 － ( d － n2 δb ) 2 ］ 4 n2 1 + n2 槡 2 + n2 3 2n1 n2 n3 ·［d2 － ( d － n1 δa ) 2 － ( d － n2 δb ) 2 － ( d － n3 δc) 2 ］ 5 n2 1 + n2 槡 2 + n2 3 2n1 n2 n3 ·［d2 － ( d － n1 δa ) 2 － ( d － n2 δb ) 2 + ( d － n1 δa － n2 δb ) 2 ］ 运用以上方法，理论上可以较为准确地计算直六 面体网格单元内的界面面积; 但对于斜六面体等其他 类型的网格单元，上述方法并不适用． 对于采用结构 化网格划分的仿真，正六面体网格往往占总网格数 的绝对多数，斜六面体等其他类型网格单元占网格 总数的比例较低并且其中多数扭曲度较小． 因此，在 采用结构化网格的仿真中，将斜六面体作为同体积、 同长宽高比的直六面体网格处理不会产生很大的计 算误差． 1. 3 计算机实现方法 本文提出的两相界面积提取方法可以利用 FLUENT 软件提供的用户编程接口———用户自定义函数 ( user-defined function，UDF) 编程实现，也可以将数值 仿真结果中网格单元的几何尺寸、中心坐标、体积、目 标流体体积分数及其梯度向量等数据导出后利用其他 编程软件编程实现． 具体的计算机实现过程如图 3 所示． 2 算法验证及应用实例分析 2. 1 铜锍 PS 转炉造铜期吹炼过程数值模拟 为了验证上节提出的相界面积计算方法并说明其 用途及意义，本文将之用于分析国内某冶炼厂一台 PS 铜锍吹炼转炉的数值模拟结果． 该转炉有 59 个 0°喷 嘴，在转炉一侧沿水平方向等距布置，炉膛截面直径 D = 3. 1 m，喷 嘴 间 距 L = 0. 152 m，喷 嘴 内 半 径 r0 = 22. 9 mm． 在造铜期初期，炉内白铜锍高度 H = 1. 3 m， 喷嘴浸入熔体深度 S = 0. 735 m． 为准确再现熔池内气 图 3 计算机实现方法流程图 Fig． 3 Flow chart of the computing method 泡的生成过程，计算区域应采用足够细的网格进行离 散，由于受到网格尺寸以及计算机性能的限制，难以对 整个转炉进行数值模拟，因此本文仅对单个喷嘴区域 进行研究，其几何结构如图 4 所示． 图 4 PS 转炉单喷嘴模型的几何结构 Fig． 4 Geometry of the single tuyere model of a PS converter 对计算区域采用分块划分六面体结构化网格，并 对喷嘴附近网格进行加密处理． 通过网格独立性分 析，选用 1 cm × 1 cm × 1 cm 尺寸网格划分计算区域． 网 格总数约 130 万，熔池内大部分区域为正六面体网格 ( 约占网格总数的 80% ) ，壁面附近有少量斜六面体网 格( 约占网格总数的 20% ，其中中间角失真大于 0. 1 仅占约 4% ) ，因此可以采用本文提出的方法计算其界 面积． 固体壁面看作无滑移边界条件，壁面附近用标准 壁面函数进行处理; 两端圆面为对称边界条件; 喷嘴入 口为质量入口边界条件，鼓入空气由于马赫数 Ma ＞ 0. 3，将其看作可压缩气体，质量流量为 0. 23 kg·s － 1，湍 · 4901 ·
宋彦坡等：基于数值模拟的气液相界面积计算方法 ·1095· 流强度为3.2%，入口温度为300K:烟气出口采用压 多相流模型选用VOF模型：湍流模型选用RNGk- 力出口边界条件，出口压力为-10Pa ε模型：动量方程的离散用二阶迎风格式：压力的离散 铜锍吹炼过程分为造渣期和造铜期两个阶段，造 用PRESTO格式：压力速度耦合方式选用PISO算法， 渣期结束后，熔池内的熔体由铜锍（主要为C,S与 这种耦合方法适用于非定常流动的计算：为确保计算 FeS的混合物)转变为白铜锍（主要为Cu,S),白铜锍 的精度以及收敛性，基本时间步长定为0.00001s,计 在造铜期吹炼过程中被继续氧化生成粗铜.本文仅针 算过程中根据收敛的实际情况适当减小或增大时间步 对造铜期吹炼初期进行数值模拟，故假设熔池中只存 长；为提高计算收敛性，对部分控制方程采用欠松弛法 在白铜锍和富氧空气两种流体.白铜锍的物性参数如 求解.计算稳定后部分时刻过喷嘴中心横截面的气相 表4所列. 体积分数分布如图5所示 表4白铜统的物性参数 Table 4 Property parameters of white metal 密度/(kgm3) 表面张力/(Nm1) 温度K 黏度/(Pa's) 比热容/(kg1K)热导率/(Wm1K) 5400 0.31 1473 0.004 540 8.9 (b) 828 图5不同时刻中心截面的气相体积分数分布.(a)8.6s:(b)8.8s:(c)9s Fig.5 Volume fraction distribution of gas in the center plane at different moments:(a)8.6s:(b)8.8s:(c)9s 2.2界面积计算方法在模拟结果分析中的应用 利用前述方法，根据数值模拟结果，计算流体运动 8.0 相对稳定的8~9s时间内多个时刻点熔池中总相界面 积，如图6所示.由图可知，由于气液相混合流体内部 的扰动及自由液面的波动，熔池内总相界面积随着时 间的变化，一般在7.53~8m范围内波动. 为进一步了解熔池中相界面的分布情况，计算 7.6 8~9s内不同熔池高度以下相界面积的时间平均值，如 15 图7所示.由图可知：在喷嘴以下区域（熔池高度 8.0 8.2 8.48.68.89.0 <0.565m)相界面积接近零，喷嘴以下存在的富氧空 时间/ 气极少：相界面主要分布在喷嘴以上、自由液面以下区 图6不同时刻熔池内总相界面积 域：自由液面以上仍存在一定的相界面积，主要是由于 Fig.6 Total interfacial area in the bath at different times 熔体的喷溅引起 由数值模拟结果图5和图7可知，气相在喷嘴以 在铜锍吹炼的造铜期，发生的主要化学反应为硫 下区域的分布极少，相应地，喷嘴以下区域发生的氧化 化亚铜的氧化反应，在近1200℃的温度下，该反应的 反应可忽略不计.为了进一步研究喷嘴以上区域的化 速率受气相中氧的传质速率所控制，即氧化反应速率 学反应情况，对喷嘴以上(0.565~1.300m)熔池内划 等于氧的传质速率： 分为48个水平层.利用仿真数据以及本文方法可获 1o.=ko,"A-co. (7) 取各层的气液相界面积、气相的平均停留时间等，根据 式中：i。为气相中氧传递速率，mol.s;k。,为氧传质系 下式可计算各层中的氧气消耗量： 数，ms1,根据文献6]，此处取0.6ms1:A为气液 lo,ko,"A"cojt (8) 相界面积，m2;co,为气相中氧浓度，molm3 式中，下标j表示第j层，是气体在第j层的平均停留
宋彦坡等: 基于数值模拟的气液相界面积计算方法 流强度为 3. 2% ，入口温度为 300 K; 烟气出口采用压 力出口边界条件，出口压力为 － 10 Pa． 铜锍吹炼过程分为造渣期和造铜期两个阶段，造 渣期结束后，熔池内的熔体由铜锍( 主要为 Cu2 S 与 FeS 的混合物) 转变为白铜锍( 主要为 Cu2 S) ，白铜锍 在造铜期吹炼过程中被继续氧化生成粗铜． 本文仅针 对造铜期吹炼初期进行数值模拟，故假设熔池中只存 在白铜锍和富氧空气两种流体． 白铜锍的物性参数如 表 4 所列． 多相流模型选用 VOF 模型; 湍流模型选用 ＲNG k- ε 模型; 动量方程的离散用二阶迎风格式; 压力的离散 用 PＲESTO 格式; 压力速度耦合方式选用 PISO 算法， 这种耦合方法适用于非定常流动的计算; 为确保计算 的精度以及收敛性，基本时间步长定为 0. 00001 s，计 算过程中根据收敛的实际情况适当减小或增大时间步 长; 为提高计算收敛性，对部分控制方程采用欠松弛法 求解． 计算稳定后部分时刻过喷嘴中心横截面的气相 体积分数分布如图 5 所示． 表 4 白铜锍的物性参数 Table 4 Property parameters of white metal 密度/( kg·m － 3 ) 表面张力/( N·m － 1 ) 温度/K 黏度/( Pa·s) 比热容/( J·kg － 1·K － 1 ) 热导率/( W·m － 1·K － 1 ) 5400 0. 31 1473 0. 004 540 8. 9 图 5 不同时刻中心截面的气相体积分数分布． ( a) 8. 6 s; ( b) 8. 8 s; ( c) 9 s Fig． 5 Volume fraction distribution of gas in the center plane at different moments: ( a) 8. 6 s; ( b) 8. 8 s; ( c) 9 s 2. 2 界面积计算方法在模拟结果分析中的应用 利用前述方法，根据数值模拟结果，计算流体运动 相对稳定的 8 ～ 9 s 时间内多个时刻点熔池中总相界面 积，如图 6 所示． 由图可知，由于气液相混合流体内部 的扰动及自由液面的波动，熔池内总相界面积随着时 间的变化，一般在 7. 53 ～ 8 m2 范围内波动． 为进一步了解熔池 中 相 界 面 的 分 布 情 况，计 算 8 ～ 9 s内不同熔池高度以下相界面积的时间平均值，如 图 7 所示． 由 图 可 知: 在 喷 嘴 以 下 区 域 ( 熔 池 高 度 ＜ 0. 565 m) 相界面积接近零，喷嘴以下存在的富氧空 气极少; 相界面主要分布在喷嘴以上、自由液面以下区 域; 自由液面以上仍存在一定的相界面积，主要是由于 熔体的喷溅引起． 在铜锍吹炼的造铜期，发生的主要化学反应为硫 化亚铜的氧化反应，在近 1200 ℃ 的温度下，该反应的 速率受气相中氧的传质速率所控制，即氧化反应速率 等于氧的传质速率［15］: l · O2 = kO2 ·A·cO2 ． ( 7) 式中: l · O2 为气相中氧传递速率，mol·s － 1 ; kO2 为氧传质系 数，m·s － 1，根据文献［16］，此处取 0. 6 m·s － 1 ; A 为气液 相界面积，m2 ; cO2 为气相中氧浓度，mol·m － 3 ． 图 6 不同时刻熔池内总相界面积 Fig． 6 Total interfacial area in the bath at different times 由数值模拟结果图 5 和图 7 可知，气相在喷嘴以 下区域的分布极少，相应地，喷嘴以下区域发生的氧化 反应可忽略不计． 为了进一步研究喷嘴以上区域的化 学反应情况，对喷嘴以上( 0. 565 ～ 1. 300 m) 熔池内划 分为 48 个水平层． 利用仿真数据以及本文方法可获 取各层的气液相界面积、气相的平均停留时间等，根据 下式可计算各层中的氧气消耗量: lO2 = kO2 ·A·j cO2 ·j tj ． ( 8) 式中，下标 j 表示第 j 层，tj 是气体在第 j 层的平均停留 · 5901 ·
·1096· 工程科学学报，第38卷，第8期 1,8 -0 > 1.4 12 3 0.8 2 ◆ 0.4 02 0 0 0.4 0.81.21.62.02.4 85 0.60.70.80.91.01.11.21.3 层底相对于炉底高度m 距离炉底的高度m 图10氧气和二氧化硫在不同高度处的摩尔流量 图7指定高度以下熔池内相界面积的时均值 Fig.10 Molar flow rate of 0,and SO,at different heights Fig.7 Mean interfacial area below the given heigh 时间.由于层高较小，各层气相中平均氧浓度可用底 反应会进一步发生，因此可以认为这一计算结果与文 面（气相入口）氧浓度近似.由于喷嘴处的氧浓度已 献6]报道和实际测量结果(95%~98%)基本一致. 知，因此利用式(8)可自下而上逐层计算氧气的消耗 另一方面，综合图5和图8~图10信息可以发现，氧化 量（如图8）、体积分数（如图9）以及摩尔流量（如 反应主要发生于0.56~1.00m高度范围的靠近风嘴 图10). 侧壁面的熔池区域（通常称为“高效反应区”），由于该 0.14 反应属于放热反应，因此该区域的壁面承受的热负荷 相对较大，理论上更易受到破坏.这一推论与工程实 0.12 际相符合 0.10 0.08 3结论 0.06 借鉴PLIC-VOF界面重构的思想，提出一种通过 对气液两相流数值模拟结果后处理，从而计算获取两 0.02 相相界面积的方法.该方法形式简单，易于计算机实 0.5 0.6 0.70.80.91.01.11.21.3 现，并可以方便地与FLUENT等流体仿真软件平台集 层底相对于炉底高度m 成，与气液相界面积的实验测定方法相比，不需搭建实 图8氧气在各层的消耗量 验平台，对两相物质的透光性等无特殊要求，具有更广 Fig.8 Consumption of 02 in each layer 的应用范围.以FLUENT软件为平台开发了相应的应 用程序，并将其用于分析铜锍吹炼PS转炉的数值模拟 0.25 结果.应用效果表明：(1)该方法能够有效地获取两相 0.20 流体系中不同时间、空间范围内的相界面积信息，为定 量分析两相流体系内的传热、传质以及物理化学反应 0.15 状态提供了依据.(2)基于获取的相界面积，分析了熔 0.10 池内主要化学反应的动力学特点，计算了氧气利用率， 识别出熔池内部高效反应区，计算和识别结果与实际 相吻合，在一定程度上证实了本文方法的准确性 85060.70.80.91.0111.21.3 相对于炉底的高度m 参 考文献 图9氧气在不同高度处的体积分数 Song HH.Wang X L,Li HH.Measurement and prediction of Fig.9 Volume fraction of 0,at different heights interfacial area on distillation tray.J Chem Ind Eng,2003,54 由图8~图10可知，气体在熔池内上升过程中氧 (8):1112 不断被氧化反应消耗，气体中氧含量逐渐降低（反应 (宋海华，王秀丽，李红海.精馏塔板上气液相界面积的测量 与预测.化工学报，2003,54(8)：1112) 产物二氧化硫含量逐渐升高)，在自由界面附近，氧气 2] Yang W G,Wang J F,Wang T F,et al.Experimental study on 体积分数约为1.17%，由此计算得到的氧气利用率约 gas-iquid interfacial area and mass transfer coefficient in three- 为94.4.考虑到在自由界面以上气液仍有接触，氧化 phase circulating fluidized beds.Chem Eng J,2001,84(3):485
工程科学学报，第 38 卷，第 8 期 图 7 指定高度以下熔池内相界面积的时均值 Fig． 7 Mean interfacial area below the given height 时间． 由于层高较小，各层气相中平均氧浓度可用底 面( 气相入口) 氧浓度近似． 由于喷嘴处的氧浓度已 知，因此利用式( 8) 可自下而上逐层计算氧气的消耗 量( 如图 8 ) 、体 积 分 数 ( 如 图 9 ) 以 及 摩 尔 流 量 ( 如 图 10) ． 图 8 氧气在各层的消耗量 Fig． 8 Consumption of O2 in each layer 图 9 氧气在不同高度处的体积分数 Fig． 9 Volume fraction of O2 at different heights 由图 8 ～ 图 10 可知，气体在熔池内上升过程中氧 不断被氧化反应消耗，气体中氧含量逐渐降低( 反应 产物二氧化硫含量逐渐升高) ，在自由界面附近，氧气 体积分数约为 1. 17% ，由此计算得到的氧气利用率约 为 94. 4． 考虑到在自由界面以上气液仍有接触，氧化 图 10 氧气和二氧化硫在不同高度处的摩尔流量 Fig． 10 Molar flow rate of O2 and SO2 at different heights 反应会进一步发生，因此可以认为这一计算结果与文 献［16］报道和实际测量结果( 95% ～ 98% ) 基本一致． 另一方面，综合图 5 和图 8 ～ 图 10 信息可以发现，氧化 反应主要发生于 0. 56 ～ 1. 00 m 高度范围的靠近风嘴 侧壁面的熔池区域( 通常称为“高效反应区”) ，由于该 反应属于放热反应，因此该区域的壁面承受的热负荷 相对较大，理论上更易受到破坏． 这一推论与工程实 际相符合． 3 结论 借鉴 PLIC--VOF 界面重构的思想，提出一种通过 对气液两相流数值模拟结果后处理，从而计算获取两 相相界面积的方法． 该方法形式简单，易于计算机实 现，并可以方便地与 FLUENT 等流体仿真软件平台集 成，与气液相界面积的实验测定方法相比，不需搭建实 验平台，对两相物质的透光性等无特殊要求，具有更广 的应用范围． 以 FLUENT 软件为平台开发了相应的应 用程序，并将其用于分析铜锍吹炼 PS 转炉的数值模拟 结果． 应用效果表明: ( 1) 该方法能够有效地获取两相 流体系中不同时间、空间范围内的相界面积信息，为定 量分析两相流体系内的传热、传质以及物理化学反应 状态提供了依据． ( 2) 基于获取的相界面积，分析了熔 池内主要化学反应的动力学特点，计算了氧气利用率， 识别出熔池内部高效反应区，计算和识别结果与实际 相吻合，在一定程度上证实了本文方法的准确性． 参 考 文 献 ［1］ Song H H，Wang X L，Li H H． Measurement and prediction of interfacial area on distillation tray． J Chem Ind Eng，2003，54 ( 8) : 1112 ( 宋海华，王秀丽，李红海． 精馏塔板上气液相界面积的测量 与预测． 化工学报，2003，54( 8) : 1112) ［2］ Yang W G，Wang J F，Wang T F，et al． Experimental study on gas-liquid interfacial area and mass transfer coefficient in threephase circulating fluidized beds． Chem Eng J，2001，84( 3) : 485 · 6901 ·
宋彦坡等：基于数值模拟的气液相界面积计算方法 ·1097· B3]Cents A HG,Brilman D W F,Versteeg G F.Gas absorption in 1013 agitated gas-iquid-liquid system.Chem Eng Sci,2001,56(3): (张勤，陈泽敬，李志信.PLC方法在浮法成型模拟中的应 1075 用.工程热物理学报，2010,31(6)：1013) 4]Liang Y C,Zhou Z,Zhang Z B.Interfacial area on valve trays. 02] Liu CL,Wu S Z,Zhang J M,et al.Numerical simulation on Chem Ind Eng,2007,58(5):1232 the surface fluctuation of molten steel in a wide slab continuous (粱银春，周政，张志炳.浮阀塔板气液相界面积.化工学报， casting mold.J Univ Sci Technol Beijing,2009,31 (2):229 2007,58(5):1232) (刘国林，吴苏州，张炯明，等.宽板坯连铸结品器内液面波 5]Danckwerts PV,Lannus A.Gas-iquid reactions.Electrochem 动的数值模拟.北京科技大学学报，2009,31(2)：229) Soc,1970,117(10):369C [13]Maurya R S,Sundararajan T,Das S K.Development of a PLIC- 6]Ma D,Liu M Y,Zu Y G,et al.Two-dimensional volume of fluid -VOF method for the dynamic simulation of entry region flow in a simulation studies on single bubble formation and dynamics in laminar falling film.Int J Comput Fluid Dyn,2009,23 (5): bubble columns.Chem Eng Sci,2012,72(16):61 391 7]Zhang Y J,Liu M Y,Xu Y G,et al.Three-dimensional volume [14]Yang J Y,Jia Y X,Ding YY,et al.Numerical simulation of of fluid simulations on bubble formation and dynamics in bubble mold filling stage in resin transfer molding processes based on columns.Chem Eng Sci,2012,73 (7):55 VOF/PLIC method.Chin J Mater Res,2009,23(6):622 [8]Xu L J,Chen G,Shao J B.Numerical simulation of bubble be- (杨俊英，贾玉玺，丁妍羽，等.基于VOF/PLIC界面追踪方 havior using VOF method.Electr Technol Cir Eng,2011,22 法的RTM工艺充模过程数值模拟.材料研究学报，2009,23 (24):2408 (6):622) 9]Mencinger J,Zun I.A PLIC-VOF method suited for adaptive [15]Ashman D W,Mckelliget J W,Brimacombe J K.Mathematical moving grids.J Comput Phys,2011,230(3):644 model of bubble formation at the tuyeres of a copper converter. [10]Lopez J,Hernandez J,Gomez P,et al.An improved PLIC-VOF Can Metall,1981,20(4):387 method for tracking thin fluid structures in incompressible two- [16]Zhu ZZ.He J Q.Modern Copper Metallurgy.Beijing:Science phase flows.J Comput Phys,2005,208 (1):51 Pres5,2003 [11]Zhang Q.ChenZJ,LiZ X.Application of PLIC method in sim- (朱祖泽，贺家齐.现代铜治金学.北京：科学出版社， ulation of float glass forming.J Eng Thermophys,2010,31(6): 2003)
宋彦坡等: 基于数值模拟的气液相界面积计算方法 ［3］ Cents A H G，Brilman D W F，Versteeg G F． Gas absorption in agitated gas-liquid-liquid system． Chem Eng Sci，2001，56 ( 3) : 1075 ［4］ Liang Y C，Zhou Z，Zhang Z B． Interfacial area on valve trays． J Chem Ind Eng，2007，58( 5) : 1232 ( 梁银春，周政，张志炳． 浮阀塔板气液相界面积． 化工学报， 2007，58( 5) : 1232) ［5］ Danckwerts P V，Lannus A． Gas-liquid reactions． J Electrochem Soc，1970，117( 10) : 369C ［6］ Ma D，Liu M Y，Zu Y G，et al． Two-dimensional volume of fluid simulation studies on single bubble formation and dynamics in bubble columns． Chem Eng Sci，2012，72( 16) : 61 ［7］ Zhang Y J，Liu M Y，Xu Y G，et al． Three-dimensional volume of fluid simulations on bubble formation and dynamics in bubble columns． Chem Eng Sci，2012，73( 7) : 55 ［8］ Xu L J，Chen G，Shao J B． Numerical simulation of bubble behavior using VOF method． Electr Technol Civ Eng，2011，22 ( 24) : 2408 ［9］ Mencinger J，Zun I． A PLIC-VOF method suited for adaptive moving grids． J Comput Phys，2011，230( 3) : 644 ［10］ Lopez J，Hernandez J，Gomez P，et al． An improved PLIC--VOF method for tracking thin fluid structures in incompressible twophase flows． J Comput Phys，2005，208( 1) : 51 ［11］ Zhang Q，Chen Z J，Li Z X． Application of PLIC method in simulation of float glass forming． J Eng Thermophys，2010，31( 6) : 1013 ( 张勤，陈泽敬，李志信． PLIC 方法在浮法成型模拟中的应 用． 工程热物理学报，2010，31( 6) : 1013) ［12］ Liu G L，Wu S Z，Zhang J M，et al． Numerical simulation on the surface fluctuation of molten steel in a wide slab continuous casting mold． J Univ Sci Technol Beijing，2009，31( 2) : 229 ( 刘国林，吴苏州，张炯明，等． 宽板坯连铸结晶器内液面波 动的数值模拟． 北京科技大学学报，2009，31( 2) : 229) ［13］ Maurya Ｒ S，Sundararajan T，Das S K． Development of a PLIC- -VOF method for the dynamic simulation of entry region flow in a laminar falling film． Int J Comput Fluid Dyn，2009，23 ( 5 ) : 391 ［14］ Yang J Y，Jia Y X，Ding Y Y，et al． Numerical simulation of mold filling stage in resin transfer molding processes based on VOF /PLIC method． Chin J Mater Ｒes，2009，23( 6) : 622 ( 杨俊英，贾玉玺，丁妍羽，等． 基于 VOF /PLIC 界面追踪方 法的 ＲTM 工艺充模过程数值模拟． 材料研究学报，2009，23 ( 6) : 622) ［15］ Ashman D W，Mckelliget J W，Brimacombe J K． Mathematical model of bubble formation at the tuyeres of a copper converter． Can Metall Q，1981，20( 4) : 387 ［16］ Zhu Z Z，He J Q． Modern Copper Metallurgy． Beijing: Science Press，2003 ( 朱祖泽，贺 家 齐． 现 代 铜 冶 金 学． 北 京: 科 学 出 版 社， 2003) · 7901 ·