一、二类边界条件下扩散方程的一个稳态近似结果

王成善^{1, 2}，穆小静³

(1. 重庆大学 材料科学与工程学院，重庆，400044；

2. 重庆大学 重庆市冶金工程重点实验室，重庆，400044；

3. 重庆大学 化学化工学院，重庆，400044)

摘要：将动力学处理非稳态问题的浓度随时间不变的稳态假设发展为浓度变化率随时间不变的稳态假设，对有限长度区间内扩散方程进行稳态近似法处理，获得一、二类边界条件下扩散方程的一个稳态近似解。并将其与精确数值解对比。研究结果表明：稳态近似法获得的结果和精确解随时间变化是同步的，利用近似解可以准确地预测达到最终稳态的时间；近似解与接近最终稳态的情形吻合程度好，与远离最终稳态的情形吻合程度较差；稳态近似法获得的结果基本上满足总体质量守恒。

关键词：冶金动力学；扩散方程；非稳态；稳态法；近似分析解

中图分类号：TF512           文献标志码：A         文章编号：1672-7207(2013)07-2712-08

Analysis of diffusion equation with the first and second boundary conditions based on steady state approximation

WANG Chengshan^{1, 2}, MU Xiaojing³

(1. School of Materials Science and Engineering, Chongqing University, Chongqing 400044, China;

2. Chongqing Key Laboratory of Metallurgical Engineering, Chongqing University, Chongqing 400044, China;

3. School of Chemistry and Chemical Engineering, Chongqing University, Chongqing 400044, China)

Abstract: An assumption of constant concentration variance ratio was made and substituted for the assumption of constant concentration frequently used in kinetics of process metallurgy. A detailed process to deal with diffusion equation based on steady state approximation was given, and the approximate solutions of the diffusive equation with the first and the second boundary conditions were obtained in the meantime. The approximate solutions were compared with the numerical solutions. The results show that the steady state approximation’s results and solutions synchromze with time. The diffusive process is considerably well predicted by the approximate solutions, and the approximate solutions accord with the situations closing to the final steady state a little better than those far from the final steady state. The results obtained by steady station approximation fairly satisfies the total mass balance.

Key words: kinetics of process metallurgy; diffusive equation; non-steady state; steady state approximation; analytic solutions

一维菲克扩散方程在无限长度区间[0, ∞]或[-∞, ∞]里有分析解；其一是误差函数解^[1-4]，其二是高斯解^[2]。在[a, b]有限长度的区间里，除了无穷级数^[1-4]解之外，它没有其他形式的分析解。无穷级数解利用起来不方便^[3]，一般结合线图使用而且最终利用结果是近似值；由于对质量守恒原理、扩散物理本质体现地不明显、不直接，导出误差解、高斯解和级数解等所用的分离变量等纯数学方法工程上很少应用。寻求有限长度区间[a, b]上扩散方程解的近似处理方法，可避免对误差函数解的依赖(即要获得物质由此位置经有限距离向彼位置扩散过程的解时，可以避免将彼位置处理在无穷远处)，不仅具有理论意义，也具有现实意义。因为工程中许多基本过程都会涉及在有限长度区间内的一维扩散。仅在冶金工程领域就可举出一些涉及这些扩散的较典型例子：确定对流传质系数时，涉及的边界层层流底层内的扩散或相界面两侧膜内扩散^[1]；钢液精炼时，喷入钢液中的粉粒内的传质^[1]；扩散系数的实验确定涉及的扩散实验^[2]等。这些例子涉及的扩散是分子扩散，扩散方程中的扩散系数D是分子扩散系数。当气体在多孔固体内扩散时，扩散机理不仅包括分子扩散，还可包括克努岑扩散和表面扩散。几种扩散机理同时作用形成的物质浓度的变化仍可用扩散方程描述，但方程中扩散系数应换为有效扩散系数D_e，因此，也将这至少存在2种机理的扩散叫有效扩散。冶金中涉及一维有效扩散的例子有：气体和致密的铁块矿反应时，气体通过固相产物层的扩散^[1-2]；气体和多孔的铁矿球团反应时，气体通过球团内部多孔固体的扩散^{[1, 5-11]}。而且这些有限长度区间内的一维扩散速率问题基本上应该沿有明确分析结果的方向来研究。铁矿团块、喷入钢液内每个粉粒等内部涉及的扩散都是在比装置更小的尺度上描述反应器内传质问题。从反应器整体来讲，这些传质可以看作发生在装置内某位置。当对整个反应器内流动与混合、热和质量传输、化学反应等进行综合数值模拟研究时，一般要求装置内的更小尺度上的传质速率、反应速率等最好有明确的数学分析式^[12-16]。若装置内各个更小尺度上的(或者说某位置上的)传质速率等也需像整个装置的传输问题处理一样，采用数值的方法迭代计算，这样势必增大问题的复杂性，降低处理问题的效率。本文作者将直接基于物质守恒原理、菲克扩散定律等，利用稳态法获得定义在[a, b]上非稳态扩散方程的近似分析解，并将分析解和有足够精确度的数值解进行对比。

