李聰新，任 成，楊星團，姜勝耀，孫艷飛

（清華大學核能與新能源技術研究院先進反應堆工程與安全教育部重點實驗室，北京 100084）

20世紀60年代日本Shindo發(fā)明以聚丙烯腈（PAN）為原料制造碳纖維，經過50多年的發(fā)展，碳纖維材料行業(yè)已發(fā)展成獨立完整的工業(yè)體系［1］，碳纖維氈作為絕熱保溫材料獲得了廣泛應用［2］。列入《國家中長期科學和技術發(fā)展規(guī)劃綱要》16個重大專項之一的高溫氣冷堆核電站項目，由于采用了石墨堆芯結構和石墨基體燃料元件，相關設計和實驗采用了大量性質兼容的碳纖維材料。高溫氣冷堆運行溫度為750℃［3］，極限事故工況下燃料元件的安全溫度限值為1 600℃［4］，相應的熱工、余熱排出等高溫試驗需要高性能的耐高溫保溫材料。碳氈中90%為空隙，保溫隔熱性能優(yōu)越、穩(wěn)定［5］，成為高溫氣冷堆相關實驗優(yōu)先選擇的保溫材料。

國內對碳氈材料力學［6］、熱學、電學［7］等方面的性能、參數(shù)的研究和測定不足，尤其是在其更具優(yōu)勢的高溫領域（1 000℃以上）。碳氈的導熱性能受溫度、密度、石墨化度、碳纖維取向等多種因素影響［8］，并非常量［9］，導熱系數(shù)仍以實驗測量為主。高溫條件下通常采用激光閃射法［9－13］和熱線法［14］。這兩種方法的測量精度受兩個因素影響較大，一是小尺寸試樣選取不具有代表性，二是測試環(huán)境與使用環(huán)境不一致。

現(xiàn)代工業(yè)越來越多的要求實時精確的測量熱物性，因此基于導熱反問題的材料物性參數(shù)識別方法取得了廣泛應用［15］。自60年代Stolz發(fā)表簡單形狀的導熱反問題數(shù)值計算方法［16］以來，導熱反問題的計算方法極大豐富，如正則化方法、LM方法、神經網(wǎng)絡算法等［17－19］。本文以實測的200多萬組實驗數(shù)據(jù)進行導熱系數(shù)反問題計算，詳細介紹采用非線性導熱反問題方法確定材料溫度相關導熱系數(shù)在物理模型簡化、非線性導熱正問題求解、導熱系數(shù)反問題識別等過程的具體算法和步驟，為用導熱反問題方法測量物性參數(shù)提供一個基本思路。

1 實驗

1．1 材料性能測試裝置

為了給高溫氣冷堆工程及相關實驗提供可靠材料驗證和技術支持，填補高溫氣冷堆在材料相關配套工業(yè)方面的不足，清華大學核能與新能源技術研究院建立了一套模擬高溫氣冷堆溫度、環(huán)境氛圍的材料性能測試裝置，以進行相關材料1 600℃以下的性能測試，溫度涵蓋高溫氣冷堆從正常運行到極限事故工況下的全部溫度范圍。該裝置設計為石墨電阻爐形式，中心采用高純度石墨加熱元件，上、下和四周為碳氈保溫材料。石墨發(fā)熱體和碳氈保溫材料之間形成模擬高溫氣冷堆高溫、碳還原環(huán)境的環(huán)形材料測試區(qū)，該區(qū)域可模擬真空（＜30Pa）、氮氣、氦氣環(huán)境。高溫區(qū)采用特殊設計的可在碳還原環(huán)境中耐受高溫的鎢錸熱電偶穩(wěn)定測量溫度場，精度±0．5%。測試裝置的內部結構和部分測點位置如圖1所示。

通常，保溫材料導熱系數(shù)越小越好，但考慮到承重和強度因素，在材料性能測試裝置中，上保溫層采用軟質碳氈，底部和側面保溫層采用導熱系數(shù)稍大的硬質碳氈。保溫材料的導熱系數(shù)一般隨溫度的升高而增加［20］，在設計時，選取的導熱系數(shù)為廠家給出的設計參考值（表1）。要在很寬溫度范圍內測量保溫材料的導熱系數(shù)有很大難度，這也是本文采用導熱反問題方法的原因之一。

圖1 材料測試裝置內部結構與測點布置

