張 洋,王 元,王 麗

(西安交通大學(xué)流體機械及工程系,西安 710049)

0 引 言

風(fēng)沙流中沙粒在外力作用下沿地表的滾動和滑動稱為蠕移,作為風(fēng)吹沙粒的三種基本運動方式之一,它在風(fēng)沙動力學(xué)研究中有著舉足輕重的作用[1]。一方面,蠕移是沙床中較大顆粒的主要運動方式;另一方面,蠕移是床面構(gòu)型變化的主要原因：由于蠕移和躍移粒徑及速度的差異,沙床粒配會隨風(fēng)沙遷移發(fā)生變化,沙粒的堆積和輸移交替發(fā)生,最終導(dǎo)致沙波紋等床面構(gòu)型的發(fā)展。

由于測量上的困難,目前研究沙粒蠕移的文獻(xiàn)較少,已有的工作也僅以床面構(gòu)型模擬及輸沙率推算為主。Werner[2]應(yīng)用離散元法(DEM)、王光謙等[3]通過給定蠕移動量交換規(guī)則模擬床面沙波紋;鄭曉靜等運用粒子流動模型對沙粒蠕移流量進(jìn)行預(yù)測,并提出離散粒子追蹤法(DPTM)來確定蠕移沙粒終止位置[4-5], Dong等[6]通過風(fēng)洞實驗中輸沙率隨高度分布曲線來反推蠕移平均輸沙率。這些工作都需要更為精細(xì)的實驗支持。在高速光學(xué)測量的基礎(chǔ)上[7]發(fā)展了系統(tǒng)性的沙粒蠕移軌跡計算方法,并對兩種按起動方式劃分的蠕移運動速度進(jìn)行了統(tǒng)計分析。

1 蠕移運動的高速光學(xué)測量

實驗在西安交通大學(xué)流體工程系的大氣邊界層(ABL)風(fēng)洞中進(jìn)行[7],實驗段尺寸為7.4m×0.6m× 0.5m。采用PCO-1200hs高速相機以俯視角拍攝蠕移運動,圖像分辨率1280pixel×1024pixel,采集頻率500Hz。如圖1,采取兩項措施：(1)環(huán)繞風(fēng)洞壁放置高亮度鹵素?zé)?1kW)形成無影體照明系統(tǒng),最大程度消除躍移運動顆粒在床面上的投影影響;(2)相機前加裝具備大光圈的長焦鏡頭,并在鏡頭與相機間加置數(shù)量可變的延伸環(huán)來盡可能地縮短工作距離。上述措施將拍攝景深控制在理想范圍之內(nèi),從而使蠕移層之上的運動顆粒無法被聚焦拍攝,保證了采集對象的單一性。

圖1 蠕移實驗平臺Fig.1 Sketch of experiment system

實驗采用8組沙樣：＜100μ m,100～125μ m, 125～200μ m,200～250μ m,250～400μ m,400～500μ m,500～800μ m,800～1000μ m。對每個工況重復(fù)三次試驗,每次試驗采集一個由3000個連續(xù)瞬時組成的圖像序列;所提取的蠕移軌跡總數(shù)量均在4000條以上,且統(tǒng)計計算包含所有這些軌跡,統(tǒng)計分析的可靠性得到保障。如3.1所述,統(tǒng)計分析蠕移速度分布中的“雙峰現(xiàn)象”在前4組結(jié)果中較好地復(fù)現(xiàn),但在后4組中因粒徑較大、躍移運動十分微弱而不明顯,原因見3.1;而且受篇幅限制,3.1中的討論僅以第二組100～125μ m為例。

2 圖像處理及蠕移軌跡提取

2.1 顆粒識別

床面蠕移灰度圖像如圖2(a)所示。幀間顆粒匹配是獲取蠕移軌跡和單顆粒速度的前提,相對于計算微區(qū)域平均歐拉速度的傳統(tǒng)PIV算法,拉格朗日假設(shè)下的PTV(particle tracking velocimetry,粒子追蹤測速)算法無疑更為合適。為達(dá)到PTV的計算前提,需提取單像素化的蠕移顆粒。首先,運用非負(fù)圖像減法[7]去除不動顆粒組成的圖像背景：

其中,Ai-1、Ai、Ai+1分別代表i-1、i、i+1瞬時的灰度矩陣,S()代表非負(fù)算子。由于非負(fù)圖像減法的中心差分格式,圖像序列首尾兩張不作運算;得到的僅包含運動沙粒及圖像噪點。然后利用自動閾值法[8]較為準(zhǔn)確地實現(xiàn)粒子單像素化及圖像噪點的消除,見圖2(b)。

2.2 顆粒匹配及軌跡提取

