孫 浩 劉晉浩 黃青青

(北京林業(yè)大學(xué)工學(xué)院，北京100083)

沙障地表形態(tài)衍化數(shù)值模擬方法研究

孫 浩 劉晉浩 黃青青

(北京林業(yè)大學(xué)工學(xué)院，北京100083)

提出了一個風(fēng)沙相互作用數(shù)值模型，該模型以概率統(tǒng)計理論為基礎(chǔ)，通過計算不同地表風(fēng)速下沙粒的風(fēng)蝕與沉積概率來實(shí)現(xiàn)風(fēng)-沙相互作用過程，最后通過動網(wǎng)格技術(shù)，改變下邊界節(jié)點(diǎn)坐標(biāo)，實(shí)現(xiàn)沙床表面的起伏變化。數(shù)值計算與野外沙障凹坑的測量結(jié)果對比表明，沙障內(nèi)地表形態(tài)的數(shù)值計算結(jié)果與實(shí)測地表形態(tài)較為接近，凹坑最低點(diǎn)位置向下風(fēng)向一側(cè)移動，誤差最大位置為凹坑的最低點(diǎn)處，t=4 d時，平均絕對誤差為17.86%，隨著積沙增加，誤差逐漸減小，t=8 d時，平均絕對誤差為8.52%。所提模型無論從定性還是從定量上，都與野外測量結(jié)果較為一致，可以正確模擬沙障地表的發(fā)展過程。

沙障;地表過程;動網(wǎng)格;數(shù)值方法

引言

沙障在防風(fēng)固沙工程中被大量應(yīng)用并占有重要地位。一般認(rèn)為，沙障可以有效降低近地表風(fēng)速，抑制沙粒啟動，削弱地表風(fēng)沙流強(qiáng)度。科研人員對不同類型沙障的防風(fēng)效能［1－4］、沙障周圍流場特征［5－8］進(jìn)行了大量的研究，而對沙障的衰退過程，目前研究還相對較少。一般認(rèn)為，沙障方格內(nèi)積沙是導(dǎo)致沙障失去固沙能力的主要原因。沙障的衰退過程即為沙障逐漸被沙掩埋的過程。沙障在鋪設(shè)完成后，首先在沙障方格前沿開始產(chǎn)生積沙，并逐漸向縱深發(fā)展，上風(fēng)側(cè)的方格積沙量相對較高，下風(fēng)向沙障內(nèi)積沙較少，隨著沙的淤積，沙障露出地表的高度逐漸降低，最終失去固沙能力，同時積沙飽和帶也逐漸前移，當(dāng)沙障內(nèi)達(dá)到蝕積動態(tài)平衡時，沙粒不再沉積［9］。一種情況為沙障完全被沙掩埋，固沙帶變成近乎于平坦地表;另一種情況為沙障形成穩(wěn)定凹曲面，由于沙障后的渦流氣動作用達(dá)到蝕積平衡，此時固沙帶完全喪失阻沙能力。許多國內(nèi)學(xué)者對不同尺寸的沙障內(nèi)凹曲面進(jìn)行研究，結(jié)果表明，1 m×1 m方格沙障容易形成穩(wěn)定凹曲面，以凹坑深度與方格尺寸比值表示蝕積系數(shù)，在1/10～1/8之間［10－13］。穩(wěn)定凹曲面的形成條件與風(fēng)力、沙障材料及地表坡度密切相關(guān)，穩(wěn)定凹曲面的蝕積系數(shù)也會發(fā)生相應(yīng)變化［14］。

沙障凹曲面的發(fā)展是一個動態(tài)過程，目前對沙障的凹曲面發(fā)展速度、沙障衰退時間方面仍然研究不足，主要是野外環(huán)境較為復(fù)雜，同時需要長時間的監(jiān)測工作。在風(fēng)沙相互作用的研究中，數(shù)值計算是一個重要的研究方法，大量的研究表明，數(shù)值計算可以很好地描述風(fēng)沙運(yùn)動的瞬態(tài)形為，與實(shí)驗結(jié)果較為接近，但當(dāng)前的數(shù)值計算方法只能模擬短時間的流固兩相運(yùn)動過程［15－16］，對于沙障凹坑的形成過程這種大時間尺度問題的模擬，無論是歐拉-歐拉［17－18］還是歐拉-拉格朗日［19－22］雙流體模型都無法滿足要求。

為了對沙障地表發(fā)展過程這種大時間尺度的風(fēng)沙運(yùn)動進(jìn)行模擬，本文基于概率理論對風(fēng)沙相互作用數(shù)值模型進(jìn)行研究，為沙障衰退過程及障內(nèi)凹曲面形成預(yù)測提供方法，從而為工程治沙及沙障尺寸設(shè)計提供理論支撐。

