擋風(fēng)墻周圍風(fēng)沙流運動特性SPH數(shù)值分析

胡 曉,金阿芳,熱依汗古麗·木沙

(新疆大學(xué)機械工程學(xué)院,新疆 烏魯木齊 830047)

1 引言

以風(fēng)沙過程為主要標(biāo)志的沙漠化已成為全球最突出的生態(tài)環(huán)境問題之一。我國的沙漠及沙漠化土地面積高達(dá)128萬平方公里,占據(jù)國土面積的13.3%。風(fēng)沙運動對農(nóng)業(yè)、牧業(yè)、交通、水利等造成了極大危害,沙害造成的直接經(jīng)濟損失已達(dá)到數(shù)百億人民幣[1]。風(fēng)沙兩相流作為一門研究沙粒在風(fēng)力作用下運動規(guī)律的科學(xué),奠定了近代風(fēng)沙物理學(xué)的理論基礎(chǔ),是認(rèn)識和研究風(fēng)沙運動規(guī)律的重要組成部分,為防治地表土壤風(fēng)蝕、土地荒漠化提供科學(xué)依據(jù)和理論支撐。因此研究風(fēng)沙兩相流運動對防沙、治沙工程來說具有重要的指導(dǎo)意義。

風(fēng)沙兩相流的研究分為直接研究和間接模擬,通常包括野外觀測、風(fēng)洞模擬和數(shù)值模擬等三種方法。由于前兩種方法容易受外界因素干擾導(dǎo)致實驗數(shù)據(jù)不準(zhǔn)確,加之受場地限制需要耗費大量的人力、物力,所以計算機數(shù)值模擬方法近年來逐漸受到學(xué)者們的廣泛研究和關(guān)注。Anderson 基于元胞自動機原理,模擬了沙床表面顆粒的周期性跳躍過程,并最終演化成沙波紋[2]。鄭曉靜等在求解沙粒垂向脈動方程的基礎(chǔ)上,通過數(shù)值模擬得到風(fēng)沙中隨機沙粒運動軌跡[3]。李萬清等利用隨機方法生成符合自然分布規(guī)律的沙床,進而結(jié)合離散單元法(DEM)建立一套粒-床碰撞數(shù)值模擬系統(tǒng)[4]?？毫姷仁褂妙w粒碰撞的歐拉-拉格朗日數(shù)值模擬方法,對風(fēng)沙躍移中顆粒沖擊多粒徑床面的碰撞過程進行了數(shù)值計算[5]。劉博等采用離散單元法和硬球模型對單顆沙粒碰撞規(guī)律進行數(shù)值分析,模擬出顆粒的運動規(guī)律和受力狀況[6]。何麗紅等基于點電荷場的基本理論,建立一種風(fēng)-沙-電相互作用的多場耦合模型,得到了風(fēng)沙電場沿高度的分布規(guī)律[7]?，F(xiàn)有的研究鮮有涉及用無網(wǎng)格法分析擋風(fēng)墻對風(fēng)沙兩相流運動特性的影響,而在我國西北沙化地區(qū)擋風(fēng)墻已然成為最為經(jīng)濟和應(yīng)用最廣的防風(fēng)結(jié)構(gòu),因此擋風(fēng)墻對躍移沙粒防護作用的數(shù)值模擬研究在實際工程應(yīng)用中更具有實踐意義。

光滑粒子流體動力學(xué)(Smoothed Particles Hydrodynamics,SPH)作為一種純拉格朗日粒子方法,也是最早的無網(wǎng)格方法[8,9]。該方法最初廣泛應(yīng)用于天體物理領(lǐng)域,由于SPH方法本質(zhì)上不需要借助網(wǎng)格,規(guī)避了網(wǎng)格畸變現(xiàn)象,因此能夠處理拉格朗日格式的大變形問題。同時SPH核函數(shù)光滑長度的自適應(yīng)性可以動態(tài)捕捉流場局部變化特征,通過“核函數(shù)”的近似積分將流體力學(xué)基本方程組轉(zhuǎn)換成數(shù)值計算的SPH方程[10],用算法將氣固兩相計算域的粒子離散化,并賦予這些離散粒子所有的力學(xué)量。近年來,SPH方法已拓展應(yīng)用于流體特性的數(shù)值模擬。本文引入虛粒子作為擋風(fēng)墻,在此基礎(chǔ)上對SPH搜索算法進行改進,建立了一種基于SPH方法的擋風(fēng)墻-風(fēng)沙流模型,以此研究自然界中擋風(fēng)墻對躍移沙粒的防護效應(yīng)。

