恒餾出液組成常規(guī)間歇精餾過程模擬的數(shù)值計算誤差

曾 真，王為國，羅旌崧，覃遠航，王存文，趙子傲

1.武漢工程大學機電工程學院，湖北 武漢 430205；2.武漢工程大學化工與制藥學院，教育部綠色化學過程教育部重點實驗室（武漢工程大學）；湖北省新型反應器與綠色化學工藝重點實驗室（武漢工程大學），湖北 武漢 430205

模擬是研究間歇精餾過程的重要方法。1902年，Rayleigh首先將數(shù)學引入間歇精餾過程的模擬研究，并提出了著名的Rayleigh方程。Rose等［1］，Huckaba等［2］，Distefano［3-4］，Gallun等［5］，Sadotomo和Miyahara［6］，Cho和Joseph［7］等均對間歇精餾過程的物理模型和計算方法進行了研究，并得到有價值的結果。Galindez等［8］提出了模擬間歇精餾過程的擬穩(wěn)態(tài)法。在化學工程單元操作教程中［9-10］，二元間歇精餾過程的物系及過程條件均被理想化，即，其物理模型簡化為理想模型。近年來，課題組基于理想模型的模擬（計算）結果，對常規(guī)［11-13］、全回流［14-16］和提餾［17-18］等間歇精餾過程的基礎理論進行了研究。模擬實際間歇精餾過程時，常用實驗結果來檢驗模擬結果的合理性及優(yōu)劣。模擬結果與實驗結果之間的偏差來源于物理模型誤差和數(shù)值計算誤差，基于理想模型的模擬（計算）結果，難以用實驗方法驗證其優(yōu)劣。非穩(wěn)態(tài)間歇精餾過程離散后進行數(shù)值計算，截斷誤差和舍入誤差是其固有誤差。擬穩(wěn)態(tài)法模擬理想條件下恒餾出物組成（恒頂濃）操作常規(guī)二元間歇精餾過程，計算過程總汽化量的程序為嵌套控制結構，即總汽化量的數(shù)值結果中含有中間變量計算誤差的傳播誤差，即，其數(shù)值計算結果的誤差為截斷誤差、舍入誤差和傳播誤差的耦合。由于舍入誤差和傳播誤差是隨機的，對于給定的分離任務，當理論板數(shù)有限時，根據(jù)大數(shù)定律［19］，通過對截斷誤差可以忽略且足夠多不同步長（或過程離散段數(shù)）的數(shù)值計算結果進行算術平均，可計算出與總汽化量真實值充分接近的數(shù)值，從而近似計算出總汽化量數(shù)值計算結果誤差。顯然，要近似計算出誤差，必須先找出什么條件下截斷誤差可以忽略。Richardson外推法［20］是估計數(shù)值計算結果截斷誤差的常用方法，當隨機的舍入誤差和傳播誤差在誤差中顯著呈現(xiàn)時，直接用Richardson外推法估計總汽化量數(shù)值結果的截斷誤差，估計結果的穩(wěn)定性較差，表明估計結果的準確性較差。本文將大數(shù)定律與Richardson外推法有機結合，提出通過擬合多組總汽化量的數(shù)值計算結果與對應的過程離散數(shù)，計算Richardson外推法誤差系數(shù)的擬合法，有效地降低了（或消除）隨機誤差（舍入誤差和傳播誤差）對計算誤差系數(shù)的影響，提高了Richardson外推法估計總汽化量數(shù)值結果截斷誤差的準確（穩(wěn)定）性，對完善間歇精餾過程的模擬技術有意義。

1 的數(shù)值計算方法

1.1 基本方程

假定常規(guī)二元間歇精餾過程的理想條件如下：（1）二元混合物是理想的混合物，相對揮發(fā)度視為常數(shù)。（2）塔內(nèi)汽相和液相流動為恒摩爾溢流。（3）每一板塊均為理論板。（4）再沸器相當于一板理論塔板數(shù)。（5）忽略理論板上的持液量和持汽量。（6）塔頂冷凝器中的汽相完全冷凝。（7）忽略塔頂?shù)某忠毫亢统制俊＠硐霔l件下，常規(guī)二元間歇精餾過程的基本公式如下：

