化學(xué)動力學(xué)方程的數(shù)值解法

1、增刊化學(xué)世界105化學(xué)動力學(xué)方程的數(shù)值解法陸兆仁(華東大學(xué),上海200051)在化學(xué)動力學(xué)的領(lǐng)域中,動力學(xué)方程的求解并非是簡單的事。通常提出一個反應(yīng)機理需要用實驗來驗證,時間的關(guān)系;,至沒有解析解,法來作近似處理,使用是有條件的1;,反應(yīng)機理各種各樣,求微分方程組的解往往不僅是很麻煩的事,而且需很多數(shù)學(xué)上的技巧,以制備膽甾醇異硫脲的甲苯磺酸鹽(CholesterylIsothiuroniumTosylate)的反應(yīng)為例,R.GPearson提出反應(yīng)機理:A+C,無:龍格2庫塔(Runge22卡羅(MonteCarlo)法。龍格-庫塔(Runge-Kutta)法3設(shè)微分方程為dC dt=f(c,

2、t);初始條件是t=0時,C=C0;把所需計算的時間t分為若干個格點,(最簡單可等間距)相鄰兩點時間間距為t;對應(yīng)的濃度C由于第一點(t=0;C=C0)已知,以后各點則由前一點的數(shù)值逐個根據(jù)以下公式推算出來:(第i個格點ti,時為Ci:由前一點ti-1時的Ci-1計算)6;Ci=Ci-1+(K1+2K2+2K3+K4)t式中:K1=f(Ci-1,ti-1);2,ti-1+t 2);K2=f(Ci-1+K1t2,ti-1+t 2);K3=f(Ci-1+K2tK4=f(Ci-1+K3t,ti-1+t);其動力學(xué)方程:設(shè)t=0時A=A0;B=C=D=0;dA dt=-(k1+k3)A;dB dt=d

3、C dt=k1A-k2BC=2k1A-k2B;2dD dt=k2BC+k3A=k2B+k3A對于微分方程組,在計算時,須把濃度C和函數(shù)f都視為向量;同樣可逐點計算各組分濃度3。圖1為上述制備(CholesterylIsothiuroniumTosylate)反應(yīng);-1A0=0.004425mol L;k1=0.0044min;k2=即便是這表面看來簡單的反應(yīng);雖然可得到解析解2:1 2( B=A0)k1+k32.96min-1mol-1L;k3=0.0021min-1。圖中曲J0(2i1(1()11)+iH0(2i)H0+i-(k1+k3)t其中:=e;=iJ1(2i) H1(2i);(k1+

4、k3)2;=k1k2A0J0和J1為虛宗量的貝塞爾(Bessel)函數(shù);H0和H1為虛宗量的漢克(Hankel)函數(shù);i是單位虛數(shù)。由于它不能僅用初等函數(shù)描述;(其中有含復(fù)數(shù)的貝塞爾函數(shù)和漢克函數(shù));要得到某時刻的濃度,仍需求助計算機。為此,提出了多種數(shù)值解法。此地,為開設(shè)課程,培養(yǎng)學(xué)生結(jié)合物修稿日期:2000208220線為解析解,為Runge2Kutta法計算值。由圖1可見,符合的很好。用Runge2Kutta法解微分方程運算快,程序相對有也較短,但由于涉及向量,需一些線性代數(shù)的知識,同學(xué)解微分方程時有點困難,容易出錯。2蒙特-卡羅(MonteCarlo)法蒙特2卡羅法是一種抽樣統(tǒng)計法,常

5、用于計算機模擬,從1963年L.JChaad最早引入動力學(xué)方程求解至今已近半個世紀(jì),其方法隨著計算機的發(fā)展而不斷改進,文獻4介紹的模擬方法由于把所取的隨機數(shù)與速率常數(shù)比較而決定是否反應(yīng),僅適用于一級反應(yīng)。本文作了如下的改進:將所取的在0和1之間的隨機數(shù)與含速率常數(shù)和各組分濃度的函數(shù)比較106化學(xué)世界2000年率;)如其數(shù)值大于零,說明它的數(shù)量增加,如小于零則說明它減少;因此可按下面步驟模擬反應(yīng)過程:把一個單位時間分為M(2000)等分,每一等分作為模擬中的一步,(每一步相當(dāng)1 2000單位時間)設(shè)第i步結(jié)束時A,B,D的濃度分別為a,b,d,它們分別都用N(10000)個球來代表;每一個球分別

6、對應(yīng)于a N;b N和d N;(i+1)濃度時:A,之間的隨機數(shù),統(tǒng)計1+S;(dt=1 2000)則反應(yīng)Sa N,經(jīng)過此步留下未反應(yīng)的A1-S N)a再計算B,令k=(k1A B-k2B)dt1來決定是否反應(yīng),。方法如下:(1)簡單n級反應(yīng):n-dC dt=kC;t=0時,C=a。根據(jù)速率常數(shù)中使用的時間單位(如分或秒等)確定若k>0說明經(jīng)過此步B的濃度增加;取N個0和1之間的隨機數(shù),統(tǒng)計其中小于k的個數(shù)S;B增加的濃度為Sb N;所以經(jīng)過此步后B的濃度為(1+S N)b;若k<0;說明經(jīng)過此步B的濃度減少;令k=-(k1A B-k2B)dt;同樣取N個0和1之間的隨機數(shù),統(tǒng)計其