2 SPH方法

2.1 SPH基本原理

SPH方法是近20多年來逐步發(fā)展起來的一種無網(wǎng)格方法,最初是由Lucy于1977 年提出并用來解決天體物理學(xué)中涉及三維開放空間的任意流動問題[11]。經(jīng)過Benz等學(xué)者的不斷改進和完善,逐漸發(fā)展成一種較為成熟的粒子模擬方法[12],推廣到流體動力學(xué)的各個研究領(lǐng)域。其基本思想是將連續(xù)的流場離散成可以運動的粒子,利用核近似法將臨近粒子關(guān)聯(lián)起來,并在粒子上建立離散動力學(xué)方程,然后通過求解離散控制方程來更新粒子的位置和相應(yīng)物理量信息。因其自身的無網(wǎng)格、自適應(yīng)性、光滑性以及拉格朗日公式與近似法的和諧結(jié)合,使得SPH方法在研究大變形、復(fù)雜界面運動等問題時有著獨特的優(yōu)勢。

2.2 核函數(shù)構(gòu)造

光滑核函數(shù)決定函數(shù)近似式的形式、一致性和精度,定義粒子支持域的尺寸,所以在SPH方法中占據(jù)非常重要的地位。核函數(shù)的主要特征可以歸納為:

1) 歸一性 保證連續(xù)函數(shù)積分的零階連續(xù)性

2) 緊支性W(x-x′)=0,|x-x′|>kh

3) 非負(fù)性 對支持域內(nèi)任意點x,當(dāng)x-x′≤kh時,W(x-x′)≥0,|x-x′|>kh

4) 衰減性 核函數(shù)的值隨粒子間間距增大而單調(diào)遞減

5) 對稱性 光滑核函數(shù)需為偶函數(shù)

6) 光滑性 對x支持域內(nèi)任一點x′,W(x-x′,h)均無窮可導(dǎo)

2.3 構(gòu)建風(fēng)沙流SPH基本方程

SPH方程的構(gòu)造通常分為兩個主要步驟來進行,首先為積分表示法,其次為粒子近似法。對任意函數(shù)和光滑函數(shù)逐步積分,然后通過最近相鄰粒子的值進行累加求和近似。

在SPH方法中,f(x)的積分表示式由式(1)定義

(1)

式中,Ω是包含x的積分體積,δ(x-x′)是狄拉克函數(shù),由下式表示

(2)

若用光滑函數(shù)W(x-x′,h)取代f(x)積分表達(dá)式中的狄拉克函數(shù)δ(x-x′),則有

(3)

因為W不是狄拉克函數(shù),故上式的積分表達(dá)式只能近似表示,式中的角括號<>在SPH的習(xí)慣用法中稱為核近似算子標(biāo)記;W(x-x′)為光滑函數(shù);h是定義光滑函數(shù)影響區(qū)域的光滑長度。

將式(3)用粒子近似法轉(zhuǎn)換為支持域中相鄰粒子疊加求和的離散值。則粒子i處函數(shù)的粒子近似形式可表示為

(4)

式中,N是粒子i支持域中所有粒子的總數(shù);xi,xj分別為粒子i和j的坐標(biāo)向量;mj,ρj為粒子j的質(zhì)量和密度。

2.4 控制方程的SPH離散

2.4.1 連續(xù)密度法

對速度散度應(yīng)用SPH近似

(5)

2.4.2 動量方程的粒子近似法

對動量方程等號右端的梯度項直接應(yīng)用SPH粒子近似法變換后可得到

(6)

2.5 雙相耦合的全配對搜索法