（1）相平衡方程

式（1）中，x為液相中輕組分摩爾分數(shù)，y為汽相中輕組分摩爾分數(shù)，α為相對揮發(fā)度。

（2）微分物料衡算方程

式（2）和式（3）中，nD為瞬間的餾出物總量（kmol），nW為再沸器內(nèi)瞬時液相量（kmol），xD為餾出物中輕組分的瞬時摩爾分數(shù)，xW為再沸器內(nèi)液相中輕組分的瞬時摩爾分數(shù)。

（3）操作線方程

式（4）中，j為個理論塔板序號，右下角標j為第j理論塔板的物理量，R為回流比。

（4）一批物料的物料衡算方程

在式（5）-式（7）之中，nF為一批物料的原料總量（kmol），nDe為一批物料的塔頂產(chǎn)品量（kmol），nWe為一批物料的塔底產(chǎn)品量（kmol），xDP為塔頂產(chǎn)品量中輕組分的摩爾分數(shù)，xWe為塔底的產(chǎn)品量中輕組分的摩爾分數(shù)。η為塔頂產(chǎn)品量中輕組分的收率。

（5）一批物料的總汽化量

式（8）中，nVT為一批物料的總汽化量（kmol）

1.2 離散方程

常規(guī)二元間歇精餾過程是非穩(wěn)態(tài)過程。選擇xW為自變量，將由xF變?yōu)閤We的xW等段長離散為k段得離散方程如下：

式（9）-式（11）中，k數(shù)值計算時過程的離散段數(shù)，ΔxW為xW的離散步長，i為離散段序號，右下角標i為第i離散段的物理量，上標“—”為一個離散段物理量的平均。

對每一離散段物料衡算得

式（12）和式（13）中，nD,i為從開始到第i個離散段結束的餾出物總量（kmol），ΔnD,i為第i個離散段的餾出液量（kmol）。對于恒餾出物組成操作，式（12）中的=xDP(i=1,2,…,k)。

由式（8）得

1.3 數(shù)值計算的編程步驟

擬穩(wěn)態(tài)法模擬理想條件下恒餾出物組成操作常規(guī)二元間歇精餾過程，數(shù)值計算的編程步驟如下：

Step1：輸入：α、nF、xF、η、xDP、N+1、k和e；

Step3：由 式（9）、（10）和（11）計 算ΔxW，xW,i(i=0,1,…,k)和(i=1,2,…,k)；

Step4：由子程序計算Ri(i=1,2,…,k)；

Step7：輸出計算結果。

主程序中N+1為理論塔板數(shù)（包括再沸器），e為數(shù)值計算時Ri(i=1,2,…k)規(guī)定的收斂精度。

本文將由上述步驟編制的計算程序稱為主程序。主程序步驟Step4是使用子程序計算Ri。在子程序中，Ri是采用二分法［21］迭代計算。子程序編程步驟如下：

Step1：由主程序輸入：α、N+1、i、xDP和e；

Step2：t=0，賦Ri的初始搜索區(qū)間[R(0)i,a,R(0)i,b]；

①R(t)i達到收斂精度（ε(t)i≤e），轉到Step6；

②R(it)未達到收斂精度（ε(it)＞e），當時，和，當時，和，轉到Step3；

Step6：Ri=R(t)i，回主程序。

子程序中，t為數(shù)值計算過程中子程序的迭代次數(shù)，右上角標(t)為數(shù)值計算過程中子程序第t次迭代的物理量，εi為Ri(i=1,2,…k)數(shù)值計算結果的相對殘差。

1.4 計算程序驗證

由文獻［22］可知，在理想條件下，當N+1→∞時可由式（15）解析計算。

式（15）中，右上角標(∞)為N+1→∞時的物理量。

對于給定的分離任務，隨著N+1增大，向漸近。為了驗證計算程序，給定的分離任務Ⅰ（α=2.500 000 ，xF=0.400000 0，η=0.900 000 0，xDP=0.960 000 0）和分離任務Ⅱ（α=2.500 000 ，xF=0.700 000 0，η=0.800 000 0，xDP=0.9000000）。對于分離任務Ⅰ，由式（5）、（6）和（7）得xWe=0.064 000 0，由 式（15）得=1.486 198［右上角標(a)為解析計算結果］。對于分離任務Ⅱ，由式（5）、（6）和（7）得xWe=0.370 5882，由式（15）得1.013255。

