基于陣列式靜電傳感器的密相氣力輸送煤粉顆粒運(yùn)動(dòng)特性分析

付飛飛 許傳龍 王式民 閆 勇

(1 東南大學(xué)能源熱轉(zhuǎn)換及其過程測控教育部重點(diǎn)實(shí)驗(yàn)室，南京210096)

(2 華北電力大學(xué)控制與計(jì)算機(jī)工程學(xué)院，北京102206)

超濃相氣力輸送系統(tǒng)具有復(fù)雜的非線性、非平衡動(dòng)力學(xué)特性.由于缺乏對濃相氣固兩相流流動(dòng)動(dòng)力學(xué)的了解，給高壓濃相氣力輸送系統(tǒng)的設(shè)計(jì)及優(yōu)化運(yùn)行帶來了較大的不確定性［1］.濃相氣固兩相流中顆粒運(yùn)動(dòng)特性的獲取有2 種方法:①通過理論分析、數(shù)值模擬方法獲取;②利用氣力輸送系統(tǒng)試驗(yàn)得到的顆粒流動(dòng)信息并結(jié)合非線性理論獲取.近似熵(approximate entropy，ApEn)作為一種非線性分析方法，從統(tǒng)計(jì)學(xué)的角度來判斷時(shí)間過程的復(fù)雜性，可以用來表征動(dòng)力學(xué)系統(tǒng)的差異或變化，近些年來，已經(jīng)廣泛應(yīng)用于流化床流型識(shí)別、氣液兩相流電導(dǎo)波動(dòng)信號(hào)復(fù)雜性測度分析以及故障檢測等［2-4］.

粉體氣力輸送過程中，由于顆粒之間以及顆粒與管壁之間的碰撞、分離，顆粒上將產(chǎn)生靜電.粉體顆粒靜電載有大量的動(dòng)態(tài)信息，是流動(dòng)特性、輸送形態(tài)、輸送管道幾何特性和粉體顆粒等因素內(nèi)在作用的綜合體現(xiàn)［5-7］.近些年來發(fā)展起來的環(huán)狀靜電傳感器在氣固兩相流檢測與控制中引起了極大關(guān)注［8-11］.對于氣固兩相流中顆粒分布及流型的測量，環(huán)狀靜電傳感器卻存在不足，這是由于環(huán)狀靜電傳感器無法獲得顆粒在管道截面的位置及空間局部區(qū)域的流動(dòng)信息.通過在氣固兩相流系統(tǒng)外部空間位置上布置陣列式靜電傳感器，獲取其臨近區(qū)域的顆粒信息［12］，這些顆粒靜電信息從不同的時(shí)間和空間反應(yīng)了稠密氣固流動(dòng)特性.結(jié)合適當(dāng)?shù)男畔⑻幚硎侄?，可以全面地獲得顆粒的運(yùn)動(dòng)及分布信息.電容層析成像(ECT)技術(shù)是一種多相流檢測技術(shù)，可用于輸送管道中的氣固、氣液兩相流等連續(xù)相為電介質(zhì)的工業(yè)過程成像.ECT 技術(shù)在兩相流流動(dòng)狀態(tài)檢測和參數(shù)檢測等方面的應(yīng)用研究已取得一定的成果［13-14］.將ECT 技術(shù)與陣列式靜電傳感器系統(tǒng)二者結(jié)合，可深入地了解加壓稠密氣固兩相流顆粒流動(dòng)狀態(tài)及其變化規(guī)律.

本文首先系統(tǒng)分析了陣列式靜電傳感器的局部靈敏度分布特性，并在滿管流和1/4 流2 種流型下進(jìn)行了實(shí)驗(yàn)驗(yàn)證.之后在高固氣比煤粉輸送裝置上，結(jié)合陣列式靜電傳感器和電容層析成像技術(shù)獲得了3 種表觀氣速(5.7，6.9，8.1 m/s)工況下的靜電輸出信號(hào)和管道截面煤粉分布ECT 成像圖.結(jié)合ECT 成像圖，利用FFT 變換和近似熵方法對傳感器輸出信號(hào)進(jìn)行了分析.

1 陣列式靜電傳感器局部傳感特性分析

