王為國,羅旌崧,曾真,王存文,吳元欣,覃遠航
(1武漢工程大學化工與制藥學院,湖北 武漢430074;2武漢工程大學機電工程學院,湖北 武漢430074;3武漢工程大學綠色化工過程教育部重點實驗室,湖北 武漢430074)
二元提餾式間歇精餾的優(yōu)化操作與最小汽化總量
王為國1,羅旌崧1,曾真2,王存文3,吳元欣3,覃遠航3
(1武漢工程大學化工與制藥學院,湖北 武漢430074;2武漢工程大學機電工程學院,湖北 武漢430074;3武漢工程大學綠色化工過程教育部重點實驗室,湖北 武漢430074)
理想操作條件下二元提餾式間歇精餾優(yōu)化操作的汽化總量與最小汽化總量的計算是約束函數(shù)優(yōu)化問題。本文采用罰函數(shù)法,將此約束函數(shù)優(yōu)化轉(zhuǎn)變?yōu)闊o約束函數(shù)優(yōu)化,并采用固定雙步長因子梯度法數(shù)值求解該函數(shù)的極值。計算表明:固定雙步長因子梯度法具有良好的收斂性,同時,降低分段數(shù)較多時,數(shù)值截斷誤差積累對計算結(jié)果的影響。二元提餾式間歇精餾優(yōu)化操作較恒殘液組成操作的能耗低的原因如下:在理論板數(shù)相對較少(接近二元提餾式間歇精餾恒殘液組成操作所需的最少理論板)時,優(yōu)化操作通過控制再沸比提高了能耗效率;在理論板數(shù)相對較多時,優(yōu)化操作通過控制再沸比,在保證過程的能耗效率較高的同時,可盡可能快地將物料移出系統(tǒng),減少了精餾過程中塔頂貯槽內(nèi)液體的混合熵產(chǎn)。通過對計算結(jié)果的歸納與外推,得到了理想操作條件下理論板數(shù)為無窮多時二元提餾式間歇精餾優(yōu)化操作再沸比的變化方式以及最小汽化總量的計算公式。
提餾式間歇精餾;優(yōu)化操作;汽化總量;模擬;能耗
間歇精餾是經(jīng)常用于小規(guī)模生產(chǎn)的一個重要的單元操作。間歇精餾分離二元混合物,根據(jù)出料方式,具有代表意義的操作方法有常規(guī)間歇精餾(主要特征:間歇精餾過程中塔頂連續(xù)出料、塔底存料,精餾終了塔底集中出料)、塔頂累積全回流間歇精餾(主要特征:間歇精餾過程中塔頂、塔底均存料,精餾終了塔頂、塔底集中出料)與提餾式間歇精餾(主要特征:間歇精餾過程中塔頂存料,塔底連續(xù)出料,精餾終了塔頂集中出料)等[1]。間歇精餾提餾式操作方法最早由Robinson和Gilliland[2]于1950年提出。Hasebe等[3]及Mujtaba等[4-5]都對提餾式間歇精餾塔進行了不同程度的研究。Chiotti等[6-7]針對二元物系的間歇精餾過程分別建立了基于亞穩(wěn)態(tài)過程的提餾式間歇精餾塔和精餾式間歇精餾塔的數(shù)學模型。通過模擬研究發(fā)現(xiàn),與常規(guī)間歇精餾相比,在分離原料中含有大量難揮發(fā)組分的精餾過程中,提餾式間歇精餾塔更加經(jīng)濟適用。Sorensen等[8]對常規(guī)間歇精餾塔和提餾間歇精餾塔的動態(tài)特性及最優(yōu)操作時的能耗進行了比較,得到當原料中輕組分含量很低的時候,使用提餾式間歇精餾操作明顯優(yōu)于常規(guī)間歇精餾操作。許松林等[9-10]對提餾式間歇精餾塔的設(shè)計方法及控制方法進行了研究。張雪梅等[11]對二元恒再沸比提餾式間歇精餾過程的最小再沸比進行了理論推導(dǎo)。近來作者對間歇精餾基礎(chǔ)理論進行了研究[12-20],文獻[19]將提餾式間歇精餾分為恒再沸比與恒殘液組成兩種典型操作, 完善了二元恒再沸比提餾式間歇精餾最小再沸比計算方法,推導(dǎo)出二元提餾式間歇精餾恒殘液組成操作,在理想操作條件下的最小汽化總量的計算公式,分析了二元提餾式間歇精餾,相對于恒殘液組成操作,要求殘液中重組分的收率高和輕組平均濃度低時,采用恒再沸比操作能耗較高的原因。文獻[20]由理想操作條件下,恒殘液組成提餾式間歇精餾二元理想混合物,理論板數(shù)與汽化總量對應(yīng)關(guān)系的數(shù)值計算結(jié)果,關(guān)聯(lián)出理論板數(shù)與汽化總量的近似關(guān)系式(捷算法)。文獻[13]采用罰函數(shù)法,將理想操作條件下二元常規(guī)間歇精餾優(yōu)化操作的汽化總量與最小汽化總量轉(zhuǎn)變?yōu)闊o約束函數(shù)優(yōu)化求解,并采用固定單步長因子梯度法數(shù)值求解該函數(shù)的極值,在分段數(shù)較多時,優(yōu)化操作的汽化總量與最小汽化總量取不同方法賦回流比初值的收斂結(jié)果中最小者。明顯這種有限窮舉賦回流比初值方法獲得計算結(jié)果,在許多情況下,與其真數(shù)值解有較大的偏差。本文以過程能耗最小為目標,采用罰函數(shù)法,將理想操作條件下二元提餾式間歇精餾優(yōu)化操作的汽化總量與最小汽化總量轉(zhuǎn)變?yōu)闊o約束函數(shù)優(yōu)化求解,研究函數(shù)優(yōu)化求解的數(shù)值計算方法,并由數(shù)值計算結(jié)果對過程能耗進行分析,最后,由理想操作條件下理論板數(shù)為無窮多時二元提餾式間歇精餾優(yōu)化操作的數(shù)值計算結(jié)果,歸納與外推二元提餾式間歇精餾優(yōu)化操作的再沸比變化方式以及最小汽化總量計算公式。
1.1 數(shù)值計算方法
二元提餾式間歇精餾流程、操作條件及基本方程均同文獻[19],本文所謂優(yōu)化操作是指對于給定的分離任務(wù),通過控制過程再沸比,使塔釜汽化總量達到最小。二元提餾式間歇精餾一批料的塔釜汽化總量
對于給定的分離任務(wù)和理論板數(shù),采用提餾式間歇精餾求優(yōu)化操作時的汽化總量nVT,opt,在數(shù)學上是約束函數(shù)優(yōu)化問題。用罰函數(shù)法將其轉(zhuǎn)化成無約束函數(shù)優(yōu)化,其表達式為
式中,μ為罰因子(正常數(shù))。
1.1.1 工作方程 將塔頂貯槽內(nèi)液體中輕組分濃度xD由xF變?yōu)閤De的精餾過程進行等距離散,取離散段數(shù)為k,則
式(11)中,ωin與ωout為步長因子(小正數(shù))。
1.1.2 計算步驟N+1為有限時,采用固定雙步長因子梯度法求解Rb,i,opt(i=1,2,…,k)與nVT,opt的具體步驟如下:
(1)輸入:α、nF、xF、ηh、xWm、N+1、k、ωin和ωout。
(2)計算nWe、nDe、xDe、ΔxD和xD,i(i=0,1,2,…,k)。
(3)s=0,ss=0,賦初值。
(10)計算nVT,opt(由式(1)得,判斷是否達到收斂要求
①如果達不到,轉(zhuǎn)到(11);
②如果達到,轉(zhuǎn)到(12)。
(12)輸出計算結(jié)果。
對于采用上述步驟編程時,補充說明如下:
計算實踐顯示,在本文實例計算時,δ較為適合取值為δ=0.0001。理論上δ取值愈小,計算出愈接近真值,但δ取值過?。ㄈ绂?0.00001),數(shù)值截斷誤差相對較大,導(dǎo)致的計算結(jié)果不穩(wěn)定。
(3)離散段數(shù)k的取值。理論上k取值愈大,計算結(jié)果愈接近真值,但計算量愈大。對于本文分離任務(wù)(α=2.5、Fx=0.5、Wmx=0.06和hη=0.9,在后文中所說的本文分離任務(wù)均指此分離任務(wù)),由圖17可見,在N+1→∞與k=100時,數(shù)值計算結(jié)果nVT,min,opt/nF=2.0435,由式(31)得到解析計算結(jié)果nVT,min,opt/nF=2.0430,兩者相對偏差為0.0245%。說明分段數(shù)k=100時,對于本文分離任務(wù),數(shù)值計算已有較好的精度,因此,如果不做特殊說明,后文所有計算結(jié)果均是在k=100條件下算出。
(4)罰因子μ的取值與收斂判據(jù)。對于本文分離任務(wù),由構(gòu)造的罰函數(shù)[式(2)],數(shù)值計算結(jié)果中約束條件最小殘差,為了計算結(jié)果盡可能符合約束條件,后文所有計算均取μ=107。
采用固定雙步長因此梯度法求解nVT,opt的計算程序是由兩循環(huán)構(gòu)成,因此,收斂判據(jù)有兩個,分別稱之為內(nèi)循環(huán)判據(jù)與終點判據(jù)。由于內(nèi)循環(huán)是對Rb,i(i=1,2,…,k)進行滿足約束條件的迭代修正,因此,為了使內(nèi)循環(huán)收斂結(jié)果盡可能符合約束條件,內(nèi)循環(huán)判據(jù)直接為約束條件的殘差,或者經(jīng)多次內(nèi)循環(huán),如果多次迭代不滿足上述收斂判據(jù),表明ωin取值過大。終點判據(jù)即判斷nVT是否收斂至nVT,opt,本文未能很好地解決,計算終點判斷有較強主觀性。當nVT收斂至nVT,opt后,再進行迭代計算,計算出的nVT應(yīng)不變,因此,終點判據(jù)為:計算出的nVT在一較小范圍內(nèi)振蕩。由于本文計算過程愈近終點,收斂速度愈慢,在接近終點時,會出現(xiàn)nVT,在較小范圍內(nèi)振蕩接近nVT,opt,這為終點判斷帶來偏差。對本文分離任務(wù),ωin與ωout選取適宜條件下,當Rb,i=Rb,i,opt(i=1,2,…,k)時,內(nèi)循環(huán)的(?F/?Rb,i)s,ss呈現(xiàn)近似周期性變化,最顯著特征是(?F/?Rb,1)s,ss≈ (?F/?Rb,1)s,s+1。以此作為輔助終了判據(jù),可減少主觀判斷的偏差。
1.2 固定雙步長因子梯度法的提出與驗證
梯度法]22[,又稱最速下降法,是一種求解無約束多元函數(shù)極值的數(shù)值方法,具有良好的收斂性。本文采用罰函數(shù)法,將求解理想操作條件下,二元提餾式間歇精餾優(yōu)化操作的汽化總量nVT,opt與最小汽化總量nVT,min,opt轉(zhuǎn)變求解無約束函數(shù)極值。本文首先采用固定單步長因子梯度法迭代求解本文分離任務(wù)的nVT,opt。在理論板數(shù)相對較少時(接近恒殘液組成提餾式間歇精餾所需的最少理論板),選取適宜的步長因子,可收斂至穩(wěn)定結(jié)果。對于本文分離任務(wù),二元提餾式間歇精餾恒殘液組成操作所需的最少理論板為6.45。圖1為N+1=7與k=100時,采用固定單步長因子梯度法數(shù)值求解本文分離任務(wù)nVT,opt的收斂過程與結(jié)果。由圖1可見。在計算終點,三者的|ωs(?F/?Rb,i)/Rb,i|<10?7(i=1,2,…,k),因此,采用固定單步長因子梯度法,nVT,opt最終收斂結(jié)果受數(shù)值截斷誤差積累影響,ωs取值愈小,數(shù)值截斷誤差積累對計算結(jié)果影響愈大。當ωs=0.004時,在迭代230萬次后,nVT,opt/nF在4.446875~ 4.448421范圍內(nèi)波動,在 10?6~10?5范圍內(nèi)波動。即采用固定單步長因子梯度法數(shù)值求解,ωs取值過大,不能得到穩(wěn)定的收斂結(jié)果。此外,圖4中ωin,c也是固定單步長因子梯度法度法迭代求解nVT,opt,使Rb,i(i=1,2,…,k)收斂到滿足約束條件(收斂精度為10?7)的ωs最大取值(N+1=7,ωin,c≈0.0031)。
圖1 固定單步長因子梯度法的收斂過程(N+1=7)Fig.1 Convergence process of fixed single step length factor gradient method (N+1=7)
對于本文分離任務(wù),N+1較多時(N+1≥12),采用固定單步長因子梯度法迭代求解nVT,opt得到:在k=100與步長因子取值滿足計算過程收斂條件下,如果給定的Rb,i(i=1,2,…,k)初值滿足約束條件(如:以二元提餾式間歇精餾恒殘液組成操作的計算結(jié)果為初值),出現(xiàn)了“初值即結(jié)果”現(xiàn)象,如果給定的初值不滿足約束條件,出現(xiàn)了“初值決定結(jié)果”現(xiàn)象。采用精確線性搜索最優(yōu)步長因子梯度法迭代求解nVT,opt,也得到類似的現(xiàn)象。
對于本文分離任務(wù),采用固定單步長因子梯度法迭代求解nVT,opt,計算結(jié)果顯示,在N+1較多(N+1≥12)、k=100與步長因子取值滿足計算過程收斂條件下,nVT收斂至穩(wěn)定結(jié)果時,|ωs(?F/?Rb,i)/Rb,i|(i=1,2,…,k)中絕大部分的數(shù)值小于10?7,只極少部分(或沒有)的數(shù)值在一定范圍內(nèi)振蕩。再根據(jù)采用固定雙步長因子梯度法迭代求解的本文分離任務(wù)的nVT,opt結(jié)果,本文認為N+1較多時(對于本文分離任務(wù),N+1≥12),采用固定單步長因子梯度法迭代求解nVT,opt,出現(xiàn)了“初值即結(jié)果”或“初值決定結(jié)果”現(xiàn)象原因為:在離散段數(shù)較多時,Rb,1,opt與Rb,2,opt不呈漸變[參見圖8(c)與(d)]。在接近收斂結(jié)果時,|?F/?Rb,1|?|?F/?Rb,i|(i=1,2,…,k)(參見圖2),相對較大|?F/?Rb,1|使Rb,i(i=1,2,…,k)收斂到滿足約束條件(收斂精度為10?7)的步長因子允許取值較?。▍⒁妶D4中ωin,c)。而步長因子取值較小時,數(shù)值截斷誤差積累對計算結(jié)果影響較大。
圖2 優(yōu)化操作時?F/?Rb,i,opt隨D,i的變化(N+1=15)Fig.2 Changing of ?F/?Rb,i,optwithD,iunder optimal operation (N+1=15)
對于本文分離任務(wù),采用精確線性搜索法確定最優(yōu)步長因子的梯度法迭代求解nVT,opt,計算結(jié)果顯示,在滿足約束條件后,再繼續(xù)迭代,Rb,i與nVT不再變化。此時,精確線性搜索法確定的最優(yōu)步長因子近似為0。由式(7)可見,?F/?Rb,i由兩部分構(gòu)成,其一為汽化總量貢獻,記為Ai,其二為約束條件貢獻,記為。圖3表示的采用固定雙步長因子梯度法,以二元提餾式間歇精餾恒殘液組成操作的計算結(jié)果,Rb,i(i=1,2,…,k)為初值,編程計算本文分離任務(wù)的nVT,opt,在數(shù)值計算的起始與終了時,Ai與Ci的變化。圖3中,線1、2分別表示在數(shù)值計算的起始與終了時,Ai的數(shù)值(點,A1=0.29435未在線2中示出),線3、4分別表示在數(shù)值計算的起始與終了時,Ci的數(shù)值(點,C1=0.28314未在線4中示出)。由圖3可見,Ai與Ci在計算過程中變化很小且很接近,因此,本文認為|B|變化是引起?F/?Rb,i變化的主要因素,對不滿足約束條件的Rb,i(i=1,2,…,k),由于|B|較大,B×Ci在?F/?Rb,i占比例較大,梯度法迭代修正Rb,i(i=1,2,…,k)時,B×Ci起作用較大,即:梯度法迭代修正Rb,(ii=1,2,…,k)優(yōu)先滿足約束條件。采用精確線性搜索法確定最優(yōu)步長因子的梯度法,線性搜索的目標是||?F/?Rb,i||達到最小,對滿足約束條件的Rb,i(i=1,2,…,k)進行修正,如果步長因子不為0,修正后必然會使其不再滿足約束條件,|?F/?Rb,i|會大幅增大,因此,在滿足約束條件后,再繼續(xù)迭代,精確線性搜索法確定的最優(yōu)步長因子近似為0。
圖3 數(shù)值計算過程中Ai與Ci的變化(N+1=15)Fig.3 Changing ofAiand Ciin process of numerical calculation (N+1=15)
由N+1=7與k=100時,采用固定單步長因子梯度法數(shù)值求解本文分離任務(wù)nVT,opt的收斂結(jié)果(參見圖1)可知,,即:采用固定單步長因子梯度法,ωs取值愈大,數(shù)值截斷誤差積累對計算結(jié)果影響愈小。為了解決采用固定單步長因子梯度法迭代求解nVT,opt,離散段數(shù)較多時,步長因子取值滿足Rb,i(i=1,2,…,k)可收斂至約束條件,收斂到的穩(wěn)定結(jié)果,在許多情況下,由于數(shù)值截斷誤差積累對計算結(jié)果影響較大,數(shù)值計算結(jié)果與真數(shù)值解有較大的偏差;而步長因子取值較大,數(shù)值計算過程不收斂的矛盾,本文提出一迭代求解nVT,opt的計算方法,具體計算過程為:先采用較小步因子長的梯度法對Rb,i(i=1,2,…,k)進行迭代修正,當Rb,i(i=1,2,…,k)滿足約束條件后,再采用較大步長因子的梯度法對Rb,i(i=1,2,…,k)進行修正一次,再轉(zhuǎn)入Rb,i(i=1,2,…,k)滿足約束條件的小步因子長梯度法迭代修正,此過程重復(fù)進行(詳見1.1.2節(jié)),本文稱此計算方法為固定雙步長因子梯度法。
對于本文分離任務(wù),17N+=時,采用固定雙步長因子梯度法求解nVT,opt的最終結(jié)果為:,較采用固定單步長因子梯度法計出的最小略小,其原因為:采用固定雙步長因子梯度法的進行數(shù)值計算時,ωout可取較大數(shù)值,在計算終點,|ωout(?F/?Rb,i)/Rb,i|(i=1,2,…,k)中絕大部分的數(shù)值大于10?6,即:在迭代過程中,外循環(huán)的修正量ωout(?F/?Rb,i) (i=1,2,…,k)始終對Rb,i(i=1,2,…,k)起到修正作用。因此,固定雙步長因子梯度法,ωout可取較大數(shù)值,有效降低了數(shù)值截斷誤差積累對計算結(jié)果影響,數(shù)值計算結(jié)果與實際的偏差更小。采用固定雙步長因子梯度法求解nVT,opt,計算結(jié)果顯示,對本文分離任務(wù),Rb,i(i=1,2,…,k)初值適宜,ωin與ωout取值合理,均可計算出合理的結(jié)果。
1.3 實例的數(shù)值計算結(jié)果與討論
1.3.1 內(nèi)循環(huán)步長因子的取值 采用固定雙步長因子梯度法編程計算本文分離任務(wù)的Rb,i,opt(i=1,2,…,k)與nVT,opt,ωin對內(nèi)循環(huán)的次數(shù)有較大的影響。圖4中虛線ωin,c是采用收斂結(jié)果Rb,i,opt(i=1,2,…,k)的反算出來,主觀性較強,由收斂結(jié)果Rb,i,opt(i=1,2,…,k)經(jīng)一次外循環(huán)后,保證程序內(nèi)循環(huán)收斂的最大步長因子。在Rb,i(i=1,2,…,k)接近Rb,i,opt(i=1,2,…,k)與ωin=(0.4~0.6)ωin,c時,數(shù)次內(nèi)循環(huán)即可滿足約束條件。圖4中實線為本文舉例計算過程中內(nèi)循環(huán)步長因子ωin實際取值,ωin≈0.5ωin,c。此外,當Rb,i(i=1,2,…,k)與Rb,i,opt(i=1,2,…,k)差別較大時,取ωin≈0.5ωin,c,每一次外循環(huán)后,要再收斂至滿足約束條件,需內(nèi)循環(huán)次數(shù)較多。如N+1=15時,以恒殘液組成操作的計算結(jié)果為初值,ωout=0.011和ωin=0.000025時,計算起始時,每一次外循環(huán)后,要再收斂至滿足約束條件,需經(jīng)300余次內(nèi)循環(huán)。
圖4ωin與N+1的關(guān)系Fig.4 Relationship betweenωinandN+1
1.3.2 再沸比的初值對計算結(jié)果的影響 采用固定雙步長因子梯度法編程計算本文分離任務(wù)的Rb,i,opt(i=1,2,…,k)與nVT,opt,計算結(jié)果顯示,當N+1≤12時,計算結(jié)果僅受ωout影響,而當N+1≥13時Rb,i(i=1,2,…,k)初值對計算結(jié)果有細微的影響,以1N+=15為例:
①以N+1=14時計算出的為初值、ωout=0.011和ωin=0.000025時,外循環(huán)10余萬次,計算出的nVT,opt/nF=2.056478;
②以N+1=16時計算出的為初值、ωout=0.011和ωin=0.000025時,外循環(huán)10余萬次,計算出的nVT,opt/nF=2.056417;
④以Rb,1=0.1Rbmc,1,Rb,i=Rbmc,i(i=1,2,…,k)(Rbmc用本式(25)計算)為初值,ωout=0.011和ωin=0.000025外循環(huán)10萬余次,計算出的nVT,opt/nF=2.056405。
四者的nVT,opt差別出現(xiàn)在小數(shù)點后第5位,最大相對偏差為0.005%。由此可見:采用固定雙步長因子梯度法編程計算本文分離任務(wù)的Rb,i,opt(i=1,2,…,k)與nVT,opt,ωin與ωout取值合理時,初值對計算結(jié)果影響可忽略,但影響計算工作量。此外,在計算終點,四者的內(nèi)循環(huán)均呈現(xiàn)近似周期性變化,因此,本文認為四者的偏差是由計算過程中數(shù)值截斷誤差積累引起的,即:采用固定雙步長因子梯度法編程計算本文分離任務(wù)的Rb,i,opt(i=1,2,…,k)與nVT,opt,可降低分段數(shù)較多時,數(shù)值截斷誤差積累對計算結(jié)果的影響,但不能完全消除。
1.3.3 外循環(huán)步長因子取值 圖5中ωout是采用固定雙步長因子梯度法編程計算本文分離任務(wù)的Rb,i,opt(i=1,2,…,k)與nVT,opt的外循環(huán)步長因子取值。外循環(huán)步長因子取值在(1±0.5)ωout范圍內(nèi),計算出的nVT,opt/nF數(shù)值差異為±10?6,ωout取值超出2ωout時,未觀察到計算終點的近似周期性,當N+1≥13時, 外循環(huán)步長因子取值超出5ωout時,直接以Rb,i,opt(i=1,2,…,k)為初值,計算過程明顯發(fā)散。1.3.4 實例的數(shù)值計算結(jié)果 對于本文分離任務(wù),采用固定雙步長因子梯度法的計算結(jié)果見圖6~圖8。圖6為數(shù)值計算出隨N+1變化,由圖6可見,優(yōu)化操作的能耗均低于恒殘液組成操作,其原因?qū)⒃诒疚牡?節(jié)討論。
圖5ωout與N+1的關(guān)系Fig.5 Relationship betweenωoutandN+1
圖6 優(yōu)化操作與恒殘液組成操作的汽化總量Fig.6 Total evaporation of optimize operation and constant residual fraction composition operation
圖7(a)~(d)分別為理論板數(shù)7、10、12和20塊時,優(yōu)化操作及恒殘液組成操作再沸比的變化規(guī)律。由圖6可見,優(yōu)化操作與恒殘液組成操作的再沸比均隨著xD的增大而增大。由圖7(c)、(d)可見,理論板數(shù)較大時,優(yōu)化操作的再沸比,在起始階段不呈漸變。
圖7N+1對Rb,i,opt隨變化規(guī)律的影響Fig.7 Effect ofN+1 on rule ofRb,i,optwith change of
圖8W,i,opt(i=1,2,…,k)隨N+1的變化規(guī)律Fig.8 Rule ofW,i,opt(i=1,2,…,k) with change ofN+1
由圖6可見,優(yōu)化操作的能耗均低于恒殘液組成操作,下文將對其原因進行分析。圖9表示隨N+1的變化規(guī)律,數(shù)據(jù)來源于圖6,由圖9可見,存在極小值,表明至少有兩個因素對其有影響,下文從精餾過程中的能耗效率與塔頂貯槽內(nèi)液體的混合熵產(chǎn)兩方面進行討論。
2.1 能耗效率的計算方法與結(jié)果
文獻[19]提出
圖9(nVT,xW=C?nVT,opt)nF隨N+1的變化規(guī)律Fig.9 Rule of (nVT,xW=C?nVT,opt)nFwith change ofN+1
用φ數(shù)值大小的表示其對應(yīng)的能耗對分離效果的影響程度,0≤φ≤1,對應(yīng)φ值較小的這一部分能耗為低效能耗,本文直接稱φ為能耗效率??偲骄芎男师誸為能耗效率對總能耗平均,由式(1)得
將式(14)離散得
式(15)中,φa,i為第i段的平均能耗效率,第i段的再沸比為Rb,i,此段的能耗中能包括0<Rb<Rb,i所對應(yīng)的能耗, 由文獻[19]圖5可見,φ隨Rb變化,因此,第i段的平均能耗效率為Rb在區(qū)間{0,Rb,i}內(nèi)對應(yīng)能耗效率平均,即
式(13)代入式(16)積分整理得
對給定分離任務(wù),φt的計算方法為:先將過程離散進行數(shù)值計算,應(yīng)用其計算結(jié)果,由式(17)計算出第i段的φa,i(i=1,2,…,k),再式(15)計算出φt。本文分離任務(wù)φt的計算結(jié)果見圖10,計算φt所需的變量數(shù)值源于獲得圖6的計算結(jié)果。
圖10φt隨N+1的變化規(guī)律Fig.10 Rule ofφtwith change ofN+1
2.2 塔頂貯槽內(nèi)液體混合熵產(chǎn)的計算方法與結(jié)果
混合是分離的逆過程,二元提餾式間歇精餾過程中,塔頂貯槽內(nèi)不同濃度液體混合必然有熵產(chǎn),減少塔頂貯槽內(nèi)液體的混合熵產(chǎn)意味著降低過程能耗。對塔頂貯槽內(nèi)液體混合熵衡算得:
式(18)中
將式(19)、式(20)離散得
對給定分離任務(wù),塔頂貯槽內(nèi)流動液體混合熵產(chǎn)Sm,p的計算方法為:先將過程離散進行數(shù)值計算,應(yīng)用其計算結(jié)果,由式(23)、式(22)與式(21)分別計算流出系統(tǒng)混合熵Sm,o、流入系統(tǒng)混合熵Sm,i與系統(tǒng)混合熵積累Sm,a,再通過式(18)計算出系統(tǒng)混合熵產(chǎn)Sm,p。對于本文分離任務(wù),Sm,p的計算結(jié)果見圖11,計算Sm,p所需的變量數(shù)值源于獲得圖6的計算結(jié)果。
圖11Sm,p/nF的隨N+1的變化規(guī)律Fig.11 Rule ofSm,p/nFwith change ofN+1
2.3 優(yōu)化操作的能耗分析
圖12中nW,1,opt/nF源于獲得圖6的計算結(jié)果,表示優(yōu)化操作時第1段的塔底出料量nW,1,opt/nF與N+1的關(guān)系。由圖12可見,nW,1,opt/nF隨N+1的增大而增加,由圖11可見,Sm,p/nF隨N+1的增大而減小。因此,本文認為,二元提餾式間歇精餾優(yōu)化操作,N+1較多時,精餾起始階段出料量較大,減小過程后期的塔頂貯槽的存液量,從而減小塔頂貯槽內(nèi)液體的混合熵產(chǎn)。
圖12 ΔnW,1,opt/nF隨N+1的變化規(guī)律Fig.12 Rule of ΔnW,1,opt/nFwith change ofN+1
結(jié)合圖10~圖12可得,對于本文分離任務(wù),N+1=7,8時,。再由圖7(a)可見,優(yōu)化操作的Rb與恒殘液組成操作有明顯差異。因此,在理論板數(shù)相對較少時,二元提餾式間歇精餾優(yōu)化操作的能耗較恒殘液組成操作低的原因是優(yōu)化操作通過控制再沸比,提高過程的能耗效率。N+1≥ 12時,但均接近1,。因此,在理論板數(shù)相對較多時,二元提餾式間歇精餾優(yōu)化操作的能耗較恒殘液組成操作低的原因是優(yōu)化操作通過控制再沸比,在保證過程的能耗效率較高前提下,盡可能快地將物料移出系統(tǒng),減少了過程中塔頂貯槽內(nèi)液體的混合熵產(chǎn)。此外,N+1=10時,。再由圖7(b)可見,優(yōu)化操作與恒殘液組成操作的Rb,i(i=1,2,…,k)差異很小,因此,兩者能耗差達到最小。
對于給定的分離任務(wù),在理想操作條件下,二元提餾式間歇精餾優(yōu)化操作,N+1→∞,nVT,opt→nVT,min,opt。
3.1 最小汽化總量的數(shù)值計算方法
由文獻[19]圖4可見,N+1→∞時,塔底殘液中輕組分瞬時濃度與塔頂貯槽液體中輕組分的瞬時濃度及再沸比之間的關(guān)系可表達成為
式(24)中,Rbmc為N+1→∞時塔頂和塔底同時為恒濃區(qū)的再沸比,由文獻[19]知
由式(24)得
將1.1節(jié)所介紹的nVT,opt計算過程中,步驟(4)的計算方法,取塔頂貯槽液體中輕組分平均濃度,改用式(24)計算;步驟(6)用數(shù)值微分法計算(i=1,2,…,k),取塔頂貯槽液體中輕組分平均濃度,改用式(26)計算,即可得到最小汽化總量的數(shù)值計算方法。
采上述方法編程對本文分離任務(wù)的nVT,min,opt/nF進行計算。Rb,i,min(i=1,2,…,k)的初值對計算結(jié)果有細微的影響,k=100、ωin=10?6與ωout=0.001時,Rb,i,min(i=1,2,…,k)初值不同,nVT,min,opt/nF計算結(jié)果在2.043471~2.0434786內(nèi)波動,本文取nVT,min,opt/nF=2.0435是本文分離任務(wù)優(yōu)化操作時所需的最小汽化總量。
圖13中,實線表示N+1與nVT,opt/nF的關(guān)系,數(shù)據(jù)同圖6,虛線是N+1→∞時,nVT,min,opt/nF的數(shù)值計算結(jié)果。由圖13可見,隨N+1增大,nVT,opt/nF的nVT,min,opt與直接數(shù)值計算出的nVT,min,opt是相同的。
圖13nVT,min,opt/nF與nVT,opt/nF的比較Fig13 Comparison ofnVT,min,opt/nFwithnVT,opt/nF向nVT,min,opt/nF趨近,表明nVT,opt的計算結(jié)果外推出
3.2 最小汽化總量的解析計算方法
由式(1)可見,要得到nVT,min,opt的計算公式,必須得到N+1→∞時精餾過程的Rb,min,opt變化規(guī)律。由圖14可見,除第一點Rb,min,1,opt外,數(shù)值計算出的Rb,min,i,opt與Rbmc,i是重合的,直接對兩者數(shù)值進行比較得到|Rb,min,i,opt?Rbmc,i|<10?5。的數(shù)值計算在0~10?6范圍內(nèi)波動,結(jié)合圖8的(i=1,2,…,k)隨N+1的變化趨勢可得:二元提餾式間歇精餾優(yōu)化操作,N+1→ ∞時,除起始階段,Rb,min,opt=Rbmc,xW,opt=0 。
圖14N+1→∞時Rb,min,i,opt與關(guān)系Fig.14 Relationship betweenRb,min,i,optandunderN+1→∞
圖15、圖16表示過程的離散段數(shù)k對Rb,min,1,opt與的影響,由圖15可見,隨著分段數(shù)的增加,Rb,min,1,opt向0趨近,即:,由圖16可見,隨著分段數(shù)的增加,向xF趨近。即:。由此可得:二元提餾式間歇精餾優(yōu)化操作,N+1→∞時,起始階段,Rb,min,opt=0,xW,opt=xF。
綜上所述,N+1→∞時,二元提餾式間歇精餾優(yōu)化操作,Rb,min,opt的變化方式為由起始Rb,min,opt=0(xW,opt=xF)突變?yōu)镽b,min,opt=Rbmc(xW,opt=0),并保持到過程結(jié)束。用nWC表示精餾起始至突變點前的釜底出料量,Rb,min,opt的變化方式可表達為
式(27)代入式(1)得
圖15k對Rb,min,1,opt計算結(jié)果的影響Fig.15 Effect ofkon numerical results ofRb,min,1,opt
圖16k對W,1,opt計算結(jié)果的影響Fig.16 Effect ofkon numerical results ofW,1,opt
二元提餾式間歇精餾優(yōu)化操作,塔底出料量由0到nWC精餾過程中,xW,opt=xF;釜底出料量由nWC到nWe精餾過程中,xW,opt=0。對輕組分物料衡算得
將式(30)微分得
將式(29)與式(31)代入式(28)積分整理得
對于本文分離任務(wù),由式(29)計算得:nWC/nF= 0.0574,由式(32)計算得:nVT,min,opt/nF=2.0430。
圖17中,實線表示數(shù)值計算結(jié)果nVT,min,opt/nF與離散段數(shù)k關(guān)系,虛線表示nVT,min,opt/nF的解析計算結(jié)果。由圖17可見,隨k增大,數(shù)值計算結(jié)果向解析計算結(jié)果趨近,表明由不同離散段數(shù)數(shù)值計算結(jié)果外推出的nVT,min,opt/nF與解析計算結(jié)果是相同的。
圖17k對nVT,min,opt/nF計算結(jié)果的影響Fig.17 Effect ofkon numerical results ofnVT,min,opt/nF
3.3 最小汽化總量的解析計算方法的合理性說明
由數(shù)值計算結(jié)果歸納出nVT,min,opt解析計算方法的合理性需用數(shù)學分析法證明,本文未能完成。下文通過對給定分離任務(wù)的nVT,min,opt和解析計算結(jié)果比較及熱力學分析,對nVT,min,opt解析計算方法的合理性進行說明。
表1中9個分離任務(wù)是理想操作條件下,常用范圍內(nèi),對確定二元提餾式間歇精餾分離任務(wù)的4個獨立變量(α,xF,xWm,ηh)采用3水平正交設(shè)計方案取值擬給定的,二元恒殘液組成提餾式間歇精餾的塔釜最小汽化總量nVT,min,xD=CnF按文獻[19]式(25)計算。由表1可見,對于所有的分離任務(wù)nVT,min,opt/nF均小于nVT,min,xD=CnF,因此,從具體分離任務(wù)計算結(jié)果看,nVT,min,opt解析計算方法是合理的。
由圖12可見,ΔnW,1,opt/nF隨N+1的增大而增加,并向ΔnW,1,opt,N+1→∞/nF趨近,即:對本文分離任務(wù),k=100時,ΔnW,1,opt,max/nF=ΔnW,1,opt,N+1→∞/nF= 0.0650。對于本文分離任務(wù),由式(29)得:nWC/nF=0.0574。由圖18可見,隨著分段數(shù)增大,ΔnW,1,opt,N+1→∞/nF向nWC/nF=0.0574趨近,即,nWC/nF=0.0574是優(yōu)化操作時起始階段最大出料量,結(jié)合圖10、圖11可得,此時,塔頂貯槽內(nèi)液體的混合熵產(chǎn)達最小。
圖18k對ΔnW,1,opt,N+1→∞/nF計算結(jié)果的影響Fig.18 Effect ofkon numerical results of ΔnW,1,opt,N+1→∞/nF
綜上所述,歸納與外推得到,N+1→∞時,在理想操作條件下,二元提餾式間歇精餾優(yōu)化操作,Rb,min,opt的變化方式由起始Rb,min,opt=0(xW,opt=xF)突變?yōu)镽b,min,opt=Rbmc(xW,opt=0),并保持到過程結(jié)束,實現(xiàn)了過程能耗效率(φt=1)最高,塔頂貯槽內(nèi)液體的混合熵產(chǎn)最小,即過程能耗最小。因此,從熱力學上說,nVT,min,opt解析計算方法是合理的。
表1nVT,min,opt/nF和nVT,min,xD=C解析計算結(jié)果的比較Table 1 Analytical results ofnVT,min,opt/nFcompared withnVT,min,xD=C
(1)采用罰函數(shù)法,將求解理想操作條件下,二元提餾式間歇精餾優(yōu)化操作的汽化總量與最小汽化總量轉(zhuǎn)變求解無約束函數(shù)極值。通分析該函數(shù)的結(jié)構(gòu)與具體計算數(shù)據(jù),設(shè)計一數(shù)值法求解極值的計算過程,先采用較小步因子長的梯度法對Rb,i,(i=1,2,…,k)進行迭代修正,當Rb,i,(i=1,2,…,k)滿足約束條件后,再采用較大步長因子的梯度法對Rb,i,(i=1,2,…,k)進行修正一次,再轉(zhuǎn)入Rb,i,(i=1,2,…,k)滿足約束條件的小步因子長梯度法迭代修正,此過程重復(fù)進行(本文稱之為固定雙步因子長梯度法)。模擬的計算實踐顯示,固定雙步長因子梯度法求解二元提餾式間歇精餾優(yōu)化操作的nVT,opt與nVT,min,opt具有良好的收斂性,同時,降低分段數(shù)較多時,數(shù)值截斷誤差積累對計算結(jié)果的影響。
(2)由模擬計算結(jié)果對過程能耗分析得到提餾式間歇精餾優(yōu)化操作的能耗較恒殘液組成操作低的原因。在理論板數(shù)相對較少(接近恒殘液組成提餾式間歇精餾所需的最少理論板)時,優(yōu)化操作通過控制過程再沸比,提高了過程的能耗效率;在理論板數(shù)相對較多時,優(yōu)化操作通過控制過程再沸比,保證過程的能耗效率較高同時,盡可能快地將物料移出系統(tǒng),減少了過程中塔頂貯槽內(nèi)液體的混合熵產(chǎn)。
(3)由模擬計算結(jié)果歸納與外推得到,N+1→∞時,在理想操作條件下,二元提餾式間歇精餾優(yōu)化操作再沸比的變化方式。Rb,min,opt變化方式為由起
始Rb,min,opt=0(xW,opt=xF)突變?yōu)镽b,min,opt=Rbmc
(xW,opt=0),并保持到過程結(jié)束。由Rb,min,opt的此變化方式推導(dǎo)出最小汽化總量的計算公式。Rb,min,opt的此變化方式,實現(xiàn)了過程能耗效率(φt=1)最高,塔頂貯槽內(nèi)液體的混合熵產(chǎn)最小,即過程能耗最小。
符 號 說 明
k——數(shù)值計算時,全過程離散的段數(shù)
N+1 ——理論板數(shù)(包括塔釜)
nD、nF、nDe——分別為塔頂貯槽中瞬時存液量、原料量與塔頂產(chǎn)品量,kmol
nVT——精餾一批料的汽化總量,kmol
nW、nWe——分別為精餾起始至某瞬時塔底出料液總量與塔底產(chǎn)品量,kmol
nWC——N+1→∞時,起始Rb,min,opt=0的塔底出料液總量,kmol
ΔnW,i、ΔnV,i——精餾過程第i段塔底出料液總量與汽化量,kmol
R——通用氣體常數(shù),R=8.314472 kJ·K?1·kmol?1
Rb——再沸比
Rbmc——N+1→∞時,塔頂和塔底同時為恒濃區(qū)的再沸比
Sm,a、Sm,i、Sm,o、Sm,p——分別為塔頂貯槽內(nèi)液體混合熵積累、流入塔頂貯槽混合熵、流出塔頂貯槽混合熵與塔頂貯槽內(nèi)液體混合熵產(chǎn),kJ·K?1·kmol?1
xD、xDe——分別為精餾過程與終了時,塔頂貯槽液體中輕組分的摩爾分數(shù)
ΔxD——xD由xF變?yōu)閤De的精餾過程的離散距離
xD,i、——分別為第i段后塔頂貯槽液體中輕組分的摩爾分數(shù)與第i段塔頂貯槽液體中輕組分的平均摩爾分數(shù)
xF——原料中輕組分的摩爾分數(shù)
xW、、xWm——分別為底出料中輕組分的瞬時摩爾分數(shù)、第i段塔底出料中輕組分的平均摩爾分數(shù)與塔底產(chǎn)品中輕組分的平均摩爾分數(shù)
yD——與xD平衡汽相中輕組分的摩爾分數(shù)
α——相對揮發(fā)度
ηh——塔底產(chǎn)品中重組分的收率
ε——收斂精度
φ、φa,i、φt——分別為瞬時、第i段平均與總平均能耗效率
ωs——固定單步長因子梯度法的步長因子
ωin、ωout——分別為固定雙步長因子梯度法的內(nèi)循環(huán)與外循環(huán)的步長因子
ωin,c——可收斂至約束條件的內(nèi)循環(huán)最大步長因子
上角標
0,0——數(shù)值計算時,各物理量的初值
s,ss——計算程序的第ss次外循環(huán)中s次內(nèi)循環(huán)時各物理量的值
下角標
i、j——數(shù)值計算時,第i段與第j段各物理量的值
opt ——優(yōu)化操作時,各物理量的值
min ——N+1→∞時,各物理量的值
xW=C——恒殘液組成操作時,各物理量的值
[1] Editorial Committee of Handbook Chemical Engineering (《化學工程手冊》編輯委員會). Handbook Chemical Engineering (化學工程手冊) [M]. Beijing: Chemical Industry Press, 1989.
[2] Robinson C S, Gilliland E R. Elements of Fractional Distillation [M]. 4th ed. New York: McGraw-Hill, 1950.
[3] Hasebe S, Abdul-Aziz B B, Hashimoto I, Watanabe T. Optimal design and operation of complex batch distillation column//Proceedings of the IFAC Workshop on Interaction between Process Design and Process Control [C]. London, 1992: 177-182.
[4] Mujtaba I M, Macchietto S. Optimal operation of reactive batch distillation//AIChE Annual Meeting [C]. Miami Beach, USA, 1992.
[5] Mujtaba I M, Macchietto S. Optimal operation of multicomponent batch distillation—a comparative study using conventional and unconventional columns//Rreprints IFAC symposium ADCHEM’95 [C]. Kyoto, Japan, 1994: 415-420.
[6] Chiotti O J, Iribarren O A. Simplified models for binary batch distillation [J].Computers & Chemical Engineering, 1991, 15 (1): 1-5.
[7] Chiotti O J, Salomone H E, Iribarren O A. Selection of multicomponent batch distillation sequences [J].Chemical Engineering Communications, 1993, 119 (1): 1-21.
[8] Sorensen E, Skogestad S. Comparison of regular and inverted batch distillation [J].Chem.Eng.Sci., 1996, 51 (22): 4949-4952.
[9] Xu Songlin, Salomone H E, Iribarren O A. Shortcut procedure for inverted batch distillation column (Ⅰ): Multicomponent ideal system [J].Chinese Journal of Chemical Engineering, 2001, 9 (1): 28-33.
[10] Xu Songlin, Espinosa J, Salomone H E, Iribarren O A. Operation of a batch stripping distillation column [J].Chinese Journal of Chemical Engineering, 2001, 9 (2): 141-144.
[11] Zhang Xuemei (張雪梅), Guo Jintang (郭錦棠), Zhang Weijiang (張衛(wèi)江). Binary stripping batch distillation [J].Journal of Tianjin University(天津大學學報), 2005, 38 (11): 950-954.
[12] Wang Weiguo (王為國), Wu Yuanxin (吳元欣), Wang Cunwen (王存文), Zhen Zhen (曾真). Algorithm of minimum reflux ratio of batch distillation under constant reflux ration and its energy consumption analysis [J].Journal of Chemical Industry and Engineering(China) (化工學報), 2004, 55 (8): 1285-1290.
[13] Wang Weiguo (王為國), Wang Cunwen (王存文), Wu Yuanxin (吳元欣) , Zhen Zhen (曾真). Minimum total vapor capacity of binary regular batch distillation [J].Journal of Chemical Industry and Engineering(China) (化工學報), 2006, 57 (11): 2647-2651.
[14] Wang Weiguo (王為國), Zeng Zhen (曾真), Bi Yafan (畢亞凡). Minimum total vapor capacity of batch distillation of binary mixture under total reflux and top accumulation [J]. Journal of Chemical Industry and Engineering(China) (化工學報), 2001, 52 (5): 460-463.
[15] Wang Weiguo (王為國), Zeng Zhen (曾真), Bi Yafan (畢亞凡), Sun Wei (孫煒). Mixing of flowing liquid in top vessel with batch distillation under total reflux [J].Journal of Chemical Industry and Engineering(China) (化工學報), 2002, 53 (8): 857-861.
[16] Wang Weiguo (王為國), Wang Cunwen (王存文), Wu Yuanxin (吳元欣), Zeng Zhen (曾真). Energy consumption of batch distillation of binary mixture under “total reflux” [J].Journal of Chemical Industry and Engineering(China) (化工學報), 2004, 55 (9): 1474-1480.
[17] Wang Weiguo, Wang Cunwen, Zhou Fanglei, Qin Yuanhang, Wu Yuanxin, Zeng Zhen. Energy consumption in binary batch distillation of total reflux and top accumulation [J].Advanced Material Research, 2012, 396: 699-705.
[18] Wang Weiguo (王為國), Zeng Zhen (曾真), Wang Cunwen (王存文), Wu Yuanxin (吳元欣). Minimum theoretical trays for binary batch distillation under total reflux and top accumulation [J]. Petrochemcal Technology(China) (石油化工), 2011, 40 (11): 1205-1210.
[19] Wang Weiguo (王為國), Zeng Zhen (曾真), Qin Yuanhang (覃遠航), Wang Cunwen (王存文), Wu Yuanxin (吳元欣). Energy consumption analysis and minimum reboil ratio of stripping batch distillation under constant reboil ratio [J].CIESC Journal(化工學報), 2012, 63 (7): 2106-2112.
[20] Wang Weiguo (王為國), Duan Xiaoling (段曉玲), Zeng Zhen (曾真), Wang Cunwen (王存文), Wu Yuanxin (吳元欣). Shortcut method of stripping binary batch distillation of under constant residual fraction composition [J].Journal of Wuhan Institute of Technology(武漢工程大學學報), 2014, 36 (1): 9-13.
[21] Li Qingyang (李慶楊), Wang Chaoneng (王超能), Yi Dayi (易大義). Numeric Analysis (數(shù)值分析) [M]. Wuhan: Huazhong University of Technology Press, 1982.
[22] Xue Yi (薛毅). Optimization Theory and Methods(最優(yōu)化原理與方法) [M]. Beijing: Beijing University of Technology Press, 2004.
Optimization operation of binary stripping batch distillation and its minimum total evaporation
WANG Weiguo1, LUO Jingsong1, ZENG Zhen2, WANG Cunwen3, WU Yuanxin3, QIN Yuanhang3
(1School of Chemical Engineering and Pharmacy,Wuhan Institute of Technology,Wuhan430074,Hubei,China;2School of Mechanical and Electrical Engineering,Wuhan Institute of Technology,Wuhan430074,Hubei,China;3Key Laboratory for Green Chemical Process of Ministry of Education,Wuhan Institute of Technology,Wuhan430074,Hubei,China)
Under ideal operation conditions, the calculation for total evaporation and the minimum total evaporation of binary stripping batch distillation under optimal operations is a constrained function optimization problem. In this work, the constrained function optimization problem is transformed to an unconstrained function optimization problem based on the penalty function method and the fixed double step length factor gradient method is used to solve the unconstrained optimization problem. The calculations show that the fixed double step length factor gradient method has good convergence, and the effect of numerical truncation error accumulation on the calculation results is reduced when the number of segments is large. The energy consumption of binary stripping batch distillation under optimal operation is lower than that under constant residual fraction composition operation. It may be attributed to the following reasons: when the theoretical plate number is relatively small (close to the minimum theoretical plate number required for binary stripping batch distillation under constant residual fraction composition operation), the energy efficiency is improved because of controllable reboil ratio under the optimal operation; when the theoretical plate number is relatively large, the entropy generation in liquidmixing in the top tank is reduced because the material is removed out of the distillation system as soon as possible in the premise of relatively high energy efficiency under the optimal operation because of the controllable reboil ratio. Under ideal operating conditions and infinite number of theoretical plates, the change strategy of the reboil ratio and the calculation formula of the minimum total evaporation of binary stripping batch distillation under optimal operation are obtained by summarizing and extrapolating the above calculation results.
stripping batch distillation; optimal operation; total evaporation; simulation; energy consumption
ZENG Zhen, zengzhen415@126.com
10.11949/j.issn.0438-1157.20141687
TQ 015
:A
:0438—1157(2015)10—4047—14
2014-11-12收到初稿,2015-04-30收到修改稿。
聯(lián)系人:曾真,副教授。
:王為國(1963—),碩士,副教授。
Received date: 2014-11-12.