1 數(shù)值方法

本文數(shù)值模型計算結(jié)構(gòu)與傳統(tǒng)兩相流數(shù)值計算結(jié)構(gòu)相同，由流體計算和沙體計算兩部分構(gòu)成，流體計算采用傳統(tǒng)計算流體力學(xué)方法，沙體計算由本研究提出的模型進(jìn)行計算。

1.1 沙相網(wǎng)格劃分

在現(xiàn)有的兩相流模擬中，氣固兩相之間的相互作用行為表示為氣固兩相之間的動量交換行為，在網(wǎng)格數(shù)量與顆粒數(shù)量都較多時，需要大量計算資源。而本研究為了實(shí)現(xiàn)大空間尺度的模擬，在模擬沙相時提出采用雙重網(wǎng)格思想，原理如圖1所示。

圖1 沙相網(wǎng)格劃分原理Fig．1 Principle of sand phasemeshing

將整個沙障地表看作二維平面，劃分網(wǎng)格時首先按草方格尺寸劃分網(wǎng)格，作為母網(wǎng)格，網(wǎng)格單元數(shù)量與沙障網(wǎng)格數(shù)量相等。再對母網(wǎng)格進(jìn)行網(wǎng)格劃分，生成子網(wǎng)格。首先將沙床表面進(jìn)行粗糙網(wǎng)格劃分，計算粗糙網(wǎng)格總沙量變化，然后對粗糙網(wǎng)格再次細(xì)化，計算母網(wǎng)格單元的積沙量在子網(wǎng)格單元的分布情況，最終計算子網(wǎng)格單元厚度，實(shí)現(xiàn)地表動態(tài)變化。

1.2 風(fēng)沙兩相耦合計算過程

風(fēng)沙兩相流計算主要包含氣相計算和沙相計算，具體步驟為:①施加邊界條件。對計算模型的氣相和沙相的計算域施加邊界條件，氣相計算需要設(shè)置壁面類型、速度入口、時間步長等參數(shù)，沙相計算需要對沙粒物理參數(shù)、地表輸沙量、時間步長及沙相模型相關(guān)參數(shù)等進(jìn)行初始化。②風(fēng)場計算。應(yīng)用傳統(tǒng)的流體計算理論對空氣相進(jìn)行計算，并獲得下邊界壁面剪力數(shù)據(jù)。③沙相計算。通過壁面剪力計算沙床表面每個子網(wǎng)格單元輸沙率，進(jìn)而計算整個草沙障輸沙率，計算沙床沉降量分布和侵蝕量分布。④基于沙障間沙粒沉降概率模型計算每個子網(wǎng)格單元的沙粒沉降量，進(jìn)而確定子單元網(wǎng)格厚度。⑤循環(huán)檢測相鄰子網(wǎng)格單元厚度，超過沙體坍塌臨界角度時，子單元發(fā)生坍塌行為，修改子單元網(wǎng)格厚度。⑥采用動網(wǎng)格技術(shù)根據(jù)子單元厚度，調(diào)整下邊界網(wǎng)格節(jié)點(diǎn)坐標(biāo)，實(shí)現(xiàn)沙床表面動態(tài)變化，然后根據(jù)邊界網(wǎng)格坐標(biāo)變化，調(diào)整內(nèi)部網(wǎng)格節(jié)點(diǎn)坐標(biāo)。⑦判斷程序終止時間，如果到達(dá)終止時間程序結(jié)束，否則返回。計算流程如圖2所示。

1.3 風(fēng)沙相互作用數(shù)值模型

將每個草方格內(nèi)部看做一個母單元，草方格內(nèi)部積沙和風(fēng)蝕總量確定之后，計算沙障內(nèi)部具體的風(fēng)沙蝕積行為。草方格內(nèi)風(fēng)沙相互作用原理如圖3所示。圖中Qin為進(jìn)入方格上空的總輸沙率;a為進(jìn)入方格區(qū)域內(nèi)的來流風(fēng)沙的進(jìn)入概率;aQin為來流風(fēng)沙影響方格內(nèi)部沙床蝕積變化的輸沙率;Qerosion/n為方格內(nèi)子單元由于風(fēng)蝕產(chǎn)生的平均起沙輸沙率。

圖2 程序流程圖Fig．2 Flow chart of program

圖3 草方格風(fēng)沙作用原理圖Fig．3 Schematic of interaction between wind and sand in straw checkerboard

