戴金東， 艾佳莉， 孫 巍

（北京化工大學(xué)化學(xué)工程學(xué)院，北京 100029）

Belousov-Zhabotinsky（BZ）反應(yīng)是一個(gè)非常經(jīng)典的非線性化學(xué)動力學(xué)系統(tǒng)。由于反應(yīng)與擴(kuò)散的耦合，BZ 反應(yīng)有著豐富的化學(xué)自組織現(xiàn)象，例如時(shí)間維度的化學(xué)振蕩和時(shí)空維度的化學(xué)斑圖等。對于反應(yīng)擴(kuò)散系統(tǒng)中化學(xué)斑圖的研究，艾佳莉等[1]、李才偉等[2]采用偏微分（PDE）模型與元胞自動機(jī)模型對BZ 反應(yīng)進(jìn)行了模擬，得到了BZ 反應(yīng)形成的螺旋波等化學(xué)波圖像。

在豐富的化學(xué)斑圖中，有一種與動態(tài)的化學(xué)波有所不同的特殊圖像?圖靈斑圖。圖靈斑圖是一種逐步產(chǎn)生的定態(tài)條紋，前人已經(jīng)證明了圖靈斑圖在化學(xué)系統(tǒng)中的存在性[3]，并討論了Brusselator 這一理想化學(xué)反應(yīng)系統(tǒng)產(chǎn)生圖靈斑圖的數(shù)學(xué)機(jī)理[4]，但是對于實(shí)際存在并且在非線性化學(xué)動力學(xué)中最為經(jīng)典的BZ 反應(yīng)，鮮有人討論圖靈斑圖產(chǎn)生的條件。

理論上，在一個(gè)體系中只要同時(shí)存在反應(yīng)與擴(kuò)散兩種機(jī)制，就可以通過反應(yīng)速率方程和擴(kuò)散系數(shù)寫出該體系對應(yīng)的數(shù)學(xué)模型，而圖靈斑圖就是對應(yīng)的反應(yīng)擴(kuò)散方程的平衡解經(jīng)歷了“由擴(kuò)散引起的不穩(wěn)定性”之后，產(chǎn)生的穩(wěn)定的、非一致的空間結(jié)構(gòu)。圖靈斑圖的現(xiàn)象存在于化學(xué)與生物體系中，比如鋁在酸性電解液中作為陽極被氧化時(shí)形成的多孔氧化鋁，就是一種圖靈斑圖[5]，其過程實(shí)質(zhì)上是一個(gè)反應(yīng)擴(kuò)散過程，其形成機(jī)理可以用斑圖動力學(xué)來解釋；再比如浮游生物在水中的分布，也可以用反應(yīng)擴(kuò)散模型來表示，進(jìn)而以斑圖的形式呈現(xiàn)出浮游生物在水中的分布規(guī)律[6]。本文對BZ 反應(yīng)進(jìn)行了圖靈不穩(wěn)定分析，用圖像展示了BZ 反應(yīng)中圖靈斑圖的可能形態(tài)。BZ 反應(yīng)是一種簡單的反應(yīng)擴(kuò)散體系，對BZ 反應(yīng)中圖靈不穩(wěn)定的研究可以為其他復(fù)雜的反應(yīng)擴(kuò)散體系提供參考。

1 BZ 反應(yīng)及其圖靈不穩(wěn)定分析

1.1 BZ 反應(yīng)模型的建立

BZ 反應(yīng)的機(jī)理非常復(fù)雜，經(jīng)過不斷簡化最終得到了五步三變量的俄勒岡機(jī)理模型[7]。俄勒岡機(jī)理模型認(rèn)為BZ 反應(yīng)豐富的時(shí)空有序現(xiàn)象是由HBrO2、Br?、Ce4+這三種關(guān)鍵物質(zhì)的不斷相互轉(zhuǎn)化實(shí)現(xiàn)的，通過寫出這三種物質(zhì)的反應(yīng)速率方程，以及經(jīng)過量綱為一處理可以得到現(xiàn)在普遍使用的Tyson模型：

其中 ?a 與 ?b 為擴(kuò)散項(xiàng)：

