沙障地表形態衍化數值模擬方法研究

孫 浩 劉晉浩 黃青青

(北京林業大學工學院，北京100083)

沙障地表形態衍化數值模擬方法研究

孫 浩 劉晉浩 黃青青

(北京林業大學工學院，北京100083)

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

沙障;地表過程;動網格;數值方法

引言

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

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

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

1 數值方法

本文數值模型計算結構與傳統兩相流數值計算結構相同，由流體計算和沙體計算兩部分構成，流體計算采用傳統計算流體力學方法，沙體計算由本研究提出的模型進行計算。

1.1 沙相網格劃分

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

圖1 沙相網格劃分原理Fig．1 Principle of sand phasemeshing

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

1.2 風沙兩相耦合計算過程

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

1.3 風沙相互作用數值模型

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

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

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

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

1.3.1 進入概率a數學模型

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

式中 t——時間，s m——高度調整系數

a(t)——t時刻的風沙流進入概率

a0——初始風沙流進入概率

Ht——t時刻沙障內部地表平均高度，m

Hbarrier——沙障草頭高度，m

1.3.2 逃逸概率b數學模型

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

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

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

n——沙障內子單元數量

1.4 草方格障間風蝕量計算

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

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

Rt——沙粒啟動時臨界剪切風速和當前剪切風速的比值

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

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

u*——剪切風速，m/s

d——沙粒平均直徑，mm

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

草方格內總風蝕率為

式中 i——草方格內子網格單元編號

Qelement(i)——單元i的風蝕率

1.5 子網格單元厚度計算

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

1.5.1 風沙蝕積作用產生的地表沙量變化

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

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

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

圖4中H為子單元厚度，L為子單元順風向長度，W為子單元垂直風向寬度。研究表明，沙床表面從靜止狀態到形成飽和風沙流狀態需要1～2 s的時間，故由于沙源的供給量和承載力不均衡導致的沙量變化也需要一定的時間，本模型中應用沙粒交換速率c對此結果進行修正。子單元網格沙粒變化率計算式為

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

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

子單元高度計算式為

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

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

1.5.2 沙障阻擋及地表起伏產生的沙量沉降

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

圖5 沙障內地表小地形Fig．5 Sub landform in sand barrier checkerboard

沙障內子單元的坡度角變化量計算式為

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

單元i判斷為坡頂地形時，θdiff(i)=0。沙障方格內子單元的概率計算式為

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

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

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

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

2 數值計算

本文以開源流體力學計算軟件OpenFOAM為基礎，編寫沙相數值模型，并與流場耦合進行數值計算。

2.1 幾何模型及網格劃分

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

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

計算域由空氣相計算域和沙障計算域構成，計算域高度為1 m，沙障高度為0.1 m，沙障間距為1m，入口邊界到第1個沙障的距離為1 m，計算域右側沙障到計算域出口的距離為8 m，為了后面便于動網格調整，避免出現低質量網格而增加計算誤差，模型對沙障部分采用孔隙域模型，設置其孔隙為零。

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

圖8 計算域網格劃分Fig．8 Mesh grid of computational domain

2.2 邊界條件及模型參數初始化

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

表1 模型參數設置Tab．1 Model parameters

3 數值計算結果

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

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

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

圖10 仿真結果與實驗結果對比Fig．10 Comparison of simulated results and field measurement results

4 討論

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

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

5 結束語

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

參 考 文 獻

1 龐營軍，屈建軍，謝勝波，等．高立式格狀沙障防風效益［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 屈建軍，凌裕泉，俎瑞平，等．半隱蔽格狀沙障的綜合防護效益觀測研究［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 黃強，雷加強，王雪芹．塔里木沙漠公路不同地貌部位的高立式沙障阻沙特征［J］．干旱區地理，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沙障與傳統草方格沙障防風效益［J］．北京林業大學學報，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］．中國沙漠，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 常兆豐，仲生年，韓福桂，等．粘土沙障及麥草沙障合理間距的調查研究［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 馬學喜，李生宇，王海峰，等．固沙網沙障積沙凹曲面特征及其固沙效益分析［J］．干旱區研究，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 張登山，吳汪洋，田慧麗，等．青海湖沙地麥草方格沙障的蝕積效應與規格選?。跩］．地理科學，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耦合的水力旋流器水沙運動三維數值模擬［J/OL］．農業機械學報，2016，47(1):126－132．http:∥www．j-csam．org/jcsam/ch/reader/view_abstract．aspx?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耦合的迷宮流道水沙運動數值模擬［J/OL］．農業機械學報，2016，47(8): 65－71．http:∥www．j-csam．org/jcsam/ch/reader/view_abstract．aspx?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—)，男，博士生，主要從事環境流體力學和風沙物理學研究，E-mail:251045257@qq．com

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

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