王茂枚，李志鵬，趙 鋼*，魯程鵬，孫洪廣

(1.江蘇省水利科學(xué)研究院，南京 210017； 2.河海大學(xué) 水文水資源學(xué)院，南京 210098； 3.河海大學(xué) 力學(xué)與材料學(xué)院，南京 211100)

拋石護(hù)岸是河道治理的傳統(tǒng)方法之一，作為一種在國內(nèi)外廣泛使用的防護(hù)技術(shù)，其工程造價小，實施方法簡便，被廣泛的應(yīng)用于各類河道治理工程中[1]。然而，不同地域的各類河流自然條件差異巨大，在工程上廣泛使用的經(jīng)驗公式往往與實際觀測產(chǎn)生較大的出入[2]，因此研究拋石的輸運行為在河道治理、工程應(yīng)用上成為了迫切需要解決的問題。目前國內(nèi)外對拋石輸運行為的研究方法主要以實驗和傳統(tǒng)水沙數(shù)值模型為主[3-4]。然而研究發(fā)現(xiàn)，河流中泥沙及拋石的輸運行為往往是隨機反常的，即使同樣的水沙條件，河床中礫石的啟動行為也是一個概率事件[5]。在拋石的輸運過程中，由于受到河床及流場演變的影響，使得拋石的輸運行為往往是非馬爾科夫性質(zhì)的，因此基于菲克定律的經(jīng)典方程在描述該類輸運行為時往往達(dá)不到預(yù)期效果。這類具有非馬爾可夫性質(zhì)的輸運行為被定義為“反常輸運”[6]。

近年來，“反常輸運”的研究受到國內(nèi)外的廣泛關(guān)注?！胺闯］斶\”在金融、醫(yī)學(xué)、半導(dǎo)體以及地質(zhì)領(lǐng)域上都得到了成功的應(yīng)用[7-8]。目前該方向研究以兩類方法為主，一類是以確定性方程為基礎(chǔ)的，例如：“分?jǐn)?shù)階導(dǎo)數(shù)”、“Hausdorff導(dǎo)數(shù)”[9]；而另一類方法以統(tǒng)計理論為基礎(chǔ)，基于事件發(fā)生的統(tǒng)計性質(zhì)得到輸運行為的隨機運動方程[10]。本文所應(yīng)用的連續(xù)時間隨機行走理論就是基于隨機運動方程描述拋石的輸運行為的。

連續(xù)時間隨機行走模型是研究“反常輸運”行為最直觀的工具，在連續(xù)時間隨機行走理論框架中，一次成功的行走被分解為兩部分：(1)隨機的空間跳躍間隔；(2)兩次成功跳躍之間的等待時間；兩者分布由各自的概率分布密度支配[6，11]。本文應(yīng)用連續(xù)時間隨機行走模型，通過長尾的等待時間及跳躍步長分布模擬拋石在河道中的“反常輸運”行為，并與實驗室尺度下的實驗觀測數(shù)據(jù)對比進(jìn)一步驗證了模型對于拋石輸運行為預(yù)測的有效性。

1 連續(xù)時間隨機行走理論

連續(xù)時間隨機行走模型中粒子的跳躍過程由空間跳躍距離和兩次成功跳躍之間的等待時間組成，兩者由共同的聯(lián)合概率密度分布決定[6，8]。在實際研究中發(fā)現(xiàn)跳躍步長與等待時間往往是不相關(guān)的，即跳躍步長和等待時間分別由各自的分布密度決定。研究中λ(x)，ω(t)表示跳躍步長與等待時間的概率密度函數(shù)，則聯(lián)合概率密度函數(shù)可表示為

ψ(x,t)=λ(x)·ω(t)

(1)

連續(xù)時間隨機行走模型的廣義主方程可描述為

(2)

式中：η(x,t)表示粒子t時刻到達(dá)x處的概率密度，粒子的初始分布為P0(x)。

粒子的概率密度函數(shù)P(x,t)滿足方程

(3)

將廣義主方程(2)代入(3)式,最終得到連續(xù)時間隨機行走的概率方程為

(4)

在連續(xù)時間隨機行走框架中，空間和時間都被定義為連續(xù)變量，因此對方程(4)進(jìn)行Fourier-Laplace變換，得

(5)

傅里葉—拉普拉斯域內(nèi)的方程(5)被稱為Montroll-Weiss方程，該方程中通過跳躍步長和等待時間的概率密度函數(shù)的性質(zhì)可以區(qū)分不同輸運行為，當(dāng)跳躍步長的概率密度分布二階矩發(fā)散而等待時間的概率密度分布一階矩收斂時，輸運行為定義為超擴散輸運；跳躍步長的概率密度分布二階矩收斂而等待時間概率密度分布一階矩發(fā)散時，輸運行為定義為次擴散輸運；當(dāng)兩者皆發(fā)散，即同時具有超擴散輸運行為和次擴散輸運行為，兩者皆收斂則為正常輸運行為。

