基于熱輻射成像建模求解的管式爐爐管溫度檢測

張向宇，鄭樹，周懷春，徐宏杰

?

基于熱輻射成像建模求解的管式爐爐管溫度檢測

張向宇1，鄭樹2，周懷春2，徐宏杰1

（1西安熱工研究院有限公司電站鍋爐煤清潔燃燒國家工程研究中心，陜西西安 710032；2清華大學(xué)熱能動力仿真與控制研究所，北京 100084）

準(zhǔn)確在線檢測管壁表面分布式溫度是優(yōu)化管式爐加熱工藝的關(guān)鍵所在。以一臺工業(yè)管式裂解爐為試驗(yàn)對象，結(jié)合輻射圖像處理方法，開展了管式爐爐管表面溫度可視化檢測研究。采用基于Monte Carlo的DRESOR法求解具有復(fù)雜邊界條件的管式爐輻射成像模型，實(shí)現(xiàn)了爐管輻射與火焰輻射、爐壁輻射的解耦計(jì)算，對爐管表面溫度與熱通量分布進(jìn)行了在線監(jiān)測，并研究了二者隨工質(zhì)流動方向的變化趨勢。經(jīng)過驗(yàn)證，溫度測量誤差小于2%，測量誤差主要出現(xiàn)在最高和最低溫區(qū)域。該項(xiàng)研究將有助于指導(dǎo)管式爐燃燒調(diào)整，改進(jìn)加熱工藝，提高爐管表面受熱均勻性，延長爐管工作壽命。

管式爐；輻射；爐管；建模；熱通量

引 言

管式加熱爐種類繁多，廣泛應(yīng)用于制氫、裂解、重整等石化工藝中，是石油化工行業(yè)的關(guān)鍵設(shè)備，也是主要能耗設(shè)備，其設(shè)計(jì)和操作直接影響整套制烴裝置的經(jīng)濟(jì)性。以管式裂解爐為例，其主要通過間壁加熱的方式來為裂解過程提供能量，熱量經(jīng)由輻射管段的高溫管壁傳遞給裂解物料，裂解溫度受到爐管的內(nèi)外壁溫度限制；此外，由于管內(nèi)裂解時(shí)的結(jié)焦過程使得熱阻和管壁溫度都會增加，造成流動阻力降增加[1-2]。因此，準(zhǔn)確檢測管壁表面分布式溫度是改進(jìn)管式爐裂解技術(shù)的關(guān)鍵所在。管式爐內(nèi)的熱量傳遞主要通過輻射傳熱進(jìn)行，為優(yōu)化加熱工藝，保證爐管表面受熱均勻，需要深入分析管式爐內(nèi)輻射換熱過程，精確定量輻射傳遞過程中的溫度場、熱通量等熱物理量場[3-4]。

管式爐作為給定幾何形狀的受限空間，求解爐內(nèi)輻射換熱方程時(shí)必須給定輻射在空間邊界上的傳輸特征，但數(shù)量眾多、排列密集、材料特殊、結(jié)構(gòu)異型爐管的存在使得爐內(nèi)輻射傳遞邊界條件變得極其復(fù)雜，再加上高溫爐壁輻射、火焰輻射、高溫爐管表面輻射的相互耦合，精確求解管式爐內(nèi)輻射傳遞過程將非常困難。

相關(guān)研究主要通過實(shí)驗(yàn)分析和數(shù)值模擬來建立管式爐內(nèi)化學(xué)反應(yīng)動力學(xué)模型，控制裂解產(chǎn)物生成速率[5]。Froment等[6-8]利用計(jì)算流體力學(xué)軟件（CFD）對丙烷裂解爐內(nèi)的質(zhì)量、動量和能量傳遞過程進(jìn)行了模擬研究。Mehdi等[9]研究了裂解爐內(nèi)氣體溫度對產(chǎn)物收率的影響。Masoumi等[10]建立了裂解爐動態(tài)操作模型，并結(jié)合小型化試驗(yàn)臺對反應(yīng)器的最佳溫度分布進(jìn)行了推導(dǎo)，提出了最佳反應(yīng)溫度區(qū)間。韋劉軻等[11-12]采用數(shù)值模擬方法對乙烯裂解爐輻射段爐管溫度的分布和均勻性問題進(jìn)行了分析。