圖1顯示，對于2個分離任務，隨著N+1增大［右上角標(n)為數(shù)值計算結果］均向=1.013255漸近，表明的計算程序是正確的。本文中稱上述驗證計算程序的方法為漸近法。

圖1 理論塔板數(shù)對總汽化量的影響：（a）分離任務Ⅰ，（b）分離任務ⅡFig.1 Effect of number of theoretical plates on total vaporization：（a）separation task I，（b）separation task II

1.5 的數(shù)值計算方法

文獻［12］中圖1顯示，對于常規(guī)二元間歇精餾過程，當N+1→∞時，x(∞)D可表示為：

在 式（16）中，R(∞)C=1/[(α-1)xW]，為N+1→∞時，塔頂和塔底同時為恒濃度區(qū)的回流比。

在恒餾出物組成操作常規(guī)二元間歇精餾過程中，xD=xDP≤1，R≤R(∞)C。因此，當N+1→∞時，其瞬時回流比為

過程離散化后，由式（17）得到每一離散段的回流比為：

圖2 過程的離散段數(shù)與總汽化量數(shù)值計算結果的關系：（a）分離任務Ⅰ，（b）分離任務ⅡFig.2 Relationship between number of discrete segments of process and numerical calculation results of total vaporization：（a）separation task I，（b）separation task II

2 [](n)的計算誤差分析

非穩(wěn)態(tài)間歇精餾過程離散后進行數(shù)值計算，截斷誤差Et和舍入誤差Er是其固有誤差。擬穩(wěn)態(tài)法模擬理想條件下恒餾出物組成操作常規(guī)二元間歇精餾過程，數(shù)值計算的程序為嵌套控制結構，在的誤差E中，應包括來自R(n)i(i=1,2,…k)計算誤差的傳播誤差Ep，則

式（19）中，E為的 誤 差，Er為的舍入誤差，Et為的截斷誤差，Ep為來源于R(n)i(i=1,2,…k)數(shù)值計算誤差傳播引起的計算誤差（傳播誤差）。

2.1 [](n)的誤差E

對于一定分離任務，當N+1→∞時，可由式（15）解析計算，誤差E可由式（19）計算。當N+1有限時，不能準確計算出的真實值，因此，不能直接由式（19）計算準確計算出的誤差E(n)。由于Er與Ep是隨機的，當截斷誤差Et很小，可以忽略時，根據(jù)大數(shù)定律［19］，通過對足夠多不同步長（或過程離散段數(shù)）的進行算術平均，可數(shù)值計算出與充分接近的數(shù)值從而近似計算出與總汽化量數(shù)值計算結果的E(n)真實值充分接近的數(shù)值。舉例驗證和應用如下：

（1）對于分離任務I，N+1→∞時，由圖2（a）可見，k≥1 000，圍繞波動，表示截斷誤差很小，可以忽略。由圖4（a）可見，C(∞),(f)k=0.389 4［右上角標(f)為擬合計算結 果］。k=1 000時，由 式（24）可 得，E(∞)t≈E(∞),(f)t=3.894×10-7。k=1 000-2 000時 ，，本文使用的計算軟件是32 bit系統(tǒng)，數(shù)值計算結果為7位有效數(shù)），當N+1→∞和k=100時，數(shù)值計算 得，=1.486159，則，由 式（19）得則 其 相 對 誤 差

（2）對于分離任務Ⅱ，N+1→∞時，由圖2（b）可見，k≥1 000，圍繞波動，表示截斷誤差很小，可以忽略。由圖4（b）可見，C(∞),(f)k=-0.111 1。k=1 000時，由式（24）可得，E(∞)t≈E(∞),(f)t=-1.111×10-7。由 圖2（b）可見 ，k=1 000-2 000時 ，1.013255=。當N+1→∞和k=100時，=1.013266，則，由 式（19）得=-1.1×10-5。則 其 相 對 誤 差

