周晅毅，劉長卿，顧 明

（同濟(jì)大學(xué) 土木工程防災(zāi)國家重點實驗室，上海 200092）

風(fēng)致積雪運動的基本過程按雪粒離開地面的程度分為蠕移、躍移及懸移運動，許多學(xué)者認(rèn)為，躍移運動對雪顆粒的遷移有很大的貢獻(xiàn).研究風(fēng)致積雪運動的方法大致分為4類：理論分析、實地觀測、風(fēng)洞試驗和數(shù)值模擬.數(shù)值模擬方法因其成本低，周期短，可以方便地改變各種參數(shù)來研究其對結(jié)果的影響規(guī)律，越來越受到雪工程研究人員的關(guān)注.

近年來，計算流體力學(xué)方法在氣固兩相流的數(shù)值模擬中應(yīng)用越來越廣，受到不同行業(yè)學(xué)者的關(guān)注［1－2］.Uematsu，Naaim 和 Beyers等人都曾致力于風(fēng)雪運動的數(shù)值模擬研究.Uematsu等將雪顆粒躍移和懸移考慮在模型中，提出了躍移層風(fēng)雪流輸運質(zhì)量方程［3］；Naaim等構(gòu)建了雪粒堆積和風(fēng)蝕模型來描述雪表面和風(fēng)雪流之間的質(zhì)量交換［4］；Gordon等通過實地觀測，得到了雪顆粒形狀、粒徑等分布情況［5］；Beyers等人結(jié)合實測數(shù)據(jù)，討論了雪顆粒沉積和侵蝕對雪運動的影響［6］；X.Y.Zhou等采用兩相流理論模擬風(fēng)吹雪運動，在計算流體動力學(xué)軟件FLUENT的平臺上進(jìn)行了二次開發(fā)，對北京首都機場T3航站樓屋面的雪荷載分布規(guī)律進(jìn)行了分析，為結(jié)構(gòu)設(shè)計提供了依據(jù)［7］；周晅毅等在考察雪顆粒運動機理的基礎(chǔ)上，采用RANS方法對無建筑開闊場地的風(fēng)雪流進(jìn)行了模擬，并與經(jīng)典文獻(xiàn)中的實測結(jié)果進(jìn)行了對比［8］.以上學(xué)者的研究都是采用歐拉-歐拉方法，即把雪相和空氣相都假定為連續(xù)相處理；黃寧等采用歐拉－拉格朗日方法對風(fēng)沙（雪）運動進(jìn)行研究［9－11］，即把固相認(rèn)為是離散介質(zhì)，其方法對風(fēng)雪流的研究有著借鑒意義.

本文基于雪顆粒躍移質(zhì)量傳輸率等相關(guān)的經(jīng)驗公式，建立了雪面侵蝕質(zhì)量通量與粒子躍移數(shù)目之間的關(guān)系，結(jié)合拉格朗日方法，分析了研究粒子起跳速度、粒子直徑和摩擦速度等參數(shù)對運動軌跡的影響.在此基礎(chǔ)上，對雪顆粒質(zhì)量傳輸率進(jìn)行了計算，并與已有的經(jīng)驗公式進(jìn)行比對分析，驗證了本文提出的假設(shè)和計算方法的合理性和有效性.

1 理論分析

1.1 躍移質(zhì)量傳輸率定義

雪躍移質(zhì)量傳輸率，是指單位時間躍移層內(nèi)單位寬度通過的雪的質(zhì)量（kg·m－1·s－1），即

式中：hsal表示雪顆粒的躍移高度；usal（y）表示y高度處雪顆粒躍移水平平均速度；φsal（y）表示y高度處雪躍移的濃度，即躍移層中y高度處單位體積內(nèi)雪的質(zhì)量.usal（y），φsal（y）都隨高度變化，需將兩者乘積沿高度積分才能獲得Qsal.

1.2 雪躍移量經(jīng)驗公式

質(zhì)量傳輸率作為描述雪顆粒躍移的重要參數(shù)，很多學(xué)者都曾經(jīng)做過大量的研究工作，提出了相關(guān)的經(jīng)驗公式.其中，對于空曠環(huán)境下的雪顆粒傳輸率計算，Iversen與Pomeroy的經(jīng)驗公式應(yīng)用十分廣泛.

Iversen假定躍移雪顆粒的起動速度與（u＊－u＊t）相關(guān)，給出躍移雪質(zhì)量傳輸率的經(jīng)驗公式［12］為

