擴(kuò)散過(guò)程的元胞自動(dòng)機(jī)模擬

李延升，侯珂珂，張保林

(1.許昌學(xué)院 化學(xué)化工學(xué)院，河南 許昌 461000； 2.鄭州大學(xué) 化工與能源學(xué)院，河南 鄭州 450007)

擴(kuò)散過(guò)程是重要的質(zhì)量傳遞方式之一，元胞自動(dòng)機(jī)則是近年來(lái)新興的仿真模擬方法.但是，將兩者結(jié)合，即用元胞自動(dòng)機(jī)研究擴(kuò)散過(guò)程的文獻(xiàn)卻相對(duì)較少.有關(guān)此類的報(bào)道，多是和反應(yīng)過(guò)程相關(guān)聯(lián)且研究的側(cè)重點(diǎn)放在后者[1]，很少見(jiàn)到用元胞自動(dòng)機(jī)專門研究擴(kuò)散過(guò)程的報(bào)道.另外，在有關(guān)的報(bào)道中，元胞自動(dòng)機(jī)模型的參數(shù)往往過(guò)多且相互之間的關(guān)系復(fù)雜，分析時(shí)需要綜合多門學(xué)科的理論，這限制了它的實(shí)際應(yīng)用.作為一種新興的研究手段，元胞自動(dòng)機(jī)的意圖是以極其簡(jiǎn)單的規(guī)則解釋或模擬復(fù)雜的現(xiàn)象，而有關(guān)的報(bào)道多數(shù)違背了這一意圖，無(wú)法體現(xiàn)元胞自動(dòng)機(jī)的優(yōu)越性.本文即從此方面入手，將元胞自動(dòng)機(jī)方法與數(shù)學(xué)分析方法相對(duì)比，給出了生動(dòng)形象的擴(kuò)散過(guò)程動(dòng)態(tài)畫(huà)面，為進(jìn)一步的研究奠定了基礎(chǔ).

1 元胞自動(dòng)機(jī)簡(jiǎn)介

元胞自動(dòng)機(jī)有時(shí)也被稱為細(xì)胞自動(dòng)機(jī)、點(diǎn)格自動(dòng)機(jī)、分子自動(dòng)機(jī)或單元自動(dòng)機(jī)，它是現(xiàn)代計(jì)算機(jī)之父Neumann及其追隨者提出的想法.20世紀(jì)末21世紀(jì)初，Stephen將這種帶有強(qiáng)烈的純游戲色彩的原始想法從學(xué)術(shù)上加以分類整理，最終使之上升到了科學(xué)方法論[2].

元胞自動(dòng)機(jī)是一時(shí)間和空間都離散的動(dòng)力系統(tǒng)，散布在規(guī)則格網(wǎng)(Lattice Grid)中的每一元胞(Cell)取有限的離散狀態(tài)，遵循同樣的演化規(guī)則，依據(jù)確定的局部規(guī)則同步更新，大量元胞通過(guò)簡(jiǎn)單的相互作用而構(gòu)成動(dòng)態(tài)系統(tǒng)的演化.不同于一般的動(dòng)力學(xué)模型，元胞自動(dòng)機(jī)不是由嚴(yán)格定義的物理方程或函數(shù)確定，而是由一系列模型構(gòu)造的規(guī)則構(gòu)成，凡是滿足這些規(guī)則的模型都可以算作元胞自動(dòng)機(jī)模型.因此，元胞自動(dòng)機(jī)是一類模型的總稱，或者說(shuō)是一個(gè)方法框架[1-2].元胞自動(dòng)機(jī)最基本的組成為元胞、元胞空間、鄰居及規(guī)則這4部分.簡(jiǎn)單來(lái)講，元胞自動(dòng)機(jī)可以視為由一個(gè)元胞空間和定義于該空間的變換函數(shù)所組成[1-3].

2 擴(kuò)散過(guò)程元胞自動(dòng)機(jī)模型的提出

同普通的元胞自動(dòng)機(jī)一樣，本文的元胞自動(dòng)機(jī)也是采用等間隔的點(diǎn)作為元胞.不失一般性，以二維的圓形區(qū)域?yàn)殡x散域，以處于中心的某個(gè)元胞為圓心，取適當(dāng)長(zhǎng)度的半徑，在圓域內(nèi)的元胞即為該元胞自動(dòng)機(jī)的元胞，正好在圓上的交叉點(diǎn)也看作圓域內(nèi)的元胞，見(jiàn)圖1.以縱橫方向上均勻分布的點(diǎn)為元胞，圓域內(nèi)外的元胞區(qū)分非常明確.為了讓畫(huà)面清晰，該圖的元胞較少，在實(shí)際應(yīng)用時(shí)，元胞數(shù)量要比該圖多得多.

