左路,夏蘭君

(湖北大學(xué)化學(xué)化工學(xué)院,湖北 武漢 430062)

基于estk-Berggren公式的動(dòng)態(tài)固化動(dòng)力學(xué)模型的非線性回歸分析

左路,夏蘭君

(湖北大學(xué)化學(xué)化工學(xué)院,湖北 武漢 430062)

動(dòng)態(tài)自催化固化反應(yīng)動(dòng)力學(xué)研究中的機(jī)理函數(shù)基于SB公式,動(dòng)力學(xué)模型需要擬合4個(gè)參數(shù):指前因子A、活化能Ea、反應(yīng)階數(shù)m、n.本文中采用非線性最小二乘回歸方法,擬合結(jié)果表明4個(gè)參數(shù)的非線性回歸結(jié)果不穩(wěn)定,尤其表現(xiàn)在參數(shù)A與m的估計(jì)結(jié)果上.筆者對(duì)造成上述結(jié)果的原因進(jìn)行了分析,指出由于模型參數(shù)效應(yīng)曲率相當(dāng)大,增加了模型的非線性強(qiáng)度,使得最小二乘回歸分析結(jié)果不可靠.根據(jù)實(shí)驗(yàn)數(shù)據(jù)特征,本文中在不確定機(jī)理函數(shù)形式前提下,運(yùn)用Model-free方法,先獲得活化能平均意義下的估計(jì)值,再進(jìn)行參數(shù)A、m、n的非線性最小二乘回歸分析.結(jié)果表明減少回歸函數(shù)參數(shù)個(gè)數(shù),可極大程度降低模型的非線性強(qiáng)度,從而獲得有效的動(dòng)力學(xué)參數(shù)估計(jì).

化學(xué)動(dòng)力學(xué);非線性回歸分析;非線性強(qiáng)度;曲率

0 引言

化學(xué)動(dòng)力學(xué)從動(dòng)態(tài)角度觀察化學(xué)反應(yīng),研究反應(yīng)對(duì)時(shí)間與溫度等因素的依賴性,建立解釋這種依賴的經(jīng)驗(yàn)?zāi)Ｐ?通過(guò)實(shí)驗(yàn)數(shù)據(jù)獲得模型參數(shù)來(lái)描述動(dòng)力學(xué)特性.這種方法快速便捷,利用適當(dāng)?shù)膭?dòng)力學(xué)模型從現(xiàn)象角度量化反應(yīng)特征,化學(xué)反應(yīng)機(jī)理研究則退居其次.實(shí)際應(yīng)用中,通過(guò)差式掃描量熱分析法(differential scanning calorimeter,以下簡(jiǎn)稱DSC),獲得動(dòng)態(tài)(非等溫)條件,或者等溫條件下自催化固化反應(yīng)實(shí)驗(yàn)數(shù)據(jù),對(duì)反應(yīng)動(dòng)力學(xué)模型進(jìn)行回歸分析,從而得到描述反應(yīng)動(dòng)力學(xué)特征的動(dòng)力學(xué)參數(shù).實(shí)驗(yàn)中掃描不同升溫速率下(動(dòng)態(tài)反應(yīng)時(shí))或者不同溫度時(shí)(等溫反應(yīng)時(shí))反應(yīng)熱流值,獲得熱流-時(shí)間曲線,熱流曲線與基線間面積為整個(gè)反應(yīng)熱量值ΔHtotal.將時(shí)刻t掃描的熱流值與整個(gè)反應(yīng)熱流值求比值得到該時(shí)刻的反應(yīng)固化程度α(=ΔHt/ΔHtotal),α是時(shí)間t的函數(shù)α(t),取值于0≤α(t)≤1,在反應(yīng)初始階段α(t0)=0.

