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

彭芳麟，鞠筱林，劉 暢

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

?

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

彭芳麟，鞠筱林，劉 暢

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

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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

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

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

當(dāng)α<(β+d)2/4d時(shí)，圖靈斑圖開始生長(zhǎng).

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

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

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

1 一維情形

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

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

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

2 二維情形

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

圖2 實(shí)驗(yàn)中獲取的斑圖照片

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

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

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

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

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

3) 初始的擾動(dòng)是一條直線.

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

4) 初始擾動(dòng)是隨機(jī)的濃度分布.

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

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

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

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

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

3 討論

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

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

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

圖7 文獻(xiàn)[2]中表示活化子與阻滯子擴(kuò)散速度的示意圖

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

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

斑圖動(dòng)力學(xué)仍然是目前非線性物理中熱門的研究方向，有許多模型在進(jìn)行理論、實(shí)驗(yàn)和數(shù)值模擬的研究.本文的數(shù)值模擬提供了一些有趣的初步結(jié)果，限于所選取的符合分岔?xiàng)l件的參數(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é)與斑圖動(dòng)力學(xué)導(dǎo)論[M].北京：北京大學(xué)出版社，2010:156，138-164.

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

[4] 劉迎東.圖靈斑圖動(dòng)力學(xué)的數(shù)學(xué)機(jī)制[J].北方交通大學(xué)學(xué)報(bào)，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ìn)展[J].淮陰工學(xué)院學(xué)報(bào)， 2006，15(1):54-57.

[8] 鄭佳喻，周華喜，孫金芳，等. B-Z反應(yīng)中斑圖的數(shù)值模擬及實(shí)驗(yàn)研究[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é)系教授，主要從事計(jì)算物理教學(xué)和研究工作.

教學(xué)研究

O 411.3

A

1000- 0712(2016)12- 0001- 06