圖1 模擬擴(kuò)散過(guò)程的元胞自動(dòng)機(jī)離散域及元胞Fig.1 The discrete region and cells of cellular automata used for simulate the diffusion process

現(xiàn)以元胞自動(dòng)機(jī)分析某一典型的擴(kuò)散過(guò)程.該擴(kuò)散過(guò)程的假設(shè)是：(1)參與擴(kuò)散的物質(zhì)為A和溶劑M.(2)擴(kuò)散過(guò)程在某一圓形區(qū)域及其外界進(jìn)行.(3)初始時(shí)刻，物質(zhì)A以溶液的形式，均勻地分布在該圓形區(qū)域內(nèi).隨后，A即開(kāi)始由內(nèi)到外擴(kuò)散.顯然，在圓形區(qū)域中，越靠近圓心的位置，A的濃度越大.(4)圓形區(qū)域內(nèi)同時(shí)還有物質(zhì)B，該物質(zhì)不參與擴(kuò)散過(guò)程，其濃度均一且不發(fā)生變化.(5)圓域外的區(qū)域視為無(wú)窮大.根據(jù)這一假設(shè)，如果設(shè)該圓形區(qū)域的半徑為a，某一元胞所在的位置處半徑為r，那么在r≥a處，A的濃度為0.

根據(jù)擴(kuò)散過(guò)程的特點(diǎn)，元胞的設(shè)定及演化規(guī)則是：

(1)以元胞0、1、2分別代表物質(zhì)M、A、B，再用k代表演化次數(shù).顯然,k和時(shí)間t對(duì)應(yīng)，也是變量.

(2)用元胞1的隨機(jī)行走代表物質(zhì)A的擴(kuò)散.具體規(guī)定是：第k次演化時(shí)，某元胞1的上、下、前、后、左上、左下、右上、右下的8個(gè)鄰居中，若有n個(gè)為元胞0(n≥1)，則在第k+1次演化時(shí)，該元胞1與這n個(gè)元胞0的任意一個(gè)交換位置.若該交換完畢后，元胞1的位置已經(jīng)在圓域的邊界處，則該元胞1演化為元胞0.

(3)元胞2不發(fā)生演化.

(4)根據(jù)前文的假設(shè)(5)，圓域外界只有元胞0，無(wú)元胞1和元胞2.

3 結(jié)果及討論

3．1 演化結(jié)果

用一個(gè)實(shí)例分析.設(shè)初始時(shí)刻，元胞自動(dòng)機(jī)中，元胞0、1、2所占的量分別為40%、20%、40%，且均勻分布，元胞自動(dòng)機(jī)的離散域是直徑為100的圓，用Matlab 編寫程序得到的演化圖見(jiàn)圖2，該圖可模擬擴(kuò)散過(guò)程且較為生動(dòng)形象.

圖2 模擬擴(kuò)散過(guò)程的元胞自動(dòng)機(jī)演化示意圖Fig.2 Sketch maps simulating the diffusion process with cellular automata

定義：(1)演化次數(shù)為k時(shí)，圓域內(nèi)已經(jīng)擴(kuò)散到外界的元胞1的數(shù)量占初始時(shí)刻元胞1的數(shù)量的比例為累計(jì)擴(kuò)散率.(2)累積演化的元胞1的數(shù)量占初始時(shí)刻數(shù)量的99%時(shí)，對(duì)應(yīng)的演化次數(shù)k為ke，ke即為擴(kuò)散的終點(diǎn).元胞1的累計(jì)擴(kuò)散率曲線ke/t-k見(jiàn)圖3.從圖2與圖3可以看出，隨著k的增大，累計(jì)擴(kuò)散率曲線由陡變緩，說(shuō)明擴(kuò)散速率由大逐漸變小.

圖3 元胞自動(dòng)機(jī)模擬的累積擴(kuò)散量曲線以及用微積分求得的累積擴(kuò)散量率曲線Fig.3 The accumulate diffuse amount curves obtained with cellular automata and the calculous respectively

改變?cè)詣?dòng)機(jī)的參數(shù)，如元胞的數(shù)量、比例等，運(yùn)行程序，得到的累計(jì)擴(kuò)散率曲線基本不變，說(shuō)明該曲線只受演化規(guī)則的影響，與元胞自動(dòng)機(jī)的參數(shù)無(wú)關(guān).

3.2 微積分分析

再用微積分的方法分析這一擴(kuò)散過(guò)程，并與元胞自動(dòng)機(jī)方法對(duì)比.該擴(kuò)散過(guò)程實(shí)際上是球形擴(kuò)散體系在球的大圓截面上的擴(kuò)散過(guò)程.為了更符合實(shí)際情況，以球坐標(biāo)分析.取球心為坐標(biāo)原點(diǎn)，建立球坐標(biāo).設(shè)球體內(nèi)，物質(zhì)A的濃度為cA，則cA是空間半徑變量r和時(shí)間變量t的函數(shù)，角度方向上無(wú)擴(kuò)散分量.物質(zhì)A的傳遞方程為