7、中小于k的個數(shù)S;經(jīng)過此步后B的濃度等于(1-S N)b;如k=0則經(jīng)單位時間,如1分(或1秒)定為1個單位時間;將1個單位時間分為M等分(如2000等分)計算每個間隔點的濃度C;如t=i 2000時的濃度為Ci;設(shè)t=(i-1) 2000時(前一點)的濃度為Ci-1;將Ci-1用N個完全相同的球(即樣本)代表;每個球代表的濃度為Ci-1 N;(在實際計算中取N應(yīng)足夠大;如N=10000);由于每個球完全相同,其反應(yīng)的概率當(dāng)然也相同,就等于整體的反應(yīng)概率。將微分方n程-dC dt=kC兩邊同除C再乘dt:n-1-dC C=kCdt;-dC為濃度C減少數(shù)量;-dC C即減少的百分?jǐn)?shù),也就是反應(yīng)體

8、系整體的反應(yīng)概率;從而也是每個球的反應(yīng)概率;故可模擬如下:第i步計算時Ci過此步B的濃度不變。仍為b。然后再計算D;因D是最終產(chǎn)物,隨時間增加而增加;令k=(k3A D+k2B2 D)dt同樣取N個隨機數(shù),統(tǒng)計其中小于k的個數(shù)S,經(jīng)過此步后D的濃度為(1+S N)d;如此一步步計算直至計算到所需要計算的時間為止。圖1為上述制備CholesterylIsothiuronium取N=10000個在0和1之間的隨機數(shù);統(tǒng)計其中小于kCn-1dt的隨機數(shù)個數(shù)S;(式中dt即相鄰兩點時間的間隔dt=1 2000;kCn-1dt當(dāng)n=1就等于kdt;此時和文獻4一致)留下未反應(yīng)的球為N-S;Ci=(N-S

9、)Ci-1 N=(1-S N)Ci-1;再同法進行下一步模擬計算。直至最后一點。圖2為計算結(jié)果,初濃度為1.5mol L。(2)復(fù)雜反應(yīng)以上述制備CholesterylIsothiuroniumTosy2late的反應(yīng)為例;將動力學(xué)方程組改寫為:A=-(k1+k3)dtdAB=(k1A B-k2B)dtdBD=(k3A D+k2B2 D)dtdDA,dB B,dD D分別代表A,dAB,D的變化百分?jǐn)?shù);(絕對值即該物質(zhì)反應(yīng)的概圖2圖中圓點為模擬值,從上而下依次為0,1,2,3級反應(yīng)。k=0.5(下轉(zhuǎn)第184頁)184化學(xué)世界ua2000年=O的吸收峰面積;AA1600為脲C=O的吸收峰面積;為

10、苯環(huán)C=C的吸收峰面積。計算結(jié)果表明隨著降解時間的延長,(Xb)CO和Ch都在逐漸減小,也就是說硬段之間的氫鍵在降解氨酯軟段和硬段的相分離。在經(jīng)過微生物降解之后,隨著氫鍵的減弱,這種軟硬段的相分離程度也有所減弱。而這則會使微生物對聚氨酯的作用變得更為容易。(下轉(zhuǎn)第87頁)之后減弱了。有關(guān)資料顯示,在硬段含量達到一定數(shù)值后,大部分氫鍵存在于硬段內(nèi)部,并因此產(chǎn)生了聚(上接第106頁)Tosylate的反應(yīng),初始條件t=0時A=0.mol L;B=C=D=0。2(MonteCarlo)前2塔(Runge2Kutta),和解析解(曲線)符合很好圖3為連續(xù)二級反應(yīng):A+Bk1,;運行也更快一些。而用蒙特

11、,(所以稱其為模擬法)但運算時間稍長一些。對于化學(xué)類學(xué)生來講,用蒙特2卡羅法來解方程組更妥當(dāng);至少在結(jié)合專業(yè)學(xué)習(xí)程序編制方面而言不失為較好內(nèi)容。(2)計算的誤差:對于龍格2庫塔法,格點分割得越細,誤差就越小,運行時間也長。對于蒙特2卡羅法代表組分濃度的球N越多,把單位時間分割的越細(M越大),誤差也越小,當(dāng)然運行時間也越長。為了便于與文獻值比較,本文N取10000,M取2000,同樣在龍格2庫塔法中也把一單位時間分為2000個格點,計算結(jié)果所有點和解析解比較,相對誤差都小于2%;已滿足一般動力學(xué)需要;如需更精確,可格點分割的更細一些。(3)用蒙特2卡羅法時,如速率常數(shù)數(shù)值太小;有可能N個隨機數(shù)全部比反應(yīng)的概率小(因k很小,dt=1 M也很小),在這種情況下,只要把時間的單位換成大一點就行了(如把分換成小時)否則就有可能運行很長時間,濃度不改變。參考文獻:1穆爾JW,等著.孫承鍔,等譯.化學(xué)動力學(xué)和歷程2勻+Dk2E初始條件:t=0時:A=B=1.5mol L;C=D=E=0;-1-1-1k1=7.5molLmin;k2=5.2molLmin-1;(解析解可根據(jù)模擬結(jié)果和解析解符合也很好。文獻5方法解得)相化學(xué)反應(yīng)研究(第二版).北京:科學(xué)出版社,1987.361,369.2ChaadLJ,AmJ.Chem.Soc,1963,(85):3588.圖3