圖1為陣列式靜電傳感器探頭結(jié)構(gòu)，由8 個(gè)電極片、2 個(gè)接地圓環(huán)保護(hù)極板、絕緣管和屏蔽罩等部分組成.其中，8 個(gè)電極組成的環(huán)狀陣列及圓環(huán)保護(hù)極板緊貼在絕緣管的外壁.

圖1 靜電傳感器探頭結(jié)構(gòu)

陣列式傳感器的數(shù)學(xué)模型在文獻(xiàn)［12］中已有詳細(xì)介紹.陣列式靜電傳感器結(jié)構(gòu)為三維形式，較為復(fù)雜，目前尚難以獲得模型中靜電場的解析解，通常采用數(shù)值法進(jìn)行求解.近些年發(fā)展起來的各種電磁場計(jì)算軟件包使得復(fù)雜的靜電場計(jì)算問題成為可能，本文利用軟件Ansoft 完成陣列式靜電傳感器內(nèi)部靜電場和極片上感應(yīng)電量的計(jì)算.

陣列式靜電傳感器的靈敏度定義為在敏感空間內(nèi)，單位點(diǎn)電荷作用下電極上感應(yīng)電量的絕對值.從陣列式靜電傳感器的有限元模型可知，靈敏度與點(diǎn)電荷所在的敏感空間位置(x，y，z)有關(guān)，即

式中，s(x，y，z)為靈敏度;q 為點(diǎn)電荷帶電量;Q為點(diǎn)電荷電量為q 時(shí)電極上的感應(yīng)電量.考慮靜電傳感器陣列電極的軸對稱性及三維有限元仿真計(jì)算量大等原因，在本實(shí)驗(yàn)中模擬5 條較有代表性的點(diǎn)電荷流線(記為a，b，c，d，e，見圖2)上的靈敏度分布，并以此來表征陣列式靜電傳感器的軸向靈敏場分布特性.模型中的管道內(nèi)徑設(shè)為10 mm，5 條流線在徑向位置上等間隔(間距2 mm)分布，其中c 過圓點(diǎn).8 個(gè)極片的軸向?qū)挾葹?0 mm，長為6 mm.已知帶電顆粒在流線上沿軸向移動(dòng)時(shí)，8 個(gè)極片會(huì)同時(shí)產(chǎn)生感應(yīng)電量.但為便于分析，只計(jì)算電極片1 上的靈敏度分布來表征陣列式靜電傳感器的軸向靈敏場分布特性，其余7 個(gè)電極片都可以旋轉(zhuǎn)到極片1 的位置而得到同樣的靈敏度分布.

圖2 8 極片陣列式靜電傳感器探頭橫截面

圖3為陣列式靜電傳感器敏感空間內(nèi)靈敏度沿軸向的分布.圖4為陣列式靜電傳感器在中心截面(z =0)上的靈敏度分布.從圖3可看出，沿軸向位置的增大，同一流線上的靈敏度值逐漸減小，且流線位置離極片越近，其上的靈敏度變化越快.從圖4可看出，靠近極片壁面處紅色區(qū)域代表的靈敏度值最大，離極片越遠(yuǎn)靈敏度值越小.藍(lán)色代表的最小值只有紅色代表的最大值的0.1 倍.結(jié)合靈敏度的軸向分布和截面分布，陣列式靜電傳感器單電極片靈敏場的空間分布具有局部敏感特性.

圖3 陣列式靜電傳感器靈敏度沿軸向的分布

圖4 陣列式靜電傳感器中心截面靈敏度分布

2 近似熵理論

近似熵方法是一種度量時(shí)間序列隨機(jī)性的非線性統(tǒng)計(jì)方法［15］.對數(shù)據(jù)長度為N 的時(shí)間序列{x(i)}，重構(gòu)m 維相空間，其中的矢量y(i)為

式中，m 又稱為窗口長度.定義矢量y(i)與y(j)之間的距離d［y(i)，y(j)］為各矢量之間的最大距離，即

式中，Φm(r)表示矢量序列{y(i)}的平均自相關(guān)程度.m 越大，相空間中的狀態(tài)點(diǎn)越稀，關(guān)聯(lián)機(jī)會(huì)越小，因而(r)和Φm(r)都隨m 增大而變小，對于Φm(r)，其偏差可用下式表示:

對于不同信號(hào)，式(6)的結(jié)果會(huì)有很大差別.白噪聲均勻充滿整個(gè)相空間各維不大，但Φm極小，因而隨著m 增大，Φm(r)－Φm+1(r)→∞.對于規(guī)則信號(hào)，當(dāng)m 增大到一定程度后，吸引子全部被打開，Φm便不再隨m 增大而增大，即接近于0.因此，可以用Φm－Φm+1刻畫系統(tǒng)隨機(jī)性的程度，近似熵值Eap可定義為

熵值越大，表明時(shí)間序列越具有隨機(jī)性或不規(guī)則性，其非周期性越強(qiáng)，復(fù)雜度越高;熵值越小，表明數(shù)據(jù)周期性越強(qiáng)，復(fù)雜度越小.

3 實(shí)驗(yàn)裝置

圖5為高壓濃相氣力輸送實(shí)驗(yàn)系統(tǒng)，主要由氣源系統(tǒng)、儲(chǔ)料罐系統(tǒng)、輸送管道系統(tǒng)、測量傳感器及數(shù)據(jù)采集系統(tǒng)4 部分組成.高壓氣體經(jīng)過緩沖罐分3 路經(jīng)流量計(jì)向輸送系統(tǒng)提供沖壓風(fēng)、流化風(fēng)和補(bǔ)充風(fēng).發(fā)料罐和收料罐容積均為0.648 m3，輸送管為內(nèi)徑10 mm 的不銹鋼管，總長53 m.接收罐的支座上布置了CHJL2 型荷重變送器.流化風(fēng)、充壓風(fēng)和補(bǔ)充風(fēng)采用AM-1521Q 型金屬管浮子流量計(jì)測量.緩沖罐氣源維持在4 MPa 左右.輸送載氣為CO2，輸送煤粉為云南褐煤.

ECT 系統(tǒng)由傳感器探頭、數(shù)據(jù)采集系統(tǒng)和圖像重建程序3 個(gè)基本部分組成，最高成像速度為74 幀/s，分辨率為32 ×32 像素.ECT 系統(tǒng)與陣列式靜電傳感器并列安裝在同一水平輸送管道中.靜電傳感器的結(jié)構(gòu)尺寸與上述數(shù)值模型的結(jié)構(gòu)尺寸一致.靜電傳感器系統(tǒng)的數(shù)據(jù)采集部分由電流放大電路、PCI9112 采集卡和計(jì)算機(jī)構(gòu)成.在軟件的控制下，可以實(shí)現(xiàn)8 通道同時(shí)采集，本實(shí)驗(yàn)中靜電傳感器輸出信號(hào)的頻率范圍在500 Hz 以內(nèi)，因此每通道的采樣頻率設(shè)定為1 000 Hz.

圖5 高壓濃相氣力輸送實(shí)驗(yàn)臺(tái)系統(tǒng)圖

4 實(shí)驗(yàn)結(jié)果與討論

4.1 陣列式靜電傳感器局部傳感特性驗(yàn)證

為了驗(yàn)證陣列式靜電傳感器的局部傳感特性，設(shè)計(jì)了簡易的重力輸送試驗(yàn)裝置.顆粒在自身重力的作用下由料斗進(jìn)入管道，再進(jìn)入靜電傳感器和電容傳感器測試探頭.輸送管段可垂直或者傾斜放置，可實(shí)現(xiàn)1/4 流和滿管流2 種流型.實(shí)驗(yàn)用的固體顆粒是直徑為1 mm 的玻璃珠，管道半徑為25 mm.圖6為1/4 流ECT 成像及能量分布(第2，3 和4 s 的靜電數(shù)據(jù))，其中極片相對位置如圖2所示.由圖6(a)可見，截面煤粉所占面積約為管道截面面積的1/4，且極片3 和4 附近顆粒濃度最大(深色代表該處顆粒密集).從圖6(b)可看出，極片3 和4 上的靜電信號(hào)能量大于其他6 個(gè)極片上的信號(hào)能量.這證明了陣列傳感器極片的靈敏度具有局域性，靠近極片的區(qū)域具有較高的靈敏度.圖7為滿管流的ECT 成像圖及8 極片的靜電信號(hào)能量分布.與圖7對比發(fā)現(xiàn)，靜電信號(hào)能量值波動(dòng)平緩并在能量均值(值為1.94)處上下波動(dòng)，盡管極片7 附近的顆粒濃度相對較低(見圖7(b))導(dǎo)致了其上的靜電信號(hào)能量值較小.因此，陣列式靜電傳感器的輸出信號(hào)的能量分布可有效地反應(yīng)管道截面上顆粒的分布特性.