1  处理非稳态问题的稳态法

本文获得的近似解是基于稳态近似法。它是处理连串反应动力学问题的一种方法^[17]。对连串反应A→R→C，为了推得A生成最终产物C的速率方程，假设中间产物R的浓度c_R=const (或)。在多相反应中，对由外膜传质、内扩散、界面反应等环节组成的串联式的动力学方程，也可由稳定近似法导出。例如，铁矿石未反应核模型，就是在产物层还原物质浓度c=const的稳态假设下推出的^[1-2]；金属氧化反应速率就是在假设氧化膜内氧化气体浓度不随时间变化的前提下得到的^[18]。

稳态假设就是在非稳态问题处理过程中，附加的某个(些)量随时间不变的有限约束。既然是假设，就会使结果出现偏差，所以，又要求它是有限的约束，不致结果偏差大。为了保证是有限约束，稳态假设的对象应该是诸多串行环节中的1个或几个环节。另外，为了减小稳态假设的约束，将动力学中常用的浓度随时间不变的稳态假设(c=const假设)发展为浓度变化率随时间不变的稳态假设(const或的假设)。c=const假设是const的假设的特例，换言之，c=const假设比const的假设施加的约束强，与应用c=const假设获得的结果相比，应用const的假设获得的结果更接近真实情况。

利用const的假设待处理的非稳态扩散问题如下：

初始条件：t=0：，c=0。

边界条件：t＞0：x=0，c=c₀；x=l，。

2  利用稳态法获得扩散方程近似解的过程

上述非稳态问题的解就是(0, l)内任一位置x处的浓度c_x随时间变化的关系，c_x只是时间t的函数，它的变化率为。

2.1  划分网格和离散方程

首先将位于x左边的区间长度x均分为n份，则每一份的长度是x/n；位于x右边的区间长度l-x均分为m份，则每一份的长度是(l-x)/m，如图1所示。这样，(0, l)内有(m-1)+(n-1)+1个间隔位置。将位于x左边第1个间隔位置处浓度记为c₁，位置x左边第2个间隔位置处浓度记为c₂，依次类推，位于x左边第n-1间隔位置处浓度记为c_n-1，位于x左边第n间隔位置是区间[0, l]的左边界；按照同样做法，将位于x右边第1个间隔位置处浓度记为c_l,1，位于x右边第2个间隔位置处浓度记为c_l,2，依次类推，位于x右边第m-1间隔位置处浓度记为c_l,m-1，位于x右边第m个间隔位置是区间[0, l]的右边界。

图1  非稳态扩散方程离散用图

Fig.1  Discretization grid for diffusion equation

与一般的数值方法相同^[19]，每个间隔位置处的浓度均有控制容积，浓度c₁, c₂, …, c_n-1的控制容积为x/n；浓度c_l,1, c_l,2, …, c_l,m-1的控制容积为；c_x的控制容积为。间隔位置处的浓度可以用控制容积内的平均浓度代替，间隔位置处的浓度随时间的变化率可以用控制容积内的平均浓度随时间的变化率代替；若控制容积越小，则间隔位置处的浓度、浓度的变化率分别和控制容积内的平均浓度、平均浓度的变化率越接近。将扩散方程分别针对每个控制容积离散。对c_x的控制容积，扩散方程离散为

    (1)

