彭芳麟，鞠筱林，劉 暢

(北京師范大學(xué) 物理學(xué)系，北京 100875)

?

斑圖動力學(xué)中Duffet-Boissonade 方程的數(shù)值模擬

彭芳麟，鞠筱林，劉 暢

(北京師范大學(xué) 物理學(xué)系，北京 100875)

對斑圖動力學(xué)中的Duffet-Boissonade反應(yīng)擴散方程組作了一維和二維的數(shù)值模擬.設(shè)置了4種不同的初始條件并按照圖靈分岔的條件來選取方程中的參數(shù)，模擬得到的六邊形斑圖和條紋狀斑圖與實驗照片很相似.

斑圖動力學(xué)；Duffet-Boissonade 反應(yīng)擴散方程組；數(shù)值模擬

斑圖是在空間或時間上具有某種規(guī)律性的非均勻宏觀結(jié)構(gòu).自然界的斑圖可分為兩類：一類是在熱力學(xué)平衡態(tài)條件下的斑圖，如無機化學(xué)中的晶體結(jié)構(gòu)、有機聚合物中自組織形成的斑圖，它的形成機理可以用平衡態(tài)熱力學(xué)和統(tǒng)計物理原理來解釋；另一類是在偏離熱力學(xué)平衡態(tài)的條件下產(chǎn)生的斑圖，如天上的條狀云、水面上的波浪、動物體皮膚表面的花紋等，它們不能用平衡態(tài)熱力學(xué)原理來解釋，但可以從動力學(xué)角度來進行研究.斑圖動力學(xué)就是研究這類斑圖形成的原因及其規(guī)律的學(xué)科.

1952年英國數(shù)學(xué)家圖靈在論文《形態(tài)形成的化學(xué)基礎(chǔ)》[1]中，首次提出一個數(shù)學(xué)模型來描述生物體表面圖案(如斑馬身上的斑圖)的生成過程.他設(shè)想，在生物體內(nèi)存在一種“形態(tài)子”(比如是某種大分子)，在生物發(fā)育過程中，“形態(tài)子”與其他反應(yīng)物發(fā)生了生物化學(xué)反應(yīng)并在體內(nèi)擴散.當 “形態(tài)子”濃度在空間分布變得不均勻時就會形成斑圖.擴散本是一種常見的物理現(xiàn)象，如不同濃度物質(zhì)之間的擴散，它們并不會形成斑圖.為何圖靈的反應(yīng)擴散模型能產(chǎn)生斑圖呢？因為他提出了一種新的假設(shè)：在擴散過程中存在兩種相反的作用子，“活化子”與“阻滯子”，活化子催化反應(yīng)，阻滯子阻礙反應(yīng)，二者的空間傳播速度不同，互相影響，造成空間各個區(qū)域的活化子與阻滯子濃度分布不同，影響到形態(tài)子的濃度在空間的分布也不同，于是形成了相對穩(wěn)定的斑圖.這種由體系內(nèi)在的反應(yīng)擴散特性所引起的空間均勻態(tài)被破壞的過程稱為圖靈失穩(wěn)(圖靈分岔)，所形成的圖紋稱為圖靈斑圖.

圖靈斑圖的研究直到上世紀60 年代末才開始引起重視，當時諾貝爾化學(xué)獎獲得者普里高金(L.Prigogine) 為首的比利時布魯塞爾熱力學(xué)小組從熱力學(xué)角度研究圖靈斑圖[6].他們證明，在遠離熱力學(xué)平衡態(tài)的條件下，體系可以通過自組織行為形成斑圖，并稱之為耗散結(jié)構(gòu)，這是自然界不同系統(tǒng)中斑圖形成的共性.同時，他們還提出了一個簡單的、不違反任何化學(xué)反應(yīng)動力學(xué)常識的反應(yīng)模型——布魯塞爾子(brusselator).此后科學(xué)家們對這個模型從理論、實驗以及計算機模擬方面進行了大量的研究，大大推動了耗散理論以及斑圖動力學(xué)的發(fā)展.

利用耗散結(jié)構(gòu)理論研究化學(xué)振蕩反應(yīng)的另一個著名的例子是B-Z反應(yīng).B-Z反應(yīng)是指1951年蘇聯(lián)生物學(xué)家Belousov 發(fā)現(xiàn)的一種在溶液中顏色會出現(xiàn)周期變化的化學(xué)反應(yīng)，由于這個現(xiàn)象是違背熱力學(xué)第二定律的，所以不被當時的科學(xué)界所接受.到了1961年研究生Zhabotinsky對該反應(yīng)體系作了進一步的研究與改進，使溶液的顏色在紅色和藍色之間來回顯著地變化，由于當時非線性理論已經(jīng)有了突破，該反應(yīng)體系也就逐步被人們所接受.該實驗的發(fā)現(xiàn)是化學(xué)體系中的非單調(diào)行為概念的一個實驗突破，它的誕生豐富了非線性動力學(xué)理論的內(nèi)容，人們?yōu)榱思o念該反應(yīng)的發(fā)現(xiàn)者就把它稱為B-Z反應(yīng).