考慮進(jìn)入任意沙障方格單元內(nèi)的風(fēng)沙流微元的運(yùn)動過程，從風(fēng)沙流微元進(jìn)入沙障方格范圍后，與地表發(fā)生沙粒交換，風(fēng)沙流率產(chǎn)生變化，一直到流出沙障范圍這一過程。由于沙床表面的地表是否發(fā)生侵蝕還取決于此位置上方的輸沙率是否達(dá)到飽和，因此需要對比地表當(dāng)前位置由風(fēng)蝕產(chǎn)生的飽和輸沙率與地表上方此時的輸沙率，本文應(yīng)用參數(shù)c對此差值進(jìn)行修正，定義為風(fēng)沙流與地表相互作用產(chǎn)生的沙粒交換速率。c(Qerosion/n－aQin)即為地表的沙量變化速率，地表上方的風(fēng)沙流與地表沙床相互作用后，沙障地表內(nèi)輸沙率將發(fā)生變化，Qair為變化后的沙床表面輸沙率，Qair=c(Qerosion/n－aQin)+aQin。沙障內(nèi)地表附近的風(fēng)沙流微元會因為繼續(xù)運(yùn)動一部分流出沙障上空范圍，與一部分的入口風(fēng)沙流共同組成方格的出口輸沙率，即Qout，一部分則會發(fā)生沉降，沉積在沙床表面。此模型中需要確定3個參數(shù)，入口風(fēng)沙流參與地表相互作用的這部分風(fēng)沙流的比例參數(shù)a，定義為風(fēng)沙流進(jìn)入概率;沙障方格內(nèi)風(fēng)沙流飛出沙障的概率b，定義為逃逸概率;沙粒交換速率c。

1.3.1 進(jìn)入概率a數(shù)學(xué)模型

草方格入口風(fēng)沙流進(jìn)入沙障地表附近直接參與地表沙粒交換的概率，受沙障內(nèi)部沙床表面高度的影響，本文假定a只與沙床表面高度有關(guān)，因此a是變量為高度h的函數(shù)。本文假定進(jìn)入概率a隨沙障內(nèi)部平均高度變化服從指數(shù)函數(shù)分布，當(dāng)沙障內(nèi)部地表高度與草頭高度相等時，進(jìn)入沙障范圍內(nèi)的風(fēng)沙流將全部對地表沙粒運(yùn)動產(chǎn)生影響，此時進(jìn)入概率為1。進(jìn)入概率a的表達(dá)式為

式中 t——時間，s m——高度調(diào)整系數(shù)

a(t)——t時刻的風(fēng)沙流進(jìn)入概率

a0——初始風(fēng)沙流進(jìn)入概率

Ht——t時刻沙障內(nèi)部地表平均高度，m

Hbarrier——沙障草頭高度，m

1.3.2 逃逸概率b數(shù)學(xué)模型

沙障內(nèi)風(fēng)沙流流出沙障的概率與沙障內(nèi)部地表與空中沙粒交換速率、沙障內(nèi)地表平均高度、沙障內(nèi)風(fēng)沙流速度有關(guān)，與進(jìn)入概率a的數(shù)學(xué)模型相同，本文定義逃逸概率b的數(shù)學(xué)模型為指數(shù)函數(shù)模型。b的初值b0由初始化的沙障地表整體輸沙率和沙粒交換速率c確定得出，c值假定僅與沙粒物理參數(shù)有關(guān)，是一個常數(shù)。初始時刻t=0時，沙障方格的出口風(fēng)沙流率和入口風(fēng)沙流率由沙相邊界條件進(jìn)行初始化。逃逸概率b的表達(dá)式為

式中 b(t)——t時刻的風(fēng)沙流逃逸概率

b0——t=0時的逃逸概率

n——沙障內(nèi)子單元數(shù)量

1.4 草方格障間風(fēng)蝕量計算

當(dāng)沙床表面剪切風(fēng)速超過臨界起沙風(fēng)速時，沙粒啟動形成風(fēng)沙流，并在1～2 s內(nèi)逐漸保持穩(wěn)定。有關(guān)剪切風(fēng)速和穩(wěn)定輸沙率的數(shù)學(xué)模型，前人已經(jīng)對其進(jìn)行了大量的研究，主要分為2類，一種為O’brien-Rindlaub型方程，將地表以上一定高度的風(fēng)速作為輸沙率函數(shù)變量。另一種為Bagnonld型方程，將地表剪切風(fēng)速作為輸沙率函數(shù)變量。由于沙障地表并非平坦地面，沙障上方風(fēng)速廓線不完全服從對數(shù)分布，因此，對于鋪設(shè)沙障的地表，以地表剪切風(fēng)速作為計算輸沙率的依據(jù)是合理的。本文使用的輸沙率計算公式為

式中 Qelement——子單元網(wǎng)格的輸沙率，kg/(m·s)

Rt——沙粒啟動時臨界剪切風(fēng)速和當(dāng)前剪切風(fēng)速的比值

ρs——沙粒密度，kg/m3

g——重力加速度，m/s2