風(fēng)沙流運動是氣體粒子和沙粒兩種不同相粒子間的耦合運動?？紤]到不同粒子對周圍粒子的支持域不同,故采用SPH方法里的光滑長度來處理風(fēng)沙兩相流之間的耦合問題。由于目標(biāo)粒子對臨近粒子的搜索隨時間變化而變化,因而需要尋求一種精確、有效的方法進行鄰域搜索。本文針對風(fēng)沙流模型提出一種新的臨域粒子搜索算法即雙向耦合的全配對搜索法。

圖1中的大球體代表粒子的支持域,位于球心處的粒子稱為目標(biāo)粒子i。對于目標(biāo)粒子i1(橘色大球體球心粒子)、i2(藍(lán)色大球體球心粒子)需要根據(jù)不同的支持域半徑kh,來計算它到整個問題域內(nèi)每一個粒子j(j=1, 2,…, N, N是問題域內(nèi)的粒子總數(shù),圖中用黃色實心小球表示沙粒,藍(lán)色實心小球表示氣體粒子)之間的距離。

圖1 支持域和光滑長度的設(shè)定

本文采用了對稱光滑長度,使粒子i和j互為一組相互作用的粒子對,即粒子i也在粒子j的支持域內(nèi),這樣在搜索問題域內(nèi)的鄰域粒子時不需要再次考慮粒子i對j的影響,故縮減了一般搜索的重復(fù)計算量,提高了搜索效率。支持域半徑為光滑長度h和光滑因子k的乘積,考慮到光滑長度h的選擇與粒子直徑相同,而雙相的不同目標(biāo)粒子對應(yīng)不同的k因子取值,因此光滑因子k的選擇決定了目標(biāo)粒子對鄰域粒子的影響范圍。

分3種情況討論,如圖1所示(黃色小球代表沙粒,藍(lán)色小球代表氣體粒子):1)若目標(biāo)粒子與作用的鄰域粒子同為固體,即沙粒。這兩種固體粒子只有在相互實質(zhì)性的接觸碰撞時方可互相作用,所以光滑因子k3=1.1*h。2)目標(biāo)粒子是氣相粒子,鄰域粒子是沙粒,因為風(fēng)沙流場中氣流的運動狀態(tài)和速度對沙粒的運動影響較大,所以光滑因子取k1=2.5*h。同理,粒子間的作用是相互的,若目標(biāo)粒子是沙粒,鄰域粒子是氣相粒子,則光滑因子為k1=2.5*h。3)當(dāng)目標(biāo)粒子和鄰域粒子同是氣相粒子時,考慮到目標(biāo)粒子的影響半徑超過為3.0h以后對周圍粒子的作用十分微弱,可忽略不記,光滑因子取k2=3.0*h。

2.6 邊界處理方法

SPH數(shù)值模擬中,邊界或鄰近邊界的粒子存在著積分時邊界截斷的固有缺陷,所以SPH方法并不能完全適用于整個模型區(qū)域。且位于邊界或臨近邊界處的粒子在邊界以外不存在粒子,故只受到邊界內(nèi)粒子的影響作用。這種單邊影響作用會導(dǎo)致求解結(jié)果錯誤。為保證積分完整性,分別在計算域上下固定邊界處各設(shè)置一層虛粒子,虛粒子與內(nèi)部氣固流體粒子參與控制方程的求解,并實時更新密度、壓強等物理信息,通過對鄰近實體粒子施加沿軸向的中心強排斥力,來阻止臨近邊界粒子的非物理穿透邊界[13]。

用虛粒子處理邊界的方法提高了SPH近似法在邊界區(qū)域處積分求解的精度。在此基礎(chǔ)上進一步改進虛粒子物理屬性使之符合流固耦合碰撞方程,以此來代替擋風(fēng)墻模型,可模擬、研究流體粒子與擋風(fēng)墻作用的風(fēng)沙流運動規(guī)律。

2.7 時間積分

離散近似后的SPH流體控制方程已經(jīng)轉(zhuǎn)化為常微分方程,可以使用數(shù)值分析中標(biāo)準(zhǔn)的計算方法,如leap-frog格式的時間積分方法,對其進行數(shù)值積分,以此來求解每一個粒子的物理變化。本文采用leap-frog格式積分,其優(yōu)點是計算效率高,計算時所占用的存儲低。粒子的密度、速度、內(nèi)能和位置可由下面的公式遞推。

t=t+Δt