在對布魯塞爾子和B-Z反應(yīng)等一系列模型的研究中，人們發(fā)現(xiàn)了大量的非線性現(xiàn)象.例如在時間序列上出現(xiàn)的多峰的周期振蕩，遲豫振蕩到陣發(fā)振蕩的過渡，直到后來確認化學(xué)混沌現(xiàn)象和分岔現(xiàn)象的存在.在時空上也發(fā)現(xiàn)了多姿多彩的斑圖，如靶形波，單螺旋波，雙螺旋波，破碎螺旋波等[7].在《大學(xué)物理》雜志上也曾介紹過有關(guān)B-Z反應(yīng)的實驗[8].這些研究的進展使斑圖動力學(xué)成了非線性動力學(xué)研究的重要方向.

為了更直觀地了解斑圖動力學(xué)，我們下面討論一下《非線性科學(xué)與斑圖動力學(xué)導(dǎo)論》[2]一書中介紹的一個簡單的反應(yīng)擴散模型即Duffet-Boissonade方程組，這是V.Duffet和J.J.Boissonade 在1992年發(fā)表在Chamical Physics上的論文《常規(guī)的與非常規(guī)的圖靈斑圖》[5]中首次建立的一個無量綱化的方程：

(1)

在方程(1)中u、v分別代表活化子與阻滯子的濃度，其中α、β、γ>0代表阻滯子與活化子的擴散系數(shù)比，d是阻滯子的擴散系數(shù).

將這個方程組與擴散系數(shù)為a的物質(zhì)擴散方程：

(2)

作個對比可知Duffet-Boissonade方程組中活化子的擴散系數(shù)是1，因此按照圖靈斑圖生成的條件，阻滯子的擴散系數(shù)d必須大于1.此外，它與物質(zhì)擴散方程還有一個顯著不同，就是在方程組中出現(xiàn)了u、v的交叉項，其中方程的三次項系數(shù)為負數(shù)，保證系統(tǒng)變量在失穩(wěn)后不會因微擾而無限增長．二階項的存在除去了方程在(u，v)→(-u，-v)變換下的不變性，這種對稱性在化學(xué)系統(tǒng)中—般不存在；當γ=0時，可以回到這種對稱的情形.這個模型的優(yōu)點是可以對系統(tǒng)的線性行為與非線性行為分別處理．系統(tǒng)的均勻定態(tài)解的穩(wěn)定條件由α、β、γ決定.在失穩(wěn)后系統(tǒng)的行為由非線性項控制.

根據(jù)文獻[2]中的介紹，方程有3個均勻定態(tài)解，當γ很小時，如果有1<β<α，系統(tǒng)在定態(tài)解u0=0，v0=0附近可能出現(xiàn)圖靈分岔，要出現(xiàn)圖靈分岔還要滿足：

(d-β)2-4d(α-β)=0，d>β

當α<(β+d)2/4d時，圖靈斑圖開始生長.

文獻[2]的分析對定性理解Duffet-Boissonade方程組很有幫助，但是由于沒有給出這個方程組的數(shù)值模擬，很難想象方程組所呈現(xiàn)的復(fù)雜圖像.本文利用文獻[2]給出的圖靈斑圖的生成條件，用MATLAB計算了方程組在一維與二維情形的數(shù)值解，并將結(jié)果可視化，顯示了這種反應(yīng)擴散方程描述的斑圖出現(xiàn)的過程與圖案，這些斑圖與實驗出現(xiàn)的某些圖像有一定的相似性.我們的模擬計算表明，文獻[2]中對控制參數(shù)的分析基本正確，必須根據(jù)圖靈分岔的條件去選擇參數(shù)才能生成斑圖，另外我們發(fā)現(xiàn)，初始條件決定了斑圖的具體形狀，不同的初始條件會生成不同的斑圖，斑圖穩(wěn)定以后能維持很長時間，一般不會轉(zhuǎn)化成別的斑圖類型.這些結(jié)果證實Duffet-Boissonade模型確實能夠描述反應(yīng)擴散中的基本現(xiàn)象.

斑圖動力學(xué)是非線性物理的重要分支，適當?shù)亓私獍邎D動力學(xué)可以幫助學(xué)生開闊視野，提高興趣.本文的數(shù)值模擬所需要的知識僅僅是非線性偏微分方程組，對于學(xué)過計算物理學(xué)的學(xué)生而言，并不難掌握.這種學(xué)習(xí)既能幫助形象地理解與學(xué)習(xí)斑圖動力學(xué)，也是計算物理在教學(xué)科研中的一個應(yīng)用.

