風沙兩相流數值模擬研究進展

王萍，鄭曉靜

蘭州大學 湍流-顆粒研究中心 西部災害與環境力學教育部重點實驗室，蘭州 730000

風沙運動是沙質地表上的沙塵顆粒在風力作用下以地表蠕移、近地表躍移、遠離地表懸移而進行輸運的一種自然現象。常見的沙質地表有占全球陸地面積1/3的沙漠和荒漠化地區[1]，海邊的沙灘和裸露的農田和建筑工地等。沙塵顆粒的輸運，在近地表1～2 m以下的高度范圍內以粒徑大于約100 μm的顆粒躍移運動為主，其輸沙率占總輸沙率的75%以上[2]，是風沙運動的主導；在離地表1～2 m乃至數百米以上高度范圍內以粒徑小于100 μm甚至數μm(如PM10)的顆粒懸移運動為主，其輸沙率不及躍移風沙流，僅占總數不到5%[2]，但這些沙塵粒子可以長期漂浮于對流層上部和平流層并遠距離輸運。當高超聲速飛行器再入行星大氣層時，這些粒子可能與飛行器高速撞擊，使其表面遭到破壞，甚至與飛行器頭部激波層相互干擾，導致飛行器表面加熱率改變[3]。躍移風沙流往往導致土壤風蝕、農林作物受損、交通道路掩埋等危害，而風沙運動的極端情形沙塵暴則往往導致大范圍空氣質量下降、無線通訊信號中斷、汽車火車翻車甚至人員傷亡。例如，2010年4月24—25日在甘肅發生的沙塵暴共造成甘肅19個縣約132.56萬人受災；農作物受災面積達20.68萬公頃。沙塵暴不僅在中國西北地區頻發，在非洲的撒哈拉地區、沙特阿拉伯首都利雅得、澳大利亞悉尼、美國亞利桑那州鳳凰城等也時有發生，而且在地外行星如火星上也經常發生。風沙運動引發的風沙流和沙塵暴已經成為影響人類經濟社會發展的一個重要環境問題，防治風沙災害被聯合國《2030年可持續發展議程》確立為17個可持續發展目標之一[4]，也被列入《全國重要生態系統保護和修復重大工程總體規劃(2021—2035年)》[5]。

風沙運動這種自然現象在本質上是大氣邊界層可蝕地表上摩阻雷諾數高達106～107的氣-固兩相流。大氣邊界層特別是近地表湍流的強非線性、巨量的沙塵顆粒大小和形狀以及起動條件的隨機性、可蝕地表上顆粒沖擊床面和碰撞后的反彈和濺起的復雜性、顆粒運動與湍流相互作用的耦合效應以及沙塵顆粒帶電形成風沙電場[6]與湍流流場的多場耦合、由瞬間啟動和碰撞的微米大小的顆粒導致的在數秒鐘即可達到飽和的數米高度范圍的風沙流、進而造就長達數百年形成演化出的數千平方公里范圍沙漠的這種時空多尺度和跨尺度性等，使得風沙運動具有當前許多學科面臨的多重非線性、多場耦合性、多尺度交織的共性，其定量模擬和準確預測十分困難。著名流體力學家Balachandar在其綜述文章中指出：“湍流和多相流是流體力學中最具挑戰性的2個問題，湍流所固有的隨機性由于顆粒的隨機分布而變得更為復雜。顆粒相的存在使得多相流的試驗觀測和數值模擬比單相流復雜得多……即使是在稀疏擴散相的條件下，當兩者聯合起來也會成為更加難以對付的挑戰”[7]。而《Science》列出的125個未解決的重要科學問題[8]中唯一一個直接與力學相關的問題就是“能否發展湍流動力學與顆粒運動學的統一模型”。

風沙運動的科學研究始于20世紀40年代，早期以野外觀測和風洞(包括固定式風洞和便攜式風洞)試驗為主要研究手段。這一時期最為重要的著作是英國工程師Bagnold所著的《The physics of blown sand and desert dunes》[2]。他運用流體力學當時的最新研究成果，即普朗特邊界層理論,結合野外考察與風洞試驗，為風沙物理學奠定了基于試驗觀測的經驗性定量分析的基礎。在他的研究中，近地表的躍移風沙流被劃分為可蝕床面的沙粒在風場作用下的流體起動，起躍沙粒在風場作用下的躍移運動，躍移沙粒與可蝕床面的沖擊碰撞、反彈及床面沙粒的濺起，躍移沙粒對風場的反作用這4個子過程。但他只是對這幾個過程進行了定性的描述。20世紀60～80年代，由于人類對太空探測的興趣，研究者們開始了風成地貌與風沙運動的定量模擬以對無法直接觀測的地外行星，如火星、金星上的風沙運動進行預測。這一時期的開創性工作是Owen[9]發表的《Saltation of uniform grains in air》。他通過假定所有從床面起跳的沙粒運動服從同一簡單形狀的軌道，在歐拉-拉格朗日框架下建立了考慮重力和氣流對沙粒拖曳力的沙粒運動方程、考慮沙粒反作用力與剪應力梯度平衡的流體速度方程，理論推導出受沙粒運動影響的流場速度分布、躍移風沙流的輸沙率等。模型中，Owen通過假設拖曳力正比于沙粒流體相對速度VR，將二維沙粒運動方程組解耦為2個二階常系數非齊次微分方程，因而獲得沙粒速度的理論解。但這一假設僅在沙粒雷諾數Rep(Rep=VRdp/ν, 其中dp為沙粒粒徑，ν為空氣運動黏度)小于1、阻力位于Stokes區的情況下近似成立，而實際躍移沙粒的雷諾數應為1