ρi(t+Δt/2)=ρi(t-Δt/2)+Δt.Dρi(t)

vi(t+Δt/2)=vi(t-Δt/2)+Δt.Dvi(t)

ui(t+Δt/2)=ui(t-Δt/2)+Δt.Dui(t)

ri(t+Δt)=ri(t)+Δt.vi(t+Δt/2)

(7)

應(yīng)注意的是,leap-frog格式需要滿足穩(wěn)定性準(zhǔn)則CFL條件,此條件保證時間步長與光滑長度成比例,本文的時間步長取作

Δt=min(ξhi/[hiΔ.vj+cj+1.2(αcj+

βmaxj(vi-vj))])

(8)

式中,ξ為Courant系數(shù),一般取0.3[14]。

2.8 SPH方法的程序?qū)崿F(xiàn)

針對風(fēng)沙顆粒兩相流的耦合特性,在SPH方法基本計算流程的基礎(chǔ)上,基于Linux,C/C++語言開發(fā)環(huán)境,編制了可用于風(fēng)沙防護的擋風(fēng)墻SPH數(shù)值模擬程序,其基本流程如圖2所示。

圖2 SPH計算程序流程圖

3 數(shù)值計算

3.1 計算模型和計算參數(shù)

本文使用離散的方法來模擬擋風(fēng)墻對躍移沙粒運動的影響。為了提高模型的精確度,在算例的計算區(qū)域內(nèi)布置19940個粒子進行仿真,規(guī)模上已達(dá)到統(tǒng)計分析的數(shù)量級,由于大粒子群具有較高的計算時間復(fù)雜度,故仿真算力依托基于Linux編譯環(huán)境的中國科技云·超算平臺。計算域及初始布點如圖3所示,模型寬0.02m,高0.01m;其中沙粒子個數(shù)為1990個,氣體粒子個數(shù)為17950個,擋風(fēng)墻采用與邊界不同物理屬性的虛粒子替代,粒子個數(shù)為60個。根據(jù)SPH方法中設(shè)定的光滑半徑范圍,選用五次樣條函數(shù)的光滑函數(shù)來計算。計算中采用雙向耦合的全配對法搜索鄰域粒子,粒子間光滑長度的選取與粒子直徑相同;時間積分選用“蛙跳法”,步長為10-6s,計算區(qū)域進出口處采用周期邊界;為防止邊界粒子在積分時發(fā)生邊界截斷導(dǎo)致計算崩潰,分別在計算區(qū)域的上邊界和下邊界各設(shè)置一層虛粒子來固定邊界。這些邊壁虛粒子會對靠近邊界的內(nèi)部實粒子產(chǎn)生一個強排斥力,以此來阻止實粒子對邊壁的非物理穿透。初始計算模型如圖3所示,入口風(fēng)速隨高度變化遵循指數(shù)分布規(guī)律,氣體屬性和自然界大氣屬性相同。數(shù)值計算參數(shù)如表1所示。

表1 計算參數(shù)

圖3 初始平坦沙床擋風(fēng)墻模型

4 數(shù)值計算結(jié)果及分析

4.1 擋風(fēng)墻周圍顆粒運動的時空變化

圖4顯示了時間步為(a) step=100,(b) step=6000,(c) step=10000,(d) step=14000,(e) step=25000,(f) step=30000時擋風(fēng)墻周圍躍移沙粒瞬時速度場的發(fā)展過程。從圖中可以看出風(fēng)沙流初始階段的動力學(xué)特性與Zhang等[15]描述的基本一致。

圖4 擋風(fēng)墻周圍風(fēng)沙兩相流運動的時空變化

如圖a,b所示,平坦沙床最外層沙粒最先從氣流中獲得動量起跳,隨著時間推移,大量起跳的沙粒在風(fēng)場中持續(xù)獲能加速,沙粒躍移高度逐漸增大,原始沙床開始松動。由于擋風(fēng)墻對氣流的阻擋作用,使得靠近擋風(fēng)墻周圍的風(fēng)速呈斷崖式遞減,因此擋風(fēng)墻附近的沙粒無法從流場中獲得風(fēng)能起跳。