本文的計算都假定在區(qū)域的邊界上存在無通量邊界條件[4]，即函數(shù)在邊界上法向?qū)?shù)為零，表示物質(zhì)沒有流出與流入.初始條件則是在區(qū)域存在一個小的微擾，微擾的傳播形成圖靈斑圖.在將偏微分方程離散化時，對時間變量的一階偏微分取向前差分公式，對空間變量的二階差分取中心差分公式[3].計算中的參數(shù)按照前面介紹的圖靈分岔的條件取值，為了對比，二維與一維的計算取相同的參數(shù)，一律取α=1.12，β=1.1，γ=0.01，d=1.5.

1 一維情形

計算區(qū)域為0≤x≤100，x的步長為1，時間t步長也取1，為了表現(xiàn)斑圖生成的過程，步數(shù)t取了從小到大的6個值，分別為t=1， 15， 100，200，1000，15000.計算表明，步數(shù)超過15000以后，斑圖形狀基本穩(wěn)定.初始條件為區(qū)域中心點存在一個小微擾.為了對比我們將u(實線)和v(點線)的圖形畫在同一張圖中.

可以看到t=15時，原來為正值的微擾點兩側(cè)各出現(xiàn)了一個較小的負值的微擾；隨著步數(shù)t增加，振蕩以波動方式傳開，振幅逐漸增大，振幅量級由最初的10-6增大到10-1；當步數(shù)t足夠大時(例如到15 000以后)，振幅趨于平穩(wěn)，形成穩(wěn)定的條紋.圖中u、v的波型基本重合，只是v波的振幅都要小于u波.這表明在相同位置上的活化子濃度小于阻滯子濃度，這正好符合圖靈給出的生成圖靈斑圖的假設(shè).

圖1u，v的濃度分布曲線.

2 二維情形

文獻[2]中登載有幾張實驗圖片如圖2所示，其中標志性的特征是存在六邊形的點狀斑圖與直線或曲線型的條紋斑圖.上面計算的一維情形雖然簡單，但是不便與實驗對比.為此，下面再作二維的數(shù)值模擬.

圖2 實驗中獲取的斑圖照片

下面計算中二維平面區(qū)域取為80×80.取時間步數(shù)t=10，200，1 000，15 000.計算考慮了4種不同的初始條件：1) 即中心點的小擾動，2) 在整個平面區(qū)域中均勻分布的點擾動，3) 直線形的擾動，4) 平行線狀的擾動.由于u、v的圖形很相似，所以只畫出了u的圖形.

1) 初始條件為平面區(qū)域的中心點有一小擾動.

開始出現(xiàn)近似于同心圓的斑圖，隨后圖形輪廓越來越像六邊形.這個結(jié)果對應(yīng)于實驗上所觀測到的六邊形圖案.

2) 初始條件是在區(qū)域中均勻分布的點狀的小擾動.

在斑圖形成的過程中，初始擾動從各個擾動點中心對稱地向外擴散，初期整個區(qū)域仍然有點狀的分布，經(jīng)過足夠長的時間后有些點聯(lián)成了條紋.注意它們始終不會成為均勻的分布.這些點狀分布的圖形4(a)、(b)與實驗圖片2(a)、(c)相似.

3) 初始的擾動是一條直線.

經(jīng)過擴散形成的斑圖是一系列平行線狀條紋，這種圖案形成后便保持基本穩(wěn)定.它與實驗圖片2(b)有些相似.

4) 初始擾動是隨機的濃度分布.

這個結(jié)果究竟如何很難想象.計算模擬表明，開始是呈現(xiàn)不規(guī)則的點狀的分布，經(jīng)過足夠長的時間，也能形成條紋，但是其形狀無法預(yù)測.下圖顯示的條紋是彎曲的，與實驗圖片2(d)有些相似.在數(shù)值模擬中也出現(xiàn)過比較直的條紋.

圖3 圖(a)、(b)、(c)、(d)是初始條件為區(qū)域中心點有一點狀擾動在步數(shù)t=10，200，1000，15000時所生成的u的圖形

圖4 圖(a)、(b)、(c)、(d)是初始條件為區(qū)域中心有均勻分布的點狀擾動在時間步數(shù)t=10，200，1000，15000時所生成的u的圖像

圖5 圖(a)、(b)、(c)、(d)是初始條件為區(qū)域中部有一直線形擾動時在時間步數(shù)t=10，200，1000，15000時所生成的u的圖像