20世紀80年代以后，隨著對兩相流問題研究的不斷深入，逐漸發展起風沙兩相流的三流體模型[13](即將向上和向下運動的沙粒用2種擬流體表征)和懸移粉塵平衡模型[12](即粉塵的沉降與擴散平衡)等。與此同時，一些學者，如Ungar 和Haff[14]以及Werner[15]，逐漸開始在風沙流的模擬時考慮躍移風沙流的第3個子過程，即躍移沙粒與可蝕床面沖擊碰撞后的反彈/濺起，同時更加完整地考慮躍移風沙流的第4個子過程，即躍移顆粒對風場的反作用。這方面的代表性工作為Anderson 和Haff[16]于1988年發表的《Simulation of eolian saltation》。他們把由試驗測得的躍移顆粒起跳初速度分布和由單顆顆粒單次沖擊可蝕床面的沖擊速度與濺起顆粒數/速度的擬合關系(擊濺函數)作為顆粒運動的初值條件，并通過顆粒對流場的反作用將顆粒運動方程與流場方程耦合。不僅如此，顆粒運動方程中表征拖曳力的風場速度也不再僅由平均速度廓線確定，而是由雷諾平均Navier-Stokes(RANS)方程求解后代入。通過聯立求解顆粒運動方程和流場方程，進而實現對躍移風沙流形成和演化直至飽和的全過程的數值模擬。由此把20世紀70年代的流體力學在兩相流理論和Navier-Stokes方程的RANS求解用于風沙流的模擬，為風沙物理學奠定了基于數值模擬的統計性定量分析的基礎。然而，基于RANS僅能模擬湍流平均風場及其與顆粒的相互作用，所得到的輸沙率往往隨時間、空間不變，這與實際觀測發現的輸沙存在顯著時空變化[17-20]不符，而這種時空變化被認為與湍流結構[21-23]，尤其是外區的大尺度/超大尺度運動(Large Scale Motions/Very Large Scale Motions, LSMs/VLSMs)[24-26]密切相關[27-29]。同時，風洞試驗與野外觀測均發現了低于風沙物理學中所給出的理論起動風速條件下的輸沙[30-32]，意味著湍流的近壁結構對起沙具有重要作用[33-34]。因此，在對風沙兩相流預測中僅用RANS描述風場是遠遠不夠的，需要更多考慮湍流脈動與結構影響。

本文將回顧風沙運動數值模擬近30年的研究進展。注意到，盡管兩相流數值模擬方法大致可分類為歐拉-歐拉、歐拉-拉格朗日、拉格朗日-拉格朗日3類，但已有風沙運動數值模擬仍以前2類為主，因此本綜述也僅側重前二者。在第1節將介紹基于RANS的風沙兩相流數值模擬，并作為應用，給出風成地貌沙波紋和沙丘場形成演化過程的模擬。第2節介紹壁解析大渦模擬(Wall Resolved Large Eddy Simulation，WRLES)方法在風沙兩相流模擬研究中應用的現狀及前景。第3節介紹基于壁模型大渦模擬(Wall Modeled Large Eddy Simulation，WMLES)的風沙運動數值模擬研究，第4節給出風沙運動數值模擬未來發展的趨勢以及值得進一步研究的若干關鍵問題。

1 基于RANS的風沙運動數值模擬

RANS方法由于僅求解湍流平均運動方程，因此，能夠以較小的計算量得到工程上所關心的平均統計特性。假定大氣表面層氣流(Atmospheric Surface Layer，ASL)為不可壓縮流動，其所滿足的雷諾平均動量方程和連續性方程為

(1)

(2)

式中：xi為坐標，i=1,2,3分別代表流向x、垂向y、展向z；Ui為流體平均速度，流向、垂向、展向速度分量分別記為U、V、W；t為時間；p為壓力；ρ、μ分別為空氣的密度和動力黏度；Fi為顆粒對不同高度上流體的反作用力；τij=μt·(?Ui/?xj+?Uj/?xi)為雷諾應力，μt為渦黏系數。風沙運動RANS模擬中一般采用混合長度模型和k-ε模型等進行方程封閉。沙粒運動采用拉格朗日點力模型計算時，當考慮拖曳力FDi、重力FG和電場力FE時，所滿足的運動方程為

(3)

式中：upi為沙粒速度，流向、垂向、展向速度分量分別記為up、vp、wp；ρp、mp分別為沙粒的密度和質量；g為重力加速度；Ei為電場強度；qE為沙粒所帶電荷量；CD為拖曳力系數；u@i為沙粒所在位置處流場速度；δi2為Kronecker記號。|u@p-up|為流場與沙粒速度矢量差的模。需要指出的是：附加質量力、Basset力、Saffman升力、沙粒旋轉產生的Magnus力以及顆粒空中碰撞產生的碰撞力等也對沙粒運動以及輸沙率等有定量影響，具體可參見Zheng的專著《Mechanics of wind-blown sand movement》[1]。

Ungar和Haff[14]運用混合長度理論封閉無限大平坦沙床表面的一維氣流動量方程，首先建立了風場-沙粒相互耦合的穩態風沙流數學模型，即將式(1)簡化為

(4)

式中：Fx為顆粒對流體反作用力的流向分量；κ為von Kármán常數。假定所有從床面起跳的沙顆粒具有相同的起跳速度和彈道形式的軌跡，其軌跡方程僅考慮拖曳力和重力。沙顆粒在一次飛行過程中2次經過同一高度y，則其對單位體積風的阻力可表示為