式中：ρ為空氣密度；g為重力加速度，取為9.8m·s－2；u＊為壁面摩擦速度；u＊t為閾值摩擦速度.

Pomeroy結(jié)合平坦流域雪躍移的實測數(shù)據(jù)，給出雪躍移質(zhì)量傳輸率的經(jīng)驗公式［13］為

1.3 雪侵蝕量經(jīng)驗公式

雪是否發(fā)生侵蝕或沉積由近壁面的摩擦速度（風(fēng)速）決定.當(dāng)摩擦速度u＊超過閾值摩擦速度u＊t，壁面上的雪被風(fēng)吹起，即發(fā)生侵蝕.雪侵蝕量（即質(zhì)量通量，kg·m－2·s－1）可用下式表達(dá)［6］：

式中：Aero為常系數(shù)，取Aero＝7×10－4.

此外，還有其他學(xué)者，如Uemastu提出經(jīng)驗公式［3］為

式中：wf為粒子下落速度；usal為粒子躍移速度；usal＝2.8u＊t；hsal為粒子平均躍移高度，可用下式表達(dá)［13］：

2 躍移雪顆粒運動模型

2.1 雪顆粒受力分析

雪顆粒在氣流中受力較多，例如氣動升力、Saffman力、Magnuss力等，但并不是每一種力在運動分析時都需要考慮.同時，考慮到重力、阻力和浮力是雪顆粒運動的主要驅(qū)動力，后續(xù)的分析都近似忽略其他外力的作用.

式（7）～（9）中：mp為雪顆粒質(zhì)量；ρp為雪顆粒密度；Vp為雪顆粒體積；dp為雪顆粒直徑；ur為雪顆粒與空氣的相對速度；A為雪顆粒迎風(fēng)面橫截面面積；CD為阻力系數(shù)，采用經(jīng)驗公式［14］Re），Re為顆粒雷諾數(shù).

前文所述的粒子下落速度wf按照文獻(xiàn)經(jīng)驗公式［15］取為，將CD與Re的表達(dá)式代入其中，下落速度可變換為

式中：ν為空氣的運動黏性系數(shù).

可見，下落速度是雪顆粒密度及直徑的函數(shù)，如圖1所示，隨著雪顆粒粒徑和密度的增大，下落速度也隨之增大.

圖1 雪顆粒下落速度變化規(guī)律Fig.1 Settling velocity as a function of snow particle diameter

2.2 雪顆粒豎向起跳速度分布函數(shù)

假定在躍移運動中，雪顆粒全部為垂直起跳.當(dāng)顆粒以不同速度起跳時，記豎向初速度分布函數(shù)為f（v0），即單位時間從單位面積上起跳的顆粒中，以速度v0起跳的粒子數(shù)目與顆?？倲?shù)之比為f（v0）.采用Anderson and Hallet提出的如下指數(shù)函數(shù)形式的初速度分布［16］，見圖2.

觀察圖2，在同一摩擦速度下，以較小速度起跳的雪粒的概率比較大，并且隨著雪粒起跳速度的增大，相應(yīng)的概率逐漸減小.隨著摩擦速度增大，以較小速度起跳的雪粒的概率較低，以較大速度起跳的雪粒的概率變大.

2.3 雪顆粒運動方程

經(jīng)過上述受力分析，把粒子運動簡化為二維平面運動，通過坐標(biāo)與時間（x，y，t）表示，運動軌跡方程為

圖2 顆粒豎向起跳速度概率密度分布Fig.2 Probability distribution of liftoff velocity of particles

3 雪顆粒運動軌跡計算分析

3.1 躍移軌跡與雪顆粒粒徑關(guān)系

由上述運動方程（12），（13）可知，影響雪顆粒運動軌跡的參數(shù)因素包括粒子的直徑、密度和摩擦速度等因素.表1給出了我國天山地區(qū)積雪參數(shù)的大致范圍［17］.可見，雪顆粒的直徑大約為毫米量級，而雪顆粒的密度大約為40～300kg·m－3不等.

表1 我國新疆天山地區(qū)積雪參數(shù)Tab.1 Snow parameters of Tianshan region in Xinjiang Uygur Autonomous Region of China

