一種多組分分離級聯的數值計算方法

叢藝坤

(核工業理化工程研究院，天津 300180)

穩定同位素被廣泛應用于基礎科學、地質學、生物學、醫療診斷以及半導體技術等領域。元素周期表中穩定的化學元素有83種，其中39種含有3種以上同位素組分，這類元素被稱為多組分同位素元素。已有許多關于多組分同位素分離的理論和實驗研究[1-7]。

離心法分離穩定同位素通常在各級流量相等的矩形級聯中進行。研究多組分矩形級聯的分離特性，必須求解描述級聯內部物質輸運規律的級聯方程。由于多組分矩形級聯方程是復雜的非線性方程組，只能用數值方法來求解，因此快速、穩定、精確地求解級聯方程的數值計算方法是級聯設計和優化的必要工具。關于多組分矩形級聯方程求解的計算方法，在許多文獻中進行了探討，其中準線性法[8-9]、q迭代法[10-11]以及非穩態法[12-13]比較快速、穩定。

本研究將提供一種基于迭代序列差的能穩定收斂的多組分矩形級聯的數值計算方法，并研究其準確性、穩定性和收斂特性。

1 級聯數學模型

1.1 級聯結構

本研究的級聯模型為各級流量相等的逆流型矩形級聯,結構示意圖示于圖1。級聯總級數為n，供料級為f，各級供料流量為G。級聯供、精、貧流量及第i組分的豐度分別用F、P、W，CFi，CWi，表示。F、P、W三個參量之間存在守恒關系，給定其中兩個量，第三個量能唯一確定。

圖1 逆流型矩形級聯示意圖Fig.1 Countercurrent square cascade

1.2 級的守恒及分離方程

假設第i組分在級聯第s級的供料、精料和貧料流中的豐度分別為Cfi,s、Cpi,s、Cwi,s，第s級的分流比為θs，根據級聯分離理論，存在以下組分物質守恒方程

θsCpi,s+(1-θs)Cwi,s，i=1，……，m

(1)

m為同位素混合物的組分數。

第s級的分離方程可表示為

(2)

γij為第i組分對第j組分的相對全分離系數。在通常情況下，相對全分離系數隨級的供料豐度和分流比的變化而改變，但是對于大部分基于分子動力學的分離方法，相對全分離系數不隨上述參量變化，并可表示為

γij=γ0Mj-Mi

(3)

γ0為基本分離系數，Mi為第i組分的摩爾質量。

另外，各組分豐度之和滿足歸一化條件，即

∑Cpi,s=1，i=1，……，m

(4)

對級聯任意s級，在Cpi,s、γ0、Mi、θs已知的條件下，要求解Cpi,s、Cwi,s，存在m個守恒方程(1)(因為i=1，……，m)，m-1個分離方程(2)(因為在i=j時，γij=1，沒有分離效應)，1個豐度歸一化方程(4)，共2m個方程。作為未知量，有m個Cpi,s和m個Cwi,s，即共2m個未知量。因此利用方程(1)～(4)可求出Cpi,s和Cwi,s。

以消去Cwi,s為目的，把方程(1)代入(2)，化簡可得

Cfj,s·Cpi,s-γijCfi,s·Cpi,s+θs(γij-1)·Cpi,s·Cpj,s=0，i=1，……，m，j=1，……，m

(5)

由(5)式求出Cpi,s之后，再利用(1)式可求出Cwi,s。

1.3 級的分流比及供料豐度的確定

要求解非線性方程組(5)，必須知道級聯各級的分流比和供料豐度，各級分流比由具體的分離任務和分離工藝確定，一般是給定某一級的分流比(比如級聯精料級分流比θn)，其他各級的分流比θs由第s級到級聯端部間的物質守恒方程求出，即

在級聯貧化段

(6)

在級聯濃縮段

(7)

供料豐度Cfi,s實際上表示的是級聯中各組分的豐度分布，其分布與級聯的結構和管道連接方式密切相關。對如圖1所示的逆流型矩形級聯，Cfi,s可表示為

(8)

對于不同的級聯形式，級間的豐度混合條件不同，(8)式也將隨之改變。

2 計算方法

2.1 算法的基本思想

在Cfi,s、γij、θs已知的情況下，方程(5)是以Cpi,s為變量的非線性方程組，無法得到其解析解，只能通過數值方法求解。下面將具體介紹非線性方程組(5)的求解方法。由于各級的非線性方程組的求解過程完全相同，與s無關，因此，下面對方程組(5)的求解過程的闡述，將隱去角標s，即Cpi,s直接用Cpi表示。

(9)

也就是說，下一次的迭代值可由本次迭代值加上相鄰兩次迭代的差值(迭代序列差)來求得，將(9)式代入(5)式中得