輻射圖像處理方法被嘗試用于在線測量一臺小型管式試驗(yàn)爐三維溫度場[13]，其屬于輻射逆問題的研究范疇。研究人員對試驗(yàn)爐內(nèi)布置的風(fēng)管表面的溫度分布和輻射率進(jìn)行了測量研究，初步分析了圓形爐管表面輻射傳遞計(jì)算的差異性[14]。

本文以一臺工業(yè)管式裂解爐上為試驗(yàn)對象，采用Monte Carlo方法精確求解爐內(nèi)復(fù)雜邊界條件下的輻射傳遞方程，建立包含爐膛壁面、爐管和火焰的輻射成像模型，并結(jié)合輻射圖像處理方法在線檢測爐管表面溫度分布。

1 原理簡述

管式爐熱輻射成像模型如圖1所示，布置在爐膛壁面上的電荷耦合器件（CCD）攝像機(jī)接收到的其實(shí)是高方向分辨率的邊界輻射強(qiáng)度圖像。

圖1 管式爐熱輻射成像模型

作為封閉空腔內(nèi)的輻射傳遞問題，輻射成像模型的本質(zhì)是求解輻射在空間邊界上的傳輸特征，從而建立邊界輻射強(qiáng)度分布與空間溫度分布之間的對應(yīng)關(guān)系，而這需要精確求解管式爐內(nèi)輻射傳遞方程[15]。管式爐內(nèi)輻射傳遞過程的復(fù)雜性在于同時(shí)存在平面爐壁和圓形爐管兩種壁面，而這兩者在輻射傳遞方式上卻不盡相同。

Monte Carlo方法的優(yōu)勢在于處理復(fù)雜的熱輻射問題，這里采用基于Monte Carlo的DRESOR法[16]來求解復(fù)雜邊界條件下的管式爐輻射成像模型。根據(jù)圖1，視線方向上的輻射強(qiáng)度由式(1)確定

式(1)中一共包含6項(xiàng)，其中第1項(xiàng)表示空間介質(zhì)單元的直接發(fā)射，第2項(xiàng)表示爐膛壁面直接發(fā)射的輻射能量，第3項(xiàng)表示爐管表面的直接發(fā)射，第4項(xiàng)表示爐壁發(fā)射分別被爐壁、管壁反射和介質(zhì)單元散射的能量，第5項(xiàng)表示爐管發(fā)射分別被爐壁、管壁反射和介質(zhì)單元散射的能量，第6項(xiàng)表示介質(zhì)單元發(fā)射分別被爐壁、管壁反射和介質(zhì)單元散射的能量。

對式(1)中的積分項(xiàng)進(jìn)行離散化，將介質(zhì)單元?jiǎng)澐譃閭€(gè)網(wǎng)格，爐壁劃分為個(gè)網(wǎng)格，所有爐管劃分為個(gè)網(wǎng)格，CCD靶面劃分為個(gè)像素單元。根據(jù)維恩輻射定律，將源項(xiàng)4換成單色輻射源項(xiàng)15exp(2)得到式(2)，計(jì)算波長取值為610 nm[17]，對應(yīng)所用CCD攝像機(jī)的R通道特征波長。由于采用可見光CCD攝像機(jī)，其光譜響應(yīng)的波長范圍為500～800 nm，主要接收的是火焰中碳黑顆粒的輻射，三原子氣體輻射只占很小份額。另外，管式爐內(nèi)氣體燃料火焰中顆粒物濃度較低，散射對爐內(nèi)輻射傳遞過程的影響很小，可忽略其各向異性。因此假定介質(zhì)單元為灰性吸收、各向同性散射，爐壁及爐管均為灰性發(fā)射、吸收、漫反射表面[18]。

式(2)用矩陣形式表示如下

式(3)～式(5)即為管式爐內(nèi)輻射成像模型。從相機(jī)拍攝的爐內(nèi)火焰圖像中提取邊界單色輻射強(qiáng)度分布后，需要對式(3)進(jìn)行逆求解得到源項(xiàng)分布，這屬于輻射傳遞逆問題。輻射逆問題求解的難點(diǎn)在于未知量大于檢測量所引起的方程病態(tài)。

一種修正的Tikhonov正則化方法[19]被用來求解該類病態(tài)方程，其基本原理是找到一個(gè)()使得式(6)極小化

()由式(7)確定

爐膛溫度分布可以由()計(jì)算得到

爐管表面離散單元的能量方程由式(9)給出