給定雪顆粒運動的初始條件：假定粒子的初始坐標(biāo)為（0，0.005），初始速度只有豎向速度，水平速度為零，具體見表2.運動方程（12），（13）便可通過龍格－庫塔方法進(jìn)行求解，得到雪顆粒的運行軌跡.如無特殊說明，本文的計算按表2中取值.

在同等摩擦速度（u＊＝0.5m·s－1）條件下，粒子起跳速度分別為0.5m·s－1和2.0m·s－1時，分別繪出粒子直徑為0.30，0.45，0.60與0.75mm 的顆粒運行軌跡.從圖3可得，在相同的起跳速度時，隨著粒子直徑的增大，粒子的躍移高度h逐漸變大；并且，由圖3a與3b對比可知，隨著粒子初始速度增大，其所能達(dá)到的最大躍移高度也隨之增大.例如，當(dāng)起跳速度由0.5m·s－1增加到2.0m·s－1時，粒徑為0.45mm的顆粒的最大躍移高度由0.014m增加到0.053m.另外，觀察圖3，當(dāng)顆粒起跳速度不同時，水平躍移距離l的變化規(guī)律是相反的：圖3a中，粒子起跳速度較小時，隨著粒子直徑的增大，粒子的躍移長度逐漸變??；而圖3b中，粒子起跳速度較大時，隨著粒子直徑的增大，粒子的躍移長度逐漸變大.阻力占豎向力中的比例相對較小，故躍移高度較大.因此，在物理意義上，圖3中的現(xiàn)象便得到了合理的解釋.而且，這與文獻(xiàn)［18］的計算結(jié)論是類似的，在一定程度上也驗證了本文計算結(jié)果的準(zhǔn)確性.

表2 數(shù)值計算的相關(guān)參數(shù)Tab.2 Numerical calculation related parameters

圖3 雪顆粒躍移軌跡隨直徑的變化規(guī)律Fig.3 Saltation trajectory as a function of snow particle diameter

圖4 雪顆粒豎向受力與重力之比隨躍移高度的變化曲線Fig.4 Variation curves of the ratio of vertical force to gravity with saltation height

圖4中箭頭↑、↓分別表示粒子處于上升和下降階段.通過前文的受力分析可知，粒子在上升過程中，重力和空氣阻力豎向分量方向均向下.在粒子達(dá)到最大躍移高度處，粒子只有水平方向的速度，空氣阻力的豎向分量為零，雪顆粒在豎直方向只受到重力的作用（浮力相比較?。?，此時豎向受力與重力之比為1.不難發(fā)現(xiàn)，在上升過程中，當(dāng)粒子以相同的速度起跳時，由于空氣阻力起減速作用，較小粒徑的雪顆粒受到的空氣阻力占豎向力中的比例較大，因而最大躍移高度較小；而較大粒徑雪顆粒受到的空氣

圖5所示為雪顆粒水平受力與重力之比隨躍移高度的變化規(guī)律.圖中橫坐標(biāo)軸的正負(fù)號表示水平受力的方向，“＋”表示起加速作用，“－”表示起減速作用.圖5a中，粒子的粒徑越小，其水平方向受力與重力之比越大，且圖5a中水平力幾乎都為正值，即當(dāng)粒子起跳速度較小時，在運動過程中，粒子速度一直小于風(fēng)速，水平作用力起加速作用，故而小粒徑的雪顆粒躍移長度較大；圖5b中，下降階段，水平力為負(fù)值，即當(dāng)粒子起跳速度較大時，在下降過程中，粒子速度一直大于風(fēng)速，水平作用力起減速作用.且粒徑越小，受到的阻力越大，故小粒徑的雪顆粒的躍移距離可能會小一些.當(dāng)然，由于本文對雪顆粒做了球體的假設(shè)，這樣的結(jié)論尚需進(jìn)一步的驗證.

3.2 躍移軌跡與摩擦速度關(guān)系

圖5 雪顆粒水平受力與重力之比隨躍移高度的變化曲線Fig.5 Variation curves of the ratio of horizontal force to gravity with saltation height