(5)

式中：N1為單位時間內從單位面積床面上起跳的沙粒數；vp↑及vp↓分別為上升沙粒與下降沙粒在高度y處的垂向速度分量；FDx↑及FDx↓分別為相應沙粒上所受拖曳力FDi的水平分量。在此基礎上，Ungar和Haff[14]提出了一個簡化的δ函數形式的擊濺模型，描述沙粒以速度Vimp沖擊床面后，以V0起跳的沙粒數：

(6)

(7)

以風速U(y)廓線和目標剪切風速u*作為收斂條件，通過雙重迭代可獲得風沙流穩定狀態。這一方法也被Zheng等[35]進一步應用，揭示出輸沙率沿高度的分層現象，即：輸沙量在貼近沙床附近處隨高度線性增加、達到極值進入飽和層、然后沿高度按負指數規律衰減。這為“沙割”現象的解釋提供了理論依據。需要指出的是，此類方法只通過顆粒-流場的耦合計算獲得穩態自平衡的風沙流模擬結果，無法再現風沙流的形成發展過程。

Anderson和Haff[16,36]率先完整實現了包含風沙運動全部4個子過程的、從沙粒起動到風沙流自平衡狀態的時間發展全過程的模擬，隨后McEwan和Willetts[37-38]又對風沙流的這一發展過程進行了更為細致的模擬。在他們的模擬中，流體起動顆粒數Na被假設為正比于流體剪切應力與床面顆粒起動的臨界剪切應力之差

Na=Ca(τw-τc)

(8)

式中:τw為床面流體切應力；τc為風力直接吹起沙粒時的臨界應力，可由理論分析結合風洞試驗得到；Ca為比例常數。擊濺過程則通過離散單元法(Discrete Element Method，DEM)模擬或采用試驗觀測給定的經驗函數[1,36,39-40]來計算。Anderson和Haff[36]基于其二維DEM(沙粒碰撞簡化為2個均質圓盤)模擬結果提出

Pr=0.95(1-exp(-2.0Vimp))

(9)

(10)

圖1 單寬輸沙率隨時間發展過程[1]Fig.1 Time evolution of sand transport rate[1]

風沙流/沙塵暴中由于顆粒間、顆粒與床面間的接觸摩擦等會導致顆粒帶電并形成風沙電場，其強度可到數十至上百kV/m[47-49]。帶電沙粒所受靜電力的大小與沙粒自身的重力相當，因而風沙電場直接影響沙粒運動。同時風沙電場大小與沙粒帶電量以及空中沙粒數密切相關。因此，由沙粒帶電所產生的風沙電場與沙粒運動和風場之間也存在著相互耦合作用。在均勻穩態風沙流中只存在垂向電場，Zheng等[50]結合風洞試驗研究所獲得的沙粒帶電和風沙電場規律，建立了考慮風-沙-電多場耦合的風沙流模型，如圖2[1]所示。根據試驗觀測，假定空中沙粒在平均意義下帶負電荷，沙床表面沙粒即蠕移沙粒在平均意義上帶正電荷，空中某高度y處的垂向電場E(y)視為由晴天大氣電場E0、空中帶負電荷的沙粒在點y處所產生的電場Es、沙床表面帶正電荷的沙粒在點y處所產生的電場Ec這3部分組成：

圖2 點電荷產生電場示意圖[1]Fig.2 Schematic diagram of electric field produced by a point charge[1]

E(y)=E0(y)+Es(y)+Ec(y)=-0.1+

(11)

式中：N為單位體積內空中運動的沙粒數，在穩態風沙流中只隨高度變化；ε0為空氣介電常數；cq=qE/mp為荷質比，即單位質量沙粒所帶電荷量，試驗觀測給出其均值約在±O(100 μC/kg)范圍內[47-48]；a和b為電荷影響空間范圍，實際計算發現，a和b的取值大于100 m后電場計算結果不受影響。數值模擬結果表明，考慮垂向風沙電場，帶正電(荷質比為60 μC/kg)的顆粒運動軌跡的長度比不帶電的增加了163%，帶負電(荷質比為-60 μC/kg)的卻減小了57%[47]。而當顆粒帶正電(荷質比為60 μC/kg)時，模擬的多場耦合風沙流隨時間演化過程能夠與風洞試驗觀測結果更吻合(見圖1)。但是式(11)中沙粒平均帶電量(或荷質比)為常數，不能反映沙粒碰撞帶電物理過程。為此，Kok和Renno[51]基于硬球模型和接觸勢差帶電機制，模擬沙粒碰撞帶電以及多場耦合風沙流。該模型能夠解釋“大顆粒帶正電小顆粒帶負電”現象，但其中涉及的物理參數為經驗參數，需試驗測定。Hu等[52]進一步采用軟球模型計算碰撞顆粒間電荷轉移，不同于硬球模型不考慮碰撞過程中的細節，軟球模型能夠反映顆粒尺寸和碰撞速度對于單次碰撞帶電量的影響。帶有電量qEi、qEj的2個顆粒碰撞后，各自的帶電量分別為

(12)

其中：ρH,0為高能態束縛空穴面密度；e為單位元電荷；CE為常參數；Si、Sj為碰撞接觸面積。將該模型應用至風沙流中可以更好地模擬荷質比隨著風速的變化規律。此外還發現，顆粒荷質比服從正態分布且其均值隨著空中碰撞效應的增強而逐漸趨向于0。