式(9)中，左邊第1項(xiàng)是爐管單元輻射出的能量，第2項(xiàng)是進(jìn)入此爐管單元的輻射凈熱量，即局部壁面輻射熱流。由于假定了平衡條件，爐管單元流入的能量等于該單元流出的能量。右邊第1項(xiàng)對應(yīng)于被爐管單元所吸收的來自所有介質(zhì)單元的輻射能，第2項(xiàng)是被爐管單元所吸收的來自所有爐壁單元的輻射能，第3項(xiàng)是被爐管單元所吸收的來自其他所有爐管單元的輻射能。獲得全爐膛溫度分布后，即可根據(jù)式(9)獲得爐管壁面輻射熱流分布。

2 計(jì)算模型

考察一臺工業(yè)管式裂解爐，其工藝為首先在對流段將反應(yīng)管內(nèi)烴和水蒸氣的混合物預(yù)熱到開始裂解的溫度，接著將混合物輸送到高溫的輻射管段繼續(xù)加熱升溫使它們發(fā)生裂解反應(yīng)。對每個(gè)輻射段，物料基本走向?yàn)椤吧线M(jìn)上出”。在第1管程內(nèi)，氣態(tài)物料從橫跨管出來自輻射室的頂部進(jìn)入輻射段，接著向下，經(jīng)過U形彎管從而進(jìn)入第2管程。在第2管程內(nèi)，物料以自下向上的方向流動。

本文主要研究裂解爐輻射段內(nèi)的爐管溫度檢測。計(jì)算模型如圖2所示。

圖2 溫度場可視化檢測模型

模型尺寸嚴(yán)格參照實(shí)際工業(yè)裂解爐確定，輻射段近似為長方體型，尺寸為14000 mm×3470 mm× 10688 mm。16支CCD攝像機(jī)被安裝在爐膛壁面上拍攝爐內(nèi)火焰圖像。攝像機(jī)分兩層布置，布置高度分別為5200 mm和11000 mm，每層8支，其中每個(gè)側(cè)墻上各安裝2支，以保證視場角覆蓋整個(gè)爐膛空間及所有爐管。輻射段內(nèi)排布了8組構(gòu)型相同的爐管，接收火焰輻射以提供管內(nèi)裂解原料反應(yīng)所需的熱量。爐管垂直懸吊在爐膛中部，單排布置，爐管根數(shù)96根，分為8組，其中每組12根，管間距70 mm。

對計(jì)算模型進(jìn)行網(wǎng)格劃分，如圖3所示。為簡化計(jì)算，忽略橫跨段爐管。首先劃分空間單元網(wǎng)格，（寬度）方向劃分為15個(gè)網(wǎng)格，并保證每組爐管包含在一個(gè)空間網(wǎng)格內(nèi)；（深度）方向劃分為8個(gè)網(wǎng)格，（高度）方向劃分為12個(gè)網(wǎng)格。這樣介質(zhì)單元網(wǎng)格總數(shù)1440個(gè)，爐壁單元網(wǎng)格總數(shù)792個(gè)。

圖3 網(wǎng)格劃分示意圖

在每組爐管內(nèi)再對每根爐管進(jìn)行網(wǎng)格細(xì)分，對每根爐管沿方向劃分1個(gè)網(wǎng)格，方向劃分2個(gè)網(wǎng)格，方向劃分11個(gè)網(wǎng)格，這樣每根爐管劃分22個(gè)網(wǎng)格，爐管網(wǎng)格總數(shù)為2112個(gè)。

3 試驗(yàn)結(jié)果與分析

在一臺工業(yè)管式裂解爐上開展了爐管表面溫度可視化檢測試驗(yàn)研究。首先利用CCD攝像機(jī)拍攝爐內(nèi)火焰圖像，經(jīng)視頻分割器合成后送入工控機(jī)進(jìn)行處理。火焰圖像如圖4所示，從上到下、從左到右編號為1～16，對應(yīng)的相機(jī)快門速度為CCD1～CCD11為1/250 s，CCD12～CCD16為1/500 s。

圖4 火焰圖像