表1 保溫材料熱物性參考值

為使材料測試區(qū)近似柱坐標一維傳熱，上下保溫層厚度遠大于側保溫層厚度。環(huán)形區(qū)高度為1m，測點T1布置在環(huán)形區(qū)軸向和徑向中心，測點T9、T8距上下保溫層分別為125mm，三者溫差很小。連同上下保溫層中測點，在1 400℃和1 600℃保溫階段，軸向溫度分布如圖2所示，可認為爐膛區(qū)域軸向溫度恒定，用一維徑向傳熱模型近似誤差可忽略。

側保溫層厚度為200mm，內半徑R1為250mm，外半徑R3為450mm，在其中心和外邊界分別布置測點T2、T3。側保溫層內側無測點，考慮到環(huán)形區(qū)溫度分布較為均勻，可用測點T1溫度作近似代替保溫層內壁面T′1處溫度來進行導熱計算。近似替代引入的誤差和修正方法將在后文介紹。

圖2 軸向溫度分布

1．2 實驗過程及結果

一次典型的材料測試實驗包括以下過程：1）對爐體小功率預熱，同時抽真空，防止碳氈和石墨材料被氧化；2）當爐膛內壓力降到100Pa以下，調節(jié)功率升溫；3）當材料測試區(qū)溫度上升到1 400℃后，保溫24h；4）再次緩慢升溫到材料測試區(qū)溫度為1 600℃，保溫22h；5）停止加熱、保持冷卻水流量不變，緩慢降溫；6）降溫到165℃時，停止抽真空，充入氮氣冷卻；7）溫度降到120℃時，停止冷卻水循環(huán)，實驗完成。

整個實驗過程持續(xù)時間約為76h，數(shù)據(jù)采集系統(tǒng)對各測點的溫度采集速率為8次／s，共得到219萬組數(shù)據(jù)，3個用于碳氈導熱系數(shù)計算的測點T1、T2、T3溫度如圖3所示。足夠多的測量數(shù)據(jù)對采用靈活的反問題計算方法是有利的。

圖3 溫度歷史曲線

2 正問題

導熱正問題求解是反問題計算的基礎，正問題求解的精度和速度對反問題計算影響很大，須首先建立合理的正問題的計算模型。

2．1 數(shù)學描述與簡化

由于軸向導熱可忽略，正問題可近似為一維柱坐標下變物性非穩(wěn)態(tài)導熱，將求解區(qū)域確定為R1～R3之間的碳氈保溫層區(qū)域，其數(shù)學描述如下。

控制方程：

邊界條件：

T（R3，τ）＝Y3（τ）

（3）

補充條件：

T（R2，τ）＝Y2（τ）

（4）

目標函數(shù)：

其中：r為空間徑向坐標；τ為時間；τi為測量時間點；T（r，τ）為r處τ時刻溫度計算值；Y1（τ）、Y2（τ）、Y3（τ）為測點T1、T2、T3在τ時刻的溫度測量值。待識別的導熱系數(shù)在滿足式（1）～（4）的前提下使得由式（5）定義的誤差函數(shù)值最小。

為了簡化正問題計算，用測點T1的溫度近似替代內壁面T′1處溫度，將內壁面的輻射邊界條件轉化為第一類邊界條件，式（2）變?yōu)椋?/p>

T（R1，τ）＝Y1（τ）

（6）

假設初始條件不存在，主要考慮一長時間的測量過程，可選取任一時間范圍內的數(shù)據(jù)進行反問題計算，其起始時刻并不一定是溫度為室溫的絕對實驗開始時刻?？捎肦1、R2、R3處溫度測量值的拋物插值作為初始條件，其與此時真實初始條件之間的差異將隨計算時間的推移而消除。

至此，由方程（1）、（3）、（6）構成的正問題已滿足定解條件。

物性參數(shù)ρT為密度，cT為比熱容，λT為導熱系數(shù)，下標T表示其為溫度的函數(shù)。對于工程中常見的大多數(shù)材料，從嚴格意義上講，物性參數(shù)均非定值，但對于固體，可認為密度不隨溫度變化，即方程（1）等價于如下方程：

比熱容一般隨溫度變化，但由于兩個邊界條件均為一類邊界條件，不能同時反求導熱系數(shù)和比熱容［21］，為確定導熱系數(shù)，必須假定比熱容為定值，取為廠家給出的參考值。此時控制方程簡化為：