最近，一些學者結合RANS流場還發展了基于二維DEM的風沙流模擬方法，直接在顆粒尺度上求解沙粒-床面碰撞過程。如Carneiro等[53]采用此方法模擬了500顆具有高斯分布粒徑的球形沙粒床面上的風沙流，發現輸沙率在低風速下存在亞穩定狀態。Durán等[54]數值模擬15 000顆粒沙粒組成的床面上的沙粒運動，分析不同密度比(風沙與水沙)條件下輸沙規律的異同；在Durán等[54]僅包含重力、浮力和拖曳力的顆粒運動模型基礎上，P?htz和Durán[55-56]進一步在模擬中考慮了附加質量力，分析床面顆粒平均速度及其與躍移層平均顆粒速度關系、風沙流起動與輸沙停止的臨界風速條件、流體顆粒密度比以及Galileo數的影響等。雖然DEM模擬可以避免采用經驗擊濺函數帶來的誤差，但由于計算量限制，已有的數值模擬工作或者計算顆粒數較少、或者計算域很小，遠未達到實際風沙運動的顆粒數量和空間尺度范圍。

以上基于RANS流場的穩態風沙流數值模擬對于理解風沙運動起到了極大的促進作用，但RANS僅模擬平均流動，平均掉了所有的湍流脈動。事實上，由于顆粒所受拖曳力正比于其迎風橫截面積，粒徑越小、受力越小、軌跡越低，因此即便粒徑相當小的懸浮顆粒(<10～20 μm)，在平均流場中采用點力方法模擬，其飛行高度也僅限于離地表約1～2 m以下，與沙塵天氣中沙塵顆粒能夠被傳輸至對流層上部進行遠距離輸運的實際不符。在RANS模擬中計及風場的湍流脈動效應，目前主要有2類途徑，一種是在穩態風沙運動模型的基礎上直接施以非定常風速，如Spies等[57]、作者團隊[58]通過在計算域上邊界施加正弦變化的周期性風速脈動，模擬非平穩風沙流，獲得風沙流輸沙率Q(t)隨來流風場參數(周期與幅值)的變化規律，如圖3[58]所示，其發現脈動風速增加風沙流輸沙率，平均風速越小，增加越顯著。Zhang等[59]進一步考慮顆粒帶電的非穩態風沙流模擬發現，除了垂向電場之外，由于輸沙率在水平方向的不均勻性，風沙流發展過程中還會形成水平方向的風沙電場。另一種途徑是借鑒描述重顆粒擴散的拉格朗日隨機模型，通過在平均速度廓線上疊加采用郎之萬方程描述的流場速度脈動[60-62]或直接對顆粒運動采用修正郎之萬方程[63]描述的方法求解顆粒軌跡。以前者為例，該方法即在式(3)中對u@i的脈動部分u′@i采用以下表達式計算:

圖3 在正弦變化風場中的輸沙率Q(t)演化[58]Fig.3 Development of sand transport rate Q(t) in sinusoidal wind variations [58]

(13)

圖4 典型沙塵顆粒隨機運動軌跡[67]Fig.4 Typical stochastic trajectories of sand/dust particles[67]

對于小粒徑懸移粉塵的遠距離、大尺度輸運的計算，除采用上述拉格朗日隨機模型方法以外，還可以采用歐拉-歐拉模型求解。實際上，現有的沙塵數值預報模式[71-75]中多采用歐拉-歐拉模型計算沙塵濃度分布、提供預報數據。這主要是因為如果采用拉格朗日隨機模型跟蹤每一顆粒的運動軌跡，為獲得統計可靠的結果，需要跟蹤的顆粒數量巨大，使該方法難以在數值預報中實際應用。在歐拉-歐拉模擬中，沙塵濃度c服從有源/匯項的、對流-擴散類型的守恒方程[76]，即

(14)

風成地貌如沙波紋、沙丘和沙丘場，目前尚未見按連續介質理論成功模擬的例子。這主要是因為涉及沙粒尺度(10-4m)、風沙流尺度(10-1～100m)、沙丘和沙丘場尺度(100～104m)至少3個長度尺度，以及沙粒起跳和躍移的時間尺度(10-2～10-1s)、風沙流形成和持續的時間尺度(100～103s)、沙丘和沙丘場形成和發展的時間尺度(104～108s)至少3個時間尺度。同時，在不同時空尺度上的現象和行為有著不同的結構層次和不同的演化速率以及不同的物理機制，因而很難建立涉及地表侵蝕、風沙流沉積等典型過程的本構模型。結合沙粒的散體性和基于離散建模的計算機模擬風成地貌的方法(如元胞自動機)，雖然在不同程度上也能模擬出沙波紋和沙丘, 但由于均需要賦以一些人為的規則而使得模擬結果與實際風成地貌不符。基于RANS模擬給出的沙粒沖擊地表速度和輸沙率以及平均躍移長度等統計量，Zheng等[77]提出的離散粒子追蹤法(Discrete Particle Tracing Method， DPTM) 對風成沙波紋的發展過程實現了計算機模擬。該方法同時考慮了與真實沙波紋形成相關的3個主要因素, 即不同粒徑的眾多沙粒在沙床上方所形成的風沙流、風沙流中沙粒對沙床的碰撞以及碰撞后沙粒的反彈或濺起、沙粒在床面上方的躍移和在床面的蠕移。其中，風沙流模擬給出的顆粒沖擊地表作為由混合粒徑顆粒組成的離散床面模擬系統的能量輸入。其模擬結果(見圖5[1])不僅再現了沙波紋從平坦沙床出現—長大—分叉—合并—飽和的全過程，而且模擬得到的沙波紋尺寸和移動速度、沙波紋的自修復過程、沙波紋頂部沙粒粒徑大底部粒徑小的倒粒序結構、以及沙波紋出現和消失的臨界摩阻風速等特性與觀測結果符合。