u*——剪切風(fēng)速，m/s

d——沙粒平均直徑，mm

D——參考沙粒直徑，取0.25mm

草方格內(nèi)總風(fēng)蝕率為

式中 i——草方格內(nèi)子網(wǎng)格單元編號

Qelement(i)——單元i的風(fēng)蝕率

1.5 子網(wǎng)格單元厚度計算

沙床表面的沙粒沉降分為兩部分，一是由于沙床表面的風(fēng)沙流和地表上沙粒相互作用產(chǎn)生的地表沙量變化，二是由于沙障內(nèi)部風(fēng)沙流在運(yùn)動過程中由于沙障的阻礙和地表起伏而引起的沙粒沉降。計算網(wǎng)格厚度變化時，首先要計算沙障方格內(nèi)部各子單元網(wǎng)格的沙粒變化量。

1.5.1 風(fēng)沙蝕積作用產(chǎn)生的地表沙量變化

在氣流的作用下，沙障地表會有不同程度的風(fēng)沙流形成，而當(dāng)前位置是否發(fā)生風(fēng)蝕或者堆積取決于當(dāng)前地表的摩阻風(fēng)速和風(fēng)沙流的供給量，即當(dāng)前位置上方的風(fēng)沙流濃度。沙床表面可以產(chǎn)生的最強(qiáng)風(fēng)沙流與當(dāng)前位置的摩阻風(fēng)速有關(guān)，本文將在一定摩阻風(fēng)速下產(chǎn)生的最大風(fēng)沙流率定義為當(dāng)前位置的風(fēng)沙流承載力。當(dāng)沙源供給量大于承載能力時，沙粒在此位置發(fā)生沉降，當(dāng)沙源供給量小于承載力時，此位置發(fā)生風(fēng)蝕。

沙床表面的輸沙率分布取決于沙障內(nèi)摩阻風(fēng)速分布，顯然每個單元的沙源供給量不是均勻分布，為了對模型進(jìn)行簡化，本文假設(shè)每個子網(wǎng)格單元的沙源供給量相同，如圖4所示，即Qsupply=aQin。每個子單元的沙量承載力可以由式(3)計算得出，定義為 Qcapacity。

圖4 沙床子單元厚度計算原理圖Fig．4 Schematic of thickness calculation of sand element

圖4中H為子單元厚度，L為子單元順風(fēng)向長度，W為子單元垂直風(fēng)向?qū)挾?。研究表明，沙床表面從靜止?fàn)顟B(tài)到形成飽和風(fēng)沙流狀態(tài)需要1～2 s的時間，故由于沙源的供給量和承載力不均衡導(dǎo)致的沙量變化也需要一定的時間，本模型中應(yīng)用沙粒交換速率c對此結(jié)果進(jìn)行修正。子單元網(wǎng)格沙粒變化率計算式為

式中 Qdiff(i)——單元i的沙粒變化率，

當(dāng)Qdiff(i)＞0時，風(fēng)沙流中的沙粒發(fā)生沉降，當(dāng)Qdiff(i)＜0時，子單元發(fā)生風(fēng)蝕。

子單元高度計算式為

式中 Hi——單元i厚度，m

ΔTs——沙相計算時間步長，s

1.5.2 沙障阻擋及地表起伏產(chǎn)生的沙量沉降

風(fēng)沙流在進(jìn)入沙障范圍后，會受到沙障的阻擋而發(fā)生明顯的堆積，而地表的起伏變化同樣會引起沙粒的不均勻沉降，研究表明對于沙丘地貌，風(fēng)沙流在到達(dá)沙丘位置時，沙粒容易在沙丘坡腳處發(fā)生堆積，一方面由于沙丘坡腳處風(fēng)速變?nèi)?，沙粒的承載能力變?nèi)?，另一方面也是由于地表起伏發(fā)生變化，沙粒從水平運(yùn)動轉(zhuǎn)變?yōu)檠赜L(fēng)坡向上運(yùn)動，需要克服重力向上爬升。本模型中將沙障內(nèi)部地表風(fēng)沙流的不均勻沉降表達(dá)為概率模型，通過計算沙障內(nèi)部不同子網(wǎng)格的沉降概率，來計算子網(wǎng)格的沉降量。結(jié)合以往的研究，本文假定沙粒沉降概率與地表坡度變化率有關(guān)。坡度變化率越高，沙粒沉降概率越高，地表平坦時，沙粒沉降概率為0。本文模型將沙障內(nèi)的地表分為3種子地形，如圖5所示，可簡化為坡中、坡頂和坡谷。

圖5 沙障內(nèi)地表小地形Fig．5 Sub landform in sand barrier checkerboard

沙障內(nèi)子單元的坡度角變化量計算式為

