分形結(jié)構(gòu)及尺度對粒子擴散行為的影響

曹馨予, 孫洪廣*, 靳 垚, 李志鵬

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室， 南京 211100; 2.河海大學(xué)力學(xué)與材料學(xué)院， 南京 211100)

隨著社會經(jīng)濟的發(fā)展，人類生產(chǎn)和生活過程中產(chǎn)生的污染物通過各種途徑進入土壤。如果污染物濃度超過土壤自凈能力，則土壤質(zhì)量和功能會發(fā)生改變，最后對人類及其他生物的生存和發(fā)展造成危害。調(diào)查結(jié)果顯示，中國部分地區(qū)土壤問題嚴(yán)峻，耕地土壤環(huán)境質(zhì)量堪憂，土壤污染防治迫在眉睫[1]。

1965年Mandelbrot[2]首次提出了分形的概念，并認(rèn)為自然界中存在大量的自相似分形結(jié)構(gòu)。隨后經(jīng)過一系列研究分析表明,土壤是一種具有分形特征的多孔結(jié)構(gòu)介質(zhì)[3-5]。污染物在土壤中的輸運行為受到介質(zhì)結(jié)構(gòu)和尺度的影響[6-7]。對于實際的土壤等多孔介質(zhì)，Katz等[8]也通過電鏡實驗證實，它們的結(jié)構(gòu)往往具有不同尺度的非均勻性,并具有分形特征。

20世紀(jì)初期，Einstein[9]建立了布朗隨機運動理論，提出了溶質(zhì)擴散的位移二階矩表達(dá)式。布朗運動可以視為與時間有關(guān)的隨機分形。在歐式空間，根據(jù)隨機行走理論，粒子布朗運動的空間分布函數(shù)為高斯分布[10]。溶質(zhì)擴散的傳統(tǒng)研究方法主要基于Fick擴散定律。但是由于土壤等復(fù)雜多孔介質(zhì)不滿足均勻介質(zhì)假設(shè)，F(xiàn)ick定律不再適用。因而，從20世紀(jì) 70 年代開始，考慮介質(zhì)非均勻性對粒子隨機運動影響的隨機行走反常擴散模型得以建立和發(fā)展[11-12]。復(fù)雜介質(zhì)中粒子隨機運動的平均擴散位移r2(t) 的一般形式可以表述為〈r2(t)〉=t2H[13-14],其中H為Hurst數(shù)。H=0.5表示正常擴散；H<0.5代表次擴散,等待時間發(fā)散而跳躍步長收斂的隨機行走過程可視為此類運動；H>0.5時的擴散稱為超擴散，其特征為等待時間收斂而跳躍步長發(fā)散，如湍流中粒子的隨機運動[15]；H=1.0表示彈道擴散[16-17]。

隨著分形理論的發(fā)展，中外研究者以分形理論為基礎(chǔ)，構(gòu)建了多種多孔介質(zhì)分形模型。例如，Turcio等[18]構(gòu)造了一種分形模型，得到多孔介質(zhì)中有效滲透率表達(dá)式；施明恒等[19]通過構(gòu)建分形導(dǎo)熱模型,導(dǎo)出了多孔泡沫材料中分形導(dǎo)熱系數(shù)的計算公式；周祥春等[20]建立了多孔介質(zhì)分形模型，研究了分形維數(shù)等因素對粒子擴散行為的影響；Cai等[21]提出一種新型變孔徑三維分形模型，對頁巖地層中流體流動機理進行了分析；Li等[22]通過構(gòu)造分形多孔介質(zhì)模型探究了多孔介質(zhì)內(nèi)部結(jié)構(gòu)與有效應(yīng)力之間的關(guān)系。已有研究表明[23-24]，構(gòu)造人工的分形多孔介質(zhì)模型可以為探索復(fù)雜介質(zhì)溶質(zhì)擴散機理和理解污染物輸運機制提供有效幫助。

現(xiàn)首先構(gòu)造兩種結(jié)構(gòu)和尺度各異的Sierpinski地毯分形模型模擬真實的土壤多孔介質(zhì)結(jié)構(gòu)[25]。進而通過粒子在分形結(jié)構(gòu)中的隨機行走過程模擬探索污染物在土壤中的輸運過程[26-27]。最后通過粒子在兩種分形結(jié)構(gòu)中擴散行為的對比分析，探究分形結(jié)構(gòu)及尺度差異對粒子擴散行為的影響。