（3）對于分離任務I，當N+1=13和e=10-7時，由圖3（a）可見，k≥1000，圍繞波動，表示截斷誤差很小，可以忽略。由圖5（a）可見，C(f)k=0.3930。則，k=1000時，由 式（24）可 得，E(n)t≈E(f)t=3.930×10-7。k=1000-2 000 時 ，1.496 075。當N+1=13、e=10-7和k=100時，數(shù)值得，=1.496 035，由式（19）得，E(n)≈=4.0×10-5。相對誤差2.7×10-5。

（4）對于分離任務Ⅱ，當N+1=8和e=10-7時，由圖3（b）可見，k≥1000，圍繞波動，表示截斷誤差很小，可以忽略。由圖5（b）可見，C(f)k=-0.112 6。則，k=1000時，由式（24）可得，E(n)t≈E(f)t=-1.126×10-7。k=1000-2000時，1.021 851。當N+1=8、e=10-7和k=100時，數(shù)值計算得，=1.021862，由式（19）得，E(n)≈=-1.1×10-5。則其相對誤差

圖3 總汽化量數(shù)值計算結果與過程離散段數(shù)的關系（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）為應用，因此，在截斷誤差可以忽略的條件下，可根據(jù)大數(shù)定律，對總汽化量數(shù)值計算結果的誤差進行定量描述，可為“模擬間歇精餾過程的數(shù)值計算誤差非常小”提供的定量證據(jù)。

2.2 [n VT n F](n)的傳播誤差Ep

由式（14）可見，R(n)i(i=1,2,…k)的計算誤差可 直 接 傳 播 到。 當εi=0(i=1,2,……,k)時，Ep=0。設R(ε=0)i(i=1,2,…,k)和為 對 應 于εi=0(i=1,2,……,k)的數(shù)值計算結果。由式（14）和（19）可得

式（21）代入（20）得

2.3 [](n)的舍入誤差E r

當k不是非常大時，Er很小，可以忽略。在數(shù)學上，舍入誤差量化估計是仍然需要研究的問題，基于本文數(shù)值計算結果，和的舍入誤差量化估計舉例如下：

（1）當N+1→∞時，R(∞)i由式（18）解析計算，因 此，的 總 計 算 誤 差E(∞),(n)中，E(∞),(n)p=0。對于分離任務I，由圖2（a）可見，k=1000-2 000范圍內(nèi)，k=1 900時，|E(∞),(n)|最大。當k=1 900時，數(shù) 值 計 算 得1.486 201，由式（19）得，|E(∞),(n)|=3×10-6。由圖4（a）可 見，=0.389 4，由 式（24）可 得，由 式 （19） 可 得。假定舍入誤差正比于離散段數(shù)，則，k=100時， ||E(∞),(n)rmax在10-7數(shù)量級。

（2）對于分離任務I，當N+1=13和e=10-7時，k=1000-2 000范圍內(nèi)，k=1000時，||最大。當k=1000時，數(shù)值計算得1.496 072，由式（19）得。由圖5（a）可 見 ，C(f)k=0.3930，由 式（24）可 得 ，=3.930×10-7，由式（22）得，|Ep|max=1.121×10-7，由式（19）可得≈4×10-6。假定舍入誤差正比于離散段數(shù)，則，k=100時在10-7數(shù)量級。

2.4 [n VT n F](n)的截斷誤差Et

由圖2和圖3可見，對于兩個分離任務，當k較小時和是的單調(diào)函數(shù)，表明截斷誤差Et在計算誤差E中占主導地位，當k較大時和不再是的單調(diào)函數(shù)，表明截斷誤差不再占主導地位，隨機的舍入誤差和傳播誤差在誤差中顯著呈現(xiàn)，當k很大時，圍繞波動圍繞波動，表示截斷誤差不顯著呈現(xiàn)。

2.4.1 擬合法計算Richardson外推法誤差系數(shù)的提出Richardson外推法是一種常用的截斷誤差估計方法［20］，Richardson外推法是利用Taylor級數(shù)進行分析，用于截斷誤差估計時，在漸近線范圍內(nèi)常常只保留冪級數(shù)的首項。即