根據(jù)實(shí)驗(yàn)溫度條件,動(dòng)力學(xué)研究分為等溫與動(dòng)態(tài)方式,等溫法相對(duì)于動(dòng)態(tài)法具有可靠性強(qiáng)的優(yōu)點(diǎn),但是實(shí)驗(yàn)耗時(shí)長(zhǎng),且存在零時(shí)間效應(yīng),從上世紀(jì)60年代開(kāi)始在動(dòng)態(tài)過(guò)程中開(kāi)展快速動(dòng)力學(xué)研究.動(dòng)態(tài)實(shí)驗(yàn)中,溫度按照線性、雙曲型、似等溫型、任意型等控制程序以一定速率上升,對(duì)整個(gè)反應(yīng)溫度范圍進(jìn)行熱流掃描,可以獲得大量不同溫度下的動(dòng)力學(xué)信息.動(dòng)態(tài)法獲得反應(yīng)動(dòng)力學(xué)參數(shù)雖經(jīng)過(guò)長(zhǎng)期發(fā)展和研究,但是大量文獻(xiàn)給出的各種環(huán)境、各種反應(yīng)機(jī)理獲得的動(dòng)力學(xué)參數(shù)有著顯著不同[1,3],存在參數(shù)不穩(wěn)定、不可靠現(xiàn)象[4].造成這種現(xiàn)象的原因即有來(lái)自實(shí)驗(yàn)的誤差因素、模型假定差異,也有數(shù)據(jù)處理方法等因素.本文中將對(duì)動(dòng)態(tài)下自催化固化實(shí)驗(yàn)DSC數(shù)據(jù)進(jìn)行非線性回歸分析,從動(dòng)力學(xué)模型特征上探討造成上述差異的原因.

1 模型描述與數(shù)據(jù)預(yù)處理

動(dòng)態(tài)動(dòng)力學(xué)分析過(guò)程建立于基本速率方程dα/dt=k(T)f(α)[5]之上,dα/dt表示固化速率,k(T)為依賴絕對(duì)溫度T的Arrhenius速率常數(shù),f(α)是依賴于固化度α的反應(yīng)機(jī)理函數(shù)(微分形式).不同反應(yīng)機(jī)理下,f(α)具有不同的形式(表1).1971年estk與Berggren(以下簡(jiǎn)稱SB)將反應(yīng)機(jī)理函數(shù)描述為三參數(shù)形式f(α)=αa(1-α)b(-ln(1-α))p[6],SB公式源自于S形人口增長(zhǎng)模型的Logistic函數(shù)[7],適用于廣泛的動(dòng)力學(xué)分析.1980年Gorbachev將公式簡(jiǎn)化為自催化反應(yīng)的兩參數(shù)形式f(α)=αm(1-α)n[8],αm描述反應(yīng)結(jié)束或者產(chǎn)物形成規(guī)律,另一部分(1-α)n描述反應(yīng)中的自催化效應(yīng). 這兩部分作用相互抵消又相互依賴,指數(shù)m、n取值非負(fù)值,在m≤1時(shí)模型才有實(shí)際意義,指數(shù)n(>1)值與反應(yīng)復(fù)雜程度成正比.

根據(jù)Arrhenius方程[9],反應(yīng)速率常數(shù)為具有負(fù)指數(shù)的指數(shù)函數(shù)k(T)=Ae-Ea/RT,其中參數(shù)A為指前因子(s-1),Ea表示活化能(J/(mol-1·K-1)),R表示氣體常數(shù)(=8.314J·K-1·mol-1). 將Arrhenius方程引入速率方程,得到動(dòng)態(tài)動(dòng)力學(xué)模型為

dα/dt=Ae-Ea/RTαm(1-α)n

(1)

表1 反應(yīng)機(jī)理函數(shù)