圖6 圖(a)、(b)、(c)、(d)是初始條件為區(qū)域內(nèi)隨機分布濃度在步數(shù)t=10，200，1000，15000時所生成的u的圖像

3 討論

通過數(shù)值模擬，我們發(fā)現(xiàn)有以下幾點值得注意.

1) Duffet-Boissonade方程組雖然簡單也能比較好地描述反應(yīng)擴散現(xiàn)象的基本特性，按照穩(wěn)定性分析的條件選取參數(shù)就能模擬出豐富的反應(yīng)擴散現(xiàn)象，如同我們前面所展示的那樣.但是，我們在數(shù)值計算中也發(fā)現(xiàn)文獻[2]給出條件(d-β)2-4d(α-β)=0 并不嚴格.根據(jù)這個條件和我們選取的α、β值，似乎d的取值只能是1.439 33，而我們的數(shù)值模擬表明，d的取值從1.3到1.8都是可行的.

2) 文獻[2]在解釋活化子與阻滯子對斑圖形成的作用機制也與我們數(shù)值模擬的結(jié)果不一致.該文獻使用了如下圖形描述活化子與阻滯子的傳播過程.

圖7 文獻[2]中表示活化子與阻滯子擴散速度的示意圖

文獻[2]認為由于阻滯子的擴散速度遠大于活化子的，阻滯子將會更快地經(jīng)擴散離開微擾形成上述圖形，注意這里有間斷的波形.但在我們的數(shù)值模擬中，無論u還是v的波形，在傳播中都是連續(xù)的，不會發(fā)生間斷.更重要的是二者的波形基本相似，只是振幅不同.為此我們設(shè)想了斑圖形成機制的另一種解釋：在傳播過程中，二者的波形其實是相同的(見圖1)，但由于振幅不同造成二者在各點的濃度不同.活化子是催化反應(yīng)，而阻滯子是減緩反應(yīng).二者的濃度差決定了在各處的反應(yīng)生成物的濃度也不同，濃度不同就會產(chǎn)生斑圖.

3) 初始條件會決定斑圖的形狀.不同初始條件會生成不同斑圖，不改變控制參數(shù)則它們不會互相轉(zhuǎn)化.這個結(jié)果給我們印象很深.

斑圖動力學(xué)仍然是目前非線性物理中熱門的研究方向，有許多模型在進行理論、實驗和數(shù)值模擬的研究.本文的數(shù)值模擬提供了一些有趣的初步結(jié)果，限于所選取的符合分岔條件的參數(shù)值非常有限，所以還有許多工作值得深入去做，這可以作為今后研究的內(nèi)容.

[1] Turing A M.The Chemical basis of morphogenesis[J]. Phil， Trans. R. Soc. London Ser B， 327， 37(1952).

[2] 歐陽頎.非線性科學(xué)與斑圖動力學(xué)導(dǎo)論[M].北京：北京大學(xué)出版社，2010:156，138-164.

[3] 彭芳麟.計算物理基礎(chǔ)[M].北京：高等教育出版社，2010:330-344.

[4] 劉迎東.圖靈斑圖動力學(xué)的數(shù)學(xué)機制[J].北方交通大學(xué)學(xué)報，2004，28(3):1-3.

[5] Duffet V， Boissonade J J. Conventional and unconventional Turing patterns[J]. Chem Phys，1992，96:664.

[6] Nicolis G ，Prigongine 1. Self- Organizalion in Nonequilibriwn Chemical Syslems[M] . New York: Wiley,1977.

[7] 宗春燕，王玉梅，高慶宇.Belousov-Zhabotinsky反應(yīng)研究進展[J].淮陰工學(xué)院學(xué)報， 2006，15(1):54-57.

[8] 鄭佳喻，周華喜，孫金芳，等. B-Z反應(yīng)中斑圖的數(shù)值模擬及實驗研究[J].大學(xué)物理，2008, 27 (10)：50-54.

Numerical simulation of the Duffet-Boissonade equations in pattern dynamics

PENG Fang-lin，JU Xiao-lin, LIU Chang

(Department of Physics, Beijing Normal University, Beijing 100875, China)

The one- and two-dimensional numerical simulations of Duffet-Boissonade reaction-diffusion equations in pattern dynamics are performed. Setting up four kinds of different initial conditions and selecting the parameters in the equation according to Turing bifurcation condition, the hexagon pattern and fringe pattern obtained from the simulation are very similar to the experimental results.

pattern dynamics；Duffet-Boissonade reaction-diffusion equations ; numerical simulation

2016-03-28；

2016-05-26

彭芳麟(1947—)，男，江西泰和人，北京師范大學(xué)物理學(xué)系教授，主要從事計算物理教學(xué)和研究工作.

教學(xué)研究

O 411.3

A

1000- 0712(2016)12- 0001- 06