与一般数值方法不同的是^[19]：浓度变化率(浓度对时间的导数)不离散，仅是关于对位置的导数的项离散；在c_x的控制容积内，用代替。对c₁, c₂, …, c_n-1的控制容积，扩散方程分别离散为：

       (2)

对c_l,1, c _l,2, …, c _l,m-2的控制容积，扩散方程分别离散为：

 (3)

x=l处的边界条件离散为

             (4)

扩散方程在c_l,m-1的控制容积离散为

   (5)

结合右边界条件离散式，式(5)化为

        (6)

2.2  稳态假设

由于随时间变化的浓度c_x所在的位置x位于(0, l)内，即0＜x＜l。x左边浓度c₁的位置为：(n≥2)。x右边浓度c_l,1的位置为：(m≥2)。因此，只要m≥2，n≥2，浓度c₁和c_l,1所在的位置也都位于(0, l)内。位置相邻的c₁，c_x，c_l,1作为一个环节串联在整体之中，对它们可以进行稳态假设。假设：

，即；

，即；

，即。

2.3  c₁，c_x和c_l,1满足的关系式

由式(1)变形得：

        (7)

将c₁对时间求导，并结合假设得：

     (8)

由式(2)中第1式变形得：

          (9)

将和c₁的表达式代入c₂，整理得：

        (10)

将c₂对时间求导，得：

    (11)

由式(2)中第2式变形得：

        (12)

将，c₂和c₁的表达式代入c₃，整理得：

       (13)

综合c₂和c₃表达式的形式，令：

        (14)

则：

  (15)

将对时间求导，并结合假设得：

 (16)

由式(2)中第(i-1)式变形得：

       (17)

将，和的表达式代入，合并同类项后得：

 (18)

将此式与式(14)比照，得A_i，B_i和C_i的递推公式   (i≥4)为：

，，；

，，；

，，。

由A_i，B_i和C_i的递推公式及相应前2项可推得：；；。

在式(14)中，令i=n，得c_n的表达式；因c_n为左边界端点，令c_n=c₀，得到c_x，c_l,1满足的关系式：

 (19)

同理可得：

其中：；； (4≤i≤m-1)。

令i=m-2得：

    (20)

令i=m-1得：

     (21)

将式(21)对时间t求导，得：

 (22)

将，和的表达式代入，整理得c_x，c₁满足的关系式：

            (23)

2.4  关于c_x的微分方程

将c₁，的表达式(式(7)和式(8))代入式(23)，整理之后得到：

    (24)

将式(24)两边同时对时间求导，并结合假设得：

              (25)

将式(25)代入式(19)得：

          (26)

将等式代入式(9)得：

       (27)

将式(26) +式(27)，并整理之后得：

         (28)

在此步骤中，只有n和m趋于无穷大，每一个控制容积才会趋近于0，每个间隔位置处的浓度和浓度变化率分别与控制容积内的平均浓度、平均浓度变化率才相等。对式(28)两边同时取极限：

即

     (29)

经整理得到：

       (30)

2.5  c_x与时间t的函数关系式

式(30)为关于时间t的一阶线性微分方程，它的通解是：

         (31)

加入初始条件(t=0时，c_x=0)限制，确定积分常数C后，最终得：

         (32)

式(32)可看作前述非稳态扩散问题的一个近似解。

3  所求近似解和数值解的对比

对前述非稳态扩散问题，取c₀=10 mol/m³，l=100 m，D=0.25 m²/s；按照通常的数值方法^[19]获得它的数值解，计算数值解时，均匀网格数取1 000 (即△x=0.1 m)，时间步长(△t)取1 s。若网格再加密，时间步长再缩短，数值解均不变，则该数值解就是此例非稳态问题的精确解。非稳态问题的数值解和按式(14)求得的近似解对比如图2所示。