圖6 1/4 流ECT 成像及能量分布

圖7 滿管流ECT 成像及能量分布

4.2 氣力輸送實(shí)驗(yàn)

在表觀氣速為5.7，6.9 和8.1 m/s 的3 個(gè)工況下進(jìn)行了氣力輸送實(shí)驗(yàn).圖8為3 個(gè)工況下管道截面煤粉分布的ECT 成像圖.在水平管氣力輸送過程中，因受重力影響，煤粉的濃度分布呈不對稱性，管下部濃度較高，上部濃度較低.從成像圖中可看出，3 個(gè)工況下煤粉的輸送連續(xù)性均較好，煤粉在管道截面的分布隨時(shí)間變化較小.工況3 下濃相區(qū)與稀相區(qū)幾乎沒有分界，此時(shí)煤粉的懸浮性及分布均勻性最好.但是ECT 成像圖只能反映管道截面煤粉的分布特性以及顆粒的宏觀運(yùn)動(dòng)特性，卻無法反映顆粒隨機(jī)運(yùn)動(dòng)的本質(zhì)特征.

圖8 3 種工況下的ECT 成像圖

靜電傳感器輸出信號(hào)頻域特征與其敏感區(qū)內(nèi)顆粒帶電量、帶電顆粒軸向平均速度和顆粒徑向分布的變化密切相關(guān).顆粒帶電量增加導(dǎo)致頻域幅值增加，帶電顆粒軸向速度增加導(dǎo)致主峰頻率值增大，而顆粒的徑向運(yùn)動(dòng)也會(huì)對頻率特性產(chǎn)生非線性的影響［11-12］.顆粒在管道中運(yùn)動(dòng)時(shí)，受到氣流的曳引作用、顆粒與顆粒及顆粒與管壁間的碰撞作用力以及自身重力等，這些力的共同作用使得顆粒的運(yùn)動(dòng)非常復(fù)雜，近似為隨機(jī)運(yùn)動(dòng).近似熵值可以定量分析隨機(jī)運(yùn)動(dòng)的程度.近似熵值越大，顆粒運(yùn)動(dòng)的復(fù)雜性越大，隨機(jī)性越強(qiáng).表1為8 極片上近似熵值、主峰值頻率及主峰值頻率對應(yīng)幅值隨表觀氣速的變化規(guī)律.可看出，隨表觀氣速的增大，主峰值頻率對應(yīng)的幅值大致有減小的趨勢.這是因?yàn)楸碛^氣速增大時(shí)，粉體的質(zhì)量流量減小，碰撞產(chǎn)生的靜電總量減小，導(dǎo)致傳感器輸出信號(hào)的能量減小.但近似熵值和主峰值頻率隨表觀氣速的增大沒有一致的變化規(guī)律.

表1 陣列傳感器8 極片上靜電輸出信號(hào)隨表觀氣速增大的規(guī)律

從圖8中看到，極片2 和3 處于濃相區(qū)，極片6和7 處于稀相區(qū)，極片1，4，5，8 處于波動(dòng)區(qū)域.由于波動(dòng)區(qū)顆粒分布隨機(jī)變化，極片1，4，5，8 上靜電信號(hào)的分析結(jié)果不確定，對比性差(見表1).因此本文中采用濃相區(qū)極片2，3 和稀相區(qū)極片6，7 上靜電信號(hào)的頻域及時(shí)域特征來描述濃相區(qū)和稀相區(qū)顆粒運(yùn)動(dòng)特征.