在獲取單像素化粒子群位置坐標(biāo)的基礎(chǔ)上,運用雙向DT-PTV算法[9](Bidirectional Delaunay Tes-sellation-based PTV,基于三角分割的雙向PTV)進(jìn)行幀間粒子匹配,粒子群Delaunay網(wǎng)格以及用矢量顯示的匹配結(jié)果分別見圖2(c)、(d)(圖2數(shù)字坐標(biāo)的單位為pixel)。一方面,從假設(shè)的簡潔性、運算速度、準(zhǔn)確率及處理對象為離散相等多角度看,DTPTV優(yōu)于一些經(jīng)典PTV算法[10](諸如文獻(xiàn)[7]所采用的BICC-PTV);另一方面,經(jīng)過“雙向算法”升級的DT-PTV與其單向算法相比具備更高的準(zhǔn)確率[9]。故經(jīng)過非負(fù)圖像減法、自動閾值法以及雙向DTPTV的依次運算,可得到較為準(zhǔn)確的幀間匹配關(guān)系。蠕移軌跡從連續(xù)幀匹配關(guān)系中提取。圖3顯示摩阻風(fēng)速U＊=0.7m/s下第三組沙在采集時間T=6s內(nèi)的部分蠕移軌跡。

圖3 蠕移軌跡Fig.3 Trajectories of creeping sand

3 實驗結(jié)果及分析

3.1 蠕移速度分布

一條步長為nT的軌跡可通過差分方法提取nT -1個速度[11]。這些速度均為蠕移速度v,軌跡第一步和最后一步的速度分別稱為蠕移起動速度v1和蠕移終止速度v2。蠕移運動具備一個特性：①不同于躍移軌跡的明顯加速過程,蠕移運動軌跡速度變化較小(最后時刻驟停)。若一條軌跡滿足0.8≤|v2/v1|≤1.2,則稱其歷程穩(wěn)定;歷程穩(wěn)定的軌跡數(shù)量與軌跡總數(shù)量的比值稱為軌跡穩(wěn)定率ηV。以第二組沙(100～125μ m)為例,所有摩阻風(fēng)速U＊對應(yīng)的ηV值均不低于88%,如圖4所示。設(shè)某一工況下所有軌跡v和v1的絕對值的集合分別為V和V1,顯然V1?V。這樣處理的目的在于將速度值與軌跡步長綜合考慮從而避免給統(tǒng)計分析帶來偏差。對全體沙樣進(jìn)行V和V1的統(tǒng)計分析,發(fā)現(xiàn)其分布均存在雙峰現(xiàn)象,見圖5。如1節(jié)所述,以第二組沙為例。

圖5(a)為V的數(shù)量分布。ns是速度為V的顆粒的數(shù)量。分布中的雙峰實際上代表按起動方式劃分的兩種蠕移顆粒：風(fēng)吹起動顆粒和碰撞起動顆粒(以m1、m2表示),它們分別指在流體作用下沿床面的滾動,以及在顆粒(以躍移為主)的高速碰撞下擺脫靜止?fàn)顟B(tài)、但因能量或起動角不足以躍移而只能沿床表面運動。一般地,m2的速度大于m1速度。如圖5 (a),隨著U＊上升,沙床獲得的能量增大且躍移運動規(guī)模增大,故m1、m2數(shù)量均開始增大;當(dāng)U＊超過一定閾值后,越來越多的顆粒將獲得較大能量并以躍移方式脫離地面,m1數(shù)量開始降低。由于蠕移的運動特性①,V1與V的數(shù)量分布類似。由圖5得到蠕移的另一個運動特性：②特征粒徑確定的蠕移顆粒速度范圍變化極小,受摩阻風(fēng)速影響甚微。

圖5(b)表示V的比重分布,η是速度為V的瞬時顆粒在全體瞬時顆粒中所占比重;5(c)表示V1的相對比重分布,ˉη是速度為V1的起動顆粒在全體起動顆粒中所占比重;二者的相似性從另一側(cè)面反映出蠕移運動的特性①。在①②的前提下,可通過臨界速度VC將V(或V1)劃為分別代表m1、m2的速度Vm1和(或)的計算方法可簡述為：相同統(tǒng)計區(qū)間下各摩阻風(fēng)速對應(yīng)的比重值所組成的數(shù)列具有一個表征大小順序的度數(shù),則存在唯一的VC將全體數(shù)列分為兩個順序度相反的部分。圖5(d)表示蠕移起動速度V1的絕對比重分布,η在此表示速度為的起動顆粒在全體瞬時顆粒中所占比重?？梢钥闯?隨著U＊上升,比重基本穩(wěn)定,而比重上升,這意味著非起動速度總比重的整體下降。

3.2 蠕移運動比重隨U＊的變化

在得到VC之后,可進(jìn)一步定義若干代表m1和m2在整體運動中所占比重的量,從而分析兩種運動方式的權(quán)重隨摩阻風(fēng)速變化的規(guī)律：