从图2可以看出：两者吻合程度较好。在扩散前期(0~9 000 s)，两者曲线形状吻合稍差，但9 000 s后，两者曲线形状吻合较好。从扩散开始到最终达到稳态的过程中(图2所示为算例条件下的结果，如选实际物质固态金属中的扩散来计算，扩散时间更长)，近似解和数值解是同步的，所以，用此近似解预测达到最终稳态的时间，准确度较高；从近似解和数值解对比来看，两者没有出现错位。在图2中，两者总有2~3个交点，即使经过漫长的扩散发展过程，同一时刻的近似解和数值解也未出现完全的错位。总体来说，近似解是可以接受的。扩散方程分析解的稳态近似法有效。

若在获得前述非稳态问题近似解的过程中，稳态假设改为c₁=const，即；c_x=const，即；，即。按照前述的方法和步骤，最终可得到近似解为：c_x=c₀，这说明在各时刻c_x的近似值就是最终稳态的值。这个近似解和实际情形差别较大。但工程上得到应用的金属氧化模型、未反应核模型等都是在假设产物层氧化或还原气体浓度c=const ()时得到的，所以，c_x 也被称为近似解。

浓度不随时间变化稳态假设下的近似解就是最终稳态的解，意味着稳态近似法处理的结果可能首先从近似最终稳态开始，然后近似接近最终稳态的状态，再近似远离最终稳态的状态。这和图2中出现的近似解和数值解吻合现象比较符合。

另外，从图2还可以看出：近似解基本上是遵循总体质量守恒的，也就是说，同一时刻近似解曲线和数值解曲线下的面积基本相同，即近似解曲线位于数值解曲线上的面积与它位于数值解曲线下的面积基本相同。但在扩散前期，两者在50，3 000和6 000 s时有大差别。这是由稳态近似法的特点决定的，稳态近似法处理的结果首先从近似最终稳态开始。

图2  数值解和近似解的对比

Fig. 2 Comparison of numerical solutions and approximate solutions

浓度随时间不变的假设应看作稳态假设。在浓度变化率随时间不变的假设中，浓度随时间是变化的，但它也是稳态假设。因为它是在处理非稳态问题中施加的相应参量随时间不变的约束。一般地，稳态假设不能同时应用在全部环节；稳态假设没有应用的环节还是呈非稳态变化。应用稳态假设并综合全部环节而获得的结果也是随时间变化。由浓度随时间不变的稳态假设获得的结果(例如未反应核模型)中，浓度不是随时间不变的；由浓度变化率随时间不变的稳态假设获得的结果(例如，本文结果)中，浓度变化率也不是常数。

4  结论

(1) 首先将动力学中常用的浓度随时间不变的稳态假设发展为浓度变化率随时间不变的稳态假设。给出了在有限长度区间内应用扩散方程进行稳态近似法处理的过程。

(2) 对文中给出的非稳态问题获得的近似解是：。对此近似解与精确数值解对比得出：从扩散开始一直到最终稳态的过程中，稳态近似法获得的结果和精确解是同步的；稳态近似法获得的结果与接近最终稳态的情形吻合程度好，与远离最终稳态的情形吻合程度较差；稳态近似法获得的结果基本上满足总体质量守恒。

参考文献：

[1] 黄希祜. 钢铁冶金原理[M]. 北京: 冶金工业出版社, 2005: 67-69, 83-87, 292-297, 420-423, 430-431.

HUANG Xihu. Principle of ferrous metallurgy[M]. Beijing: Metallurgical Industry Press, 2005: 67-69, 83-87, 292-297, 420-423, 430-431

[2] 韩其勇. 冶金过程动力学[M]. 北京: 冶金工业出版社, 1983: 24-28, 58-65, 156-160.

HAN Qiyong. Kinetics of process metallurgy[M]. Beijing: Metallurgical Industry Press, 1983: 24-28, 58-65, 156-160.

[3] 沈颐身, 李保卫, 吴懋林. 冶金传输原理基础[M]. 北京: 冶金工业出版社, 2003: 200-205, 405-407.