將式（8）右端偏導項展開，得到偏導數(shù)可直接用差分格式代替的形式：

2．2 數(shù)值計算方法

雖然目前主流的商業(yè)CFD軟件如FLUENT、CFX可精確求解三維固體熱傳導方程，但是導熱反問題計算過程中，需根據(jù)反問題算法不斷修改待識別參數(shù)，進行導熱正問題計算，所以高效的做法是根據(jù)具體問題統(tǒng)一編寫正、反問題計算程序。

如圖4所示，將側保溫層沿徑向等分為20個網(wǎng)格，Δr＝（R3－R1）／20＝10mm。

圖4 離散格式

將式（9）中的偏導項用如下差分形式代替：

則可得到：

其中：Tn為節(jié)點n的當前溫度；T0n為節(jié)點n前一時刻溫度；Tn－1、Tn＋1為該節(jié)點相鄰節(jié)點溫度；Δτ為時間步長；Δr為空間步長；rn為當前節(jié)點半徑。將式（11）整理成計算形式為：

需注意，此離散格式有別于常物性一維柱坐標的離散格式：

寫成Ax＝d的矩陣形式為：

其中：

系數(shù)矩陣為標準的三對角矩陣。由于導熱系數(shù)為非定值，系數(shù)矩陣中元素為溫度的函數(shù)，可用追趕法等三對角常系數(shù)矩陣算法迭代求解。

2．3 輻射邊界處理

在上述計算過程中，左邊界條件用環(huán)形材料測試區(qū)中的測點溫度Y1代替，該測點溫度不同于真實碳氈保溫層內壁面溫度。測試區(qū)為真空氛圍，輻射是熱傳遞的主要方式，兩者差值與溫度和熱流密度有關，可采用如下“預估－校正”方法。

1）預估過程

用測點溫度Y1作為計算的邊界條件，用上述方法計算得到整個實驗段內的溫度場，將節(jié)點2的溫度記為T02。假設壁面與溫度為T預的灰體處于輻射平衡狀態(tài)，則有：

其中：ε為灰體表面發(fā)射率，實測值為0．8左右；σ為斯特潘－波爾茲曼常量。圖5為表面發(fā)射率分別取0．6、0．8、1．0時的預估溫差比較，可見，相對于1 600℃的溫度變化范圍，預估溫度對表面發(fā)射率變化并不十分敏感。由式（16）得到測點與壁面溫度的預估溫度差值為：

2）校正過程

將預估差值作為測點和壁面溫度差值，得到以后計算所用內壁面溫度：

用新壁面溫度作為以后計算的邊界條件。該過程并不需要迭代，只需在反問題計算過程中，對導熱系數(shù)值變化較大的試算過程，重新執(zhí)行一次校正過程。

用“預估－校正”方法是為了克服直接處理輻射邊界條件帶來的計算量增大和發(fā)散。預估差值與測點溫度Y1的關系示于圖5a，由于該差值不是溫度的單值函數(shù)，不能用測點溫度做簡單補償。與直接用測點溫度Y1作為內壁面溫度相比，使用該方法更接近實際傳熱過程，可使由式（5）定義的誤差值減小，尤其是實驗剛開始的升溫過程中，其原因可從圖5b看出，在剛開始的升溫過程中，壁溫和測點溫差較以后過程大很多。

圖5 預估溫度差值

3 導熱系數(shù)反問題計算

3．1 一般反問題計算方法

對導熱系數(shù)進行識別，預先對其滿足的函數(shù)關系式進行假設可簡化反問題計算。由于大多數(shù)保溫材料的導熱系數(shù)可認為是溫度的線性函數(shù)［20］，可假設其滿足如下關系式：

λ＝a＋bT

（19）

對導熱系數(shù)的函數(shù)識別變?yōu)閷?shù)a、b的識別。反問題計算即通過不同的參數(shù)a、b組合求解導熱正問題，使由式（5）定義的誤差函數(shù)最小，這是一個二維最優(yōu)化問題。

J（a，b）＝min（Y－T）T（Y－T）

（20）

由多元函數(shù)極值條件知，上式取最小值的條件為：

其中，系數(shù)矩陣稱為敏感度矩陣，敏感度系數(shù)為未知參數(shù)的函數(shù)，可采用如下方法計算：