方程組（1）包含2 個(gè)變量與5 個(gè)參數(shù)，且變量與參數(shù)均在化簡過程中經(jīng)過了量綱為一化，其中a 代表HBrO2的濃度，b 代表Ce4+的濃度（但a、b 并不等同于濃度）。ε、f、q 是反應(yīng)動力學(xué)參數(shù)，其中ε 是與初始濃度以及溫度相關(guān)的參數(shù)，其取值一般在0～1 之間；q 是只與溫度相關(guān)的參數(shù)，其取值一般也在0～1 之間，并根據(jù)經(jīng)驗(yàn)其值一般遠(yuǎn)小于1；f 是一個(gè)可調(diào)參數(shù)，其取值一般在0～5 之間。d1、d2為兩種物質(zhì)的擴(kuò)散系數(shù)。t 為時(shí)間，x、y 分別表示橫、縱坐標(biāo)。

這里只考慮物質(zhì)在二維平面上的擴(kuò)散，因此可以規(guī)定 x ∈[0,L] 和 y ∈[0,L] ，其中L 為x 和y 所能取得的最大長度。L 的長度在模擬時(shí)取400 個(gè)空間步長（h），可以模擬出完整的圖靈斑圖或化學(xué)波圖案，并通過多次模擬試驗(yàn)得到結(jié)論：當(dāng) h ≤1 時(shí)可以得到較為連續(xù)、不失真的模擬圖像。空間步長也就是單位步長，使用有限差分法進(jìn)行模擬時(shí)會對每單位步長采用拉格朗日公式處理以近似替代對應(yīng)點(diǎn)的導(dǎo)數(shù)值。為方便計(jì)算，后文在驗(yàn)證計(jì)算結(jié)果時(shí)取h=1，L=400進(jìn)行模擬。

可以看出，當(dāng)Tyson 模型不考慮擴(kuò)散項(xiàng)時(shí)，系統(tǒng)是一個(gè)常微分（ODE）系統(tǒng)，物質(zhì)濃度僅隨時(shí)間的變化而變化；當(dāng)Tyson 模型考慮擴(kuò)散項(xiàng)時(shí)，系統(tǒng)是一個(gè)偏微分（PDE）系統(tǒng)，物質(zhì)濃度同時(shí)隨時(shí)間和位置的變化而變化。

圖靈不穩(wěn)定是指平衡解在不考慮擴(kuò)散項(xiàng)的ODE 系統(tǒng)中穩(wěn)定，在考慮擴(kuò)散項(xiàng)的PDE 系統(tǒng)中變得不穩(wěn)定的一種情況[8]，因此圖靈不穩(wěn)定分析要分別在兩個(gè)系統(tǒng)下討論平衡解的穩(wěn)定性。

1.2 BZ 反應(yīng)ODE 模型及其平衡解的穩(wěn)定性分析

首先寫出不考慮擴(kuò)散項(xiàng)的Tyson 模型：

系統(tǒng)在平衡點(diǎn)附近受到的擾動同樣可以寫為式(3)的形式，只不過需要再將平衡解代入到系數(shù)中，以表示是在平衡點(diǎn)附近受到的擾動：

其中系數(shù)矩陣為：

至此就將對微分方程的分析轉(zhuǎn)化為對系數(shù)矩陣的分析，系數(shù)矩陣的特征方程為

式（9）可以表示為：

其中T 為矩陣的跡，即矩陣對角線之和，D 為矩陣行列式的值，分別如式（11）和式（12）所示。

當(dāng)系數(shù)矩陣的特征值均為負(fù)數(shù)時(shí)，對應(yīng)的平衡解在系統(tǒng)中是穩(wěn)定的[6]，而通過判斷矩陣的跡與行列式值的正負(fù)，就可以得到特征方程對應(yīng)的二次函數(shù)曲線圖，將3 組平衡解分別代入到T 與D 之中，觀察對應(yīng)的特征方程可以得到結(jié)論：第1 組平衡解（a1*，b1*)的特征值為一正一負(fù)，不管參數(shù)取何值，該平衡解均不能在ODE 系統(tǒng)中保持穩(wěn)定；第2 組平衡解（a2*，b2*)對應(yīng)的物質(zhì)濃度出現(xiàn)了負(fù)值，不符合常理；而第3 組平衡解（a3*，b3*)的穩(wěn)定性與參數(shù)的取值相關(guān)，接下來詳細(xì)討論第3 組平衡解。

根據(jù)特征方程，使特征值均為負(fù)數(shù)的條件可以等價(jià)于使矩陣的跡小于零、矩陣的行列式的值大于零，即只需要讓各參數(shù)滿足 T <0 、 D>0 這兩個(gè)條件，就可以說明此平衡解是穩(wěn)定的。

為計(jì)算使第3 組平衡解在ODE 系統(tǒng)中保持穩(wěn)定的參數(shù)范圍，首先計(jì)算使 D>0 的參數(shù)范圍，將第3 組解代入到D 的表達(dá)式中可以得到：