2 蒙特卡羅模擬方法

蒙特卡羅方法，又稱隨機模擬方法，利用重復(fù)隨機的抽樣方法，求得統(tǒng)計量的數(shù)字特征。當(dāng)重復(fù)次數(shù)足夠多時，蒙特卡羅方法得到的數(shù)值解逼近真實解。本研究中由于方程(4)在轉(zhuǎn)換到時域上具有較大的困難，因此采用蒙特卡羅方法求得其數(shù)值結(jié)果，該方法是典型的拉格朗日法，可以追蹤粒子運動過程中任意時刻的位置，完整呈現(xiàn)粒子運動軌跡。

經(jīng)典的對流擴散方程(ADE)表征流體的質(zhì)量傳輸規(guī)律，并能描述滿足菲克定律的流體中粒子輸運行為[9]。但實際觀測中發(fā)現(xiàn)很多情況下菲克定律會被打破，例如含水層的異質(zhì)性、河床結(jié)構(gòu)的非均勻性都會導(dǎo)致ADE方程無法準(zhǔn)確描述該類環(huán)境中粒子的輸運行為。這類打破經(jīng)典菲克定律的輸運行為被稱為“反常輸運”。連續(xù)時間隨機行走理論描述輸運行為時，當(dāng)概率密度分布滿足中心極限定理時，即跳躍步長分布為正態(tài)分布，等待時間分布為指數(shù)分布時，方程(5)將回歸到ADE方程。在本研究中，為了反映河床結(jié)構(gòu)的異質(zhì)性導(dǎo)致的拋石“反常輸運”行為，采用 穩(wěn)定分布作為跳躍步長分布，Mittag-Leffler分布作為等待時間分布。

本研究中我們將拋石的運動行為劃分為沿水流方向和垂直于水流方向的運動。因此，運動過程的時間計量由下式確定

tn+1=tn+τβ
t0=0

(6)

式中：τβ為滿足Mittag-Leffler分布的隨機數(shù)。

沿水流方向的位移由流場產(chǎn)生的確定位移和粒子相互碰撞的噪聲項確定，表述為

(7)

式中：v表示粒子由對流產(chǎn)生的速度，可由實驗獲得。

垂直于水流方向的位移僅僅由粒子的相互碰撞產(chǎn)生，由下式確定

(8)

式中：ξα1和ξα2是滿足α穩(wěn)定分布的隨機數(shù)。研究中隨機數(shù)的生成分別通過以下公式生成

其中滿足跳躍步長的α穩(wěn)定分布隨機數(shù)采用式(9)生成[12]

(9)

滿足Mittag-Leffler分布的等待時間隨機數(shù)采用式(10)生成[13]

(10)

式中：u1,u2∈(0，1)為均勻分布的隨機數(shù)，γx，γt是隨機數(shù)產(chǎn)生的尺度系數(shù)。

圖1驗證了產(chǎn)生隨機數(shù)的概率密度分布以及不同參數(shù)下α穩(wěn)定分布和Mittag-Leffler分布的變化特征。通過圖1-a與1-c隨機數(shù)概率密度與實際分布對比，隨機數(shù)的生成很好的滿足了所需概率密度分布。進(jìn)一步分析不同參數(shù)下分布的概率密度特性，發(fā)現(xiàn)α穩(wěn)定分布在參數(shù)α=2時回歸到正態(tài)分布，并隨著α的減小分布逐漸變“寬”，說明隨著參數(shù)α的減小，粒子的跳躍步長愈加發(fā)散，超擴散輸運行為越強烈；Mittag-Leffler分布具有和α穩(wěn)定分布相同的性質(zhì)，當(dāng)參數(shù)β=1時，Mittag-Leffler分布回歸到指數(shù)分布，當(dāng)參數(shù)β變小，等待時間的分布更加發(fā)散，表明其次擴散行為越強烈。

1-a 穩(wěn)定分布隨機數(shù)概率密度驗證1-b 不同參數(shù)下穩(wěn)定分布的狀態(tài)1-c Mittag-Leffler 分布隨機數(shù)概率密度驗證1-d 不同參數(shù)下Mittag-Leffler分布的狀態(tài)圖1 跳躍步長與等待時間隨機數(shù)概率密度驗證及不同參數(shù)下概率密度函數(shù)的差異Fig.1 Verification of the probability of random numbers for jump steps and waiting times, and the difference in probability density function under various parameters