其中，ε為一個小量。上述問題可采用任何一種基于梯度的迭代求解方法，如最速下降方法、共軛梯度法［22］、Levenberg－Marquardt方法［23］等。這些方法都經過嚴格的數(shù)學推導，具有一定的普適性，但物理概念并不明確，初值選取不當會使迭代過程發(fā)散。

3．2 動靜導熱過程法

根據(jù)實驗過程的特點和固體熱傳導規(guī)律，給出一種簡便的優(yōu)化方法。本次試驗涉及穩(wěn)態(tài)和非穩(wěn)態(tài)導熱過程，本文將利用這一特點分別識別a和b，簡稱為動靜導熱過程法。

假設導熱系數(shù)為常數(shù)，無論其大小，則穩(wěn)態(tài)過程中測點T1、T2、T3的溫度差值之比（Y3－Y2）／（Y2－Y1）將保持不變，即通常認為在熱流未知時，用穩(wěn)態(tài)方法不能測量導熱系數(shù)。要改變其比值，只有改變其導熱系數(shù)隨溫度的變化率，據(jù)此可確定導熱系數(shù)隨溫度的變化率。變化率確定后，改變導熱系數(shù)零點參考值，即絕對大小，使非穩(wěn)態(tài)過程中溫度計算值與測量值偏差減小。這樣便將一個二維搜索變?yōu)閮蓚€一維搜索。

為應用上述方法，將導熱系數(shù)的表達式（式（19））改寫為以下形式：

記a／b為c，分別計算溫度對a、c的敏感度系數(shù)，結果如所圖6所示，其中對應a每變化0．001，c變化0．000 01，6次等差變化的值?？梢娫诜€(wěn)態(tài)階段，a的變化對溫度場幾乎無影響，所以當確定c后用非穩(wěn)態(tài)過程優(yōu)化a不會改變穩(wěn)態(tài)過程的優(yōu)化效果。

在假設一個a值前提下，可使用任何一種一維最小值搜索算法，如二分法，來確定c。在a取不同值時，使穩(wěn)態(tài)過程誤差函數(shù)最小的c值幾乎不變，為0．000 403左右。測點T2處不同a值假設下搜索到最優(yōu)c值時的溫度歷史曲線如圖7所示。

圖6 溫度對a、c的敏感度系數(shù)

圖7 測點T2處不同a值假設下最優(yōu)c對應的溫度歷史

c值確定后，可使用同樣的一維最小值搜索算法，得到在確定c值下使誤差函數(shù)最小的a值，為0．099。

a和c確定后，為保證其組合為最優(yōu)解，可在極小范圍內調整a和c的值，也可以其為共軛梯度法等反問題算法提供迭代初始值，這種良好初值可使反問題算法迅速收斂。

最終導熱系數(shù)的表達式為：為方便工程應用，實驗得出的不同溫度下硬質碳氈的導熱系數(shù)列于表2。

λ＝0．098 2（1＋0．000 404 96T）

（24）

采用這種反問題方法進行物性參數(shù)識別是直接在碳氈材料的使用環(huán)境中對其整體導熱性能的衡量，在工程實踐中具有重要的實用價值。

表2 1 600℃以下硬質碳氈的導熱系數(shù)

3．3 結果分析

在反問題計算過程中，導熱系數(shù)的優(yōu)化目標為測點T2位置處溫度的測量值與計算值在不同時刻差值的二范數(shù)最小，最終二者繪于圖8。由圖8可見，識別的導熱系數(shù)使測點T2處溫度計算值與測量值在整個實驗時間范圍內基本一致。

圖8 測點T2處溫度計算值與測量值比較

T2位置處溫度計算值的相對誤差如圖9所示。溫度場恒定和平穩(wěn)變化的過程中，相對誤差基本控制在±2．5%以內。在功率瞬間調整和工況改變的實驗開始和結束時，計算溫度值的相對誤差較大，最大為15%。在這段時間內，傳熱過程更為復雜，與計算模型之間存在更大的誤差。例如時間坐標為80h處，裝置內溫度已降到200℃以下，輻射傳熱過程減弱，再充入氮氣，氣體導熱作用不可忽略，再將內壁面作為輻射邊界條件已不合適，必然引起更大的誤差。

圖9 計算溫度相對誤差

4 結論