(10)

而將(9)式代入方程(4)中可得

(11)

2.2 算法的收斂判據

收斂判據是指迭代計算停止的條件。不同的計算方法，其收斂判據也不同，有的以級聯組分守恒方程為收斂判據，有的以精、貧料組分豐度的歸一化條件為收斂判據。本方法所采用的收斂判據為級聯組分守恒方程，即定義迭代誤差為

(12)

當δ達到給定精度時，迭代停止，所得到的解為最終解。

2.3 算法的迭代步驟

算法的迭代步驟如下：

1) 根據(6)～(7)計算各級分流比θs。

3) 根據(8)式計算各級供料豐度Cfi。

4) 根據(9)～(11)式求新的Cpi，Cwi。

5) 驗證精度δ是否達到給定值，如果達到，計算結束，否則，利用所得到的Cpi和Cwi，從第3步重新迭代。

2.4 算法的優越性

從以上對算法的闡述可知，本算法不同于傳統的級聯算法，傳統的級聯算法一般是在級聯各級與級聯端部之間建立起組分物質輸運方程，且輸運方程中包含分離方程所引起的非線性項，因此各算法的基本任務就是設法求解這些非線性的組分物質輸運方程。如果級聯的形式比較復雜，就會增加非線性方程組的求解難度。而本方法的核心就是針對迭代序列差求解各級的分離方程，而各級豐度混合條件即級間管道的連接方式只作為求解各級分離方程的約束條件，所以級聯形式上的變化不會使分離方程的求解復雜化。因此本算法可以很方便地推廣到帶附加流級聯、層架級聯、非對稱級聯等復雜級聯的計算。

本算法與非穩態法[12]有些類似，其求解過程就是中間解向最終解動態靠近的過程。需要指出的是，本算法為迭代算法，在求解過程中級聯內部組分豐度的變化與時間不相關，與非穩態法有著本質的區別。

3 算法的特性研究

3.1 算法的準確性及其計算能力

根據上述算法，用Matlab軟件編寫了計算程序，由三個子程序組成，一個計算分流比，另一個計算由(8)式所確定的各級供料豐度，第3個用來求解級的分離方程。由于求解級的分離方程的子程序適用于所有形式的級聯，因此，在級聯形式改變的情況下，只需對前兩個子程序進行相應的修改，方便計算。利用所編寫的計算程序，進行了級聯數值計算，并對與級聯算法相關的一些特性進行了研究。

檢驗計算方法的準確性，可以通過兩種途徑，一是與實際級聯的分離結果進行比較；二是由不同的計算方法得到的計算結果之間的相互驗證。由于實際級聯的聯運行條件很復雜，不完全符合理論計算的假設條件。在理論研究階段，用不同計算方法相互比較來檢驗理論計算結果的準確性是可行的。

為了驗證本算法的準確性，對多遍分離硅-28級聯進行了計算，采用的級聯參數為：級聯總級數n= 30級，供料級f=15級，基本分離系數γ0=1.2，級聯分流比θc=0.6，供料流量與級間流量比為F/G=0.2，分離遍數為4遍。分離工質為SiHCl3，其組分的摩爾質量M和相應的天然豐度CFi列于表1。

表1 三氯氫硅(SiHCl3)各組分的摩爾質量及天然豐度

計算結果列于表2。四遍分離三氯氫硅，由于上一遍的精料將作為下一遍的供料，第一遍的精料豐度之和不滿足歸一化條件，誤差將傳遞給下一遍的級聯計算，并在每一遍分離計算中得到放大。到第四遍時，累積誤差可能超過分離工藝所能允許的范圍。在這種情況下，要求計算方法有比較高的精度，才能在多遍分離級聯計算之后，滿足實際分離工藝的精度要求。

表2 分離三氯氫硅(SiHCl3)的計算結果

從表2可見，本方法與準線性法計算所得到的結果完全一致。用文獻[12]中提到的第二種計算方法，即非穩態法計算上述兩種級聯時，計算結果一致。通過上面的級聯計算以及不同計算方法所得到的結果之間的比較，驗證了本計算方法的準確性。

檢驗多組分分離級聯計算能力的另一種情況是分離工質同位素組分多、級聯長、級聯分流比很小或接近全回流。雖然在實際分離工藝中很少有這種級聯，理論上，用這種級聯來檢驗一種計算方法的計算能力是完全合理的。為此，本研究以SiHCl3同位素分離為例進行了級聯計算。級聯總級數n= 50級，供料級f=10，基本全分離系數γ0=1.2。目標組分為最輕的組分。不同級聯分流比θc和不同F/G值下的目標S組分的精料豐度CP1結果示于圖2。