在工控機(jī)中對火焰圖像進(jìn)行處理，提取每個(gè)像素單元的紅色譜色值，并根據(jù)標(biāo)定數(shù)據(jù)轉(zhuǎn)化為紅色波長下的單色輻射強(qiáng)度分布。根據(jù)之前的研究結(jié)果，假定介質(zhì)吸收系數(shù)為0.2 m-1，散射系數(shù)s為0.1 m-1，爐壁發(fā)射率w為0.8，爐管表面發(fā)射率p為0.75[20]，這樣通過式(7)和式(8)可實(shí)現(xiàn)全爐膛溫度分布的在線檢測。圖4火焰圖像對應(yīng)工況下的爐管表面溫度分布如圖5所示。

圖5 爐管表面溫度云圖

這里根據(jù)管內(nèi)工質(zhì)的流向?qū)t管劃分為4組，每組爐管均完成一個(gè)獨(dú)立的工質(zhì)循環(huán)。從云圖中可以看出，由于燃燒器布置在爐膛底部，下部爐管的表面溫度明顯高于上部，并且第Ⅲ組爐管的加熱強(qiáng)度要大于其他幾組爐管。爐管表面溫度范圍為900～1100℃。

進(jìn)一步研究爐管表面溫度隨爐膛高度的變化趨勢，如圖6所示。在第1管程內(nèi)，隨著高度下降，爐管表面溫度逐漸上升，上升速度先慢后快；進(jìn)入第2管程后，工質(zhì)沿爐管向上運(yùn)動，爐管表面溫度逐漸降低。工質(zhì)出口處的爐管表面溫度略高于入口處。

圖6 爐管表面溫度變化曲線

獲得全爐膛溫度分布后，由式(9)計(jì)算得到爐管壁面輻射熱流分布，如圖7所示。

可以發(fā)現(xiàn)，沿工質(zhì)流動方向，4組爐管的表面熱流均先升高再降低，與表面溫度變化趨勢一致，但在爐膛底部，熱通量有一個(gè)下降過程，這可能是由于火焰中心位于爐膛高度1 m處，其對爐管的加熱能力大于根部火焰。

另外在圖6和圖7中能夠發(fā)現(xiàn)，第Ⅱ組和第Ⅲ組爐管的表面溫度與熱通量均大于第Ⅰ組和第Ⅳ組爐管，其主要原因是爐底中部燃燒器火焰會卷吸兩側(cè)燃燒器火焰中未燃燼的炭黑顆粒及高溫?zé)煔猓胁咳紵骰鹧娴娜紵龔?qiáng)度明顯大于兩側(cè)燃燒器，因此在高溫?zé)煔獾募訜嶙饔孟?，爐膛中部爐管的表面溫度與熱通量將大于兩側(cè)爐管。

采用兩種方法對可視化系統(tǒng)的檢測結(jié)果進(jìn)行驗(yàn)證：①利用重建得到的溫度場與輻射成像矩陣正向計(jì)算探測器接收的單色輻射強(qiáng)度分布，并與實(shí)際火焰圖像中提取的單色輻射強(qiáng)度進(jìn)行對比；②利用裂解爐上安裝的爐管紅外測溫裝置與可視化系統(tǒng)進(jìn)行對比測量，紅外測溫裝置的測點(diǎn)分布如圖3所示。

單色輻射強(qiáng)度的對比曲線如圖8所示，兩條曲線的吻合程度決定了重建誤差的大小。

圖8 單色輻射強(qiáng)度校核曲線

可以發(fā)現(xiàn)，重建強(qiáng)度與測量強(qiáng)度基本吻合，偏差主要出現(xiàn)在峰值區(qū)域，這主要是由于正則化方法對溫度重建的約束作用。

兩種測量裝置的測溫結(jié)果對比如表1所示，最大偏差小于20℃，誤差小于2%。與紅外測溫裝置相比，可視化系統(tǒng)能夠?qū)崿F(xiàn)所有爐管表面溫度的場掃描，從而指導(dǎo)運(yùn)行人員進(jìn)行燃燒調(diào)整。

表1 測量結(jié)果比較

4 結(jié) 論

（1）采用基于Monte Carlo的DRESOR法求解具有復(fù)雜邊界條件的管式爐輻射成像模型，實(shí)現(xiàn)了爐管輻射與火焰輻射、爐壁輻射的解耦計(jì)算，并建立了爐管單元的能量方程。

（2）結(jié)合輻射圖像處理技術(shù)，在一臺工業(yè)管式裂解爐上開展了爐管表面溫度可視化檢測試驗(yàn)研究，分析了爐管表面溫度和熱通量沿工質(zhì)流動方向的變化趨勢。