不同粒徑的顆粒最大水平躍移距離（即躍移長度）隨摩擦速度的變化規(guī)律如圖6所示.與前文3.1類似，粒子起跳速度不同時，水平躍移距離的變化規(guī)律也不同.當(dāng)顆粒起跳速度較小時，如圖6a中，v0＝0.5m·s－1，隨著摩擦速度的增大，顆粒的水平躍移距離先增大后減?。磺译S著粒子直徑的增大，水平躍移距離隨之減小.當(dāng)顆粒起跳速度較大時，如圖6b中，v0＝2.0m·s－1，隨著摩擦速度的增大，不同粒徑顆粒的水平躍移距離都隨之增大，但增大的幅度不同.較小粒徑的顆粒水平躍移距離隨著摩擦速度的增大幅度相比粒徑較大的顆粒增幅小一些.例如對于粒徑為0.15mm的顆粒，當(dāng)摩擦速度從0.2m·s－1增加到1.0m·s－1時，水平躍移距離的增大幅度為0.21mm；而對于粒徑為0.75mm的顆粒，當(dāng)摩擦速度從0.2m·s－1增加到1.0m·s－1時，水平躍移距離的增大幅度為0.44mm，約為前者的2.1倍.

上述現(xiàn)象的具體原因，可由圖7中所描述的變化規(guī)律予以解釋.由于篇幅所限，數(shù)據(jù)較多，且變化規(guī)律類似，圖中只列舉了一些代表性的摩擦速度數(shù)值下的變化規(guī)律（顆粒直徑為0.45mm）.例如圖7a中，當(dāng)粒子起跳速度較小時，在上升階段，此時粒子速度小于風(fēng)速，水平作用力起加速作用.當(dāng)摩擦速度u＊小于0.5m·s－1時，隨著摩擦速度的增大，粒子在水平方向的受力與重力之比隨之變大，在u＊＝0.5m·s－1時，水平受力與重力之比最大；當(dāng)摩擦速度u＊大于0.5m·s－1時，隨著摩擦速度的增大，粒子在水平方向的受力與重力之比隨之變小.這就解釋了圖6a中躍移長度先增大后減小，且在0.5m·s－1時其躍移長度最大的原因.圖7b中，起跳速度為2.0m·s－1與起跳速度為0.5m·s－1的躍移規(guī)律明顯不同，在上升階段，隨著摩擦速度的增大，粒子在水平方向的受力與重力之比隨之變大.因而當(dāng)豎向起跳速度較大時，摩擦速度增大，水平躍移長度增加.

圖6 水平躍移距離隨摩擦速度的變化曲線Fig.6 Variation curves of saltation length with friction velocity

由圖8可知，在粒徑不變的情況下，隨著摩擦速度的增大，本文所計算出的躍移高度并未發(fā)生明顯的變化，基本維持在某一固定值附近；而經(jīng)驗公式中，躍移高度與摩擦速度呈拋物線關(guān)系.可以看出，本文計算出的粒子躍移高度介于平均躍移高度hsal最大值與最小值的范圍之內(nèi)，并且隨著起跳速度的增大而增大.一方面說明了傳統(tǒng)經(jīng)驗公式的合理性，另一方面也進(jìn)一步驗證了本文計算的粒子躍移軌跡的有效性.

3.3 躍移軌跡與雪顆粒密度關(guān)系

假定雪顆粒粒徑單一不變，均為0.45mm，顆粒躍移軌跡與雪顆粒密度的變化關(guān)系如圖9所示.可見，隨著雪顆粒密度的增大，粒子的躍移高度隨之增大.

圖7 雪顆粒水平受力與重力之比隨躍移高度的變化曲線（粒徑0.45mm）Fig.7 Variation curves of the ratio of horizontal force to gravity with saltation height（d＝0.45mm）

對比圖4和圖10，通過前文的受力分析可知，粒子在上升過程中，空氣阻力豎向分量向下，在粒子達(dá)到最大躍移高度處，粒子只有水平方向的速度，空氣阻力的豎向分量為零，雪顆粒在豎直方向只受到重力的作用（浮力相比較?。?此時，豎向加速度近似等于重力加速度g.在上升過程中，粒子做減速運動，加速度越大，上升的高度越低.由圖10可知，豎向起跳速度相同時，較小密度的雪顆粒豎向加速度較大，因而最大躍移高度較低；而較大密度的顆粒的加速度相對較小，故躍移高度較大.這與前文3.1的變化規(guī)律類似.

4 雪顆粒躍移質(zhì)量傳輸率計算分析

4.1 本文計算模型

圖10 雪顆粒豎向加速度隨密度的變化曲線Fig.10 Variation curves of vertical acceleration with snow particle density