圖2 目標組分的精料豐度CP1與θc和F/G之間的關系Fig.2 Relation of CP1, θc and F/G

從圖2可以看出，當級聯分流比θc從0.6減小到0.02，F/G從0.3減小到0.01時，目標組分豐度CPi平穩地增加到99.824%，表示目標組分豐度值的整個曲面沒有突變(曲面局部嚴重向上凸起或向下凹陷)。結果表明，本算法在上述級聯參數下都能穩定地收斂。尤其是當θc=0.02、F/G=0.01這種接近全回流情況下能穩定收斂，進一步說明本計算方法具有很好的穩定性和在組分很多、級聯很長、級聯分流比很小等極限情況下的計算能力。

3.2 算法的收斂特性

不同的算法有不同的收斂特性，有的算法對初始豐度依賴性比較強，當初始豐度與最終解比較接近時，算法能穩定收斂，而當初始值偏離最終解比較遠時，算法不能收斂。有的算法在計算短級聯時，能穩定收斂，而當計算多組分的長級聯時，不能穩定收斂，或者只能收斂到一定精度。下面將以分離三氯氫硅為例，從三個方面來研究本算法的收斂特性，即收斂性與初始豐度的關系，收斂性與級聯長度的關系以及收斂性與基本全分離系數的關系。

在研究收斂特性與初始豐度的關系時，級聯總級數n= 30級，供料級f=6，基本全分離系數γ0=1.2，級聯分流比θ=0.6，供料流量與級間流量比為F/G=0.2。選擇三組不同的初始豐度，第一組是假設第一個組分的豐度為1，其余組分的豐度都為0，第二組為各組分豐度相等，第三組為表1所示的三氯氫硅的天然豐度。三種初始豐度下的迭代精度的對數lgδ與迭代次數V的變化關系示于圖3。

圖3 不同初值下迭代精度對數與迭代次數的關系Fig.3 Relation of accuracy and iterative times for different initial values

從圖3可以看出，本算法最終都能穩定收斂到10～14量級以上。不同的初值只對收斂速度產生一定的影響，比如本算法收斂到10-13時，三種初值下的迭代次數分別約為6 000、6 100和7 000次(見圖3)。為了比較各初值與最終結果的近似度，計算了初值與最終結果的協方差。第一、二、三組初值與最終解的協方差分別為12.531、4.537 6、0.526 18。三種初值下最終解的精度幾乎一致。所以，認為初值對最終解的精度沒有影響，只影響收斂速度，收斂速度由初值與最終解的近似程度決定的。最終解的各級供料豐度與天然豐度具有最高的近似度，把其作為迭代初始值時，算法收斂速度最快。

在研究算法收斂性與級聯長度的關系時，初始豐度值取表1所示的三氯氫硅的天然豐度。基本全分離系數γ0=1.2，級聯分流比θc=0.6，供料流量與級間流量比為F/G=0.2。供料級為第6級，不同級聯總級數下迭代精度的對數lgδ與迭代次數的變化關系如圖4所示。圖4中的曲線分別表示級聯總級數n=20、30、40、50下的收斂趨勢。

圖4中不同級聯長度下的收斂趨勢說明，級聯長度對最終的收斂精度不產生影響，也就是說，無論是在哪種級聯長度下，算法都能穩定地收斂到10-14量級以上。但是級聯長度對算法的收斂速度影響很大。比如收斂到近10-14時，級聯長度為20、30、40、50級所對應的收斂次數依次約為3 000、7 000、12 000、15 000次。

為研究算法收斂性與基本全分離系數的關系，初始豐度值取表1所示的三氯氫硅的天然豐度。計算時選用的級聯總級數n=30級，供料級為f=6級，級聯分流比θc=0.6，供料流量與級間流量比為F/G=0.2。不同基本全分離系數下迭代精度的對數lgδ與迭代次數的變化關系示于圖5。圖5中的曲線分別表示基本全分離系數γ0=1.1、1.2、1.3、1.4下的收斂趨勢。

圖5 不同基本全分離系數下迭代精度對數與迭代次數的關系Fig.5 Relation of accuracy and iterative times for different separation factor

圖5中不同基本全分離系數下的收斂趨勢說明，基本分離系數γ0對最終的收斂精度也不產生影響，無論是在哪種基本全分離系數下，算法都能穩定地收斂到10-14量級以上。同時對收斂速度有一定的影響，但其影響程度遠沒有級聯長度對收斂速度的影響大。

從以上算例分析可以看出，本算法最終的收斂精度不依賴于初始豐度，不受級聯長度的影響，與基本全分離系數的大小無關，無論在哪種情況下，都能穩定收斂到10-14量級以上。精度完全能滿足穩定同位素分離工藝對算法精度的要求。