在式（23）中，C是與ΔxW無關的常數(shù)，O是高階無窮小，p為計算方法的收斂精度階，是Taylor級數(shù)首項的冪，一般為正整數(shù)。離散步長與離散段數(shù)成反比，因此，可以將式（23）變?yōu)?/p>

式（24）中，Ck為計算截斷誤差的系數(shù)（本文簡稱誤差系數(shù)）。由式（24）可見，略去式（24）中高階無窮小，如果確定了p和Ck，可計算出Et與k的數(shù)值關系。將式（24）代入式（19）得

對于兩個分離任務，當e=10-7和k=100時，由2.3可知在10-7數(shù)量級，由2.2可知，在10-7數(shù)量級。結合由圖2、圖3與式（25）可 得，當k=100時，|E|＞＞ |Ep|和 |E|＞＞ |Er|Et≈E。由式（24）與（25）可見，忽略高階無窮小Ok、舍入誤差和傳播誤差的影響，當N+1→∞時，可 由 式（15）解 析 計 算，由 兩 組和相應的k可確定p和C(∞)k。當N+1有限時，僅能計算出與真實值充分接近的數(shù)值在未知情況下，由三組和相應的k可確定p和Ck。

當隨機的舍入誤差和傳播誤差在誤差中顯著呈現(xiàn)，忽略舍入誤差和傳播誤差的影響，直接由式（25）計算誤差系數(shù)Ck，計算結果的穩(wěn)定較差，即，直接用Richardson外推法估計總汽化量數(shù)值結果的截斷誤差，估計結果的穩(wěn)定性較差。目前，在數(shù)學上，提高估計截斷誤差準確（穩(wěn)定）性的方法是保留Taylor級數(shù)更多項數(shù)（采用更多組不同步長的數(shù)值計算結果），課題組前期研究表明，計算誤差系數(shù)方法是解線性方程組，保留項數(shù)越多，計算過程越復雜，即保留項數(shù)是有限的，而有限項保留不能顯著地提高截斷誤差估計結果的穩(wěn)定性。根據(jù)的大數(shù)定律，本文認為擬合足夠多組和相應的k來計算Richardson外推法的Ck，可以有效地降低（或消除）Er和Ep對計算Ck的影響，可較為準確地估計Et。

2.4.2 擬合法計算Richardson外推法誤差系數(shù)Ck的驗證數(shù)值計算方法的收斂精度階p為Taylor級數(shù)首項的冪，一般為整數(shù)。在擬合和對應的k前，對p進行整數(shù)預設，擬合過程和結果可在Microsoft Excel電子表格上呈現(xiàn)。如果預設的p正確，擬合估計法合理，在 |E|＞＞ |Ep|和條件下，由式（25）可見與k-p必定呈良好的線性關系（Microsoft Excel電子表格擬合的線性相關系數(shù)接近1），由式（19）和（24）可見必定成立。

由圖4可見，對于2個分離任務，當k=40-150時，采用Microsoft Excel電子表格擬合和k-2的線性相關系數(shù)均接近1，表明擬穩(wěn)態(tài)法數(shù)值計算具有二階收斂精度。

對于兩個分離任務，E(∞),(n)p=0，由2.3可知，k=100時在10-7數(shù)量級。對于分離任務I，當N+1→∞和k=100時，由圖2（a）可見，E(∞),(n)=3.9×10-5，即 ，和因此，≈3.9×10-5。由圖 4（a）可 見 ，C(∞),(f)k=0.389 4，則 ，。對 于 分 離 任 務Ⅱ，當N+1→∞和k=100時，由圖2（b）可見，E(∞),(n)=-1.1×10-5，即和因 此≈-1.1×10-5。由圖4（b）可見，C(∞),(f)k=0.111 1，則1.01。因此，當N+1→∞時，擬合法計算Richardson外推法的誤差系數(shù)C(∞),(f)k，可較為準確地定量估計E(∞)t。

圖4 采用Microsoft Excel電子表格擬合總汽化量數(shù)值計算結果與過程離散段數(shù)：（a）分離任務Ⅰ，（b）分離任務ⅡFig.4 Fitting of numerical calculation results of total vaporization and number of discrete segments of process by using Microsoft Excel spreadsheet：（a）separation task I，（b）separation task II