SHEN Yisheng, LI Baowei, WU Maolin. Principle of heat and mass transfer in metallurgy[M]. Beijing: Metallurgical Industry Press, 2003: 200-205, 405-407

[4] Coulson J M, Richardson J F, Backhurst J R, et al. Chemical engineering, 1B, Heat transfer and mass transfer[M]. 6th ed. Dalian: Dalian University of Technology Press, 2008: 395-397, 602-605

[5] Bonalde A, Henriquez A, Manrique M. Kinetic analysis of the iron oxide reduction using hydrogen-carbon monoxide mixtures as reducing agent [J]. ISIJ International, 2005, 45: 1255-1260.

[6] 肖兴国, 谢蕴国. 冶金反应工程学基础[M]. 北京: 冶金工业出版社, 1997: 31-34, 143-157.

XIAO Xingguo, XIE Yunguo. Metallurgical reaction engineering[M]. Beijing: Metallurgical Industry Press, 1997: 31-34, 143-157

[7] Sohn H Y, Perez-fontes S. Application of the law of additive reaction times to fluid–solid reactions in porous pellets with changing effective diffusivity[J]. Metallurgical and Materials Transactions B, 2010, 41: 1261-1267.

[8] Sohn H Y. The influence of chemical equilibrium on fluid-solid reaction rates and the falsification of activation energy[J]. Metallurgical and Materials Transactions B, 2004, 35: 121-131.

[9] Sohn H Y, Wadsworth M E. Rate processes of extractive metallurgy[M]. New York: Plenum Press, 1979: 1-51.

[10] Sohn H Y. Author’s reply[J]. Metallurgical and Materials Transactions B, 2005, 36: 897-901.

[11] Sohn H Y. The effects of reactant starvation and mass transfer in the rate measurement of fluid–solid reactions with small equilibrium constants[J]. Chemical Engineering Science, 2004, 59: 4361-4368.

[12] Wu S, Xu J, Yang S, et al. Basic characteristics of the shaft furnace of COREX smelting reduction process based on iron oxides reduction simulation[J]. ISIJ International, 2010, 50: 1032-1039.

[13] Yang K, Choi S, Chung J, et al. Numerical modeling of reaction and flow characteristics in a blast furnace with consideration of layered burden[J]. ISIJ International, 2010, 50: 972-980.

[14] Dong X F, Yu A B, Yagi J, et al. Modelling of multiphase flow in a blast furnace: recent developments and future work[J]. ISIJ International, 2007, 47: 1553-1570.

[15] Chuang H C, Kuo J H, Huang C C, et al. Multi-phase flow simulations in direct iron ore smelting reduction process[J]. ISIJ International, 2006, 46: 1158-1164.

[16] Chu M S, Nogami H, Yagi J. Numerical analysis on blast furnace performance under operation with top gas recycling and carbon composite agglomerates charging[J]. ISIJ International, 2004, 44: 2159-2167.

[17] 王淑兰. 物理化学[M]. 3版. 北京: 冶金工业出版社, 2008: 219-220.

WANG Shulan. Physical chemistry[M]. 3rd ed. Beijing: Metallurgical Industry Press, 2008: 219-220.

[18] 梁英教. 物理化学[M]. 北京: 冶金工业出版社, 2005: 272-274.

LIANG Yingjiao. Physical chemistry[M]. Beijing: Metallurgical Industry Press, 2005: 272-274.

[19] 陶文铨. 数值传热学[M]. 西安: 西安交通大学出版社, 1995: 79-83, 190-191.

TAO Wenquan. Numerical heat transfer[M]. Xi’an: Xi’an Jiaotong University Press, 1995: 79-83, 190-191.

(编辑  赵俊)

收稿日期：2012-08-31；修回日期：2012-12-06

基金项目：国家自然科学基金资助项目(50704040)；重庆市自然科学基金资助项目(CSTC 2009BB4197)

通信作者：王成善(1972-)，男，江苏赣榆人，博士，副教授，从事冶金传输和动力学、冶金工艺模拟等研究；电话：13983463604；E-mail: wangchengshan@sina.com