4 結論

本研究提出了一種新的求解多組分級聯方程的數值計算方法。該方法以針對迭代序列差的迭代法求解級聯各級的分離方程，以各級豐度混合條件即級間管道的連接方式作為求解各級分離方程的約束條件，以級聯組分守恒方程為收斂判據進行了級聯數值計算，通過與其他計算方法得到的結果進比較，論證了本計算方法的準確性。同時，本計算方法在計算相對分離系數很大、組分很多、級聯很長、級聯分流比很小或接近全回流等特殊情況下的級聯時，都能穩定收斂，計算能力。

本算法可以很方便地推廣到帶附加流級聯、層架級聯、非對稱級聯等復雜級聯的計算。

參考文獻：

[1] Xie Quanxin, Sulaberidze GA, Borisevich VD. Optimization of Model cascades for Multicomponent Isotope Mixtures [C]//Physical and Chemical Processes on Selection of Atoms and Molecules :Proceedings of XI All-Russia (International) Scientific Conference. December 11-15, 2006: 15-20.

[2] 謝全新，李大勇，李文泊，等. 多組份同位素分離級聯的優化[J]. 核科學與工程，2008， 28 (1)：86-91.

Xie Quanxin, Li Dayong, Li Wenbo, et al. Cascade optimization for multicomponent isotope separation [J]. Chinese Journal of Nuclear Science and Engineering, 2008,28(1): 86-91(in Chinese).

[3] 謝全新，李大勇，吳建軍. 回收鈾分離級聯計算[J]. 核技術，2009,32 (1): 61-65.

Xie Quanxin, Li Dayong, Wu Jianjun. Cascade calculations for separation of recycled uranium [J]. Nuclear Techniques, 2009,32 (1): 61-65(in Chinese).

[4] 周明勝，鐘松，程維娜，等. 高豐度硅-28同位素分離技術 [J]. 清華大學學報:自然科學版, 2003, 43 (6): 729-731.

Zhou Mingsheng, Zhong Song, Cheng Weina, et al. Enrichment of 28Si to high concentrations by gas centrifuge [J]. J Tshinghua Univ: Sci & Tech, 2003, 43 (6): 729-731(in Chinese).

[5] 周明勝，梁雄文，張永剛，等. 離心法生產氙同位素 [J]. 同位素，2006,19 (1): 1-3.

Zhou Mingsheng, Liang Xiongwen, Zhang Yonggang, et al. Experimental study of Xenon isotopes production by gas centrifuge [J]. Journal of Isotopes, 2006,19 (1): 1-3(in Chinese).

[6] 李大勇,牟宏, 李文泊, 等.74Se 同位素分離研究 [J]. 同位素, 2006,19 (1): 4-6.

Li Dayong, Mou Hong, Li Wenbo, et al. The study of74Se isotope separation [J]. Journal of Isotopes, 2006,19 (1): 4-6(in Chinese).

[7] 周明勝，楊光偉，聶玉光. 離心法分離餓同位素的試驗研究 [J]. 同位素，2004,17(1): 5-7.

Zhou Mengsheng, Yang Guangwei, Nie Yuguang. Experimental study on Osmium isotopes separation by gas centrifuge [J]. Journal of Isotopes, 2004,17 (1): 5-7(in Chinese).

[8] Xie Quanxin, Borisevich V D, Potapov D V. et.al. Multicomponent isotope mixture separations in square cascade with losses of working substances [J]. Engineering Physics, 2006, 2: 30-34 (in Russian).

[9] Xie Quanxin. Stationary molecular selective transfer in ordinary and multiflow cascades with large enrichment factors and losses at stages for the separation of multicomponent isotopic mixture [D]. Moscow: Moscow Engineering Physics Institute, 2007.

[10] Zeng Shi, Ying Chuntong. A robust and efficient calculation procedure for determining concentration distribution of multicomponent mixtures [J]. Separation Science and Technology, 2000, 35 (4): 613-622.

[11] Zeng Shi, Ying Chuntong. A method of separating a middle component in multicomponent isotope mixtures by gas centrifuge cascades [J]. Separation Science and Technology, 2000, 35 (14): 2 173-2 186.

[12] Wu Hongjiang, Ying Chuntong, Liu Guanjun. Calculation methods for determining the distribution of components in a separation cascade for multicomponent mixture [J]. Separation Science and Technology, 1998, 33 (6): 887-898.

[13] Zeng Shi, Ying Chuntong. A second-order time-accurate method for determination of concentration distribution of multicomponent mixtures in separation cascades[J]. Separation Science and Technology, 2000, 35 (5): 729-741.