2015, Vol. 38
四氯化钛是生产海绵钛及其合金的重要中间原料，目前，其生产工艺包括沸腾氯化(流态化氯化)、熔盐氯化等，其中，流态化氯化以氯气为气体相，钛渣和石油焦的混合物为固相，反应在气-固两相中进行，与其他生产工艺比较，流态化氯化具有成本低、原料和设备利用率高及热交换效率高等诸多优势[1-2]，因此被广泛应用于国外四氯化钛的生产。
流态化氯化温度一般控制在1 000~1 100 ℃。钛渣因含有Ca、Mg等元素，当氯化温度为800 ℃时，Ca和Mg等元素优先于Ti氯化生成CaCl2和MgCl2，由于CaCl2和MgCl2的熔点低(分别为772 ℃和714 ℃)，故在氯化后两者呈熔融状态，随着氯化时间的延长导致整个床层出现粘结失流，因此，目前流态化氯化工艺对原料质量要求较高：TiO2品位越高越好(一般大于92%)，Ca和Mg等杂质的总含量低(一般小于1.5%)[3]。目前，西昌地区的钛渣大部分不能满足流态化氯化的要求，许多研究者针对该钛渣特性，通过对流化床反应器进行设计改造以满足流态化的要求。徐聪等[4]采用一种将提升管和湍床串联使用的组合式流化床氯化钛渣，减小了流化床内粘结。流态化氯化粘结导致流化床产生失流的现象，床内气固两相流动性复杂，具有非线性，可通过对流化床反应器的设计和放大进行解决，目前，实验室研究存在投入高、耗时长等缺点。随着计算机技术的发展，数值模拟逐渐成为气固两相流研究的重要方法之一。通过研究钛渣流态化过程中不同临界流化速度下的流动特征，可进一步判断床内气固相转变，为反应器设计奠定基础。
笔者建立了一个钛渣流态化氯化流动特性预测模型，主要讨论了表观操作气速、气固相间的曳力等方面的内容。依据传统的最小流化速度和完全流化速度的概念和计算方法计算了钛渣的最小流化速度和完全流化速度；通过对B类颗粒的钛渣在相关临界流化速度下的流态化过程进行冷态数值模拟计算，得到了钛渣流态化过程的流体动力学特征；预测模型结合流体动力学软件Fluent在相关临界流化速度下得到流动特性与相关研究一致。
1 钛渣流态化临界流化速度Wen-yu在大量实验研究基础上，提出了雷诺数范围0.001~4 000，粒径范围为0.002~1.97 cm，床层孔隙率范围为0.136~1的最小流化速度计算公式[5]为
| $ \frac{{{d_{\text{p}}}{u_{{\text{mf}}}}{\rho _{\text{g}}}}}{{{\mu _{\text{g}}}}} = {\left[{C_1^2 + {C_2}\frac{{d_{\text{p}}^3{\rho _{\text{g}}}\left( {{\rho _{\text{p}}}-{\rho _{\text{g}}}} \right)g}}{{{\mu ^2}}}} \right]^{1/2}} -{C_1}, $ | (1) |
式中：dp为颗粒直径；umf为最小流化速度；ρg为气体密度；μg为气体黏度；d为颗粒密度；C1和C2为常数。
Grace等[6]在1982年提出：B类颗粒气速达到最小流化速度气泡在分布板位置形成长大最后到上表面破裂，针对气固两相流，当式(1)中C1=27.2, C2=0.040 8时更能满足气固系统最小流化速度的计算。
Gupta等[7]提出了基于多组分颗粒系统重颗粒和大颗粒开始流化的颗粒系统完全流化速度的计算公式，计算出相应的理论计算值为
| $ \frac{{1.75}}{{{\varphi _{\text{s}}}\varepsilon _{{\text{tf}}}^3}}Re_{{\text{tf}}}^2 + \frac{{150\left( {1-{\varepsilon _{{\text{tf}}}}} \right)}}{{\varphi _{\text{s}}^2\varepsilon _{{\text{tf}}}^3}}R{e_{{\text{tf}}}} = Ar, $ | (2) |
式中：
钛渣流态化数学模型以双流体模型为基础，通过Fluent软件计算钛渣在最小流化速度和完全流化速度下的流动特性。
2.1 双流体数学模型双流体模型是将颗粒作为拟流体，认为固相和流体相共同存在且互相渗透的连续介质模型[8]。Fluent软件提供的描述气-固两相流的流体连续和动量控制有以下方程。
连续方程：
气相为
| $ \frac{{\partial \left( {{\varepsilon _{\text{g}}}{\rho _{\text{g}}}} \right)}}{{\partial t}} + \nabla \cdot \left( {{\varepsilon _{\text{g}}}{\rho _{\text{g}}}{\boldsymbol{\mu} _{\text{g}}}} \right) = 0, $ | (3) |
固相为
| $ \frac{{\partial \left( {\left( {1-{\varepsilon _{\text{g}}}} \right){\rho _{\text{p}}}} \right)}}{{\partial t}} + \nabla \cdot \left( {\left( {1-{\varepsilon _{\text{g}}}} \right){\rho _{\text{g}}}{\boldsymbol{\mu} _{\text{p}}}} \right). $ | (4) |
动量方程：
气相为
| $ \begin{gathered} \frac{{\partial \left( {{\varepsilon _{\text{g}}}{\rho _{\text{g}}}{\mu _{\text{g}}}} \right)}}{{\partial t}} + \nabla \cdot \left( {{\varepsilon _{\text{g}}}{\rho _{\text{g}}}{\boldsymbol{\mu} _{\text{g}}}{\boldsymbol{\mu} _{\text{g}}}} \right) = \hfill \\ \nabla \cdot \left( {{\varepsilon _{\text{g}}}{\tau _{\text{g}}} + {\varepsilon _{\text{g}}}{\rho _{\text{g}}}g-{\varepsilon _{\text{g}}}\nabla P-\beta \left( {{\boldsymbol{\mu} _{\text{g}}}-{\boldsymbol{\mu} _{\text{p}}}} \right)} \right), \hfill \\ \end{gathered} $ | (5) |
固相为
| $ \begin{gathered} \frac{{\partial \left( {1-{\varepsilon _{\text{g}}}){\rho _{\text{p}}}{\mu _{\text{p}}}} \right)}}{{\partial t}} + \nabla \cdot \left( {1-{\varepsilon _{\text{g}}}){\rho _{\text{p}}}{\boldsymbol{\mu} _{\text{p}}}{\boldsymbol{\mu} _{\text{p}}}} \right) = \hfill \\ \nabla \cdot \left( {1-{\varepsilon _{\text{g}}}} \right){\tau _{\text{p}}} + \left( {1 - {\varepsilon _{\text{g}}}} \right){\rho _{\text{p}}}g - \left( {1 - {\varepsilon _{\text{g}}}} \right)\nabla P - \beta \left( {{\boldsymbol{\mu} _{\text{g}}} - {\boldsymbol{\mu} _{\text{p}}}} \right), \hfill \\ \end{gathered} $ | (6) |
式中：ε为空隙率；ρ为密度；μ为黏度；β为气固界面曳力函数；t为时间；ΔP为压强；下标p和g分别表示固体颗粒和气体。
2.2 模型方程主要参数的确定为求解上述模型方程，需要明确4个重要的参数，分别为气相黏性系数、固相黏性系数、颗粒相压力气和固相间曳力系数。
气固两相黏性系数根据牛顿粘性定律计算[9]为
| $ {\tau _m} = {u_m}\left[{\nabla {u_m} + {{\left( {{u_m}} \right)}^{\text{T}}}-\frac{2}{3}\left( {\nabla {u_m}} \right)I} \right], $ | (7) |
式中，下标m可取p和g，分别代表固相和气相。
颗粒相压力采用链式求导得
| $ \frac{{\partial {\rho _{\text{s}}}}}{{\partial {x_i}}} = \frac{{\partial {\rho _{\text{s}}}}}{{\partial {\varepsilon _{\text{s}}}}}\frac{{\partial {\varepsilon _{\text{s}}}}}{{\partial {x_i}}} = G\left( {{\varepsilon _{\text{g}}}} \right)\frac{{\partial {e_{\text{s}}}}}{{\partial {x_{\text{i}}}}} $ | (8) |
式中，G(εg)为固相碰撞系数，可以用指数函数表达[14]为
| $ G\left( {{\varepsilon _{\text{g}}}} \right) =-{G_0}\exp \left( {C\left( {{\varepsilon _0}-{\varepsilon _{\text{g}}}} \right)} \right), $ | (9) |
式中，相关参数取值为：G0=1 Pa，C=230。
2.3 曳力模型合理的曳力模型被用于封闭动量方程(3)和(4)中的气固相间阻力系数。为了能够正确模拟出气泡相，气固两相之间的Gidaspow曳力模型能够计算颗粒体积分数较小的情况，定量预测气泡相的效果比较好[10]，因此，选择该曳力模型进行模拟研究。
Gidaspow曳力模型为
| $ \begin{gathered} \beta = 150\frac{{{{\left( {1-{\varepsilon _{\rm{g}}}} \right)}^2}{\mu _{\rm{g}}}}}{{{\varepsilon _{\rm{g}}}{{\left( {{\Phi _{\rm{p}}}{d_{\rm{p}}}} \right)}^2}}} + 1.75\frac{{\left( {1-{\varepsilon _{\rm{g}}}} \right){\rho _{\rm{g}}}\left| {{\mathit{\boldsymbol{\mu }}_{\rm{g}}}-{\mathit{\boldsymbol{\mu }}_{\rm{p}}}} \right|}}{{{\Phi _{\rm{p}}}{d_{\rm{p}}}}}\left( {{\varepsilon _{\rm{g}}} < 0.8} \right), \hfill \\ \beta = \frac{3}{4}CD\frac{{{\varepsilon _{\rm{g}}}\left( {1 - {\varepsilon _{\rm{g}}}} \right){\rho _{\rm{g}}}\left| {{\mathit{\boldsymbol{\mu }}_{\rm{g}}} - {\mathit{\boldsymbol{\mu }}_{\rm{p}}}} \right|}}{{{\Phi _{\rm{p}}}{d_{\rm{p}}}}}{\varepsilon _{\rm{g}}}^{ - 2.65}\;\;\;\;\;\;\;\;\;\;\;\;\left( {{\varepsilon _{\rm{g}}} > 0.8} \right), \hfill \\ \end{gathered} $ | (10) |
式中，相间阻力系数β为曳力系数的函数。曳力系数CD表示为
| $ CD = \left\{ {\begin{array}{*{20}{c}} {\frac{{24}}{{{\text{Rep}}}}\left[{1 + 0.15{{\left( {{\varepsilon _{\text{g}}}{\text{Rep}}} \right)}^{0.687}}} \right]}&{\left( {{\text{Rep}} > 1000} \right), } \\ {0.44}&{\left( {{\text{Rep}} < 1000} \right);} \end{array}} \right. $ | (11) |
式中：dp为颗粒直径，m；ρg为气体密度；ρp为颗粒密度；Φs为球形度；umf为最小流化速度；τ为黏滞系数；β为曳力函数。
3 不同临界流化速度下流动特性数值模拟 3.1 钛渣流态化物理模型按照文献[11]的参数，同时为了便于气泡观察，数值模拟设计准三维流化床冷态模型。
模型具体尺寸(200 mm×20 mm×600 mm)，分布板上均匀布置(梅花型布孔方式)了46个直径1 mm的小孔，分3排，孔中心距为1.5 mm，开孔率为0.85%。数值模拟的几何模型如图 1所示。
|
图 1 流化床2D模型 Figure 1 2D model of fluidized-bed |
攀枝花钛渣的相关物性和计算机模拟条件参数如表 1所示。钛渣氯化流态化过程主要研究最小流化速度和完全流化速度下的流动特性。最小流化速度的理论计算应用Wen-yu[5]提出的经验计算公式和Grace等[6]对Wen-yu公式的修正公式；当床层进入稳定流化阶段时，则应用Gupta[7]提出的完全流化速度经验公式进行计算。
|
表 1 模拟计算的边界条件 Table 1 Boundary conditions for simulation calculation |
模型采用Fluent 6.3.26计算流体力学软件包的数值模拟，即综合式(1)-(11)的求解，其中Gidaspow模型为Fluent软件自带模型。模型方程采用有限容积法进行离散，计算方式为压力基隐式非稳态求解，应用Simple算法，选择粘性模型RNGk-ε，壁面处气固两相采用无滑移边界条件。动量差分格式为二阶迎风格式，为了使相间界面更加清晰和提高计算精度，体积分数差分格式采用QUICK格式。计算时连续性方程X和Y两个方向的动量方程残差精度为10-4，湍动能k和湍动能耗散率ε残差精度为10-3，其他参数取值见表 1所示。
3.4 数值模拟结果分析最小流化速度下的钛渣氯化流态化过程的数值模拟是研究在最小流化速度下床层的启动过程，研究对象主要是床层开始形成气泡的位置、大小和时间。在气泡认定时，将考虑气体孔隙率大于85%的区域认定为气泡[12]，由于气泡并非规则的圆形，一般以当量直径考虑，即以与气泡面积相等的圆的直径作为气泡直径，计算气泡当量直径的公式为
| $ {d_{\rm{p}}} = \sqrt {\frac{S}{{\frac{1}{4}{\rm{\pi }}}}}, $ | (12) |
式中：dp为气泡直径，m；S为气泡当量面积，m2。
3.4.1 最小流化速度下钛渣流态化氯化过程分析钛渣氯化过程最小流化速度下数值模拟结果主要表达了钛渣氯化的启动阶段。文中采用Fluent软件计算了Wen-yu方程和Grace对Wen-yu方程的修正式预测的最小流化速度下的床层行为，模拟结果如图 2所示。
|
图 2 最小流化速度下钛渣体积分数数值模拟结果 Figure 2 Fluidized bed of Titanium slag Solids volume at minimum fluidization velocity |
数值模拟对象为粒径500 μm、密度3 800 kg/cm3的钛渣，该实验对象符合Geldart [13]对颗粒分类的B类(100 μm≤d≤600 μm，1 400 kg/m3≤ρ≤4 000 kg/m3)，因其最小流化速度与初始鼓泡速度相等，气速一旦达到最小流化速度床层内就会出现气泡相和密相，因此B类颗粒常被作为鼓泡床的流化介质。本次数值模拟选择的钛渣为B类颗粒的范畴，按照B类颗粒的流动特性，当达到最小流化速度时，床层在分布板位置形成，随着表观流化时间延长，气泡上升、合并长大最后破裂。
根据数值模拟结果分析可知：在0~0.5 s之间，按照气泡理论分析床层开始并未形成气泡，床层形成的是乳相，随着流化时间的延长，乳相层沿床层逐渐上升和增宽，床层出现均匀膨胀。乳相的形成表明流态化初期沿着床层轴向存在速度梯度，这样就导致了在高度方向出现流态化不均匀现象，在流速较低的地方导致床层均匀膨胀，在流速较高的地方出现乳相。在0.5 s后，Wen-yu方程预测得到的最小流化速度下Fluent数值计算结果表明：床层同样是乳相，而Grace对Wen-yu修正公式预测得到最小流化速度在0.7 s时，床层中乳相中开始有气泡产生。因此，B类颗粒的钛渣进行流态化氯化时，Grace方程相对于Wen-yu方程能够更加准确地预测钛渣氯化流态化的最小流化速度。
3.4.2 完全流化速度下钛渣氯化流态化过程分析钛渣氯化过程完全流化速度下数值模拟主要了解钛渣进入充分流化后的流动特性，分析床层气固两相在床内体积分数的分布情况，文中采用Fluent软件计算了Gupta方程预测得到的完全流化速度下床层行为，模拟结果如图 3所示。
|
图 3 完全流化速度下钛渣体积分数 Figure 3 Fluidized bed of Titanium slag Solids volume at complete fluidization velocity |
根据图 3分析可知，在完全流化速度下，钛渣氯化流态化在0~0.5 s之间，未进入鼓泡床阶段直接进入节涌床，随着流化时间的延长，床层开始进入完全流化状态，床层进入节涌床阶段主要是由于设计的流化床反应器床体高径比较大，导致在0~0.5 s床层进入节涌床。在0.5 s之后，气泡在床层中心和壁面位置形成，随着流化时间延长，气泡分布板位置开始形成的气泡和原有气泡之间发生并聚，导致气泡在上升过程逐渐长大，最后破裂。
3.4.3 钛渣完全流化状态下速度矢量分析图 4是流态化氯化还原床层进入完全流化状态在0.7~1.1 s床层钛渣的体积分数和速度矢量图，从图中可以看出钛渣的体积分数和速度，进入完全流态化后整个床层在气泡的扰动下进行循环；中心区域固体颗粒在气泡上升过程中呈现上升趋势，中心区域和气泡之间的固体颗粒出现向下运动趋势，气泡和壁面之间颗粒向下运动。数值模拟得到的钛渣流态化流化在完全流化速度下的固体颗粒速度矢量图与相关文献的研究结果(见图 5)相同。
|
图 4 颗粒流场速度分布 Figure 4 Velocity maps |
|
图 5 相关研究结果 Figure 5 Ralated research results |
通过式(12)气泡大小以面积大小计量，结合图 3分析，从气泡理论认定气泡，钛渣进入完全流化状态时，前0.5 s床层处于节涌状态，从0.7 s开始，床层气泡直径大小随着流化时间的延长，气泡直径逐渐增大。该数值模拟结果与文献[17]中结果一致。
为进一步从气泡变化的角度定量分析该模型的正确性，文中在完全流化速度下的数值模拟结果中选取了一组气泡(气泡从分布板位置开始形成到气泡破裂)变化过程，图 6研究了气泡直径与分布板之间距离的关系。
|
图 6 气泡直径随分布板距离变化 Figure 6 Relationship bubble diameter with vertical distance along the distributor |
从图 6可知，钛渣流态化下床层的气泡直径随着床高的增加在不断增加，气泡尺寸与距分布板距离近似呈线性函数关系，玻璃珠床实验值拟合的线性函数关系为d=0.618 5 h-2.995，线性拟合度达到98.05%；同时，Price等[18]和Fan等[19]也做了相似的数值实验研究，文中数值实验结果与他们的研究结果一致。
4 结论文中以双流体模型为基础研究了钛渣流动特性数学模型；结合钛渣的物性, 采用传统的经验公式计算了钛渣流态化系统的最小流化速度、完全流化速度。研究钛渣在最小流化速度和完全流化速度下的流动特性。研究结果表明：
1) 以双流体模型为基础得到的钛渣流动特性数学模型，成功地模拟出了钛渣流态化氯化的流动特征，为反应器设计提供数值模拟基础。
2) B类颗粒的钛渣最小流化速度数值模拟结果显示，Wen-yu预测公式整个床层基本出现的是乳相，而Grace修正公式在0.7 s时能够成功地模拟出床层气泡行为。
3) 完全流态化时，数值模拟表明钛渣床层首先要经历一个节涌过程，然后再进入完全流化状态。
4) 完全流化状态下，文中数学模型对于钛渣流态化床层中固体矢量图和气泡尺寸大小变化规律的预测均与相关文献研究结果一致，表明该数学模型适合钛渣流态化氯化过程数值模拟。
| [1] | Jin Y. Fluidization engineering principles[M]. BeiJing: Tsinghua University Press, 2001. |
| [2] | Sun K. Titanium extraction metallurgy physical chemistry[M]. BeiJing: Metallurgical Industry Press, 2001: 38-42. |
| [3] | Guo M, Li H. Review and prospect of fluidization science and technology[J]. Journal of Chemical Industry and Engineering, 2012, 64(1): 52–62. |
| [4] | Xu C, Yuan Z, Wang X. One dimensional modeling of multiple unit pneumatic transport reactor for producing titanium tetrachloride. mathematical model[J]. The Chinese Journal of Process Engineering, 2005, 1: 18–22. |
| [5] | Wen C, Yu Y H. A generalized method for predicting the minimum fluidization velocity[J]. Aiche Journal, 1966, 12(3): 610–612. DOI:10.1002/(ISSN)1547-5905 |
| [6] | Grace J R, Taghipour F. Verification and validation of CFD models and dynamic similarity for fluidized beds[J]. Powder Technology, 2004, 139(2): 99–110. DOI:10.1016/j.powtec.2003.10.006 |
| [7] | Gupta S K, Agarwal V K, Singh S N, et al. Prediction of minimum fluidization velocity for fine tailings materials[J]. Powder Technology, 2009, 196(3): 263–271. DOI:10.1016/j.powtec.2009.08.003 |
| [8] | Rüdisüli M, Schildhauer T J, Biollaz S M A, et al. Bubble characterization in a fluidized bed by means of optical probes[J]. International Journal of Multiphase Flow, 2012: 56–67. |
| [9] | Gopalakrishnan P, Garg R, Shahnam M. CFD-DEM study of effect of bed thickness for bubbling fluidized beds[J]. Particuology, 2012, 10(5): 532–541. DOI:10.1016/j.partic.2012.02.006 |
| [10] | Gidaspow D. Multiphase flow and fluidization: continuum and kinetic theory descriptions[M]. Academic press, 1994. |
| [11] | Grace J R, Taghipour F. Verification and validation of CFD models and dynamic similarity for fluidized beds[J]. Powder Technology, 2004, 139(2): 99–110. DOI:10.1016/j.powtec.2003.10.006 |
| [12] | Annaland S, Van M. Critical comparison of hydrodynamic models for gas-solid fluidized beds -Part I: bubbling gas-solid fludized beds operated with a jet[J]. Chemical Engineering Science, 2005, 60(1): 57–72. DOI:10.1016/j.ces.2004.07.059 |
| [13] | Geldart D. Types of gas fluidization[J]. Powder Technology, 1973, 7(5): 285–292. DOI:10.1016/0032-5910(73)80037-3 |
| [14] | Goldschmidt M J V, Kuipers J A M, Swaaij W P M V. Hydrodynamic modelling of dense gas-fluidised beds using the kinetic theory of granular flow: effect of coefficient of restitution on bed dynamics[J]. Chemical Engineering Science, 2001, 56(2): 571–578. DOI:10.1016/S0009-2509(00)00262-1 |
| [15] | Zhou J, Zou Z, Yu A. Numerious simulation for flow characteristic in bubbling fluidized bed[J]. Journal of Materials and Metallurgy, 2007, 6(2): 126–129. |
| [16] |
杨太阳, 王安仁, 张锁江, 等.
气固鼓泡流化床的流动特性数值模拟[J]. 计算机与应用化学, 2005, 22(3): 207–210.
YANG Taiyang, WANG Anren, ZHANG Suojang, et al. Numerical simulation for fluid dynamics in a gas-solid bubbling fluidized bed[J]. Computer and Applied Chemistry, 2005, 22(3): 207–210. (in Chinese) |
| [17] |
林亮成, 郑忠, 陈伟, 等.
不同曳力模型对鼓泡床内气固两相流的模拟研究[J]. 化学反应工程与工艺, 2011, 26(5): 390–398.
LIN Liangcheng, ZHEN Zhong, CHEN Wei, et al. A simulation study of gas'solid two phase flow in a bubbling fluidized bed with various drag force models[J]. Chemical Reaction Engineering and Technology, 2011, 26(5): 390–398. (in Chinese) |
| [18] | Price G, Chandrasekaran B, Hulme I, et al. Comparison of 2-D and 3-D CFD simulations of bubbling fluidized beds with x-ray fluoroscopy and imaging experiments[J]. Powder Technology, 2011, 116: 142–154. |
| [19] | Fan X, Yang Z, Parker D J, et al. Prediction of bubble behaviour in fluidised beds based on solid motion and flow structure[J]. Chemical Engineering Journal, 2008, 140: 358–369. DOI:10.1016/j.cej.2007.10.013 |