隨著國內需求量日益增長，碳纖維材料已被列為國家化纖行業(yè)重點扶持的新產品，成為國內新材料行業(yè)研發(fā)的熱點。鑒于目前碳纖維材料在高溫條件下性能參數(shù)數(shù)據(jù)不足，本文提供了清華大學核能與新能源技術研究院高溫氣冷堆材料性能測試裝置對碳氈材料76h的驗證溫度數(shù)據(jù)，并據(jù)此采用非線性導熱反問題方法確定了1 600℃以下硬質碳氈材料的導熱系數(shù)，為高溫氣冷堆相關實驗和其他工程應用提供了高溫條件下導熱性能設計依據(jù)。文中詳細介紹了采用非線性導熱反問題方法確定材料溫度相關導熱系數(shù)的完整過程和具體算法，在正問題計算中提出了處理輻射邊界條件的“預估－校正”算法，減少了用測量數(shù)據(jù)作為第一類邊界條件代替實際輻射邊界條件引入的誤差；在反問題處理中提出了依據(jù)物理過程的動靜導熱過程法進行導熱反問題計算，物理概念清晰且簡便可行。該方法既可單獨使用，又可為其他反問題算法提供良好的迭代初值。

［1］ 趙稼祥．碳纖維的發(fā)展與應用［J］．纖維復合材料，1996，46（4）：46－50．

ZHAO Jiaxiang．Development and application of carbon fiber［J］．Fiber Composites，1996，46（4）：46－50（in Chinese）．

［2］ 孫志堅，孫瑋，傅加林，等．國內絕熱保溫材料現(xiàn)狀及發(fā)展趨勢［J］．能源工程，2001，26（4）：26－28．SUN Zhijian，SUN Wei，F(xiàn)U Jialin，et al．Current status and development of thermal insulating materials in China［J］．Energy Engineering，2001，26（4）：26－28（in Chinese）．

［3］ ZHANG Z，WU Z，SUN Y，et al．Design aspects of the Chinese modular high－temperature gas－cooled reactor HTR－PM［J］．Nuclear Engineering and Design，2006，236（5）：485－490．

［4］ NIESSEN H，BALL S．Heat transport and afterheat removal for gas cooled reactors under accident conditions，IAEA－TECDOC－1163［R］．Vienna：IAEA，2001．

［5］ 黎小平，張小平，王紅偉．碳纖維的發(fā)展及其應用現(xiàn)狀［J］．高科技纖維與應用，2005，30（5）：28－34．LI Xiaoping，ZHANG Xiaoping，WANG Hongwei．Modification and application of ultra－h(huán)igh molecular weight polyethylene fiber［J］．Hi－Tech Fiber ＆Application，2005，30（5）：28－34（in Chinese）．

［6］ 易法軍，孟松鶴，韓杰才，等．碳氈／碳復合材料超高溫力學性能實驗研究［J］．無機材料學報，2001，16（6）：1 229－1 234．YI Fajun，MENG Songhe，HAN Jiecai，et al．Mechanical behavior of carbon felt－carbon composites under ultra－h(huán)igh temperature［J］．Journal of Inorganic Materials，2001，16（6）：1 229－1 234（in Chinese）．

［7］ 邱廣瑋，劉平，韓金鐸，等．熱處理溫度對碳纖維氈導電性能的影響［J］．機械工程材料，2012，36（11）：66－68．QIU Guangwei，LIU Ping，HAN Jinduo，et al．Effect of heat－treatment temperature on conductivity of carbon felt［J］．Materials for Mechanical Engineering，2012，36（11）：66－68（in Chinese）．

［8］ 趙建國，李克智，李賀軍，等．碳／碳復合材料導熱性能的研究［J］．航空學報，2005，26（4）：501－504．ZHAO Jianguo，LI Kezhi，LI Hejun，et al．Research on the thermal conductivity of C／C composites［J］．Acta Aeronautica et Astronautica Sinica，2005，26（4）：501－504（in Chinese）．

［9］ 陳潔，熊翔，肖鵬．單向C／C復合材料導熱系數(shù)的計算［J］．炭素技術，2008，27（2）：1－4．CHEN Jie，XIONG Xiang，XIAO Peng．Calculation of thermal conductivity coefficient in unidirectional carbon／carbon composite［J］．Carbon Techniques，2008，27（2）：1－4（in Chinese）．