圖5 離散粒子追蹤法模擬得到的混合粒徑沙波紋(摩阻風速為0.5 m/s)[1]Fig.5 Sand ripple in varied-sized case simulated by DPTM (friction wind velocity is 0.5 m/s)[1]

如果說對于沙波紋的形成演化過程可以通過對一顆顆沙粒的運動模擬得到，那么對于沙丘與沙丘場，其計算量顯然是不允許的。為此，仍然基于對風沙流的RANS模擬，Zheng等[78]提出了一種類似于三級跳遠的“三級跳”方法(見圖6[1])，即：耦合尺度沙丘場模型——基于近地表輸沙的典型離散模型。該模型不僅需要擊濺函數和沙粒起跳速度分布這些涉及沙粒近地表運動的已有統計量，而且其還推導出床面的侵蝕因子和覆蓋因子以及沙粒輸運的傳輸因子這3個新的統計量。具體的推導可見文獻[1]。該模型成功實現了從一顆顆沙粒運動到風沙流再到“沙體元”進而到數百平方公里甚至更大范圍的沙丘場從任意初始狀況(包括平坦情況)形成和發展的長達數十年或數百年演化過程的定量模擬，實現這一具有時空尺度在108～109量級的跨越。其模擬結果除了定性上與野外觀測結果一致外，如沙漠邊緣的新月形沙丘場(見圖7(a))和沙漠腹地的線形沙丘場(見圖7(b))，在定量上，如沙丘運動速度隨高度的變化(見圖7(c)[79-80])和沙丘高度隨沙源厚度Hd的變化(見圖7(d)[81])也與野外觀測結果吻合。在此基礎上，Zheng等實現了對沙漠邊緣擴展速度的理論預測，使得沙漠擴展速度從原來的“后知后覺”到提前獲知；其還給出為遏制擴展速度所用的固沙草方格的優化鋪設方案，使得固沙成本可降低2/3。

圖6 風成沙丘場演化過程定量模擬示意圖[1]Fig.6 Schematic diagram in quantitatively simulating evolution process of aeolian dune field[1]

圖7 模擬沙丘場及其統計結果[1,78]Fig.7 Simulated sand dunes and their statistics[1,78]

2 風沙運動的WRLES

如前所述，湍流結構對風沙運動中沙塵顆粒的起動、輸運均有重要影響，因此，理論上，湍流直接數值模擬(Direct Numerical Simulation， DNS)更適于近地表風沙運動研究。事實上，基于DNS與點力方法的數值模擬已經廣泛應用于均勻各向同性湍流、自由剪切湍流、槽道、管道、湍流發展邊界層中的顆粒兩相流動模擬[82-95]。在相對簡單的固壁條件下，揭示了顆粒在壁湍流邊界層中的聚集特性、顆粒對湍流的調制規律。如Zhao 等[88-89]發現球形顆粒與槽道湍流相互作用會增加動能耗散，同時顆粒運動會使流場湍動能在空間再分配，顆粒湍流相互作用還抑制顆粒的轉動行為；Wang和Richter[92]發現慣性顆粒在外區傾向于分布在湍流高速區，同時顆粒增強湍流LSMs/VLSMs，且這種增強隨顆粒Stokes 數(顆粒慣性時間尺度τpar與湍流時間尺度τtul之比)非單調變化，顆粒在湍流內、外區都存在與當地湍流結構相關的非均勻分布。對于尺寸大于湍流最小尺度(Kolmogorov尺度)的顆粒，不能再將顆粒視作一點，顆粒需要被精細分辨并考慮顆粒表面的無滑移邊界條件[96-98]，采用浸沒邊界法計算顆粒與湍流間的相互作用，此即全分辨模擬，如圖8[99]所示。最近，湍流DNS模擬結合顆粒全分辨模擬也被應用于分析可蝕地表上的顆粒運動及其與湍流的相互作用[100-104]，如Ji等[100-101]采用全分辨方法與湍流DNS模擬可蝕床面上躍移顆粒兩相流，分析顆粒與湍流的相互作用，發現湍流下掃與顆粒的流體濺起密切相關，而顆粒運動將破壞近床面湍流結構。Vowinckel等[102]采用類似的方法研究全分辨可蝕床面上的顆粒湍流相互作用，不僅發現湍流被增強，還發現可蝕性導致高質量分數下床面形成長達12倍計算域高度的丘狀顆粒分布、低質量分數下形成顆粒聚集的現象等。作者所在團隊也開展了可蝕地表大顆粒(～500 μm)與湍流相互作用的全分辨直接數值模擬研究[105]，全分辨球體瞬時分布與流場速度分布如圖9[105]所示(該算例模擬總顆粒數為7 200)；利用DEM解析顆粒與顆粒、顆粒與床面之間的碰撞，基于該模擬結果識別流體起動、反彈、濺起事件，分析了顆粒速度分布。在此基礎上，針對全分辨直接數值模擬的大規模并行計算，提出了按顆粒數劃分的并行算法，并對顆粒信息存儲進行了優化，使得可計算顆粒數能達到百萬量級。

圖8 球形全分辨顆粒模擬的歐拉網格點和拉格朗日網格點分布[99]Fig.8 Distribution of Eulerian grid points and Lagrangian grid points over a fully resolved sphere [99]