1 分形結(jié)構(gòu)模型

1.1 土壤分形結(jié)構(gòu)

土壤是由不同形狀和大小的固體組分及孔隙以一定形式連接而成的多孔介質(zhì)[28]。已有研究證實[29-31]，許多土壤結(jié)構(gòu)具有分形結(jié)構(gòu)特征，表現(xiàn)出自相似性。因此，將采用Sierpinski地毯這一“人工構(gòu)造分形體”模型來表征自然界中實際土壤的介質(zhì)結(jié)構(gòu)。

1.2 兩種分形結(jié)構(gòu)

目前土壤分形結(jié)構(gòu)的相關(guān)研究[32-34]選取的分形結(jié)構(gòu)大多為單塊分形結(jié)構(gòu)，這類分形結(jié)構(gòu)一定程度上能夠反映土壤這類多孔介質(zhì)的異質(zhì)性。然而實際觀測中發(fā)現(xiàn),由于采樣土壤所處的環(huán)境各異，實驗室土樣與自然環(huán)境土樣的尺度差別不容忽視，這些因素會影響預(yù)測結(jié)果的準(zhǔn)確性[35-36]。為準(zhǔn)確模擬污染物在土壤中的輸運行為，建立兩種分形結(jié)構(gòu)：單塊分形結(jié)構(gòu)及四塊拼接分形結(jié)構(gòu)，目的是考察結(jié)構(gòu)及尺度對粒子擴散行為的影響。其中，四塊拼接分形結(jié)構(gòu)是一種內(nèi)含重復(fù)分形結(jié)構(gòu)單元的模型。基于這種拼接分形結(jié)構(gòu)模型，將著重探究粒子的擴散問題。表1列舉了上述兩種分形結(jié)構(gòu)的孔隙分形維數(shù)及空間尺度。其中孔隙分形維數(shù)由盒計數(shù)法計算得到。

表1 兩種分形結(jié)構(gòu)的孔隙分形維數(shù)及空間尺度Table 1 Fractal dimensions of pores and spatial scales of two types of fractal structure

1.3 拼接分形結(jié)構(gòu)的構(gòu)造

四塊拼接分形結(jié)構(gòu)的構(gòu)造方法如下：

步驟1構(gòu)造單塊分形結(jié)構(gòu)。

通過MATLAB軟件建立一級Sierpinski地毯分形結(jié)構(gòu)。應(yīng)用迭代方法，構(gòu)造四級Sierpinski地毯分形結(jié)構(gòu)，如圖1(a)所示。

步驟2構(gòu)造四級拼接分形結(jié)構(gòu)。

選取如圖1(a)所示的四級分形結(jié)構(gòu)。通過軟件進行圖像處理，生成四塊拼接分形結(jié)構(gòu)，如圖1(b)所示。

圖1 兩種分形結(jié)構(gòu)模型Fig.1 Two models of fractal structure

1.4 粒子行走過程描述

將分形結(jié)構(gòu)中心位置作為隨機行走起始點，在起始點處投放5 000個粒子。采用全瞎的螞蟻隨機行走模式模擬擴散過程。假定粒子在孔隙與基質(zhì)中均能夠擴散[37]，在基質(zhì)中的行走步長為1 pixel，在孔隙中的行走步長為5 pixel。只允許粒子在分形結(jié)構(gòu)邊界范圍內(nèi)隨機運動。若粒子行走一定步數(shù)后，位置超出邊界范圍，則該粒子停止運動，停留在邊界處，待下一步隨機跳躍時，從停留處重新隨機跳躍。若下一步的行走位置超出邊界，則重復(fù)上述過程，若未超出邊界，則按規(guī)則正常跳躍。所有粒子完成規(guī)定的行走步數(shù)記為一次有效的行走過程，記錄所有粒子隨機行走過程中每一步的位置坐標(biāo)。具體的數(shù)值模擬算法流程,如圖2所示。

圖2 數(shù)值模擬程序算法流程圖Fig.2 Algorithm flow chart of numerical simulation program

2 粒子擴散過程模擬

2.1 單塊分形結(jié)構(gòu)中粒子擴散過程模擬

在單塊分形結(jié)構(gòu)的中心位置(500，500)處投放5 000個粒子。記錄粒子行走500～10 000步的平面空間分布、密度分布及粒子擴散均方位移與時間的關(guān)系。