式中 θdiff(i)——第i個子單元的坡度角變化量

單元i判斷為坡頂?shù)匦螘r，θdiff(i)=0。沙障方格內(nèi)子單元的概率計算式為

式中 Psettlement(i)——第i個子單元的沉積概率

子單元沉降量和沙障內(nèi)輸沙率更新示意圖如圖6所示。

圖6 沙粒沉降及沙障輸沙率計算Fig．6 Calculation of sand settling and sand transport rate above sand barrier

圖6中，(1－b)QairPsettlement(i)為子單元網(wǎng)格的沉降量，子單元厚度改變量計算與式(6)相同。沙障方格的出口輸沙率更新為Qout，并作為下一個母網(wǎng)格的入口輸沙率，計算式為

2 數(shù)值計算

本文以開源流體力學(xué)計算軟件OpenFOAM為基礎(chǔ)，編寫沙相數(shù)值模型，并與流場耦合進(jìn)行數(shù)值計算。

2.1 幾何模型及網(wǎng)格劃分

本文風(fēng)沙兩相耦合計算模型采用二維幾何模型，模擬沙障在風(fēng)洞環(huán)境下的積沙行為，模型計算域如圖7所示。

圖7 沙障計算域示意圖Fig．7 Schematic diagram of computational setup for sand fence

計算域由空氣相計算域和沙障計算域構(gòu)成，計算域高度為1 m，沙障高度為0.1 m，沙障間距為1m，入口邊界到第1個沙障的距離為1 m，計算域右側(cè)沙障到計算域出口的距離為8 m，為了后面便于動網(wǎng)格調(diào)整，避免出現(xiàn)低質(zhì)量網(wǎng)格而增加計算誤差，模型對沙障部分采用孔隙域模型，設(shè)置其孔隙為零。

流場網(wǎng)格尺寸為25 mm×25 mm。在使用標(biāo)準(zhǔn)k-ε湍流時，一般要求y+在30～200之間，以保證近壁面流場計算準(zhǔn)確，壁面第 1層網(wǎng)格高度為1mm，壁面網(wǎng)格設(shè)置10層。網(wǎng)格劃分如圖8所示。

圖8 計算域網(wǎng)格劃分Fig．8 Mesh grid of computational domain

2.2 邊界條件及模型參數(shù)初始化

流場左邊界采用速度入口邊界條件，入口風(fēng)速沿高度服從對數(shù)函數(shù)分布，摩阻風(fēng)速為0.6m/s。右邊界采用自由流出邊界，上邊界采用光滑壁面邊界，下邊界采用無滑移壁面邊界，沙障計算域設(shè)置為孔隙域模型。主要參數(shù)如表1所示。

表1 模型參數(shù)設(shè)置Tab．1 Model parameters

3 數(shù)值計算結(jié)果

對上述模型進(jìn)行計算，不同時間的流場速度分布計算結(jié)果如圖9所示，時間t為沙相的計算時間。時間t=0 d時，流場已經(jīng)充分發(fā)展，然后開始耦合計算。從圖中可以看出，t=0.5 d時，沙床表面沙障兩側(cè)最先開始積沙。隨著時間的發(fā)展，沙障兩側(cè)積沙逐漸提高，沙坑表面逐漸發(fā)展為圓弧形態(tài)，由于沙坑內(nèi)各位置的地表風(fēng)速分布不均，圓弧最低點(diǎn)逐漸向下風(fēng)向移動。從以往的研究和實(shí)地調(diào)查中可以看出，從定性分析上看，本模型的沙坑發(fā)展與野外觀測結(jié)果較為吻合。

對野外草沙障凹坑沿主風(fēng)向中線進(jìn)行測量，并與數(shù)值計算結(jié)果進(jìn)行對比，結(jié)果如圖10所示。由圖10可以看出，時間為4 d時仿真結(jié)果與本文的測量結(jié)果較為一致。時間為8 d時，仿真結(jié)果與張登山等［14］的實(shí)地測量結(jié)果接近。當(dāng)仿真時間為4 d時，仿真結(jié)果與實(shí)地測量結(jié)果相比，相對誤差最大的位置為凹曲面的最低點(diǎn)處，仿真結(jié)果比實(shí)驗結(jié)果小30.80%，沙障內(nèi)部凹曲面的平均絕對誤差為17.86%。當(dāng)仿真時間為8 d時，仿真結(jié)果與實(shí)地測量結(jié)果相比，差距最大的位置仍為凹曲面的最低點(diǎn)處，仿真結(jié)果比實(shí)驗結(jié)果小16.94%，沙障內(nèi)部凹曲面的平均絕對誤差為8.52%。仿真結(jié)果小于實(shí)驗結(jié)果的原因為:沙障不具有透過性時，障后渦流強(qiáng)度更高，沙床表面的風(fēng)蝕量大于實(shí)際工況。而實(shí)際情況下，沙障存在一定的孔隙，加之草頭對氣流的擾動，導(dǎo)致沙障后渦流強(qiáng)度小于模擬工況，即便如此，仍可以從定量與定性上看出，模擬結(jié)果與實(shí)驗結(jié)果吻合較好。