圖9 全分辨直接數值模擬顆粒分布示意圖[105]Fig.9 Schematic diagram of fully resolved particle distribution in DNS[105]

圖10 沙塵暴沙墻結構與異重流模擬Fig.10 Sand wall and numerical simulation based on gravity current model

壁湍流的DNS，即便單相流動，目前所能達到的最高摩阻雷諾數Reτ(Reτ=uτH/ν，H為邊界層厚度，uτ為湍流壁面摩擦速度)也一直在O(103)徘徊[108-111]。而大部分兩相流數值模擬，考慮到顆粒帶來的額外計算量，Reτ則多限于O(102)。目前基于點力方法的最高雷諾數可能是Bernardini[90]的Reτ=1 000的單向耦合槽道兩相流模擬、Wang 和 Richter[92]的Reτ=950慣性顆粒雙向耦合兩相流模擬、Jie等[112]Reτ=1 000 的非球形顆粒兩相流模擬。而對于顆粒全分辨模擬，由于要分辨每個顆粒周圍的流場，計算量進一步急劇增大，如Kidanemariam 和 Uhlmann[103]中顆粒數105萬、雷諾數Reτ=293的單個算例需要計算500萬核時，這也導致目前全分辨顆粒兩相流的湍流最高雷諾數僅為Reτ=647[100]，遠低于大氣風沙兩相流起沙雷諾數(野外Reτ～105，風洞Reτ>2 000～3 000)。

湍流大渦模擬直接求解濾波尺度以上的“大渦”，而對亞格子尺度“小渦”建立模型(即亞格子模型)，從而能夠以更小的計算量得到接近DNS 的湍流脈動與結構信息，因而被認為是最具潛力的湍流模擬工具[113]。如果顆粒尺寸小于LES最小可解尺度或顆粒密度遠大于流體，點力大渦模擬是實現高雷諾數湍流兩相流數值模擬最具吸引力的方法[114]。大渦模擬的控制方程為濾波后的Navier Stokes方程和連續性方程。

(15)

(16)

圖11 采用歐拉-拉格朗日方法模擬湍流顆粒流的概念分類[115]Fig.11 Conceptual classification of Eulerian-Lagrangian modelling approaches for simulation of particle-laden flows [115]

在較低的雷諾數下，即便不考慮亞格子湍流對顆粒運動影響，基于WRLES流場的兩相流數值模擬也被證明可以獲得與基于DNS流場模擬一致的結果，且亞格子模型(Smagorinsky模型和動力Smagorinsky模型)對模擬結果影響較小[117-118]，因此，已逐步被用于分析顆粒的沉積和擴散以及顆粒對湍流的調制[119-127]等研究。值得一提的是，最近WRLES結合DEM方法對可蝕床面上懸浮兩相流數值模擬的湍流雷諾數與模擬的顆粒總數相比DNS都得到了大幅提升。比如在Schmeeckle[124]的研究中，最大摩阻雷諾數約為3 470，模擬最大顆粒總數為33萬，而Capecelatro和Desjardins[123]在圓管泥沙流動的WRLES結合DEM模擬中雷諾數約為2 000，顆粒數達到了1 600萬。對于風沙運動，本文作者團隊采用WRLES結合點力模型對近地表躍移風沙兩相流進行了較為精確的數值模擬[128]，雷諾數最高達Reτ=4 200，亞格子湍流對可解尺度湍流的貢獻采用動力Smagorinsky模型計算，同時忽略了顆粒與亞格子湍流之間的相互作用，即顆粒與湍流的相互作用只發生在可解尺度上。采用三維擊濺函數(將式(9)擴展為三維[129])作為沙粒運動下邊界條件，同時考慮沙粒流場雙向耦合，單時間步跟蹤總沙粒數達5 000萬。數值模擬統計量與風洞試驗結果吻合，且再現了沙粒在風洞尺度上的大尺度分布(見圖12)[128,130]，也這也是目前唯一采用WRLES精細模擬風沙運動的工作。

圖12 風沙運動WRLES模擬所給出的沙粒體積分數廓線和瞬時分布[128]Fig.12 Volume fraction profiles and instantaneous distribution of sand particles obtained by WRLES simulation[128]

3 風沙運動的WMLES

壁湍流邊界層含能渦尺度隨距離壁面高度的增加而迅速減小。WRLES求解內層(邊界層厚度的約10%以下)所需網格數幾乎與DNS同一量級，且隨Re的增加而增加。根據Piomelli[131]估算，當Re>106，超過99%的計算資源都用于求解內層，加之顆粒跟蹤，因此湍流兩相流模擬計算量仍是相當巨大的。例如本文作者團隊Reτ=4 200 的風沙運動WRLES算例在OpenMP與Message Passing Interface(MPI)混合并行的情況下需要約40萬核·時才能獲得統計穩定的結果。這使WRLES難以應用至雷諾數高達105的野外實際風沙運動模擬。一種有效的湍流流場求解方案是采用壁模型大渦模擬，即：直接求解外區含能渦，內區對動量傳輸有貢獻的渦由于網格非常粗糙而無法求解，因而從整個動力系統中移除，這將極大減小計算量。此時，需要將壁面應力與外區可解尺度速度相關聯或在雷諾平均意義下對內區動量傳輸的影響進行求解以獲得正確的速度廓線[131-132]。目前已經發展了許多的WMLES方法，大致上可分為混合RANS/LES模型和壁面應力模型2類。應用最為廣泛的壁面應力模型是平衡率壁模型，通過假定速度廓線代數求解壁面應力或通過局部求解簡化RANS方程估算壁面應力。Deardorff[133]和Schumann[134]最早應用了此類模型。此后，通過考慮浮力[135]、內區結構傾角[136]、外區湍流結構[137]、慣性項以及粗糙子層[138]、輻射[139]等影響不斷修正以包含更為復雜的近壁物理和非平衡過程。