［10］童文輝，沈峰滿，王文忠，等．高爐常用耐火材料導熱系數(shù)的測定［J］．金屬學報，2002，38（9）：983－988．TONG Wenhui，SHEN Fengman，WANG Wenzhong，et al．Measurement of heat－transfer coefficient of refractory for BF［J］．Acta Metallurgica Sinica，2002，38（9）：983－988（in Chinese）．

［11］王東，劉宗奎，張斌，等．激光閃射法測試氮化硅結合碳化硅材料導熱系數(shù)的研究［J］．耐火與石灰，2012，37（6）：3－6．WANG Dong，LIU Zongkui，ZHANG Bin，et al．Study on coefficient of thermal conductivity of silicon carbide material bonded with silicon nitride determined with laser flash method［J］．Refractories and Lime，2012，37（6）：3－6（in Chinese）．

［12］易法軍，張巍，孟松鶴，等．C／C復合材料高溫熱物理性能實驗研究［J］．宇航學報，2002，23（5）：85－88．YI Fajun，ZHANG Wei，MENG Songhe，et al．An experimental study on thermo physical properties of C／C composites at elevated temperature［J］．Journal of Astronautics，2002，23（5）：85－88（in Chinese）．

［13］鄒林華，黃伯云，黃啟忠，等．C／C復合材料的導熱系數(shù)［J］．中國有色金屬學報，1997，7（4）：135－138．ZOU Linhua，HUANG Boyun，HUANG Qizhong，et al．Thermal conductivity for C／C composites［J］．The Chinese Journal of Nonferrous Metals，1997，7（4）：135－138（in Chinese）．

［14］程遠貴，董岱峰，周勇，等．耐火纖維氈的高溫導熱系數(shù)研究［J］．化工設計，2001，11（2）：22－23．CHENG Yuangui，DONG Daifeng，ZHOU Yong，et al．A study on measuring thermal conductivities of Al2O3fibrous refractories under high temperature［J］．Chemical Engineering Design，2001，11（2）：22－23（in Chinese）．

［15］胡國?。畟鳠嶂械亩嘧诹糠囱菅芯浚跠］．大連：大連理工大學，2002．

［16］STOLZ G，Jr．Numerical solutions to an inverse problem of heat conduction for simple shapes［J］．Journal of Heat Transfer，1960，82：20－25．

［17］ZUECO J，ALHAMA F，F(xiàn)ERN A，et al．An inverse problem to estimate temperature dependent heat capacity under convection processes［J］．Heat and Mass Transfer，2003，39（7）：599－604．

［18］HUANG C，CHEN C．A three－dimensional inverse geometry problem in estimating the space and time－dependent shape of an irregular internal cavity［J］．International Journal of Heat and Mass Transfer，2009，52（7）：2 079－2 091．

［19］王一博，楊海天，鄔瑞鋒．基于時域精細積分算法的瞬態(tài)傳熱多宗量反演［J］．應用數(shù)學和力學，2005，26（5）：512－518．WANG Yibo，YANG Haitian，WU Ruifeng．Precise integral algorithm based solution for transient inverse heat conduction problems with multi－variables［J］．Applied Mathematics and Mechanics，2005，26（5）：512－518（in Chinese）．

［20］姜堯忠．工業(yè)電爐［M］．北京：清華大學出版社，1993．

［21］楊晨，Ulrich Gross．基于熱傳導逆問題方法預測材料熱物性參數(shù)［J］．化工學報，2005，56（12）：2 415－2 420．YANG Chen，Ulrich Gross．Estimation of thermal properties based on inverse heat conduction method［J］．Journal of Chinese Industry and Engineering（China），2005，56（12）：2 415－2 420（in Chinese）．

［22］張有為，李輝，姜培學，等．采用共軛梯度法的管內壁溫度導熱反問題求解［J］．工程熱物理學報，2009，30（7）：1 188－1 190．ZHANG Youwei，LI Hui，JIANG Peixue，et al．Inverse heat conduction problem of deducing inner wall temperature by using conjugate gradient method［J］．Journal of Engineering Thermophysics，2009，30（7）：1 188－1 190（in Chinese）．

［23］LUGON J，Jr，Neto SILVA A J，SANTANA C C．A hybrid approach with artificial neural networks，Levenberg－Marquardt and simulated annealing methods for the solution of gas－liquid adsorption inverse problems［J］．Inverse Problems in Science and Engineering，2009，17（1）：85－96．