圖9 流場速度分布Fig．9 Velocity distributions of wind flow field

圖10 仿真結(jié)果與實(shí)驗結(jié)果對比Fig．10 Comparison of simulated results and field measurement results

4 討論

沙障地表的衍化過程是沙障地表長期的風(fēng)沙蝕積過程，沙障在風(fēng)沙作用下，障內(nèi)地表形態(tài)發(fā)生改變，而障內(nèi)地表形態(tài)的改變又反過來影響流場結(jié)構(gòu)的變化。雖然實(shí)際問題中，沙體顆粒數(shù)量巨大，對于沙障地表形態(tài)改變這一大空間與時間尺度的問題，現(xiàn)有的模型無法滿足要求，但大量的沙粒運(yùn)動，在風(fēng)力作用下卻表現(xiàn)出了一定的概率統(tǒng)計規(guī)律，如地表的風(fēng)沙流結(jié)構(gòu)與剪切速度的關(guān)系、沙粒長期的沉降概率等。應(yīng)用概率統(tǒng)計模型雖然對于模擬風(fēng)沙瞬時運(yùn)動不夠精確，但在模擬長期的風(fēng)沙耦合過程方面卻表現(xiàn)出了特有的優(yōu)勢，不但避免了沙粒數(shù)量帶來的大計算量，而且避免了顆粒計算的小時間步長問題，使得計算速度加快。

本文的計算模型將沙障看作為非透過性薄板墻體，與實(shí)際的沙障有一定差距，導(dǎo)致在計算初期，地表衍化過程與實(shí)際情況有一定差距，凹坑最低點(diǎn)處的仿真結(jié)果比實(shí)驗結(jié)果小30.80%，但隨著積沙過程的持續(xù)，沙障墻體兩側(cè)被沙掩埋，墻體的孔隙作用對流場的影響減弱，使得沙障內(nèi)積沙一定程度后，模擬工況和實(shí)際工況差距縮小，沙障內(nèi)凹坑發(fā)展逐漸接近實(shí)際情況;另一方面，在野外環(huán)境下，沙障凹坑的發(fā)展是在多個風(fēng)向下共同作用形成的，而模擬工況更接近于風(fēng)洞環(huán)境，使得沙障地表的數(shù)值計算與野外地表發(fā)展過程有一定偏差。本文選擇主風(fēng)向明顯的沙障地表進(jìn)行沙障凹坑的測量，盡量縮小模擬工況和實(shí)際工況的差距，在數(shù)值計算結(jié)果與實(shí)際測量結(jié)果的對比中，數(shù)值模型計算結(jié)果雖然只代表實(shí)際工況下地表發(fā)展中經(jīng)歷的一個地表形態(tài)，不能精確地確定地表實(shí)際情況下發(fā)展到此形態(tài)所需要的時間，但通過數(shù)值模型，可以對比不同尺寸沙障在同一環(huán)境下地表的發(fā)展過程、衰退以及沙障可以保持固沙作用的時間，可以為沙障的尺寸設(shè)計提供參考。

5 結(jié)束語

本文提出的風(fēng)沙相互作用模型基于概率統(tǒng)計理論，通過動網(wǎng)格技術(shù)調(diào)整計算域下邊界，實(shí)現(xiàn)沙障地表的形態(tài)變化。數(shù)值計算表明，本文模型計算結(jié)果與實(shí)際沙障的地表形態(tài)發(fā)展過程較為吻合，凹坑發(fā)展過程中，初始積沙時，凹坑內(nèi)積沙速度較快，而障間形成凹曲面后，地表高度增長變緩，沙障凹坑最低點(diǎn)向下風(fēng)向移動。從沙障曲面的定量分析可以看出，沙障地表初期發(fā)展與實(shí)際環(huán)境下地表變化差距較大，相對誤差最大為30.80%，平均絕對誤差為17.86%。當(dāng)沙障內(nèi)部存在一定積沙后，數(shù)值計算與實(shí)際測量結(jié)果差距縮小，仿真結(jié)果與實(shí)驗結(jié)果相比相對誤差最大為16.94%，沙障內(nèi)部凹曲面的平均絕對誤差為8.52%，本文數(shù)值模型可以正確模擬野外沙障的凹坑發(fā)展過程。

參 考 文 獻(xiàn)