對于兩個分離任務，當e=10-7和k=100時，由2.3可知在10-7數(shù)量級，由2.2可知，在10-7數(shù)量級。對于分離任務I，當N+1=13、e=10-7和k=100時，10-5，即，和，因 此，E(n)t≈4.0×10-5。由圖5（a）可見，當N+1=13時，C(f)k=0.3930，C(f)kk-2=3.930×10-5，即，≈0.98。對于分離任務Ⅱ，當N+1=8、e=10-7和k=100時，E(n)≈=-1.1×10-5，即 ，和，因 此 ，E(n)t≈-1.1×10-5。由圖5（b）可見，當N+1=8時，C(f)k=-0.112 6，C(f)k k-2=-1.126×10-5，即，≈1.02。因此，擬合法計算截斷誤差Richardson外推法的誤差系數(shù)Ck，可較為準確地定量估計Et。

此外，對于分離任務Ⅱ，由圖5（b）可見，N+1=8時C(f)k是負值，由圖5（c）可見，N+1=3時C(f)k是正值，表明Et的方向有可以隨N+1變化。

圖5 采用Microsoft Excel電子表格擬合總汽化量數(shù)值計算結果與過程離散段數(shù)示例：（a）分離任務Ⅰ，N+1=13，（b）分離任務Ⅱ，N+1=8，（c）分離任務Ⅱ，N+1=3Fig.5 Examples of fitting of numerical calculation results of total vaporization and number of discrete segments of process by using Microsoft Excel spreadsheet：（a）N+1=13 for separation taskⅠ，（b）N+1=8 for separation taskⅡ，（c）N+1=3 for separation taskⅡ

2.4.3 擬合法計算Richardson外推法誤差系數(shù)Ck的應用舉例對于恒餾出物組成操作常規(guī)二元間歇精餾，由Fenske方程［22］知，過程結束時所需的最少理論達到其最大值，此最大值為所需的最少理論塔板數(shù)，即

由式（26），對于分離任務Ⅰ，Nmin+1=6.260899，對于分離任務Ⅱ，Nmin+1=2.976041。

圖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與相同過程離散段數(shù)的E(∞)t相當。②當N+1≥3(Nmin+1)時，C(f)k圍繞C(∞),(f)k波動，表明有限組和k擬合，隨機的Er和Ep對擬合結果仍有影響。

圖6 理論塔板數(shù)對誤差系數(shù)的影響（a）分離任務Ⅰ，（b）分離任務ⅡFig.6 Effect of number of theoretical plates on error coefficient：（a）separation taskⅠ，（b）separation taskⅡ

3 結 論

（1）根據(jù)大數(shù)定律，提出了N+1有限時，誤差的近似計算方法，并舉例驗證和應用。誤差的近似計算方法為：在的截斷誤差可以忽略的條件下，通過對足夠多不同步長（或過程離散進行算術平均，計算出與真實值充分接近的數(shù)值再由式（19）計算出與誤差真實值充分接近的數(shù)值。

（2）根據(jù)大數(shù)定律，提出了采用擬合法計算Richardson外推法的誤差系數(shù)Ck。擬合法通過擬合足夠多組和對應的k計算Ck，有效降低了（或消除）隨機誤差（舍入誤差和傳播誤差）對計算Ck的影響，提高了Richardson外推法估計截斷誤差的準確（穩(wěn)定）性。通過舉例對擬合法進行驗證和應用，并得到結論如下：①擬穩(wěn)態(tài)法模擬理想條件下，恒餾出物組成操作常規(guī)二元間歇精餾過程，具有二階收斂精度。②擬合法計算Ck時，對p進行整數(shù)預設，擬合過程和結果可在Microsoft Excel電子表格上呈現(xiàn)。

（3）通過舉例，對模擬理想條件下恒餾出物組成操作常規(guī)二元間歇精餾過程的和的數(shù)值計算結果的截斷誤差進行討論，得到結論如下：當N+1≥2(Nmin+1)時，C(f)k≈C(∞),(f)k，表明當理論板較大（大于2倍所需的最少理論板）時，Et與相同過程離散段數(shù)的E(∞)t相當。