圖c,d中,躍移的沙粒群正處于成長期,整體呈上升狀態(tài),在接近擋風(fēng)墻時受到風(fēng)場擾動,躍移沙粒出現(xiàn)隨機運動,部分運動的沙粒擊中擋風(fēng)墻后反向反彈,致使同一區(qū)域出現(xiàn)向前和向后的運動狀態(tài)。此時,背風(fēng)區(qū)躍移沙粒在重力作用下出現(xiàn)下降趨勢。

圖e,f迎風(fēng)區(qū)風(fēng)場中躍移的沙粒一部分因為重力作用影響開始下降,靠近擋風(fēng)墻附近的另一小群沙粒由于尾流風(fēng)速迅速減小,沙粒無法從周圍風(fēng)場獲得動能來維持躍移運動也逐漸下落,堆積。而擋風(fēng)墻背風(fēng)區(qū)因為存在一個反向渦[16],在背風(fēng)近墻區(qū)域產(chǎn)生一種向上的升力,上升氣流又會將背風(fēng)側(cè)堆積的沙粒從沙床噴射到空氣中,與傳統(tǒng)躍移沙粒相比,這些粒子呈拋物線運動軌跡,且移動距離更短。由于躍移沙粒與擋風(fēng)墻的碰撞顯著降低顆粒速度,最終穩(wěn)定后的風(fēng)沙流會在擋風(fēng)墻周圍堆積。

4.2 擋風(fēng)墻周圍積沙分布特性

圖5表示風(fēng)沙流穩(wěn)定后迎風(fēng)區(qū)和背風(fēng)區(qū)沙粒數(shù)目的橫向分布。從圖中可以看出擋風(fēng)墻對風(fēng)沙有明顯的阻擋效應(yīng),距離擋風(fēng)墻橫向距離越近,沙粒堆積的數(shù)目越多,以擋風(fēng)墻位置0.01處為中心呈正態(tài)分布特性,中間多,兩邊少。這正與石龍等人的相關(guān)風(fēng)沙兩相流模擬結(jié)果相符[17]。迎風(fēng)區(qū)靠近擋風(fēng)墻處的沙粒因為尾流風(fēng)速減小,獲取的風(fēng)場動能不足以支撐其起跳躍移,實際運動則以蠕移為主。進口處躍移沙粒與擋風(fēng)墻反復(fù)碰撞后逐漸失去動量,在重力作用下降落、沉積于擋風(fēng)墻附近。來流風(fēng)因為擋風(fēng)墻的阻擋作用會在背風(fēng)區(qū)形成一個反向渦,如圖6所示,渦流攜帶的大量沙粒與擋風(fēng)墻發(fā)生碰撞后也會沉積在擋風(fēng)墻周圍。以上現(xiàn)象可以看出擋風(fēng)墻的存在阻礙了沙粒群的運動,大量懸浮的躍移沙粒被攔截,起到了很好的防沙效果。

圖5 擋風(fēng)墻周圍積沙橫向分布

4.3 沙粒平均速度橫向變化

擋風(fēng)墻迎風(fēng)和背風(fēng)區(qū)域躍移沙粒平均速度橫向變化如圖7所示,通過對計算域內(nèi)的1990個沙粒速度矢量進行集合、平均后得到沙粒平均速度,本文主要關(guān)注粒子的合成速度。不難看出,迎風(fēng)區(qū)上游沙粒平均速度整體呈上升趨勢,波動較大,這是因為風(fēng)沙流剛起動時,進口處來流風(fēng)為紊流狀態(tài),各點的流速大小和方向隨時間脈動,沙粒和氣體粒子間發(fā)生強烈的能量交換,使得躍移沙粒的速度和運動軌跡變化得更為復(fù)雜。中游躍移沙粒在重力負(fù)功作用消耗下平均速度大幅降低,下游靠近擋風(fēng)墻的沙粒一部分因為尾流效應(yīng)無法從風(fēng)場獲取動量演變成“蠕移粒子”,速度較低。另一部分懸浮的沙粒與該區(qū)域擋風(fēng)墻發(fā)生碰撞后,落回沙床表面停止運動成為死亡沙粒[18],從而失去動能,最終導(dǎo)致下游擋風(fēng)墻附近沙粒的平均速度急劇下降。與之相反的是背風(fēng)區(qū)靠近擋風(fēng)墻周圍沙粒的平均速度表現(xiàn)出逐漸增大的趨勢,此現(xiàn)象進一步驗證了背風(fēng)區(qū)反向渦的存在,該區(qū)域風(fēng)向與來流風(fēng)向相反,形成對沙粒有抬升作用的回流區(qū)。橫向距離0.0135m之后,風(fēng)場不再受擋風(fēng)墻阻擋效應(yīng)的影響,氣流開始擴散,速度緩慢上升,并最終恢復(fù)至來流風(fēng)初始狀態(tài),此時,沙粒的平均速度逐漸增大,且波動較小,與迎風(fēng)區(qū)上游沙粒平均速度增加趨勢相符。