粒子行走1 500、2 500、6 000、10 000步后的空間分布，如圖3所示。由圖3可知，隨著行走步數(shù)增加，粒子的擴散區(qū)域逐漸增大。由于受到分形結(jié)構(gòu)中基質(zhì)的阻礙作用，基質(zhì)中的粒子數(shù)相對較少。

從圖3(a)中觀察到，粒子行走1 500步后，未出現(xiàn)粒子擴散至分形結(jié)構(gòu)邊界處，擴散區(qū)域較小，此時粒子的擴散過程只受到結(jié)構(gòu)內(nèi)基質(zhì)的阻礙作用。從圖3(b)中發(fā)現(xiàn)，粒子行走2 500步后，擴散區(qū)域進一步增大，少量粒子擴散至結(jié)構(gòu)邊界。此時粒子的擴散受到基質(zhì)和邊界效應(yīng)的共同影響。由圖3(c)、圖3(d)可知，粒子行走6 000、10 000步后，大量粒子擴散至結(jié)構(gòu)邊界，邊界效應(yīng)對粒子擴散過程的影響逐漸增大。

圖3 單塊分形結(jié)構(gòu)中粒子行走1 500、2 500、6 000、10 000步后的空間分布Fig.3 The spatial distribution of particles after 1 500, 2 500, 6 000, 10 000 steps in single fractal structure

粒子行走1 500、2 500、6 000、10 000步后的密度分布如圖4所示?？芍?，粒子的密度分布是與分形結(jié)構(gòu)有關(guān)的不均勻分布。在起始處粒子密度最大，隨著粒子行走步數(shù)增加，密度分布范圍逐漸增大。由于基質(zhì)對粒子擴散行為的阻礙作用，基質(zhì)中的粒子密度較小。大量粒子的平面密度分布在左右方向上大致對稱。

圖4 單塊分形結(jié)構(gòu)中粒子行走不同步數(shù)后的密度分布Fig.4 The density distribution of particles after different steps in single fractal structure

圖5為粒子行走1 500、2 500、6 000、10 000步后的均方位移與時間的關(guān)系曲線。表2列舉了粒子行走500～10 000步時對應(yīng)的Hurst參數(shù)值。

由圖5和表2分析可知，在分形結(jié)構(gòu)尺度為1 000 pixel×1 000 pixel、孔隙分形維數(shù)為1.90的單塊分形結(jié)構(gòu)中，粒子行走1 500步前，Hurst值小于0.5。此時粒子的擴散只受到基質(zhì)對它的阻礙作用，粒子的擴散行為表現(xiàn)為次擴散；粒子行走1 500～6 000步時，Hurst值增大到0.5附近波動。這是由于當(dāng)粒子行走到分形結(jié)構(gòu)邊界時，邊界效應(yīng)限制了粒子的行走方向,這相當(dāng)于縮小了粒子在孔隙與基質(zhì)中的步長比。粒子的擴散行為出現(xiàn)短暫的正常擴散現(xiàn)象；粒子行走6 000～10 000步時，Hurst值急劇減小。該階段大量粒子行走至邊界處，邊界效應(yīng)對粒子擴散過程的影響增大。由于空間尺度較小，粒子位移無法進一步增加，輸運行為受到阻滯，Hurst值減少。

圖5 單塊分形結(jié)構(gòu)中粒子擴散均方位移與時間的關(guān)系曲線Fig.5 The relationship between mean square displacement and time of particles in single fractal structure

表2 單塊分形結(jié)構(gòu)中的Hurst值

2.2 四塊拼接分形結(jié)構(gòu)中數(shù)值模擬

在四塊拼接分形結(jié)構(gòu)中心位置(1 000，1 000)處投放5 000個粒子。記錄粒子行走500～10 000步的空間分布、密度分布及粒子擴散均方位移與時間的關(guān)系。

圖6表示粒子在四塊拼接分形結(jié)構(gòu)中行走2 500、7 000、10 000步后的空間分布。從圖6可觀察到，粒子的空間分布范圍隨著行走步數(shù)的增加逐漸增大。粒子在孔隙中的擴散速度較在基質(zhì)中更快。由于此分形結(jié)構(gòu)呈對稱分布，粒子的空間分布也呈對稱狀態(tài)。