圖9為稀相和密相區(qū)靜電信號(hào)主峰值頻率與表觀氣速之間的關(guān)系.隨著表觀氣速的增加，極片2，3 上的靜電信號(hào)主峰值頻率一致增大，表明濃相區(qū)的顆粒軸向平均速度增大，顆粒徑向運(yùn)動(dòng)波動(dòng)相對不劇烈并且其對輸出信號(hào)頻譜特性造成的非線性影響不明顯.極片6，7 上的靜電信號(hào)主峰值的頻率雖然總體有增大的趨勢，但其間有減小的變化過程，說明稀相區(qū)顆粒的徑向運(yùn)動(dòng)相對劇烈，導(dǎo)致傳感器輸出信號(hào)的主峰值頻率呈非線性增長.在湍流運(yùn)動(dòng)中，顆粒濃度與顆粒的湍動(dòng)能有關(guān)［16］.當(dāng)表觀氣速為5.7 m/s 時(shí)，雷諾數(shù)(Re)大于2 ×105，已經(jīng)達(dá)到湍流狀態(tài).湍動(dòng)能越大，表明顆粒在徑向上的波動(dòng)越劇烈.盡管當(dāng)表觀氣速增大時(shí)，顆粒速度增大是主要的變化，但稀相區(qū)顆粒速度在徑向上的波動(dòng)也足夠明顯，以至于其對靜電信號(hào)主峰值頻率的增長造成了非線性的影響.

圖9 靜電信號(hào)主峰值頻率隨表觀氣速的變化規(guī)律

圖10為極片2，3 和極片6，7 上的輸出信號(hào)近似熵值與表觀氣速之間的關(guān)系.表觀氣速增大時(shí)，極片2，3 上輸出信號(hào)的近似熵值一致減小而極片6，7 上輸出信號(hào)的近似熵值一致增大，但前者始終大于后者.密相區(qū)的顆粒濃度大，顆粒碰撞頻率大于稀相區(qū).實(shí)際上，靜電傳感器輸出信號(hào)的近似熵值不完全表征顆粒運(yùn)動(dòng)的復(fù)雜性，同時(shí)它還與顆粒電荷在軸向和徑向上的分布梯度以及電荷隨時(shí)間的變化率等因素有關(guān).顯然，在濃相區(qū)，顆粒電荷分布的梯度以及隨時(shí)間的變化率更隨機(jī).盡管稀相區(qū)的顆粒的湍流會(huì)導(dǎo)致顆粒運(yùn)動(dòng)復(fù)雜，但是其靜電信號(hào)的復(fù)雜性還是低于濃相區(qū).但當(dāng)表觀氣速增大時(shí)，濃相區(qū)顆粒濃度減小，顆粒碰撞頻率降低，導(dǎo)致稀相區(qū)和濃相區(qū)的顆粒運(yùn)動(dòng)復(fù)雜性之間的差距減小.對照ECT 成像圖，還發(fā)現(xiàn)管道煤粉輸送的懸浮性與極片2，3 和極片6，7 上輸出信號(hào)的近似熵值之差有對應(yīng)關(guān)系:差值越小，煤粉的輸送懸浮性越好.

圖10 靜電信號(hào)近似熵值隨表觀氣速的變化規(guī)律

5 結(jié)語

在加壓密相氣力輸送實(shí)驗(yàn)裝置上，結(jié)合ECT技術(shù)，利用FFT 變換和近似熵方法分析了密相氣固兩相流陣列靜電波動(dòng)信號(hào).結(jié)果表明:陣列式靜電傳感器可利用自身的陣列結(jié)構(gòu)及靈敏度分布特性，獲取陣列中各個(gè)單傳感器靈敏區(qū)的顆粒流動(dòng)的靜電信息，從而全面地得到顆粒的運(yùn)動(dòng)及分布信息.當(dāng)表觀氣速增大時(shí)，濃相區(qū)靜電信號(hào)的主峰值頻率一致增大，而稀相區(qū)靜電信號(hào)的主峰值頻率非線性增長，表明表觀氣速增大導(dǎo)致顆粒運(yùn)動(dòng)速度增大，并且稀相區(qū)顆粒徑向速度波動(dòng)明顯;稀相區(qū)和濃相區(qū)的顆粒運(yùn)動(dòng)復(fù)雜性程度越接近，煤粉的輸送懸浮性越好.同時(shí)，隨表觀氣速的增大，濃相區(qū)靜電信號(hào)的近似熵值減小，稀相區(qū)靜電輸出信號(hào)的近似熵值增大，但濃相區(qū)的近似熵值一直大于稀相區(qū)的近似熵值，對應(yīng)ECT 成像圖，發(fā)現(xiàn)隨兩者差值的減小，管道內(nèi)煤粉分布越均勻，輸送懸浮性越好.

References)

［1］梁財(cái).高壓超濃相煤粉氣力輸送流動(dòng)特性研究［D］.南京:東南大學(xué)能源與環(huán)境學(xué)院，2007.