圖2 Martin[2012] 實驗結(jié)果及連續(xù)時間隨機行走模型和ADE模型預(yù)測效果比較Fig.2 Martin [2012] experimental results and comparison of predictive effects between continuous time random walk model and ADE

3 實驗結(jié)果預(yù)測

研究中為了驗證連續(xù)時間隨機行走模型在實際工況中的預(yù)測效果，本節(jié)對兩組不同源條件下實驗室尺度水槽實驗結(jié)果進(jìn)行預(yù)測。表1給出了兩組實驗數(shù)據(jù)下連續(xù)時間隨機行走模型參數(shù)的匹配情況；比較發(fā)現(xiàn)，在實際運動中拋石的輸運行為往往是次擴散的(均方位移指數(shù)小于1)，但運動過程中同時表現(xiàn)出超擴散行為(1<α<2)和次擴散行為( 0<β<1)。

3.1 瞬時源條件下實驗驗證及分析

圖3 連續(xù)時間隨機行走模型預(yù)測Martin[2012] 試驗的時-空分布及均方位移Fig.3 Time-Space evolution and MSD using continuous time random walk model

此處，我們選擇Martin于2012年在瞬時源條件下進(jìn)行的實驗[14]。實驗結(jié)果表明，礫石在河床上的分布并不為對稱形式的正態(tài)分布，反而表現(xiàn)出一定程度的空間右拖尾現(xiàn)象，這種現(xiàn)象可能是由于空間的異質(zhì)性導(dǎo)致了礫石在河床上的運動出現(xiàn)了反常輸運的現(xiàn)象。圖2比較了連續(xù)時間隨機行走模型與經(jīng)典ADE方程在實驗數(shù)據(jù)預(yù)測上的優(yōu)劣，通過比較發(fā)現(xiàn)，連續(xù)時間隨機行走模型很好的捕捉了礫石的拖尾現(xiàn)象，而ADE方程不能捕捉具有這類拖尾現(xiàn)象的行為。同時圖3粒子的時-空分布及均方位移指數(shù)說明該實驗中礫石在河床上的運動行為是次擴散的，即一部分礫石被水流及河床結(jié)構(gòu)加速，但河床結(jié)構(gòu)整體上對礫石運動有阻隔作用，故礫石的整體的運動效應(yīng)表現(xiàn)為次擴散輸運狀態(tài)。

3.2 連續(xù)源實驗驗證及分析

為了驗證連續(xù)時間隨機行走理論在連續(xù)源條件下以及短時對長時的預(yù)測效果，我們在這里選取了Chang and Yen于2002年進(jìn)行的連續(xù)源實驗[15]。圖4對比了連續(xù)時間隨機行走模型和經(jīng)典ADE方程的預(yù)測結(jié)果，此處我們采用初始時刻t=7.5 s作為模型參數(shù)的濾定參考項，之后各個時刻均采用初始時刻的參數(shù)。相比于經(jīng)典ADE方程結(jié)果，連續(xù)時間隨機行走模型在短時間濾定參數(shù)進(jìn)而預(yù)測長時結(jié)果時具有更好的效果。進(jìn)一步分析發(fā)現(xiàn)，ADE方程在短時預(yù)測長時的效果上，隨著時間的增加其預(yù)測誤差會逐漸增大，圖4-a明顯表現(xiàn)出ADE方程在描述空間尾部效應(yīng)時，其曲線尾部的衰減以指數(shù)形式衰減，而連續(xù)時間隨機行走模型尾部以冪律衰減，這是兩類模型在后期預(yù)測效果產(chǎn)生巨大差異的重要原因。

4-a t=7.5 s4-b t=22.5 s4-c t=37.5 s

4-d t=52.5 s4-e t=67.5 s4-f t=82.5 s圖4 Chang and Yen[2002] 實驗結(jié)果及連續(xù)時間隨機行走模型和ADE模型預(yù)測效果比較Fig.4 Chang and Yen[2002] experimental results and comparison of predictive effects between continuous time random walk model and ADE model

表1 連續(xù)時間隨機行走模型參數(shù)及實驗數(shù)據(jù)對比Tab.1 Continuous time random walk model parameters and experimental data

4 模擬結(jié)果分析

為了表征不同輸運行為在拋石輸運過程中的影響，本節(jié)對比了不同位置、不同拋石方式下，拋石的運動行為。分析中參數(shù)采用無量綱化處理,粒子的平均速度統(tǒng)一為v=0.5。

4.1 上游中心點源釋放