另一方面，在應用了壁模型的風沙運動WMLES工作中，其所采用的壁模型主要是通過假定速度廓線計算壁面應力，即假定近壁區存在當地可解尺度速度分布與當地剪應力之間的平衡關系，流向與展向壁面應力τw,yx和τw,yz的計算式為

(17)

(18)

(19)

其中：Δs是為考慮近壁結構傾角而引入的，以使得壁面應力與下游Δs處瞬時速度相關[136]。值得一提的是，文中采用了尺度依賴動力亞格子模型[144]，在粗網格、低精度數值格式下再現高雷諾數湍流中的LSMs/VLSMs結構，不僅模擬出與野外實際觀測一致的長度可達數十米的沙粒結構，還通過條件平均分析發現“Sand Streamers”出現在 LSMs/VLSMs 近壁面尾跡中，即“Sand Streamers”是 LSMs/VLSMs 在近壁面的“足跡”，見圖13[143]，由此揭示出 LSMs/VLSMs 對風沙兩相流中顆粒運動的影響，證明在風沙運動中湍流結構的重要作用。

圖13 風沙流中以高濃度c′>0條件平均的流向速度脈動u′[143]Fig.13 Conditionally averaged fluctuating streamwise velocity u′ for c′>0 in wind-blown sand flow[143]

以上采用壁模型進行風沙運動大渦模擬的工作中，所采用的壁模型本身也是不同的，雖然在風速較大時，沖擊濺起過程主導近床面沙粒運動，使得不同壁模型對于風沙運動及其統計量的預測差別較小，但對于較低風速、流體起動占優的情況，床面附近速度與應力預測變得十分重要。Zheng等[145]對比了式(14)和Sidebottom等[146]提出的考慮外區LSMs/VLSMs對壁面應力調制的內外尺度相互作用(Inner-Outer Scale Interaction, IOSI)模型發現，在低風速下，尤其是剪切風速低于臨界起動條件時，IOSI模型能預測更強的、與試驗更符合的單寬輸沙率，如圖14[145]所示。據此，他們認為，外區湍流結構尺度(或邊界層厚度)影響沙塵流體起動，這意味著野外條件下發生輸沙的臨界條件可能低于風洞，或者換句話說，同樣剪切風速下，野外輸沙率可能大于風洞的。

圖14 低風速條件下計算無量綱單寬輸沙率隨無量綱剪切應力變化(根據Zheng等[145]改繪)Fig.14 Dimensionless transport rate as function of dimensionless shear stress at low wind speed according to Zheng et al. [145]

此外，即便在以上采用的壁模型兩相流WMLES中，所使用的壁模型仍是基于單相流提出的，未考慮顆粒對壁模型的影響。本文作者團隊基于Reτ=4 200的 WRLES風沙流模擬結果后驗考察現有各類用于高雷諾數單相壁湍流的壁模型，發現WMLES在預測顆粒通量時存在可達 100%以上的差異。進而提出通過在積分壁模型(Integral Wall Model，IWM)[138]中引入顆粒體力項獲得含顆粒影響的兩相流大渦模擬壁模型，所得的WMLES結果與WRLES 結果的誤差降至20%以下[128]，這為高雷諾數顆粒兩相流WMLES提供了重要的途徑，但在積分壁模型中仍假定平均速度廓線為對數律與線性律的線性組合，而風沙流中、尤其是躍移層內流場速度廓線形式目前尚無統一結論，需要進一步研究。

WMLES也被應用于模擬大氣湍流邊界層中懸移粉塵地表釋放和分布。如Gu等[147]模擬了半徑50 m、高150 m圓柱形計算域內的塵卷風現象——一種發生在大氣對流邊界層內、能將沙塵或者碎屑等物體揚到高空、具有溫度較高的低壓核心和較短生命周期的旋風。他們基于拉格朗日方法計算了不同粒徑沙塵顆粒運動軌跡，發現較小的粉塵跟隨塵卷風氣流，而較大沙塵由于受到離心力作用偏離主流向外運動并在旋風外圍沉降，又被匯流層的高速氣流攜帶到旋渦中心處，往復循環直至無法被上升氣流揚起，進而解釋了塵卷上部細顆粒富集、下部大顆粒富集的機理。Klose 和Shao[148]在美國環境預測中心(NCEP)，美國國家大氣研究中心(NCAR)等機構開發的 WRF(The Weather Research and Forecasting)模式的大渦模擬模塊基礎上發展了大渦沙塵模式，并模擬了不同風速與穩定度條件下的地表粉塵釋放。Zhang等[149]數值模擬了Reτ= 3×107中性大氣邊界層中粉塵輸運，其中粉塵采用歐拉模型模擬(式(14))，證實粉塵濃度也存在超大尺度結構，如圖15所示。通過相關分析發現，濃度與流場流向速度脈動負相關而與垂向速度脈動正相關。此外，他們通過條件平均分析與象限分析揭示粉塵高低濃度區伴隨流場反向旋轉準流向渦對，沙塵的垂向輸運與流體大尺度渦卷模式密切相關。