圖7 迎、背風(fēng)區(qū)躍移沙粒平均速度橫向變化

4.4 擋風(fēng)墻對水平風(fēng)速的影響

風(fēng)沙流是一種貼近地表的風(fēng)沙運動,近地表風(fēng)速會對沙粒運動產(chǎn)生顯著的影響[19],因此分析近地表風(fēng)速變化趨勢對研究擋風(fēng)墻附近積沙形成機理具有指導(dǎo)意義。

圖8顯示了有無擋風(fēng)墻對水平風(fēng)速的影響情況。由圖8可見,擋風(fēng)墻對水平風(fēng)速有明顯的阻擋作用。風(fēng)場中無擋風(fēng)墻時水平風(fēng)速遵循風(fēng)速廓線的指數(shù)分布規(guī)律。增加擋風(fēng)墻后,則遮擋高度以下的風(fēng)速大幅降低,且阻擋效果隨擋風(fēng)墻高度增加而減小,當(dāng)超過擋風(fēng)墻高度時,水平風(fēng)速開始逐漸接近無擋風(fēng)墻時的風(fēng)速廓線。這一現(xiàn)象驗證了用改進虛粒子代替擋風(fēng)墻進行風(fēng)沙流SPH數(shù)值模擬的有效性。

5 結(jié)論

本文基于SPH方法,建立一種由改進虛粒子替代擋風(fēng)墻的平坦沙床模型,并在優(yōu)化鄰域粒子搜索算法的基礎(chǔ)上開發(fā)了可生成SPH粒子的風(fēng)沙流程序。通過上述方法與手段,模擬了擋風(fēng)墻周圍風(fēng)沙流運動的時空變化,可以得出以下幾點結(jié)論:

1)擋風(fēng)墻附近積沙量,以擋風(fēng)墻為中心,呈“正態(tài)性”分布。迎風(fēng)區(qū)由于擋風(fēng)墻的阻擋作用使得風(fēng)場中躍移沙粒與墻體發(fā)生反復(fù)碰撞而失去動量,最終變成“死亡沙?！倍逊e在墻腳迎風(fēng)側(cè)。同時下游因為尾流效應(yīng)的存在,靠近擋風(fēng)墻的部分沙粒無法從風(fēng)場中獲取動量而演變成蠕移粒子。

2)背風(fēng)區(qū)因為反向渦的影響,在擋風(fēng)墻附近產(chǎn)生向上的升力,上升氣流攜帶的大量沙粒與墻體再次發(fā)生碰撞導(dǎo)致動量損失,最終堆積在擋風(fēng)墻背風(fēng)側(cè)。

3)風(fēng)場無擋風(fēng)墻時,水平風(fēng)速遵循風(fēng)速廓線的指數(shù)分布規(guī)律。有擋風(fēng)墻時,對風(fēng)速的阻擋效果隨擋風(fēng)墻高度增加而減小,超過擋風(fēng)墻垂向高度,水平風(fēng)速逐漸恢復(fù)至無擋風(fēng)墻時的凈風(fēng)場狀態(tài)。

數(shù)值模擬結(jié)果表明,擋風(fēng)墻對躍移沙粒和來流風(fēng)具有明顯的阻擋效應(yīng)。本文提出的基于SPH方法顆粒兩相流仿真系統(tǒng)對防沙、治沙工程具有一定的借鑒意義,并能反映出沙障附近積沙的普遍規(guī)律。