在升溫速率分別為5、10、12.5、15、20 ℃/min時(shí),通過(guò)DSC法得到實(shí)驗(yàn)數(shù)據(jù)(T,α,dα/dt)(見(jiàn)圖1-圖2),圖1所示固化度-溫度關(guān)系圖呈現(xiàn)出典型的S形態(tài).非線性最小二乘法進(jìn)行回歸分析時(shí),擬合算法需要計(jì)算模型函數(shù)對(duì)估計(jì)參數(shù)的偏導(dǎo)數(shù),其中偏導(dǎo)數(shù)?f/?m,?f/?n表達(dá)式中將分別出現(xiàn)對(duì)數(shù)函數(shù)lnα與ln(1-α),因此當(dāng)α取值為0或者1時(shí),將導(dǎo)致該偏導(dǎo)數(shù)無(wú)意義,擬合過(guò)程無(wú)法進(jìn)行.因此,剔除實(shí)驗(yàn)數(shù)據(jù)中α為0或1的數(shù)據(jù)點(diǎn),再利用其余數(shù)據(jù)點(diǎn)進(jìn)行擬合. 本文中回歸分析計(jì)算環(huán)境為美國(guó)WolframResearch公司的計(jì)算軟件Mathematica,其強(qiáng)大的計(jì)算功能通過(guò)符號(hào)計(jì)算語(yǔ)言來(lái)實(shí)現(xiàn),Mathematica從1988年第一版發(fā)布便標(biāo)志著現(xiàn)代科學(xué)計(jì)算的開(kāi)端.Mathematica在非線性回歸的大樣本數(shù)據(jù)處理中具有穩(wěn)定、快速、精度高的優(yōu)勢(shì),并且能夠勝任回歸模型非線性強(qiáng)度分析的復(fù)雜計(jì)算.

圖1 不同升溫速率下α-T圖

圖2 不同升溫速率下dα/dt-T圖

2 模型的非線性回歸分析

非線性回歸的最小二乘法算法利用一階泰勒展開(kāi)式線性近似回歸函數(shù),并將此線性近似作為回歸函數(shù)進(jìn)行最小二乘回歸分析.模型(1)的泰勒展開(kāi)式略去增量Δi的二次及其以上的高次項(xiàng),得到如下近似形式

(2)

圖3 模型(1)擬合結(jié)果圖

2.1 4個(gè)參數(shù)的非線性擬合Mathematica中通過(guò)NonlinearModelFit函數(shù),并采用最小二乘Levenberg-Marquardt法進(jìn)行非線性擬合,Levenberg-Marquardt法是改良后的Gauss-Newton法,引入迭代阻尼因子修正了Gauss-Newton法可能存在的不收斂性.Levenberg-Marquardt法的迭代過(guò)程需要參數(shù)的初始值信息,并且對(duì)初始值的依賴性非常強(qiáng),最大迭代次數(shù)的限制也會(huì)影響到最終的結(jié)果.參考文獻(xiàn) [11]中A,Ea,m,n的估計(jì)值,本文中選擇分別用2.0×105(1/s)、40 000J/mol、0.2、1作為參數(shù)A,Ea,m,n的初始值,并設(shè)置迭代最大次數(shù)為500,擬合結(jié)果如表3和表4所示,曲線擬合結(jié)果比較見(jiàn)圖3.除了反應(yīng)后期階段有偏差,所有升溫速率下的回歸結(jié)果與實(shí)驗(yàn)數(shù)據(jù)很好的擬合.這是因?yàn)榉磻?yīng)到了一定階段,反應(yīng)過(guò)程還受到了擴(kuò)散作用的控制,模型(1)中沒(méi)有描述擴(kuò)散因子.

表3和表4中擬合的Adj-R2值均>0.99,與圖3的擬合效果吻合.但是觀察參數(shù)擬合效果發(fā)現(xiàn)A與m的估計(jì)值顯示出很大的波動(dòng)性,并且標(biāo)準(zhǔn)誤差與估計(jì)值的比例分別為60%,95%,95%置信區(qū)間跨越0值,且區(qū)間長(zhǎng)度達(dá)到估計(jì)值的2倍以上,估計(jì)結(jié)果失去了價(jià)值.