顯然對于該式，當(dāng)q>0 且q ?1時(shí)，f 在0～5 之間取任何數(shù)都可以保持 D>0 。接下來計(jì)算使 T <0 的參數(shù)范圍，將第三組解代入到T 的表達(dá)式中可以得到：

化簡該不等式便可以得到BZ 反應(yīng)圖靈不穩(wěn)定分析的第1 個(gè)參數(shù)限制條件：

對于第3 組平衡解，只要參數(shù)的取值滿足式(15)，則該平衡解在ODE 系統(tǒng)中是保持穩(wěn)定的。

1.3 BZ 反應(yīng)PDE 模型及其等價(jià)轉(zhuǎn)化

對于Tyson 模型，在考慮擴(kuò)散項(xiàng)之后，系統(tǒng)會變?yōu)楦訌?fù)雜的偏微分方程組。為了應(yīng)用Routh-Hurwitz判據(jù)對平衡解的穩(wěn)定性進(jìn)行分析，必須先對偏微分方程組進(jìn)行轉(zhuǎn)化。

Routh-Hurwitz 判據(jù)[9]給出了一種判斷數(shù)學(xué)模型解的穩(wěn)定性的方法，對于兩變量的常系數(shù)常微分方程組，其解的形式固定為 C1eλ1t+C2eλ2t，即系統(tǒng)受到的擾動可以寫為這種解的形式，其中C1與C2為常數(shù)，由初始條件確定，t 為時(shí)間， λ1與 λ2為方程組對應(yīng)系數(shù)矩陣的特征值。當(dāng)兩個(gè)特征值均為負(fù)實(shí)數(shù)時(shí)，系統(tǒng)受到擾動之后，只要經(jīng)過足夠長的時(shí)間（t→∞）總會使擾動對狀態(tài)的影響變?yōu)?，即系統(tǒng)總會回到原來的狀態(tài)，這時(shí)該狀態(tài)是穩(wěn)定的；若兩個(gè)特征值不全為負(fù)數(shù)，則該平衡狀態(tài)是不穩(wěn)定的。

對于考慮擴(kuò)散項(xiàng)的BZ 反應(yīng)Tyson 模型，需將其向兩變量常系數(shù)常微分方程組轉(zhuǎn)化。首先將系統(tǒng)在平衡點(diǎn)附近受到的擾動寫為微分方程組的形式，然后將平衡解代入到該方程組的系數(shù)中，使原來的方程組變?yōu)橐粋€(gè)常系數(shù)偏微分方程組，并寫出考慮擴(kuò)散項(xiàng)的Tyson 模型：

與1.2 節(jié)的處理方法相似，此處認(rèn)為a 與b 為系統(tǒng)在平衡點(diǎn)附近受到的擾動，其中 J11、 J12、 J21、J22均與ODE 系統(tǒng)中的系數(shù)相等。

為了使用Routh-Hurwitz 判據(jù)，需要繼續(xù)將方程組向兩變量常系數(shù)常微分方程組的形式進(jìn)行等價(jià)轉(zhuǎn)化，為此將方程組的解寫為傅里葉級數(shù)的形式：

其中k 與r 均為矩陣， k=kikj， r=[x y]T，其中 ki、kj稱為波數(shù)，且 ki=iπ/L ， kj= jπ/L ，這里只考慮物質(zhì)在二維平面上的擴(kuò)散，因此可以規(guī)定 x ∈[0,L] 和y ∈[0,L]。將方程組的解代入式(16)中可以得到：

其中 αij和 βij表示對應(yīng)項(xiàng)的變量。

再對 ?sin(kr) 做一個(gè)處理：

最終可以得到一個(gè)與式(16)等價(jià)的方程組：

該方法稱為分離變量法，通過傅里葉級數(shù)展開的方法將偏微分方程組轉(zhuǎn)化為多個(gè)常微分方程組求和的形式，再將該方程組對應(yīng)的新的系數(shù)矩陣提取出來：

至此就完成了從偏微分方程到常微分方程、從常微分方程到系數(shù)矩陣的轉(zhuǎn)化，之后再進(jìn)行與1.2 節(jié)相似的處理就可以得到使平衡解在考慮擴(kuò)散項(xiàng)的PDE系統(tǒng)中變得不穩(wěn)定的參數(shù)范圍，并且BZ 反應(yīng)模型用傅里葉級數(shù)變換得到的系數(shù)矩陣與文獻(xiàn)[10]中Lotka-Volterra 混合離散反應(yīng)擴(kuò)散模型使用不同的處理方法得到的系數(shù)矩陣具有相同的形式，由此可以證明該系數(shù)矩陣的正確性。