1 龐營軍，屈建軍，謝勝波，等．高立式格狀沙障防風(fēng)效益［J］．水土保持通報，2014，34(5):11－14．PANG Yingjun，QU Jianjun，XIE Shengbo，etal．Windproof efficiency ofupright checkerboard sand-barriers［J］．Bulletion of Soil and Water Conservation，2014，34(5):11－14．(in Chinese)

2 屈建軍，凌裕泉，俎瑞平，等．半隱蔽格狀沙障的綜合防護(hù)效益觀測研究［J］．中國沙漠，2005，25(3):329－335．QU Jianjun，LING Yuquan，ZU Ruiping，et al．Study on comprehensive sand-protecting efficiency of semi-buried checkerboard sand-barriers［J］．Journal of Desert Research，2005，25(3):329－335．(in Chinese)

3 黃強(qiáng)，雷加強(qiáng)，王雪芹．塔里木沙漠公路不同地貌部位的高立式沙障阻沙特征［J］．干旱區(qū)地理，2000，23(3):227－232．HUANG Qiang，LEIJiaqiang，WANG Xueqin．Sand-obstructing effectof the high sandbarriers in the differentmorhpologic sections along the traim desert highway［J］．Arid Land Geography，2000，23(3):227－232．(in Chinese)

4 黨曉宏，高永，虞毅，等．新型生物可降解PLA沙障與傳統(tǒng)草方格沙障防風(fēng)效益［J］．北京林業(yè)大學(xué)學(xué)報，2015，37(3): 118－125．DANG Xiaohong，GAO Yong，YU Yi，et al．Windproof efficiency with new biodegradable PLA sand barrier and traditional straw sand barrier［J］．Journal of Beijing Forestry University，2015，37(3):118－125．(in Chinese)

5 DONG Z B，LUOW Y，QIAN G Q，et al．A wind tunnel simulation of themean velocity fields behind upright porous fences［J］．Agricultural and ForestMeteorology，2007，146(1):82－93．

6 CORNELISW M，GABRIELSD．Optimal windbreak design for wind-erosion control［J］．Journal of Arid Environments，2005，61(2):315－332．

7 ZHANG N，KANG J，LEE S．Wind tunnel observation on the effect of a porous wind fence on shelter of saltating sand particles［J］．Geomorphology，2010，120(3):224－232．

8 DONG Z B，QIAN G Q，LUOW Y，et al．Threshold velocity for wind erosion:the effects of porous fences［J］．Environmental Geology，2006，51(3):471－475．

9 姚正毅，陳廣庭，韓志文，等．機(jī)械防沙體系防沙功能的衰退過程［J］．中國沙漠，2006，26(2):226－231．YAO Zhengyi，CHEN Guangting，HAN Zhiwen，et al．Declinemechanism and process ofmechanical defense system［J］．Journal of Desert Research，2006，26(2):226－231．(in Chinese)

10 王振亭，鄭曉靜．草方格沙障尺寸分析的簡單模型［J］．中國沙漠，2002，22(3):229－232．WANG Zhenting，ZHENG Xiaojing．A simplemodel for calculatingmeasurements of straw checkerboard barriers［J］．Journal of Desert Research，2002，22(3):229－232．(in Chinese)

11 常兆豐，仲生年，韓福桂，等．粘土沙障及麥草沙障合理間距的調(diào)查研究［J］．中國沙漠，2000，20(4):455－457．CHANG Zhaofeng，ZHONG Shengnian，HAN Fugui，et al．Research of the suitable row spacing on clay barriers and straw barriers［J］．Journal of Desert Research，2000，20(4):455－457．(in Chinese)

12 馬學(xué)喜，李生宇，王海峰，等．固沙網(wǎng)沙障積沙凹曲面特征及其固沙效益分析［J］．干旱區(qū)研究，2016，33(4):898－904．MA Xuexi，LIShengyu，WANG Haifeng，et al．Concave surface features in sand-fixing net barriers and evaluation of sand-fixing benefits［J］．Arid Zone Research，2016，33(4):898－904．(in Chinese)

13 周丹丹，虞毅，胡生榮，等．沙袋沙障凹曲面特性研究［J］．水土保持通報，2009，29(4):22－25．ZHOU Dandan，YU Yi，HU Shengrong，et al．Concave surface characteristics of sand bag sand barrier［J］．Bulletin of Soil and Water Conservation，2009，29(4):22－25．(in Chinese)

14 張登山，吳汪洋，田慧麗，等．青海湖沙地麥草方格沙障的蝕積效應(yīng)與規(guī)格選取［J］．地理科學(xué)，2014，34(5):627－634．ZHANG Dengshan，WU Wangyang，TIAN Huili，et al．effects of erosion and deposition and dimensions selection of strawcheckerboard barriers in the desert of qinghai lake［J］．Scientia Geographica Sinica，2014，34(5):627－634．(in Chinese)