2.2模型非線性強(qiáng)度的度量最小二乘Levenberg-Marquardt法擬合實(shí)際是一種線性近似方法,估計(jì)效果的優(yōu)劣與非線性模型有著密切的關(guān)系.Beale[12]在1960年提出了度量非線性模型的線性近似在統(tǒng)計(jì)推斷上優(yōu)劣程度的數(shù)量指標(biāo),即非線性強(qiáng)度的度量——曲率.Bates與Watts[13]將Beale的概念分解為固有曲率、參數(shù)效應(yīng)曲率,并建立比較標(biāo)準(zhǔn)——相對(duì)曲率值,通過(guò)固有曲率、參數(shù)效應(yīng)曲率與相對(duì)曲率的大小關(guān)系來(lái)判斷回歸模型的非線性強(qiáng)度.固有曲率大小反映了期望平面的平坦程度,期望平面越平坦,越可以用切平面去近似,近似效果越好.因此期望平面越平坦,固有曲率就越小.參數(shù)效應(yīng)曲率則反映模型參數(shù)對(duì)于非線性程度的貢獻(xiàn)值,若固有曲率、參數(shù)效應(yīng)曲率均小于相對(duì)曲率值,那么對(duì)于非線性模型的線性近似具有可行性,并對(duì)其進(jìn)行最小二乘法的回歸分析結(jié)果可以接受;固有曲率大于相對(duì)曲率值時(shí),說(shuō)明模型的固有非線性非常強(qiáng),不適合采用線性近似方法進(jìn)行擬合,在實(shí)用上,模型的固有非線性往往可以忽略不計(jì)[14];參數(shù)效應(yīng)曲率大于相對(duì)曲率值時(shí),說(shuō)明模型的非線性特征主要由參數(shù)造成,需要對(duì)參數(shù)進(jìn)行變換,修正模型再進(jìn)行回歸分析.

模型(1)的曲率如表5所示,每個(gè)升溫速率下模型的相對(duì)曲率一致,固有曲率、最大參數(shù)效應(yīng)曲率隨升溫速率增加呈上升趨勢(shì).由于實(shí)驗(yàn)中每秒掃描一次以獲得熱流量數(shù)據(jù),當(dāng)升溫速率低時(shí),反應(yīng)持續(xù)時(shí)間長(zhǎng),獲得的數(shù)據(jù)信息豐富,因此參數(shù)效應(yīng)曲率隨著升溫速率增加而上升.每個(gè)升溫速率下,固有曲率遠(yuǎn)小于相對(duì)曲率,表明模型自身的固有非線性程度可以忽略不計(jì),但是由于參數(shù)效應(yīng)曲率極大程度超過(guò)了相對(duì)曲率,說(shuō)明該模型的非線性特征主要受到參數(shù)的影響,因此希望通過(guò)參數(shù)變換或者減少模型(1)的參數(shù)個(gè)數(shù)降低參數(shù)效應(yīng)曲率.

2.3三個(gè)參數(shù)的非線性擬合觀察圖2,我們會(huì)發(fā)現(xiàn)dα/dt曲線隨著溫度T的升高,有單一峰值出現(xiàn).

表2 不同升溫速率下(da/dT)p,Tp值

圖4 Friedman法求活化能線性擬合圖

dα/dt=Ae-48 960.15/RTαm(1-α)n

(3)

在不同升溫速率下重新擬合得到參數(shù)(表4中模型(3))與模型的非線性曲率值(表5中模型(3)),Adj-R2幾乎沒(méi)有變化,但是參數(shù)A,m,n的值隨著升溫速率的升高呈現(xiàn)出下降的趨勢(shì),并且標(biāo)準(zhǔn)誤差比模型(1)的估計(jì)結(jié)果縮小10-1,置信區(qū)間更短且沒(méi)有跨越0值. 模型(3)的相對(duì)曲率略有減小,固有曲率縮小10-1,參數(shù)效應(yīng)曲率縮小10-4～10-5,并且兩者均小于相對(duì)曲率,減少參數(shù)個(gè)數(shù)極大程度降低了模型的非線性程度,此時(shí)再進(jìn)行的線性近似的最小二乘法是有效且可行的.

表3 模型(1)的非線性擬合結(jié)果

表4 模型(3)的非線性擬合結(jié)果

表5 不同升溫速率下模型(1)、(3)的曲率

3 結(jié)論