（3）對爐管表面溫度可視化檢測結(jié)果進(jìn)行了驗(yàn)證，最大偏差小于20℃，誤差小于2%，測量誤差主要出現(xiàn)在最高溫和最低溫區(qū)域。

（4）管式加熱爐爐管表面溫度分布的準(zhǔn)確檢測將有助于指導(dǎo)燃燒調(diào)整，改進(jìn)加熱工藝，提高爐管表面受熱均勻性，延長爐管工作壽命。

符 號 說 明

——單色輻射成像矩陣 C1，C2——Planck常數(shù) D——正則化矩陣 Eλ(T)——溫度T、波長λ下的黑體輻射能，W·m-3 Iλ——攝像機(jī)接收的單色輻射強(qiáng)度，W·m-3·sr-1 I(O,s)——布置在O點(diǎn)的攝像機(jī)在s視線方向上接收的輻射強(qiáng)度，W·m-3·sr-1 k——介質(zhì)單元吸收系數(shù)，m-1 l——能束行進(jìn)的路徑長度，m m，n——分別為介質(zhì)單元、爐壁網(wǎng)格數(shù) P——靶面像素單元 Qp——爐管單元輻射熱流，W·m-2 Rd(?j-i)——DRESOR數(shù)，表示第j個(gè)網(wǎng)格單元發(fā)射出的總輻射能被第i個(gè)成像單元在單位面積、單位角度所接收的份額 ΔSw，ΔSp——分別為爐壁和爐管單元面積，m2 Tv，Tw，Tp——分別為介質(zhì)、爐壁及爐管溫度，K u——爐管網(wǎng)格數(shù) ΔVv——空間單元體積，m3 α——正則化參數(shù) β——消光系數(shù)，m-1 εw，εp——分別為爐壁和爐管表面輻射率 λ——波長，m σ——Boltzmann常數(shù) σs——介質(zhì)單元散射系數(shù)，m-1 下角標(biāo) p——爐管 v——空間介質(zhì) w——爐壁

References

[1] Cai Huilin (蔡惠林), Zhang Zaiyu (張?jiān)阼?, Mao Wenxing (茅文星),. Olefins manufacture from vacuum gas oil in tubular furnace [J].() (化工學(xué)報(bào)), 1979, 30 (2): 165-175

[2] Shokrollahi Yancheshmeh M S, Seifzadeh Haghighi S, Gholipour M R, Dehghani O, Rahimpour M R, Raeissi S. Modeling of ethane pyrolysis process: a study on effects of steam and carbon dioxide on ethylene and hydrogen productions [J]., 2013, 215: 550-560

[3] Wen Liankui (文聯(lián)奎), Zheng Yuanyang (鄭遠(yuǎn)揚(yáng)). Mathematical models and calculation method of radiation heat transfer in furnace [J].(石油學(xué)報(bào)), 1986, 2 (1): 57-69

[4] Guo Shucai (郭樹才). Heat transfer model in tubular furnace [J].() (化工學(xué)報(bào)), 1980, 31 (3): 255-264

[5] Rao Libo (饒利波), Yang Guangjiong (楊光炯), Huang Zuqi (黃祖棋). Mathmatical simulation of fluid flow and heat transfer in cuboid radiation section of petro-chemical tube heater [J].() (化工學(xué)報(bào)), 1993, 44 (4): 410-416

[6] Sundaram K M, Froment G F. Kinetics of coke deposition in the thermal cracking of propane [J]., 1979, 34: 635-644

[7] Hu Guihua, Wang Honggang, Qian Feng. Numerical simulation on flow,combustion and heat transfer of ethylene cracking furnaces [J]., 2011, 66:1600-1611

[8] Lan Xingying (藍(lán)興英). Comprehensive numerical simulation of combustion, heat transfer and pyrolysis reaction process in ethylene cracking furnace [D]. Beijing: Beijing University of Petroleum and Chemical Engineering, 2004

[9] Mehdi Berreni, Wang Meihong. Modelling and dynamic optimization of thermal cracking of propane for ethylene manufacturing [J]., 2011, 35: 2876-2885

[10] Masoumi M E, Sadrameli S M, Towfighi J, Niaei A. Simulation, optimization and control of a thermal cracking furnace [J]., 2006, 31: 516-527