粒子于上游中心釋放，對比了正常輸運行為、次擴散輸運行為以及超擴散輸運行為下拋石的空間分布。圖5比較了次擴散-正常擴散-超擴散的輸運行為差異，隨著超擴散行為的增強，拋石的空間離散性逐漸增大。同時由于次擴散產(chǎn)生的長程等待時間，在次擴散狀態(tài)下相當(dāng)一部分拋石停留在離點源位置較近的區(qū)域。圖中X方向是粒子具有對流速度的方向，箭頭表示水流方向，Y方向是粒子僅具有隨機彌散位移的方向。

5-a 正常輸運行為下粒子不同時刻的空間分布5-b 次擴散輸運行為下粒子不同時刻的空間分布5-c 超擴散輸運行為下粒子不同時刻的空間分布圖5 沿中心釋放拋石Fig.5 Release of riprap along the center

4.2 上游岸邊點源釋放

將拋石的釋放位置放置于岸邊，對比正常擴散-次擴散-超擴散三種輸運行為下拋石的空間分布。比較圖6得出，在超擴散輸運狀態(tài)下，拋石會在最短時間到達(dá)對岸，這是由于跳躍步長概率密度分布的二階矩發(fā)散導(dǎo)致拋石在短時間產(chǎn)生長程的輸運行為；隨著次擴散輸運行為的增強，拋石到達(dá)對岸的時間逐漸增加。

6-a 正常輸運行為下粒子不同時刻的空間分布6-b 次擴散輸運行為下粒子不同時刻的空間分布6-c 超擴散輸運行為下粒子不同時刻的空間分布圖6 沿上游河岸釋放拋石Fig.6 Release of riprap along the upper reaches of the riverside

4.3 上游沿水流方向釋放拋石以及同時釋放不同粒徑拋石

研究中對比了以線源方式釋放拋石及同時釋放不同粒徑拋石的情況。由于河床上由不同粒徑的砂石組成的聚集體對一些拋石有阻隔作用，而另一些拋石有概率進(jìn)入河床上的裂隙結(jié)構(gòu)，這樣的裂隙結(jié)構(gòu)使進(jìn)入其中的拋石快速向下游運動，導(dǎo)致實際運動中拋石同時具有超擴散和次擴散輸運行為。為表征這類輸運行為，研究中參數(shù)設(shè)置為α=1.6，β=0.6。圖7-a展示了以線源釋放后，拋石的空間分布較為均勻，對水流的阻隔效果最好。圖7-b顯示次擴散行為較強的粗粒徑拋石更易停留在離源較近的區(qū)域，而超擴散行為更強的細(xì)粒徑拋石會更快向下游運動。

7-a 垂直流速方向上游連續(xù)釋放拋石的空間分布7-b 不同粒徑拋石同時釋放后的空間分布圖7 連續(xù)源釋放拋石Fig.7 Release of riprap under continuous source

4.4 不同輸運行為的MSD及時間-空間演化分析

物理上通常將輸運行為用均方位移來表征，其中滿足菲克定律的正常輸運行為描述為

〈x2〉∝t

(11)

然而，當(dāng)輸運行為的均方位移滿足下式時，則被定義為反常輸運行為

〈x2〉∝tγ(γ≠1)

(12)

式中：0<γ<1時為次擴散輸運行為，當(dāng)γ>1時為超擴散輸運行為。

8-a 次擴散輸運行為下時-空演化及均方位移對比8-b 正常輸運行為下時-空分布及均方位移8-c 超擴散輸運行為下時-空分布及均方位移圖8 三種擴散行為下粒子的空間演化圖及均方位移圖Fig.8 The space evolution and mean square displacement of particles under three diffusion behaviors

通過圖8-a、8-b、8-c比較可以發(fā)現(xiàn)隨著次擴散輸運狀態(tài)向正常輸運狀態(tài)再向超擴散輸運狀態(tài)轉(zhuǎn)變，粒子的空間分布逐漸變寬，隨著次擴散向超擴散現(xiàn)象的過渡，粒子由中心集中逐漸向空間均勻分布過渡。

5 結(jié)論

本研究通過連續(xù)時間隨機行走模型結(jié)合蒙特卡羅方法對拋石的輸運行為進(jìn)行研究，并與室內(nèi)水槽實驗結(jié)果進(jìn)行比較，得出了以下2點結(jié)論：

(1)連續(xù)時間隨機行走模型通過“寬”分布的等待時間和跳躍步長從微觀粒子隨機運動的角度解釋了拋石反常輸運的物理機理。

(2)通過與室內(nèi)試驗比較得出，連續(xù)時間隨機行走模型在捕捉拋石反常輸運行為時效果比經(jīng)典的對流擴散方程好；拋石在河道中的輸運行為同時存在超擴散輸運行為和次擴散輸運行為。