動(dòng)態(tài)條件下的快速動(dòng)力學(xué)研究提供了較等溫條件下溫度范圍廣泛的反應(yīng)信息,但是動(dòng)力學(xué)模型也因此引入了新的變量.基于SB公式的自催化反應(yīng)機(jī)理函數(shù)由兩個(gè)參數(shù)決定,模型的復(fù)雜程度因此增加,影響了非線性回歸中參數(shù)A、Ea、m、n的估計(jì)值.回歸模型是溫度和固化度的非線性函數(shù),其非線性強(qiáng)度程度大小將影響回歸分析結(jié)果的有效性.模型參數(shù)效應(yīng)曲率大于相對(duì)曲率時(shí),就要考慮對(duì)參數(shù)進(jìn)行變換,通過(guò)峰值等數(shù)據(jù)特殊點(diǎn)值預(yù)先確定某些參數(shù)值,減少回歸模型參數(shù)個(gè)數(shù),改善擬合模型的非線性性態(tài),利用最小二乘法得到的回歸結(jié)果提供有價(jià)值的動(dòng)力學(xué)特征參數(shù)值.

[1] Sergey V, Charles A W. Model-free and model-fitting approaches to kinetic analysis of isothermal and nonisothermal data[J]. Thermochimica Acta,1999(340/341):53-68.

[2] Omar Moussa, Anastasios P Vassilopoulos, Thomas Keller. Effects of low-temperature curing on physical behaviors of cold-curing epoxy adhesives in bridge construction[J]. International Journal of Adhesion & Adhesives,2012(32):15-22.

[3] Nicolas Sbirrazzuoli. Determination of pre-exponential factors and of the mathematical functionsf(α) ofG(α) that describe the reaction mechanism in a model-free way[J]. Thermochimica Acta,2013(564):59-69.

[4] 葛慶仁.非等溫氣固反應(yīng)動(dòng)力學(xué)曲線非線性擬合求解[J].中國(guó)核科技報(bào)告,1987(S1).

[5] Brown M E, Dollimore D, Galwey A K. Reactions in the solid state[J]. Comprehensive Chemical Kinetics,1980:22.

[8] Gorbachev V M. Some aspects ofestk’s generalized kinetic equation in thermal analysis[J]. J Therm Anal,1980,18:193-197.

[9] Arrhenius S. On the reaction rate of the inversion of non-refined sugar upon souring[J]. Z Phys Chem,1889(4):226-248.

[10] Macro Schwaab, José Carlo Pinto. Optimum reference temperature for reparameterization of the Arrhenius equation, Part 1: Problems involving one kinetic constant[J]. Chemical Engineering Science,2007(62):2750-2764.

[11] Sun L, Pang S, Sterling A, et al. Dynamic analysis of curing process of Epoxy Prepeg[J]. J Appl Polym Sci,2002,86(8):1911-23.

[12] Beale E M L. Confidence region in nonlinear estimation[J]. J R Statist Soc B,1960;22:41-88.

[13] Douglas M Bates, Donald G Watts. Nonlinear regression analysis and its applications[M].New York: John Wiley & Sons Inc.,1988:232-267.

[14] 韋博成.近代非線性回歸分析[M].南京:東南大學(xué)出版社,1989:35-93.

[15] Henry L. Friedman. New methods for evaluating kinetic parameters from thermal analysis data[J]. Journal of Polymer Science Part B: Polymer Letters,1969;7:41-46.

(責(zé)任編輯 胡小洋)

Regression analysis of non-isothermal cure reaction kinetic model based onestk-Berggren equation

ZUO Lu,XIA Lanjun

(School of Chemistry & Chemical Engineering, Hubei University, Wuhan 430062, China)

In chemical kinetic analysis of non-isothermal cure reaction, autocatalytic reaction mechanism function was based onestk-Berggren equation, parameters of regression were pre-exponential factor A, activation energyEa, reaction ordermandn. In this article, we fit rate equation with four parameters by nonlinear regression based on the least squares estimation, regression result was unreliable, especially the estimator ofAandm. This article pointed that nonlinearity of model with large parameter-effects curvature was responsible for the terrible regression result. In the next section, according to characters of experiment data, we estimatedEafirstly with model-free methods, rate equation would only have three parameters, parameter-effects curvature of model was cut short, and with this solution we could get effective estimation of kinetic parameters.

chemical kinetic; nonlinear regression analysis; nonlinearity; curvature

2013-08-05

左路(1975-),講師,E-mail:zlumail@hubu.edu.cn

1000-2375(2014)03-0224-06

TQ018;O643.11

A

10.3969/j.issn.1000-2375.2014.03.008