本文試圖建立雪面侵蝕質(zhì)量通量與粒子躍移數(shù)目之間的關(guān)系，假設(shè)單位時間單位面積上起跳的粒子數(shù)為s，mp為單個雪顆粒的質(zhì)量，則mps為單位時間單位面積上起跳粒子的質(zhì)量，物理意義上等價于壁面侵蝕量qero，表達(dá)形式見公式（4）與（5）.即

下文的計算模型參考了文獻(xiàn)［9］（具體計算過程與之類似），記豎向初速度分布函數(shù)為f（v0），則單位時間、單位面積內(nèi)從地面以初速度v0起跳的雪粒質(zhì)量為

進(jìn)一步可知，以初速度v0起跳的雪粒，在上升（↑）和下降（↓）階段，高度y處的質(zhì)量濃度表示為

通過式（16）與（17），可以得到在高度y處單位時間內(nèi)通過單位面積以初速度v0起跳的雪粒質(zhì)量為

將式（18）對雪顆粒所有起跳的速度積分，并且從起跳位置高度至最大躍移高度處積分，便可得到以下顆粒躍移質(zhì)量傳輸率的計算公式：

由公式（2）～（5）可知，qero有3種表達(dá)形式，即通過Iversen，Pomeroy和Beyers三者的經(jīng)驗公式變換，可分別記為qI，qP和qB，對應(yīng)的躍移質(zhì)量傳輸率記為QI，QP和QB.

4.2 與經(jīng)驗公式對比分析

將公式（19）與躍移質(zhì)量傳輸率的經(jīng)驗公式（2），（3）進(jìn)行對比，從而驗證本文計算結(jié)果的合理性，具體見圖11和圖12給出的躍移質(zhì)量傳輸率的變化規(guī)律.不難發(fā)現(xiàn)：

圖11 雪顆粒躍移質(zhì)量傳輸率隨（閾值）摩擦速度變化規(guī)律（粒徑0.30mm）Fig.11 Total snow mass flux as a function of friction velocity／threshold velocity（d＝0.30mm）

（1）對比圖11a和圖12a，在摩擦速度較小時（如u＊t＜0.4m·s－1），本文的計算結(jié)果與經(jīng)驗公式吻合較好；在摩擦速度較大時，即粒子處于高風(fēng)速情況下，Pomeroy與Iversen所計算出的結(jié)果相差很大，這是因為前者的計算結(jié)果與摩擦速度的平方成正比，而后者正比于摩擦速度的立方.此時計算結(jié)果QI的數(shù)值與經(jīng)驗公式最為接近.

（2）觀察圖11b和圖12b，可知隨著閾值摩擦速度的增大，Pomeroy公式計算的躍移雪質(zhì)量傳輸率近似呈直線增加；Iversen公式的躍移雪質(zhì)量傳輸率近似呈直線減小.Pomeroy與Iversen曲線趨勢相反，因為后者與0.1m·s－1關(guān)系是負(fù)相關(guān)，故Iversen公式計算值減小；而前者與0.1m·s－1關(guān)系既有正相關(guān)，還有立方的負(fù)相關(guān)，總體呈現(xiàn)正相關(guān)（u＊t＜1.0m·s－1）.可見，QI的數(shù)值基本介于兩經(jīng)驗公式結(jié)果之間，趨勢與Iversen曲線趨勢相同.

（3）對比圖11和圖12，Pomeroy與Iversen所計算出的結(jié)果并未變化，因其是統(tǒng)計意義的平均值；但是本文的計算模型，隨著粒徑的增大，躍移雪質(zhì)量傳輸率也隨之增大.

圖12 雪顆粒躍移質(zhì)量傳輸率隨（閾值）摩擦速度變化規(guī)律（粒徑0.45mm）Fig.12 Total snow mass flux as a function of friction velocity／threshold velocity（d＝0.45mm）

5 結(jié)論

本文通過變化不同的顆粒運動參數(shù)，分析其與雪顆粒運動軌跡之間的變化規(guī)律.建立了雪面侵蝕質(zhì)量通量與粒子躍移數(shù)目之間的關(guān)系，研究了粒子直徑、摩擦速度等參數(shù)對質(zhì)量傳輸率的影響，并與已有的經(jīng)驗公式比對分析，結(jié)論如下：

（1）隨著雪顆粒直徑、雪粒密度或者摩擦速度的增大，粒子的躍移高度也隨之增大；粒子起跳速度不同時，水平躍移距離的變化規(guī)律也不同.

