張 洋,王 元,王 麗

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

0 引 言

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

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

1 蠕移運(yùn)動(dòng)的高速光學(xué)測量

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

圖1 蠕移實(shí)驗(yàn)平臺(tái)Fig.1 Sketch of experiment system

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

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

2.1 顆粒識別

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

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

2.2 顆粒匹配及軌跡提取

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

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

3 實(shí)驗(yàn)結(jié)果及分析

3.1 蠕移速度分布

一條步長為nT的軌跡可通過差分方法提取nT -1個(gè)速度[11]。這些速度均為蠕移速度v,軌跡第一步和最后一步的速度分別稱為蠕移起動(dòng)速度v1和蠕移終止速度v2。蠕移運(yùn)動(dòng)具備一個(gè)特性：①不同于躍移軌跡的明顯加速過程,蠕移運(yùn)動(dòng)軌跡速度變化較小(最后時(shí)刻驟停)。若一條軌跡滿足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ì)分析帶來偏差。對全體沙樣進(jìn)行V和V1的統(tǒng)計(jì)分析,發(fā)現(xiàn)其分布均存在雙峰現(xiàn)象,見圖5。如1節(jié)所述,以第二組沙為例。

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

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

3.2 蠕移運(yùn)動(dòng)比重隨U＊的變化

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

ηm1和ηm2分別代表兩種蠕移運(yùn)動(dòng)在整體蠕移運(yùn)動(dòng)中的比重,η′m1和η′m2分別代表兩種蠕移起動(dòng)方式在整體蠕移運(yùn)動(dòng)中的比重,η″m1和η″m2分別代表全體“起動(dòng)”中兩種起動(dòng)方式的比重;η1代表“起動(dòng)”在全體瞬時(shí)運(yùn)動(dòng)中的比重。

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

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

4 結(jié) 論

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

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

(3)兩種蠕移運(yùn)動(dòng)的相對比重變化與沙樣特征粒徑以及躍移運(yùn)動(dòng)關(guān)系密切。因此發(fā)展蠕、躍同時(shí)測量技術(shù)將是研究床面運(yùn)動(dòng)形式相互作用及轉(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ì)算機(jī)模擬[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é)報(bào),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.