[11] Wei Liuke (韋劉軻). Numerical research on non-uniform heat transfer and its impact on ethylene cracking furnace [D]. Beijing: Beijing University of Chemical Technology, 2010

[12] Zhou Hanzhang (周瀚章). Three dimensional numerical simulation to firebox of ethylene cracking furnace [D]. Beijing: Beijing University of Chemical Technology, 2007

[13] Huang Xuekong (黃學(xué)孔). Three dimensional temperature reconstruction of tubular furnace [D]. Wuhan: Huazhong University of Science and Technology, 2009

[14] Zhang Xiangyu, Cheng Qiang, Lou Chun, Zhou Huaichun. An improved colorimetric method for visualization of 2-D inhomogeneous temperature distribution in a gas fired industrial furnace by radiation image processing [J]., 2011, 33 (1): 2755-2762

[15] Liu D, Huang Q X, Wang F, Chi Y, Cen K F,. Simultaneous measurement of three-dimensional soot temperature and volume fraction fields in axisymmetric or asymmetric small unconfined flames with CCD cameras [J]., 2010, 132: 061202-1

[16] Cheng Q, Zhou H C. The DRESOR method for a collimated irradiation on an isotropically scattering layer [J]., 2007, 129 (5): 634-645

[17] Sun Yipeng, Lou Chun, Zhou Huaichun. A simple judgement method of gray property of flames based on spectral analysis and the two-color method for measurements of temperatures and emissivity [J]., 2011, 33 (1): 735-741

[18] Lou Chun, Li Wenhao, Zhou Huaichun, Carlos T Salinas. Experimental investigation on simultaneous measurement of temperature distributions and radiative properties in an oil-fired tunnel furnace by radiation analysis [J]., 2011, 54:1-8

[19] Zhou H C, Han S D. Simultaneous reconstruction of temperature distribution, absorptivity of wall surface and absorption coefficient of medium in a 2-D furnace system [J]., 2003, 46: 2645-2653

[20] Cheng Qiang, Zhang Xiangyu, Wang Zhichao, Zhou Huaichun, Shao Song. Simultaneous measurement of three-dimensional temperature distributions and radiative properties based on radiation image processing technology in a gas-fired pilot tubular furnace [J]., 2014, 35 (6/7/8): 770-779

Visualization of pipe temperature distribution in tubular furnace based on radiation imaging model solving

ZHANG Xiangyu1, ZHENG Shu2, ZHOU Huaichun2, XU Hongjie1

(National Engineering Research Center of Clean Coal CombustionXian Thermal Power Research Institute CoLtdXianShaanxiChina2Institute of Simulation & Control of Power SystemTsinghua UniversityBeijingChina

Accurate measurement of temperature distribution of pipe surface becomes the key issue for optimization of heating process in tubular furnace. Experiment research of visualization of pipe temperature distribution was conducted on an industrial ethylene cracking furnace based on radiation image processing. Inside the furnace pipe, thermal cracking reaction of naphtha occurs, whose reaction heat is provided by eight burners located on the bottom, and the fuel is natural gas. Coking of furnace pipe was ignored in the model. In order to deal with the radiative heat transfer equation with complex boundary condition in the tubular furnace, a DRESOR method based on Monte Carlo principle was used, and the radiations from tubes, flame and furnace wall were decoupled. Sixteen CCD cameras were mounted on the furnace wall in order to capture the flame image which was transferred into boundary radiation intensity distribution by blackbody calibration. A revised Tikhonov regularization method was used to solve the morbid radiation image equation, and the distributions of pipe temperature and heat flux were measured online and their variation trends with flow direction of naphtha were also discussed. Validated by two different methods, temperature reconstruction error was less than 2% and maximum deviation was within 20K, and the major error occurred on the maximum and minimum temperature areas. This study would be useful for adjustment of combustion and heating process in tubular furnace, and would improve the uniformity of pipe surface temperature in order to extend their working life.

tubular furnace; radiation; pipe; modeling; heat flux

2014-09-28.

ZHANG Xiangyu, zhangxiangyu214@163.com

10.11949/j.issn.0438-1157.20141478

TK 224

A

0438—1157（2015）03—0965—07

國家自然科學(xué)基金項(xiàng)目（51406095，51025622）

2014-09-28收到初稿，2014-11-18收到修改稿。

聯(lián)系人及第一作者：張向宇（1984—），男，博士。

supported by the National Natural Science Foundation of China (51406095, 51025622).