15 喻黎明，鄒小艷，譚弘，等．基于CFD-DEM耦合的水力旋流器水沙運(yùn)動三維數(shù)值模擬［J/OL］．農(nóng)業(yè)機(jī)械學(xué)報，2016，47(1):126－132．http:∥www．j-csam．org/jcsam/ch/reader/view_abstract．a(chǎn)spx?flag=1＆file_no=20160117＆journal_id= jcsam．DOI:10．6041/j．issn．1000-1298．2016．01．017．YU Liming，ZOU Xiaoyan，TAN Hong，et al．3D numerical simulation of water and sediment flow in hydrocyclone based on coupled CFD-DEM［J/OL］．Transactions of the Chinese Society for Agricultural Machinery，2016，47(1):126－132．(in Chinese)

16 喻黎明，譚弘，鄒小艷，等．基于CFD-DEM耦合的迷宮流道水沙運(yùn)動數(shù)值模擬［J/OL］．農(nóng)業(yè)機(jī)械學(xué)報，2016，47(8): 65－71．http:∥www．j-csam．org/jcsam/ch/reader/view_abstract．a(chǎn)spx?flag=1＆file_no=20160810＆journal_id=jcsam．DOI: 10．6041/j．issn．1000-1298．2016．08．010．YU Liming，TAN Hong，ZOU Xiaoyan，et al．Numerical simulation of water and sediment flow in labyrinth channel based on coupled CFD-DEM［J/OL］．Transactions of the Chinese Society for Agricultural Machinery，2016，47(8):65－71．(in Chinese)

17 CHEN G H，WANG W W，SUN C F，et al．3D numerical simulation of wind flow behind a new porous fence［J］．Powder Technology，2012，230:118－126．

18 LOPESA M G，OLIVEIRA L A，ALMERINDO D，et al．Numerical simulation of sand dune erosion［J］．Environmental Fluid Mechanics，2013，13(2):145－168．

19 TONG D，HUANG N．Numerical simulation of saltating particles in atmospheric boundary layer over flat bed and sand ripples［J］．Journal of Geophysical Research:Atmospheres，2012，117(D16):16205．

20 HUANG N，XIA X P，TONG D．Numerical simulation ofwind sand movement in straw checkerboard barriers［J］．The European Physical Journal，2013，36:1－7．

21 SCHMEECKLEM W．Numerical simulation of turbulence and sediment transport of medium sand［J］．Journal of Geophysical Research Earth Surface，2014，119(6):1240－1262．

22 DURAN O，CLAUDIN P，ANDREOTTIB．Direct numerical simulations of aeolian sand ripples［J］．Proceedings of the National Academy of Sciences，2014，111(44):15665－15668．

Numerical Method of Surface Morphology Evolution of Sand Barrier

SUN Hao LIU Jinhao HUANG Qingqing
(School of Technology，Beijing Forestry University，Beijing 100083，China)

A new numerical model of sand-wind interaction was proposed to predict the surface morphology evolution and provide the basis for dimension design of sand barrier．Based on the theory of probability and statistics，the model realized sand-wind interaction process by calculating wind erosion and deposition probability of sand grains at different wind speeds．Firstly，the model parameters of air phase and sand phase were initialized，and the velocity distribution can be calculated by computational fluid dynamic(CFD)method．Secondly，the wind friction velocity of sand element was calculated by using velocity field data，and the sand transport rate can be obtained．Finally，the thickness of sand elementwas calculated by calculating sand change in one time step．The moving mesh technique was used to change the fluctuation of the sand bed by changing the coordinates of lower boundary nodes．The numerical resultswere compared with the measured results，which showed that the numerical results of the surfacemorphology in the sand barrierwere close to themeasured results．The lowest point of pitwas themaximum position thatmoved along the downwind area，the average absolute errorwas17.86%when t=4 d，the error was decreased gradually with the increase of sand accumulation，the average absolute error was8.52%when t=8 d．Both qualitative and quantitative resultswere consistentwith themeasured results，and themodel can simulate the surfacemorphology evolution of sand barrier correctly．

sand barrier;surface process;dynamic mesh;numericalmethod

S727.23;S775

A

1000-1298(2017)07-0265-07

2017-03-25

2017-05-03

“十二五”國家科技支撐計劃項目(2015BAD07B00)

孫浩(1989—)，男，博士生，主要從事環(huán)境流體力學(xué)和風(fēng)沙物理學(xué)研究，E-mail:251045257@qq．com

劉晉浩(1958—)，男，教授，博士生導(dǎo)師，主要從事林業(yè)裝備自動化及智能化研究，E-mail:liujinhao@vip．163．com

10．6041/j．issn．1000-1298．2017．07．033