ηm1和ηm2分別代表兩種蠕移運動在整體蠕移運動中的比重,η′m1和η′m2分別代表兩種蠕移起動方式在整體蠕移運動中的比重,η″m1和η″m2分別代表全體“起動”中兩種起動方式的比重;η1代表“起動”在全體瞬時運動中的比重。

如圖6(a)、(b)所示,由于蠕移運動的特性①,隨著U＊增加,ηm1和和的變化趨勢分別相似。如3.1所述,隨著U＊增加,起動條件占優(yōu)的顆粒(如大的迎風(fēng)面積)將越傾向于以躍移方式起動,而其余的顆粒將在規(guī)模增大的躍移顆粒的碰撞下開始蠕移。數(shù)值相對穩(wěn)定,表明了風(fēng)吹蠕移運動在達(dá)到臨界之前是蠕移運動中的穩(wěn)定變化部分,而的趨勢則表明碰撞蠕移運動是蠕移運動中的動態(tài)變化部分。

如圖6(c),無論是在全體起動還是整體運動中,在U＊達(dá)到一個臨界摩阻風(fēng)速之后m2相對于m1的比重開始大于1,這意味著運動形式主導(dǎo)地位的轉(zhuǎn)移。另外,起動的總比重η1隨U＊增加而增加;這是因為m2代表著躍移入射的高隨機性,其比重增加導(dǎo)致沙床表面運動隨機性增強,即單沙粒的實際軌跡更易呈斷續(xù)狀的n個分段,在采集頻率較高的前提下則對應(yīng)為n條步長較短的數(shù)字軌跡;因此非起動速度的整體數(shù)量和比重下降。m2的比重可作為床面隨機性的量化標(biāo)識之一。

4 結(jié) 論

(1)蠕移運動圖像是連續(xù)型無影照明系統(tǒng)下的俯視拍攝結(jié)果,存在不動顆粒組成的背景。而非負(fù)圖像減法、自動閾值法及雙向粒子追蹤算法的有機結(jié)合能夠很好地從此類圖像中提取有效的顆粒軌跡;

(2)蠕移總速度及起動速度的統(tǒng)計結(jié)果表明,當(dāng)近床面同時存在躍移和蠕移運動時,流體起動蠕移代表整個運動的穩(wěn)定變化部分,其比重隨摩阻風(fēng)速增大而減小;碰撞起動蠕移則代表動態(tài)變化部分,其比重隨摩阻風(fēng)速增大而增大。在摩阻風(fēng)速超過臨界值之后,蠕移起動的主導(dǎo)形式將從流體起動轉(zhuǎn)化為碰撞起動;

(3)兩種蠕移運動的相對比重變化與沙樣特征粒徑以及躍移運動關(guān)系密切。因此發(fā)展蠕、躍同時測量技術(shù)將是研究床面運動形式相互作用及轉(zhuǎn)化的重要突破口。

[1] ANDERSON R S,BUNAS K L.Grain size segregation and stratigraphy in aeloian ripples modelled with a cellular automaton[J].Nature,1993,365(2)：740-743.

[2] WERNER B T.A physical model of wind-blown sand transport[D].California,US：California Institute of Technology,1987.

[3] 王光謙,冉啟華,劉綠波.風(fēng)成沙紋過程的計算機模擬[J].泥沙研究,1998(3)：1-6.

[4] WANG Z T,ZHENG X J.Heoretical prediction of creep flux in aeolian sand transport[J].Powder Technology, 2004,139(2)：123-128.

[5] 鄭曉靜,薄天利,謝 莉.風(fēng)成沙波紋的離散粒子追蹤法模擬[J].中國科學(xué)G,2007,37(4)：527-534.

[6] DONG Z B,LIU X P,WANG H T,et al.The flux profile of a blowing sand cloud：a wind tunnel investigation [J].Geomorphology,2002,49(3-4)：219-230.

[7] WANG Y,WANG D W,WANG L,et al.Measurement of sand creep on a flat sand bed using a high-speed digital camera[J].Sedimentology,2009,56(6)：1705-1712.

[8] OHMI K,LI H Y.Particle-tracking velocimetry with new algorithms[J].M eas.Sci.Technol.,2000,11 (6)：603-616.

[9] 張 洋,王 元,王大偉.DT-PTV離散相偽矢量的概率擇優(yōu)式剔除算法[J].西安交通大學(xué)學(xué)報,2009,43 (5)：119-123.

[10]SONG X Q,YAMAMOTO F,IGUCHI M,et al.A new tracking algorithm of PIV and removal of spurious vectors using Delaunay tessellation[J].Exp.Fluids, 1999,26(4)：371-380.

[11]ZHANG W,KANG J H,LEE S J.Tracking of saltating sand trajectories over a flat surface embedded in an atmospheric boundary layer[J].Geomorphology,2007, 86(3)：320-331.