(1)

初始條件：cA=c0,

(2)

邊界條件：cA(l,t) =0,

(3)

以上三式中，cA與元胞1的數(shù)量對(duì)應(yīng)；r與元胞自動(dòng)機(jī)某處的位置對(duì)應(yīng)；t為時(shí)間變量，與演化次數(shù)k對(duì)應(yīng)；DAB為物質(zhì)A在物質(zhì)B中的擴(kuò)散系數(shù)，與元胞1的隨機(jī)行走的情況對(duì)應(yīng).這三式組成的偏微分方程的解為

(4)

由此，可求得物質(zhì)A的擴(kuò)散速率NA(t)和累計(jì)擴(kuò)散率x1：

(5)

(6)

為了進(jìn)一步闡述的需要，茲列舉一個(gè)實(shí)例，該實(shí)例來(lái)自文獻(xiàn)[6],實(shí)例為球形的暈海寧膠囊的擴(kuò)散控制的釋放系統(tǒng)，其藥物擴(kuò)散體系符合本文的假設(shè).膠囊的直徑為3.26 mm，暈海寧在控釋凝膠層中的擴(kuò)散系數(shù)DAB為3.0×10-7m2/s，由此得

(6a)

用Matlab輔助計(jì)算，取n=1 000，即

(6b)

則截?cái)嗾`差r1 000的范圍為

截?cái)嗾`差已經(jīng)足夠小，得到的累積擴(kuò)散率曲線x1/te-t亦示于圖3中，便于對(duì)比.圖中，te為累積釋放率達(dá)到99%時(shí)的t值，與元胞自動(dòng)機(jī)的ke對(duì)應(yīng).從圖3可以看出，兩條曲線非常接近.實(shí)際上，任意取一DAB的值，并作x1/te-t曲線，結(jié)果都是如此，從而說(shuō)明用元胞自動(dòng)機(jī)的研究結(jié)論與用微積分研究的結(jié)論是一致的.

4 結(jié) 論

以元胞自動(dòng)機(jī)作為方法和手段，通過(guò)編寫Matlab程序，可以形象、逼真地模擬擴(kuò)散過(guò)程.元胞自動(dòng)機(jī)的核心是，各參數(shù)與擴(kuò)散現(xiàn)象緊密相關(guān)，用元胞的演化模擬擴(kuò)散的過(guò)程.元胞自動(dòng)機(jī)模擬的結(jié)果與微積分的結(jié)論一致，并由實(shí)例為佐證.

用元胞自動(dòng)機(jī)模擬擴(kuò)散過(guò)程，雖然具有規(guī)則簡(jiǎn)單、模擬效果較直觀等優(yōu)點(diǎn)，但是所用元胞自動(dòng)機(jī)的演化規(guī)則和重要參數(shù)，基本上是根據(jù)釋放的過(guò)程而作的人為規(guī)定，可以說(shuō)沒(méi)有完全脫離經(jīng)驗(yàn)?zāi)Ｐ偷目蚣芎湍Ｊ?當(dāng)然，目前絕大多數(shù)元胞自動(dòng)機(jī)模擬的情況也都是這樣，這也是元胞自動(dòng)機(jī)尚未被廣泛應(yīng)用的重要原因之一.如何通過(guò)理論分析得到演化規(guī)則和演化參數(shù)，使該模擬過(guò)程上升到較高的理論層次，屬于進(jìn)一步研究的內(nèi)容.

參考文獻(xiàn)：

[1] Bastien Chopard, Meichel Droz.物理系統(tǒng)的元胞自動(dòng)機(jī)模擬[M].祝玉學(xué)，趙學(xué)龍，譯.北京：清華大學(xué)出版社，2003：1-47.

[2] Stephen Wolfram.A new kind of science[M].Wolfram media, 2002:1-50.

[3] 曹偉.元胞自動(dòng)機(jī)與計(jì)算機(jī)模擬[J].丹東紡專學(xué)報(bào)，2005，12(2)：1-4.

[4] Crank J.The mathematics of diffusion[M].Oxford: Oxford university press, 1975:89-103.

[5] 李延升，張保林，張雪梅，等.膜控型緩控釋肥料養(yǎng)分釋放模型的研究[J].江蘇農(nóng)業(yè)學(xué)報(bào)，2009，25(5)：1033-1038.

[6] Welty J R, Wicks C E, Wilson R E.動(dòng)量、熱量和質(zhì)量傳遞原理[M].馬紫峰，吳衛(wèi)生，譯.北京：化學(xué)工業(yè)出版社，2005：338-341.