（2）本文提出的計算模型與空曠環(huán)境下的雪顆粒傳輸率經(jīng)驗公式計算結(jié)果基本吻合，證明了本文假定和方法的合理性，其中部分結(jié)果有待與實測或?qū)嶒炦M(jìn)行對比檢驗.

［1］ Tominaga Y，Mochida A.CFD prediction of flowfield and snowdrift around a building complex in a snowy region［J］.Journal of Wind Engineering and Industrial Aerodynamics，1999，81（1／3）：273.

［2］ ZHU Zhiwen，LIU Zhenqing.CFD prediction of local scour hole around bridge piers［J］.Journal of Central South University of Technology，2012，19：273.

［3］ Uematsu T，Nakata T，Takeuchi K，et al.Three-dimensional numerical simulation of snowdrift［J］.Cold Regions Science and Technology，1991，20：65.

［4］ Naaim M，Naaim-Bouvet F，Martinez H.Numerical simulation of drifting snow：erosion and deposition models［J］.Annals of Glaciology，1998，26：191.

［5］ Gordon M，Taylor P A.Measurements of blowing snow，part I：particle shape，size distribution，velocity，and number flux at Churchill，Manitoba，Canada［J］.Cold Regions Science and Technology，2009，55（1）：63.

［6］ Beyers J H M，Sundsb P A，Harms T M.Numerical simulation of three-dimensional，transient snow drifting around a cube［J］.Journal of Wind Engineering and Industrial Aerodynamics，2004，92（9）：725.

［7］ Zhou X Y，Li X F.Simulation of snow drifting on roof surface of terminal building of an airport［J］.Disaster Advances，2010，3（1）：42.

［8］ 周晅毅，李雪峰，顧明.風(fēng)致積雪運動數(shù)值模擬的兩方程模型方法［J］.空氣動力學(xué)學(xué)報，2012，30（5）：640.

ZHOU Xuanyi，LI Xuefeng，GU Ming.Two-equation model for numerical simulation on snow drifting［J］.Acta Aerodynamica Sinica，2012，30（5）：640.

［9］ 黃寧.沙粒帶電及風(fēng)沙電場對風(fēng)沙躍移運動影響的研究［D］.蘭州：蘭州大學(xué)土木工程與力學(xué)學(xué)院，2002.

HUANG Ning.Electrification in wind-blown sand flux and its influence to wind-blown sand saltation［D］.Lanzhou：School of Civil Engineering and Mechanics of Lanzhou University，2002.

［10］ Huang N，Sang J，Han K.A numerical simulation of the effects of snow particle shapes on blowing snow development［J］.Journal of Geophysical Reserch，2011，116，D22206，doi：10.1029／2011JD016657.

［11］ LüX H，Huang N，Tong D.Wind tunnel experiments on natural snow drift［J］.Science China Technological Sciences，2012，55：927.

［12］ Iversen J D. Drifting-snow similtude-transport-rate and roughness modeling［J］.Journal of Glaciology，1980，26（94）：393.

［13］ Pomeroy J W，Gray D M.Saltation of snow［J］.Water Resources Research，1990，36（7）：1583.

［14］ Carrier C F.On slow viscous flow［R］.Providence：Brown University，1953.

［15］ 李雪峰.風(fēng)致建筑屋蓋表面及其周邊積雪分布研究［D］.上海：同濟(jì)大學(xué)土木工程學(xué)院，2011.

LI Xuefeng.Research on snowdrifting on building roof and around building.［D］.Shanghai：College of Civil Engineering of Tongji University，2011.

［16］ Anderson R S，Hallet B.Sediment transport by wind：toward a general model［J］.Geological Society of America，1986，97：523.

［17］ 朱光耀.公路風(fēng)吹雪雪害形成機理與防治［M］.哈爾濱：黑龍江人民出版社，2007.

ZHU Guangyao.The formation mechanism of Snowdrift disaster on highway and its control［M］.Haerbin：Heilongjiang Publishing Group，2007.

［18］ 李歲虹.穩(wěn)態(tài)風(fēng)沙運動的物理力學(xué)特性研究［D］.蘭州：蘭州大學(xué)土木工程與力學(xué)學(xué)院，2005.

LI Suihong.Physical and mechanical characteristics of the steady state of sand movement［D］.Lanzhou：School of Civil Engineering and Mechanics of Lanzhou University，2005.