圖15 y/H=0.05高度處流向展向平面內瞬時流向脈動速度和濃度分布模擬結果(粒徑dp=10 μm)[149]Fig.15 Visualizations of instantaneous fluctuating streamwise velocity and concentration in streamwise-spanwise planes at wall-normal position of y/H= 0.05 (particle diameter is dp=10 μm)[149]

4 結 論

風沙運動及其引發的沙塵暴、沙丘移動等自然災害是廣受關注的環境問題。采用數值模擬方法研究風沙運動是揭示輸沙機理、提高預測精度的重要手段。盡管風沙運動數值模擬從一維發展至三維、從簡單軌跡計算發展至多場耦合模擬、從僅考慮平均流的雷諾平均模擬發展至考慮湍流脈動和湍流結構的大渦模擬，但仍存在需要進一步深入研究的問題和改進之處。

1) 風沙起動、擊濺與湍流模型。大尺度風沙運動模擬中，地表顆粒很難單顆分辨，采用模型提供下邊界條件是必然的選擇。但是，目前多數基于歐拉-拉格朗日框架的風沙流模擬仍沿用了早期不精確的模型，如采用由沖擊試驗或基于平均流場模擬獲得的擊濺函數，再如沙粒的流體起動則被簡單地假設為正比于平均剪切應力和臨界起動應力之間的差值。而實際上，近壁(近床面)湍流速度、應力脈動對于沙塵顆粒流體起動以及沖擊過程中地表顆粒的運動和濺起十分重要，目前尚未見考慮湍流影響的模型。實際沙粒地表為多種沙塵顆粒組成的混合粒徑地表，而已有模型大多采用平均粒徑，也是導致現有數值模擬預測結果與觀測存在誤差的重要原因之一。此外，沙粒運動與湍流的相互作用也影響湍流本身的統計量與結構，現有大渦模擬風沙流中所采用的亞格子模型和壁模型都是基于單相湍流提出的，盡管本文作者團隊嘗試考慮顆粒影響改進已有壁模型，但仍值得進一步研究。

2) 高精度、高效風沙流數值模擬方法的發展。湍流直接數值模擬或壁解析大渦模擬結合顆粒全分辨模擬或DEM方法是認識湍流輸沙物理機制、準確刻畫顆粒湍流相互作用、提煉流體起動規律、分析沙粒與可蝕床面沖擊濺起過程的重要方法。但是所需計算量巨大。一方面，高雷諾數湍流邊界層模擬既需網格足夠小以解析最小湍渦，又需計算域足夠大以再現對湍動能與物質輸運起主要貢獻的LSMs/VLSMs；另一方面，對于顆粒的精細模擬，如全分辨方法，為了解析顆粒，每一個方向顆粒上拉格朗日點數量應為～O(10)，模擬所需網格量相比點力方法至少再增加2個量級。同時，無論全分辨顆粒模擬還是DEM模擬，顆粒碰撞搜索過程耗時且均需存儲大量顆粒碰撞信息。因此，亟待發展低存儲需求、高效、高精度的湍流風沙兩相流模擬的并行優化算法，以開展更大規模的計算、獲得足夠的統計樣本。另外還需要特別指出的是，雖然近年來，拉格朗日-拉格朗日框架下如基于光滑粒子流體動力學 (Smo-othed Particles Hydrodynamics, SPH) 方法和格子Boltzmann方法流體求解的兩相流數值模擬方法[150-152]也開始發展，但對風沙運動的模擬尚處于定性一致階段。

3) 風沙運動大尺度、跨尺度模擬的模型改進。風沙運動發生時，在垂向同時存在近地表躍移高濃度區和遠離地表懸移顆粒的稀薄濃度區，本質上也是多尺度問題。已有研究主要對于躍移沙粒和懸移粉塵分別在歐拉-拉格朗日和歐拉-歐拉框架下模擬，并通過水平躍移輸沙通量和垂直粉塵釋放通量將二者聯系起來，這種方式雖然有利于沙塵數值模式的快速預報需求，但需要引入諸多人為參數。更為重要的是，實際沙塵暴發生時，強對流流動對沙塵起動和垂向輸運的影響可能比近中性條件下的剪切湍流中湍流結構本身的影響還大。在這種情況下，粉塵分布和通量預測的參數化方案中如何引入溫度或穩定度等的影響還需要進一步深入研究。

4) 風沙運動數值模擬軟件的開發。已有的風沙運動和典型風成地貌數值模擬工作多數采用研究者自行編制in-house代碼，少數在計算流體動力學CFD軟件或開源代碼基礎上進行二次開發。但是，現有軟件如ANSYS Fluent主要的流場計算以RANS為主，雖然近年來也開發了相應LES部分，但所采用的湍流模型、壁模型等對高雷諾數大氣邊界層湍流的模擬仍需要檢驗。此外，雖然CFD軟件如ANSYS Fluent也提供了多相流求解器，但由于沒有考慮復雜的顆粒-可侵蝕床面的沖擊濺起過程、風場-顆粒-電場多場耦合相互作用等，均不能模擬實際風沙運動和風成地貌形成演化，或需要科研人員進行長期、復雜的二次開發。2018年，中國啟動了國家數值風洞(NNW)工程項目，目前已經發布多款自主研發功能先進的CFD軟件系統[153]。在NNW一期基礎上，開發針對風沙運動的專門軟件，既可以使風沙環境相關從業人員方便地利用軟件進行風沙災害的預報與防治，也可以一定程度上加快中國CFD工業軟件自主可控的進程。同時也將是NNW軍民融合的一個典型范例。