從圖6(a)中觀察到，粒子行走2 500步后仍未行走至邊界，粒子的擴散過程只受結(jié)構(gòu)內(nèi)基質(zhì)影響。由圖6(b)可知，粒子行走7 000步后，少量粒子行走至分形結(jié)構(gòu)邊界，此時粒子的擴散過程不僅受到分形結(jié)構(gòu)內(nèi)基質(zhì)的阻礙作用，還受到邊界效應(yīng)的影響。由圖6(c)中觀察到，粒子行走10 000步后，行走至邊界處的粒子數(shù)逐漸增多。

圖6 四塊拼接分形結(jié)構(gòu)中粒子行走不同步數(shù)后的空間分布Fig.6 The spatial distribution of particles after different steps in four spliced fractal structure

粒子行走2 500、7 000、10 000步后的密度分布如圖7所示。從圖7中觀察到，粒子的密度分布近似對稱。起始點周圍粒子密度較大，并由中心位置向四周逐漸減小。隨著行走步數(shù)增加，粒子擴散位移及密度范圍增大。同時從圖7中可以觀察到，基質(zhì)對粒子擴散過程的阻礙作用明顯。

圖7 四塊拼接分形結(jié)構(gòu)中粒子行走不同步數(shù)后的密度分布Fig.7 The density distribution of particles after different steps in four spliced fractal structure

圖8為粒子行走2 500、7 000、10 000步后的均方位移與時間的關(guān)系曲線。表3列舉了粒子行走500～10 000步時對應(yīng)的Hurst值。

由圖8和表3分析可知，在分形結(jié)構(gòu)尺度為2 000×2 000 pixel、孔隙分形維數(shù)為1.98的四塊拼接分形結(jié)構(gòu)中，粒子行走500～10 000步的過程中，Hurst值均小于0.5且數(shù)值波動較小，擴散行為均呈現(xiàn)次擴散現(xiàn)象。

圖8 四塊拼接分形結(jié)構(gòu)中粒子擴散均方位移與時間的關(guān)系曲線Fig.8 The relationship between mean square displacement and time of particles in four spliced fractal structure

2.3 數(shù)值模擬結(jié)果的有效性

為論證上述模擬結(jié)果的有效性，將模擬條件設(shè)定為：在單塊分形結(jié)構(gòu)中心位置投放5 000個粒子，采用“全瞎”的螞蟻隨機行走模式模擬擴散過程。假定粒子在孔隙中的行走步長為1 pixel，在基質(zhì)中不能行走。擴散過程中沒有粒子行走至分形結(jié)構(gòu)邊界。

記錄粒子行走5 000、10 000、20 000步后的擴散均方位移與時間的關(guān)系。如圖9所示，通過數(shù)值模擬，計算得到Hurst參數(shù)的平均值為0.477 0。該模擬結(jié)果同與本文模擬條件相同的相關(guān)研究的數(shù)據(jù)：0.462±0.003[38]、0.477±0.004[39]接近，且處于95%的置信限內(nèi)。因此，數(shù)值模擬結(jié)果的有效性可以在一定程度上得到驗證。

圖9 粒子擴散均方位移與時間的關(guān)系曲線Fig.9 The relationship between mean square displacement and time of particles

3 結(jié)果與討論

由上述模擬結(jié)果可知，粒子在單塊分形結(jié)構(gòu)與四塊拼接分形結(jié)構(gòu)中的擴散行為存在差異。當(dāng)粒子行走2 500步時，在單塊分形結(jié)構(gòu)中，已存在粒子擴散至邊界，其擴散行為受基質(zhì)和邊界的共同影響；而在四塊拼接分形結(jié)構(gòu)中，此時粒子并未擴散至邊界，粒子未受到邊界效應(yīng)影響。行走7 000、10 000步后，雖然已出現(xiàn)擴散至邊界的粒子，但與單塊分形結(jié)構(gòu)中的行走過程相比，邊界效應(yīng)對其擴散行為的影響較小。但由于單塊分形結(jié)構(gòu)及四塊拼接分形結(jié)構(gòu)的尺度及結(jié)構(gòu)均不相同，為分別研究結(jié)構(gòu)和尺度對污染物輸運過程的影響，需對粒子在面積放大四倍的單塊分形結(jié)構(gòu)中的擴散行為進行討論。

