《武汉工程大学学报》 2021年02期
134-142
出版日期：2021-04-30
ISSN:1674-2869
CN:42-1779/TQ
恒馏出液组成常规间歇精馏过程模拟的数值计算误差
模拟是研究间歇精馏过程的重要方法。1902年，Rayleigh首先将数学引入间歇精馏过程的模拟研究，并提出了著名的Rayleigh方程。Rose 等［1］ ，Huckaba等［2］，Distefano［3-4］ ，Gallun等［5］，Sadotomo和Miyahara［6］，Cho和Joseph［7］等均对间歇精馏过程的物理模型和计算方法进行了研究，并得到有价值的结果。Galindez等［8］提出了模拟间歇精馏过程的拟稳态法。在化学工程单元操作教程中［9-10］，二元间歇精馏过程的物系及过程条件均被理想化，即，其物理模型简化为理想模型。近年来，课题组基于理想模型的模拟（计算）结果，对常规［11-13］、全回流［14-16］和提馏［17-18］等间歇精馏过程的基础理论进行了研究。模拟实际间歇精馏过程时，常用实验结果来检验模拟结果的合理性及优劣。模拟结果与实验结果之间的偏差来源于物理模型误差和数值计算误差，基于理想模型的模拟（计算）结果，难以用实验方法验证其优劣。非稳态间歇精馏过程离散后进行数值计算，截断误差和舍入误差是其固有误差。拟稳态法模拟理想条件下恒馏出物组成（恒顶浓）操作常规二元间歇精馏过程，计算过程总汽化量的程序为嵌套控制结构，即总汽化量的数值结果中含有中间变量计算误差的传播误差，即，其数值计算结果的误差为截断误差、舍入误差和传播误差的耦合。由于舍入误差和传播误差是随机的，对于给定的分离任务，当理论板数有限时，根据大数定律［19］，通过对截断误差可以忽略且足够多不同步长（或过程离散段数）的数值计算结果进行算术平均，可计算出与总汽化量真实值充分接近的数值，从而近似计算出总汽化量数值计算结果误差。显然，要近似计算出误差，必须先找出什么条件下截断误差可以忽略。Richardson外推法［20］是估计数值计算结果截断误差的常用方法，当随机的舍入误差和传播误差在误差中显著呈现时，直接用Richardson外推法估计总汽化量数值结果的截断误差，估计结果的稳定性较差，表明估计结果的准确性较差。本文将大数定律与Richardson外推法有机结合，提出通过拟合多组总汽化量的数值计算结果与对应的过程离散数，计算Richardson外推法误差系数的拟合法，有效地降低了（或消除）随机误差（舍入误差和传播误差）对计算误差系数的影响，提高了Richardson外推法估计总汽化量数值结果截断误差的准确（稳定）性，对完善间歇精馏过程的模拟技术有意义。1 [nVTnF]的数值计算方法 1.1 基本方程假定常规二元间歇精馏过程的理想条件如下：（1）二元混合物是理想的混合物，相对挥发度视为常数。（2）塔内汽相和液相流动为恒摩尔溢流。（3）每一板块均为理论板。（4）再沸器相当于一板理论塔板数。（5）忽略理论板上的持液量和持汽量。（6）塔顶冷凝器中的汽相完全冷凝。（7）忽略塔顶的持液量和持汽量。理想条件下，常规二元间歇精馏过程的基本公式如下：（1） 相平衡方程[y=αx1+(α-1)x] （1）式（1）中，[x]为液相中轻组分摩尔分数，[y]为汽相中轻组分摩尔分数，[α]为相对挥发度。（2） 微分物料衡算方程 [dnD+dnW=0] （2）[d(nWxW)+xDdnD=0] （3）式（2）和式（3）中，[nD]为瞬间的馏出物总量（[kmol]），[nW]为再沸器内瞬时液相量（[kmol]），[xD]为馏出物中轻组分的瞬时摩尔分数，[xW]为再沸器内液相中轻组分的瞬时摩尔分数。（3） 操作线方程[yj+1=RR+1xj+xDR+1] （4）式（4）中，[j]为个理论塔板序号，右下角标[j]为第j理论塔板的物理量，[R]为回流比。（4） 一批物料的物料衡算方程[nF=nDe+nWe] （5）[nFxF=nDexDP+nWexWe] （6）[η=nDexDP(nFxF)] （7）在式（5）-式（7）之中，[nF]为一批物料的原料总量（[kmol]），[nDe]为一批物料的塔顶产品量（[kmol]），[nWe]为一批物料的塔底产品量（[kmol]），[xDP]为塔顶产品量中轻组分的摩尔分数，[xWe]为塔底的产品量中轻组分的摩尔分数。[η]为塔顶产品量中轻组分的收率。 （5） 一批物料的总汽化量[nVT=0nDe(R+1)dnD] （8）式（8）中，[nVT]为一批物料的总汽化量（[kmol]）1.2 离散方程常规二元间歇精馏过程是非稳态过程。选择[xW]为自变量，将由[xF] 变为 [xWe]的[xW]等段长离散为[k]段得离散方程如下：[ΔxW=(xF-xWe)/k] （9）[xW,i=xF-i×ΔxW??(i=0, 1, ?, k)] （10）[xW,i=xW,i-1+xW,i2??(i=1, 2, ?, k)] （11）式（9）-式（11）中，[k]数值计算时过程的离散段数，[ΔxW]为[xW]的离散步长，[i]为离散段序号，右下角标[i]为第[i]离散段的物理量，上标“—”为一个离散段物理量的平均。对每一离散段物料衡算得[ΔnD,inF=(1-nD,i-1nF)ΔxWxD,i-xW,i ?(i=1, 2, ?, k)] （12）[nD,inF=m=1iΔnD,mnF??(i=1, 2, ?, k)] （13）式（12）和式（13）中，[nD,i] 为从开始到第i个离散段结束的馏出物总量（kmol），[ΔnD,i]为第i个离散段的馏出液量（kmol）。对于恒馏出物组成操作，式（12）中的[xD,i=xDP][(i=1, 2, ?, k)]。由式（8）得[nVTnF=i=1k(Ri+1)ΔnD,inF] （14）1.3 数值计算[nVTnF]的编程步骤拟稳态法模拟理想条件下恒馏出物组成操作常规二元间歇精馏过程， 数值计算[nVTnF]的编程步骤如下：Step1： 输入： [α]、[nF]、[xF]、[η]、[xDP]、[N+1]、[k]和 [e]；Step2：由式（5）、（6）和（7） 计算[nDenF] 和 [xWe]；Step3： 由式（9）、（10）和（11）计算[ΔxW]， [xW,i(i=0,1,?,k)] 和[xW,i(i=1,2,?,k)]；Step4：由子程序计算[Ri(i=1,2,?,k)]；Step5：由式（12）和（13）计算 [ΔnD,inF(i=1,2??,k)]和[nD,inF(i=1,2,?,k)]；Step6：由式（14） 计算[nVTnF]；Step7：输出计算结果。主程序中[N+1]为理论塔板数（包括再沸器），[e]为数值计算时[Ri(i=1,2,?k)]规定的收敛精度。本文将由上述步骤编制的计算程序称为主程序。主程序步骤Step4是使用子程序计算[Ri]。在子程序中，[Ri]是采用二分法［21］迭代计算。子程序编程步骤如下：Step1：由主程序输入：[α]、[N+1]、[i]、[xW,i]、[xDP] 和[e]；Step2： [t=0]，赋[Ri]的初始搜索区间 [[R(0)i,a,R(0)i,b]]；Step3：[t=t+1]，[R(t)i=(R(t-1)i,a+R(t-1)i,b)2]；Step4： 交替使用式（1）和（4），由[xD,i]（[xD,i=xDP]）计算 [xW,i(t)] （McCabe-Thiele法［9］）；Step5：计算 [ε(t)i=(R(t-1)i,b-R(t-1)i,a)R(t)i]， 判断[R(t)i]是否达到规定的收敛精度；① [R(t)i]达到收敛精度（[ε(t)i≤e]），转到Step6；②[R(t)i]未达到收敛精度（[ε(t)i>e]），当 [xW,i(t)>xW,i]时， [R(t)i,a=R(t)i]和[R(t)i,b=R(t-1)i,b]，当[xW,i(t)[5 10 15 20 25 30N+1][2.502.201.901.601.301.00][[nVT/nF](n)][[n(∞) VT/nF](a)=1.486 198][α=2.500 000xF=0.400 000 0xDP=0.960 000 0η=0.900 000 0k=100e=10-7][0 4 8 12 16 20N+1][2.201.601.200.800.40][[nVT/nF](n)][[n(∞) VT/nF](a)=1.013 255][α=2.500 000xF=0.700 000 0xDP=0.900 000 0η=0.800 000 0k=100e=10-7][ a ][ b ]图1 理论塔板数对总汽化量的影响：（a）分离任务Ⅰ，（b）分离任务ⅡFig. 1 Effect of number of theoretical plates on total vaporization： （a） separation task I， （b） separation task II1.5 [n(∞)VTnF]的数值计算方法文献［12］中图1显示，对于常规二元间歇精馏过程，当[N+1→∞]时，[x(∞)D]可表示为：[x(∞)D=R(yW-xW)+yW??R[α=2.500 000xF=0.400 000 0xDP=0.960 000 0η=0.900 000 0N+1=13e=10-7][[[nVT/nF]](n)=1.496 075(k=1 000~2 000)][ a ][ b ][1.021 8601.021 8571.021 8541.021 8511.021 8480.021 845][[nVT/nF](n)][[[nVT/nF]](n)=1.021 851(k=1 000~2 000)][0 500 1 000 1 500 2 000 2 500k][α=2.500 000xF=0.700 000 0xDP=0.900 000 0η=0.800 000 0N+1=8e=10-7]图3 总汽化量数值计算结果与过程离散段数的关系（a） 分离任务Ⅰ，（b）分离任务ⅡFig. 3 Relationship between numerical calculation results of total vaporization and number of discrete segments of process ： （a） separation task I， （b）separation task II例（1）和（2）为验证，例（3）和（4）为应用，因此，在截断误差可以忽略的条件下，可根据大数定律，对总汽化量数值计算结果的误差进行定量描述，可为“模拟间歇精馏过程的数值计算误差非常小”提供的定量证据。2.2 [[nVTnF](n)]的传播误差[Ep] 由式（14）可见，[R(n)i(i=1,2,?k)]的计算误差可直接传播到[[nVTnF](n)]。当[εi=0(i=1,2,??,k)]时，[Ep=0]。设[R(ε=0)i(i=1,2,?,k)]和[[nVTnF](ε=0)]为对应于[εi=0(i=1,2,??,k)]的数值计算结果。由式（14） 和 （19）可得[Ep=i=1k(R(ε=0)i-R(n)i)ΔnD, inF] （20）数值计算[nVTnF]的程序中，[Ri]由子程序采用二分法迭代计算，[R(t)i(i=1,2,??,k)]的收敛准则为[ε(t)i=(R(t-1)i,b-R(t-1)i,a)R(t)i≤e]。则[(R(ε=0)i-R(n)i)R(n)i≤(R(t-1)i,b-R(t-1)i,a)R(t)i≤e] （21）式（21）代入（20）得[E(n)pmax=e{[nVTnF](n)-nDenF}] （22）由式（22）可见，[E(n)pmax] 正比于[e]，因此，在本文在数值计算[nVTnF]的程序中，取[e=10-7]（本文使用的计算软件是32 bit系统，[e]的最小有效取值为[10-7]），结合图1可得，对于任务Ⅰ和任务 Ⅱ，[E(n)pmax]在[10-7]数量级。2.3 [[nVTnF](n)]的舍入误差 [Er]当[k]不是非常大时，[Er]很小，可以忽略。在数学上，舍入误差量化估计是仍然需要研究的问题，基于本文数值计算结果，[[nVTnF](n)]和[[n(∞)VTnF](n)]的舍入误差量化估计举例如下：（1）当[N+1→∞]时，[R(∞)i]由式（18）解析计算，因此，[[n(∞)VTnF](n)]的总计算误差[E(∞),(n)]中，[E(∞),(n)p=0]。对于分离任务 I，由图2（a）可见，[k=1 000-2 000]范围内，[k=1 900]时，[E(∞),(n)] 最大。当[k=1 900]时，数值计算得 [[n(∞)VTnF](n)=1.486 201]，由式（19） 得， [E(∞),(n)=3×10-6]。由图4（a）可见，[C(∞),(f)k=0.389 4]，由式（24）可得，[E(∞),(n)t≈E(∞),(f)t≈C(∞),(f)kk-2=1.079×10-7]，由式（19）可得[E(∞),(n)rmax=E(∞),(n)+E(∞),(n)t≈3×10-6]。假定舍入误差正比于离散段数，则，[k=100]时，[E(∞),(n)rmax]在[10-7]数量级。（2）对于分离任务I，当[N+1=13]和[e=10-7]时， [k=1 000-2 000]范围内，[k=1 000]时，[E(n)] 最大。当[k=1 000]时，数值计算得，[[nVTnF](n)=1.496 072]，由式（19） 得，[E(n)=3×10-6]。由图5（a）可见，[C(f)k=0.3930]，由式（24）可得，[E(n)t≈E(f)t≈C(f)kk-2=3.930×10-7]，由式（22） 得，[Epmax=1.121×10-7]，由式（19）可得[E(n)rmax=E(n)+E(n)t+Epmax≈4×10-6]。假定舍入误差正比于离散段数，则，[k=100]时，[E(n)rmax]在[10-7]数量级。2.4 [[nVTnF](n)]的截断误差 [Et] 由图2和图3可见，对于两个分离任务，当[k]较小时，[[n(∞)VTnF](n)]和[[nVTnF](n)]是的单调函数，表明截断误差[Et]在计算误差[E]中占主导地位，当[k]较大时，[[n(∞)VTnF](n)]和[[nVTnF](n)]不再是的单调函数，表明截断误差不再占主导地位，随机的舍入误差和传播误差在误差中显著呈现，当[k]很大时，[[n(∞)VTnF](n)]围绕[[n(∞)VTnF](a)]波动，[[nVTnF](n)]围绕[[nVTnF](n)]波动，表示截断误差不显著呈现。2.4.1 拟合法计算Richardson外推法误差系数的提出 Richardson外推法是一种常用的截断误差估计方法［20］，Richardson外推法是利用Taylor级数进行分析，用于截断误差估计时，在渐近线范围内常常只保留幂级数的首项。即[Et=CΔxWp+O] （23）在式（23）中，[C]是与[ΔxW]无关的常数，[O]是高阶无穷小，[p]为计算方法的收敛精度阶，是Taylor级数首项的幂，一般为正整数。离散步长与离散段数成反比，因此，可以将式（23）变为[Et=Ckk-p+Ok] （24）式（24）中，[Ck]为计算截断误差的系数（本文简称误差系数）。由式（24）可见，略去式（24）中高阶无穷小，如果确定了[p]和[Ck]，可计算出[Et]与[k]的数值关系。将式（24）代入式（19）得[E=[nVTnF](a)-[nVTnF](n)][=Ckk-p+Ok+Er+Ep] （25）对于两个分离任务，当[e=10-7]和[k=100]时，由2.3可知，[E(n)rmax]在[10-7]数量级，由2.2可知，[E(n)pmax]在[10-7]数量级。结合由图2、图3与式（25）可得，当[k=100]时，[E>>Ep]和[E>>Er] [Et≈E]。由式（24）与（25）可见，忽略高阶无穷小[Ok]、舍入误差和传播误差的影响，当[N+1→∞] 时，[[n(∞)VTnF](a)]可由式（15）解析计算，由两组[[n(∞)VTnF](n)]和相应的[k]可确定[p]和[C(∞)k]。当[N+1]有限时，仅能计算出与[nVTnF]真实值充分接近的数值[[nVTnF](n)]，在[[nVTnF](a)]未知情况下，由三组[[n(∞)VTnF](n)]和相应的[k]可确定[p]和[Ck]。当随机的舍入误差和传播误差在误差中显著呈现，忽略舍入误差和传播误差的影响，直接由式（25）计算误差系数[Ck]，计算结果的稳定较差，即，直接用Richardson外推法估计总汽化量数值结果的截断误差，估计结果的稳定性较差。目前，在数学上，提高估计截断误差准确（稳定）性的方法是保留Taylor级数更多项数（采用更多组不同步长的数值计算结果），课题组前期研究表明，计算误差系数方法是解线性方程组，保留项数越多，计算过程越复杂，即保留项数是有限的，而有限项保留不能显著地提高截断误差估计结果的稳定性。根据的大数定律，本文认为拟合足够多组[[nVTnF](n)]和相应的[k]来计算Richardson外推法的[Ck]，可以有效地降低（或消除）[Er]和[Ep]对计算[Ck]的影响，可较为准确地估计[Et]。2.4.2 拟合法计算Richardson外推法误差系数[Ck]的验证 数值计算方法的收敛精度阶[p]为Taylor级数首项的幂，一般为整数。在拟合[[nVTnF](n)]和对应的[k]前，对[p]进行整数预设，拟合过程和结果可在Microsoft Excel电子表格上呈现。如果预设的[p]正确，拟合估计法合理，在[E>>Ep]和[E>>Er]条件下，由式（25）可见，[[nVTnF](n)]与[k-p]必定呈良好的线性关系（Microsoft Excel电子表格拟合的线性相关系数接近1），由式（19）和（24）可见，[E(n)t≈E(n)≈E(f)t≈][C(f)kk-p]必定成立。由图4可见，对于2个分离任务，当[k=40-150]时，采用Microsoft Excel电子表格拟合[[n(∞)VTnF](n)]和[k-2]的线性相关系数均接近1，表明拟稳态法数值计算[n(∞)VTnF]具有二阶收敛精度。对于两个分离任务，[E(∞),(n)p=0]，由2.3可知，[k=100]时，[E(∞),(n)rmax]在[10-7]数量级。对于分离任务I，当[N+1→∞]和[k=100]时，由图2（a）可见，[E(∞),(n)=3.9×10-5]，即，[E(∞),(n)>>E(∞),(n)p]和[E(∞),(n)>>E(∞),(n)r]，因此，[E(∞),(n)t≈3.9×10-5]。由图4（a）可见，[C(∞),(f)k=0.389 4]，则，[C(∞),(f)kk-2E(∞),(n)t≈1.00]。对于分离任务Ⅱ，当[N+1→∞]和[k=100]时，由图2（b）可见，[E(∞),(n)=-1.1×10-5]，即，[E(∞),(n)>>E(∞),(n)p]和[E(∞),(n)>>][E(∞),(n)r]，因此，[E(∞),(n)t≈][-1.1×10-5]。由图4（b）可见，[C(∞),(f)k=0.111 1]，则，[C(∞),(f)kk-2E(∞),(n)t≈][1.01]。因此，当[N+1→∞]时，拟合法计算Richardson外推法的误差系数[C(∞),(f)k]，可较为准确地定量估计[E(∞)t]。 图5为采用Microsoft Excel电子表格拟合[[nVTnF](n)] 和 k示例。当[e=10-7]和[k=40-150]时，[[nVTnF](n)]与[k-2]具有良好的线性关系（Microsoft Excel电子表格拟合的线性相关系数大于0.99，本文未全部显示），表明拟稳态法数值计算[nVTnF]具有二阶收敛精度。对于两个分离任务，当[e=10-7]和[k=100]时，由2.3可知，[E(n)rmax]在[10-7]数量级，由2.2可知，[E(n)pmax]在[10-7]数量级。对于分离任务I，当[N+1=13]、[e=10-7]和[k=100]时，[E(n)≈E(n)=4.0×10-5]，即，[E(n)>>E(n)p]和[E(n)>>E(n)r]，因此，[E(n)t≈4.0×10-5]。由图5（a）可见，当[N+1=13]时，[C(f)k=0.393 0]，[C(f)kk-2=3.930×10-5]，即，[C(f)kk-2E(n)t≈0.98]。对于分离任务Ⅱ，当[N+1=8]、[e=10-7]和[k=100]时，[E(n)≈E(n)=-1.1×10-5]，即，[E(n)>>E(n)p]和[E(n)>>E(n)r]，因此，[E(n)t≈-1.1×10-5]。由图5（b）可见，当[N+1=8]时，[C(f)k=-0.112 6]，[C(f)kk-2=-1.126×10-5]，即，[C(f)kk-2E(n)t≈1.02]。因此，拟合法计算截断误差Richardson外推法的误差系数[Ck]，可较为准确地定量估计[Et]。此外，对于分离任务Ⅱ，由图5 （b）可见，[N+1=8]时[C(f)k]是负值，由图5 （c）可见，[N+1=3]时[C(f)k]是正值，表明[Et]的方向有可以随[N+1]变化。2.4.3 拟合法计算Richardson外推法误差系数[Ck]的应用举例 对于恒馏出物组成操作常规二元间歇精馏，由Fenske方程［22］知，过程结束时所需的最少理论达到其最大值，此最大值为所需的最少理论塔板数，即[Nmin+1=lnxDP(1-xWe)(1-xDP)xWelnα] （26）由式（26），对于分离任务 Ⅰ，[Nmin+1=6.260 899]，对于分离任务Ⅱ， [Nmin+1=2.976 041]。图6显示[N+1]对[C(f)k]的影响。对于分离任务Ⅰ，当[N+1=7]时，[C(f)k=3.321 2]（未示出），对于分离任务Ⅱ，当 [N+1=3]时， [C(f)k=42.049]，［未示出，见图5（c）］。由图6可见，对于两个分离任务，①当[N+1≥2(Nmin+1)]时，[C(f)k≈C(∞),(f)k]，表明当理论板较大（大于2倍所需的最少理论板）时，[Et]与相同过程离散段数的[E(∞)t]相当。②当[N+1≥3(Nmin+1)]时， [C(f)k]围绕[C(∞),(f)k]波动，表明有限组[[nVTnF](n)]和[k]拟合，随机的[Er]和[Ep]对拟合结果仍有影响。3 结 论（1）根据大数定律，提出了[N+1]有限时，[[nVTnF](n)]误差的近似计算方法，并举例验证和应用。[[nVTnF](n)]误差的近似计算方法为：在[[nVTnF](n)]的截断误差可以忽略的条件下，通过对足够多不同步长（或过程离散[[nVTnF](n)]进行算术平均，计算出与真实值[[nVTnF](a)]充分接近的数值[[nVTnF](n)]，再由式（19）计算出与误差真实值充分接近的数值。（2）根据大数定律，提出了采用拟合法计算Richardson外推法的误差系数[Ck]。拟合法通过拟合足够多组[[nVTnF](n)]和对应的[k]计算[Ck]，有效降低了（或消除）随机误差（舍入误差和传播误差）对计算[Ck]的影响，提高了Richardson外推法估计截断误差的准确（稳定）性。通过举例对拟合法进行验证和应用，并得到结论如下： ①拟稳态法模拟理想条件下，恒馏出物组成操作常规二元间歇精馏过程，具有二阶收敛精度。②拟合法计算[Ck]时，对[p]进行整数预设，拟合过程和结果可在Microsoft Excel电子表格上呈现。（3）通过举例，对模拟理想条件下恒馏出物组成操作常规二元间歇精馏过程的[nVTnF]和[n(∞)VTnF]的数值计算结果的截断误差进行讨论，得到结论如下：当[N+1≥2(Nmin+1)]时，[C(f)k≈C(∞),(f)k]，表明当理论板较大（大于2倍所需的最少理论板）时，[Et]与相同过程离散段数的[E(∞)t]相当。