1.4 平衡解在PDE 系統(tǒng)中的穩(wěn)定性分析

新的系數(shù)矩陣的跡 T1與新的系數(shù)矩陣的行列式的值 D1為：

與第1.2 節(jié)的處理方式相同，當(dāng)系數(shù)矩陣的兩個(gè)特征值均具有負(fù)實(shí)部時(shí)，系統(tǒng)是穩(wěn)定的，而目前計(jì)算的是由于擴(kuò)散項(xiàng)的加入使平衡解變得不穩(wěn)定的參數(shù)范圍，因此 T1<0 與 D1>0 這兩個(gè)條件中至少有一個(gè)條件不滿足。首先對 T1進(jìn)行處理：

保持平衡解在常微分系統(tǒng)中穩(wěn)定時(shí)，滿足T <0，因此 T1<0 依然成立。即若產(chǎn)生圖靈不穩(wěn)定，則參數(shù)的取值一定要滿足 D1≤0 ，這也是保持平衡解在常微分系統(tǒng)穩(wěn)定的基礎(chǔ)上偏微分系統(tǒng)唯一需要滿足的條件。對 D1進(jìn)行處理：

其中 J11J22?J12J21=D ，D 為1.2 節(jié)中常微分方程對應(yīng)系數(shù)矩陣的行列式的值，保持平衡解在常微分系統(tǒng)中穩(wěn)定時(shí)滿足 D>0 ，即 H(k2) 對應(yīng)的二次函數(shù)曲線的截距為正值，此時(shí)要使 H(k2) 對應(yīng)的二次函數(shù)曲線在 (0,+∞) 上存在低于x 軸的部分，需滿足兩個(gè)條件：（1）二次函數(shù)曲線的對稱軸在y 軸右邊；（2）二次函數(shù)存在兩個(gè)實(shí)根，即判別式大于或等于零，這兩個(gè)條件寫為數(shù)學(xué)表達(dá)式即為

該不等式組表明了使BZ 反應(yīng)產(chǎn)生圖靈不穩(wěn)定的5 個(gè)參數(shù)之間的限制條件，只要一組參數(shù)符合以上3 個(gè)條件，則系統(tǒng)將會發(fā)生圖靈不穩(wěn)定。在驗(yàn)證該條件時(shí)，需要取出一組滿足不等式組的參數(shù)值，可以采用固定參數(shù)的方法，即首先固定q，再在計(jì)算出的f 的取值范圍中取一個(gè)f 值，然后根據(jù)式(15)化簡出ε 的取值范圍，比如q=0.000 8，f=0.9，ε=0.77 就是一組滿足第1 個(gè)不等式的可選參數(shù)。再根據(jù)已選參數(shù)對式(26)進(jìn)行化簡，最后選擇一個(gè)d1的值，則對應(yīng)一個(gè)d2的邊界條件。為使計(jì)算結(jié)果更加直觀，可以通過作三維圖像的方法，固定q、f、ε 這3 個(gè)參數(shù)，得到由擴(kuò)散系數(shù)d1與d2組成的圖靈不穩(wěn)定參數(shù)范圍圖像，如圖（1）所示。

圖1d1 與 d2 的取值范圍Fi g.1Value range of d1 andd2

圖1 中縱坐標(biāo)表達(dá)式為z=(J11d2+J22d1)2?4d1d2(J11J22?J12J21)，即式(27)中第3 個(gè)不等式，在此三維界面上如果取一組d1與d2使函數(shù)值z >0 ，且滿足 d1< 0.989 6d2，則系統(tǒng)產(chǎn)生圖靈不穩(wěn)定，此時(shí)系統(tǒng)對應(yīng)產(chǎn)生的圖像應(yīng)為圖靈斑圖，即自發(fā)產(chǎn)生空間定態(tài)條紋。需要說明的是，這一套計(jì)算方法與計(jì)算結(jié)果有待驗(yàn)證。

2 對BZ 反應(yīng)圖靈不穩(wěn)定的模擬驗(yàn)證