面積放大四倍的單塊分形結(jié)構(gòu)的空間結(jié)構(gòu)與單塊分形結(jié)構(gòu)相同。其空間尺度為 2 000 pixel×2 000 pixel，與四塊拼接分形結(jié)構(gòu)的尺度相同。通過粒子在面積放大四倍的單塊分形結(jié)構(gòu)中的擴散模擬，可以得到粒子在此結(jié)構(gòu)中行走500步至10 000步時對應(yīng)的Hurst值，列于表4。由表2～表4所列的Hurst值可以得到粒子在上述三種結(jié)構(gòu)中隨行走步數(shù)變化的Hurst值折線圖，如圖10所示。

表4 面積放大四倍的分形結(jié)構(gòu)中的Hurst值Table 4 Hurst value in the single fractal structure with four times larger area

由圖10分析可知，在單塊分形結(jié)構(gòu)中，隨著行走步數(shù)增加，Hurst值先在小于0.5處保持較小波動，然后增大到0.5附近波動后急劇減??；在面積放大四倍的單塊分形結(jié)構(gòu)中，隨著行走步數(shù)增加，Hurst值先在小于0.5處波動，然后增大到0.5附近；在四塊拼接分形結(jié)構(gòu)中，隨著粒子行走步數(shù)增加，Hurst值穩(wěn)定在小于0.5處。由上述分析可知，粒子在三種分形結(jié)構(gòu)中的擴散行為均存在差異。分形結(jié)構(gòu)及空間尺度的差異影響粒子的擴散行為。在空間尺度較小的單塊分形結(jié)構(gòu)中，粒子的擴散行為會在短時間內(nèi)發(fā)生變化，其擴散過程更容易受到結(jié)構(gòu)內(nèi)部及結(jié)構(gòu)邊界的共同影響。

圖10 三種分形結(jié)構(gòu)中隨粒子行走步數(shù)變化的Hurst值折線圖Fig.10 The line chart of Hurst value varying with steps taken by particles in three fractal structures

4 結(jié)論

通過建立兩種土壤分形結(jié)構(gòu)：單塊分形結(jié)構(gòu)、四塊拼接分形結(jié)構(gòu)，分別考察了尺度和結(jié)構(gòu)差異明顯的多孔介質(zhì)結(jié)構(gòu)。通過粒子在這兩種分形結(jié)構(gòu)中的隨機行走過程模擬，分析了污染物在土壤結(jié)構(gòu)中的擴散行為，研究并得出以下結(jié)論。

(1)粒子在單塊分形結(jié)構(gòu)、四塊拼接分形結(jié)構(gòu)中的擴散行為受分形結(jié)構(gòu)中基質(zhì)的阻礙作用，粒子的擴散行為呈次擴散現(xiàn)象。當(dāng)粒子行走至分形結(jié)構(gòu)邊界后，其擴散行為受邊界影響，短暫向正常擴散行為過渡。

(2)分形結(jié)構(gòu)的尺度差異對粒子的擴散行為有明顯影響。當(dāng)粒子的擴散過程未受到邊界效應(yīng)影響時，小尺度與大尺度分形結(jié)構(gòu)中的粒子擴散行為均呈現(xiàn)次擴散現(xiàn)象。然而，當(dāng)粒子的擴散過程受到邊界效應(yīng)影響時，在小尺度分形結(jié)構(gòu)中，粒子擴散行為在短時間內(nèi)發(fā)生明顯變化，由次擴散現(xiàn)象變?yōu)檎U散現(xiàn)象再過渡到次擴散現(xiàn)象；而在大尺度分形結(jié)構(gòu)中，粒子的擴散行為在相當(dāng)長的時間內(nèi)保持穩(wěn)定，擴散行為保持次擴散輸運狀態(tài)。

(3)不同的分形結(jié)構(gòu)會導(dǎo)致粒子的擴散行為產(chǎn)生較大差異。若分形結(jié)構(gòu)的空間尺度一定，在單塊分形結(jié)構(gòu)中，粒子在受到邊界效應(yīng)影響時擴散行為會發(fā)生明顯變化，由次擴散現(xiàn)象過渡為正常擴散現(xiàn)象；而這種變化在內(nèi)含重復(fù)分形結(jié)構(gòu)單元的拼接分形結(jié)構(gòu)中并不明顯，次擴散行為相對穩(wěn)定。因此，由于自然環(huán)境中的土壤結(jié)構(gòu)各異，選擇合適的分形結(jié)構(gòu)模型才能夠更準(zhǔn)確的模擬污染物等溶質(zhì)在土壤中的輸運行為。