［2］楊進(jìn)，文玉梅，李平.基于相關(guān)分析和近似熵的管道泄漏聲信號(hào)特征提取及辨識(shí)方法［J］.儀器儀表學(xué)報(bào)，2009，30(2):272-278.

Yang Jin，Wen Yumei，Li Ping.Feature extraction and identification of leak acoustic signal in water distribution pipelines using correlation analysis and approximate entropy［J］.Chinese Journal of Scientific Instrument，2009，30(2):272-278.(in Chinese)

［3］栗然，陸鳳怡，徐宏銳，等.基于局域波與近似熵的負(fù)荷分析方法［J］.中國電機(jī)工程學(xué)報(bào)，2010，30(25):51-58.

Li Ran，Lu Fengyi，Xu Hongrui，et al.Novel approach for load analysis based on local wave and approximate entropy［J］.Proceedings of the CSEE，2010，30(25):51-58.(in Chinese)

［4］林萍，王曉萍，陳伯川，等.近似熵的性質(zhì)及其在氣固流化床復(fù)雜性研究中的應(yīng)用［J］.高?；瘜W(xué)工程學(xué)報(bào)，2004，18(3):281-286.

Lin Ping，Wang Xiaoping，Chen Bochuan，et al.Character of approximate entropy and its application in the complexity measurement of gas-solid fluidized bed［J］.Journal of Chemical Engineering of Chinese Universities，2004，18(3):281-286.(in Chinese)

［5］Matsusaka S，Masuda H.Electrostatics of particles［J］.Advanced Powder Technology，2003，14(2):169-179.

［6］Nifuku M，Katoh H.A study on the static electrification of powders during pneumatic transportation and the ignition of dust clouds［J］.Powder Technology，2003，135-136:234-242.

［7］Masuda H，Komatsu T，Mitsui N，et al.Electrification of gas-solid suspensions flowing in steel and insulatingcoated pipes［J］.Journal of Electrostatics，1977，2(4):341-350.

［8］Qian Xiangchen，Yan Yong，Shao Jiaqing，et al.Quantitative characterization of pulverised coal and biomass-coal blends in pneumatic conveying pipelines using electrostatic sensor arrays and data fusion techniques［J］.Measurement Science and Technology，2012，23(8):1-8.

［9］Gajewski J B.Non-contact electrostatic flow probes for measuring the flow rate and charge in the two-phase gassolid flows［J］.Chemical Engineering Science，2006，61(7):2262-2270.

［10］Gajewsk J B.Electrostatic nonintrusive method for measuring the electric charge，mass flow rate，and velocity of particulates in the two-phase gas-solid pipe flows—its only or as many as 50 years of historical evolution［J］.IEEE Transactions on Industry Applications，2008，44(5):1418-1430.

［11］Xu Chuanlong，Wang Shimin.HHT analysis of electrostatic fluctuation signals in dense-phase pneumatic conveying of pulverized coal at high pressure［J］.Chemical Engineering Science，2010，65(4):1334-1344.

［12］Xu Chuanlong，Li Jian，Wang Shimin.Investigations into sensing characteristics of electrostatic sensor arrays through computational modeling and practical experimentation［J］.Journal of Electrostatics，2012，70(1):60-71.

［13］楊道業(yè)，周賓，王式民.電容層析成像在高壓濃相煤粉氣力輸送中的應(yīng)用［J］.化工學(xué)報(bào)，2009，60(4):892-897.

Yang Daoye，Zhou Bin，Wang Shimin.Application of electrical capacitance tomography in dense-phase pneumatic conveying of pulverized coal under high pressure［J］.Journal of the Chemical Industry and Engineering Society of China，2009，60(4):892-897.(in Chinese)

［14］Cong Xingliang，Guo Xiaolei，Gong Xin，et al.Experimental research of flow patterns and pressure signals in horizontal dense phase pneumatic conveying of pulverized coal［J］.Powder Technology，2011，208(3):600-609.

［15］Pincus S.Approximate entropy(ApEn)as a complexity measure［J］.CHOAS，1995，5(1):110-117.

［16］Sinclair J L，Jackson R.Gas-particle flow in a vertical pipe with particle-particle interactions［J］.AIChE Journal，1989，35(9):1473-1486.