采用數(shù)值模擬的方法對BZ 反應(yīng)進(jìn)行模擬，首先對設(shè)置條件進(jìn)行說明。初始條件中除了區(qū)域中心點(diǎn)（中心點(diǎn) a 與 b 為任意值）外，區(qū)域中其他處a 與b 的值均為零，邊界條件采用循環(huán)的邊界條件，最后對數(shù)學(xué)模型采用區(qū)域離散化與有限差分法，區(qū)域離散化將連續(xù)的區(qū)域近似替代為離散的區(qū)域，有限差分法則通過拉格朗日公式將兩點(diǎn)間函數(shù)的導(dǎo)數(shù)近似為兩點(diǎn)間的斜率，從而將微分方程組近似替代為代數(shù)方程組，以便于求出方程組的近似解。模擬時(shí)取空間步長h=1，在400×400 的網(wǎng)格區(qū)域進(jìn)行模擬。

基于以上思路在Matlab 中編寫出BZ 反應(yīng)的模擬程序，進(jìn)而對不同參數(shù)下的BZ 反應(yīng)現(xiàn)象進(jìn)行模擬。這里只對b(r, t)的變化進(jìn)行模擬，即可呈現(xiàn)出不同參數(shù)下BZ 反應(yīng)所產(chǎn)生的圖像。因?yàn)閎 是量綱為一變量，所以在模擬時(shí)沒有對應(yīng)的單位。模擬時(shí)用顏色表示物質(zhì)濃度，濃度越高越接近黃色，濃度越低越接近藍(lán)色。

根據(jù)1.4 節(jié)的計(jì)算結(jié)果，在固定q=0.000 8，f=0.9，ε=0.77， d1=2 的情況下， d2∈(4.457,+∞) 時(shí)系統(tǒng)產(chǎn)生圖靈不穩(wěn)定。圖2 所示是在固定q、f、ε、d1這4 個(gè)參數(shù)下取 d2=5 時(shí)的模擬結(jié)果?？梢钥闯?，當(dāng)參數(shù)取值在圖靈不穩(wěn)定范圍內(nèi)時(shí)，反應(yīng)系統(tǒng)從中心開始緩慢出現(xiàn)一圈又一圈的環(huán)形定態(tài)條紋，模擬結(jié)果驗(yàn)證了計(jì)算結(jié)果的正確性。

為完善模擬結(jié)果，固定q=0.000 8，f=0.9，ε=0.77，d1=2，取 d2=7 進(jìn)行模擬，結(jié)果如圖3 所示。發(fā)現(xiàn)同樣在圖靈不穩(wěn)定參數(shù)范圍內(nèi)，當(dāng)固定的4 個(gè)參數(shù)取值相等時(shí)，兩個(gè)擴(kuò)散系數(shù)d1、d2相差越大，則系統(tǒng)出現(xiàn)斑圖的速度越快。

固定q=0.000 8，f=0.9，ε=0.77， d1=2 ，取d2=2進(jìn)行模擬，結(jié)果如圖4 所示?？梢钥闯雠c前兩組不同，系統(tǒng)并沒有產(chǎn)生空間定態(tài)斑圖，而是出現(xiàn)了物質(zhì)濃度由內(nèi)向外的動態(tài)擴(kuò)散現(xiàn)象。因此認(rèn)為只有當(dāng)參數(shù)取值在圖靈不穩(wěn)定參數(shù)范圍內(nèi)時(shí)，系統(tǒng)才會產(chǎn)生定態(tài)斑圖，反之則不會產(chǎn)生。

d1=2 d2=5Fig.2Simulation result when and圖2、 時(shí)的模擬結(jié)果d1=2d2=5

圖3d1=2 、 d2=7 時(shí)的模擬結(jié)果Fig.3Simulation result when d1=2 andd2=7

圖4d1=2 、 d2=2 時(shí)的模擬結(jié)果Fig.4Simulation result when d1=2 andd2=2

3 結(jié) 論

（1）以BZ 反應(yīng)擴(kuò)散系統(tǒng)為研究對象，對反應(yīng)的數(shù)學(xué)模型進(jìn)行了圖靈不穩(wěn)定分析，計(jì)算了使BZ 反應(yīng)產(chǎn)生圖靈不穩(wěn)定的參數(shù)限制條件，并作出三維圖像以直觀表現(xiàn)出計(jì)算結(jié)果。

（2）使用Routh-Hurwitz 判據(jù)對平衡解的穩(wěn)定性進(jìn)行判斷，使用傅里葉級數(shù)展開對偏微分方程進(jìn)行處理，最后通過數(shù)值模擬對計(jì)算結(jié)果進(jìn)行了驗(yàn)證，得到了BZ 反應(yīng)產(chǎn)生圖靈不穩(wěn)定時(shí)所呈現(xiàn)出的定態(tài)斑圖圖像。

（3）模擬結(jié)果證明了所使用的計(jì)算方法與得到的參數(shù)范圍是正確的。