硅酸鉀材料加固土遺址耐風蝕顆粒元模擬分析

張曉平 ，王思敬，李 黎，王彥兵

（1. 敦煌研究院 古代壁畫保護國家文物局重點科研基地，甘肅 敦煌 736200；2. 中國科學院地質與地球物理研究所 中國科學院工程地質力學重點實驗室，北京 100029；3. 中國文化遺產研究院，北京 100871）

1 引 言

我國是世界上文明起源最早的國家之一，保存在西北古絲綢之路沿線數千公里的戰(zhàn)國、漢代、明代長城及新疆交河故城遺址等土遺址文物都是我國珍貴文化遺產的重要組成部分，具有非常重要的歷史、科學和藝術價值，其中交河故城則是西北干旱地區(qū)土遺址的代表?，F存土建筑20余萬平方米，是研究中西文化交流，中國西部和中亞地區(qū)歷史以及城市建筑史、宗教史、美術史、民族史等多種學科不可多得的珍貴實物資料，具有極大的保存價值。但由于遺址區(qū)常年大風、沙暴，集中式強降雨等，造成遺址嚴重風化風蝕、雨水沖刷侵蝕、墻體開裂、坍塌，急需進行搶修，目前土遺址的保護加固工作僅步入初始階段[1]。其中加固所需的 PS(potassium silicate)材料[1]及其加固機制研究是最基礎的研究內容。

土遺址 PS加固機制研究是土遺址加固的理論部分，在進行大量土遺址搶修加固應用[2－6]的同時，也已經開展了相關的研究探討[7－8]。通過PS加固遺址土過程的熱分析、加固前后物相分析、可溶鹽反應、X射線光電子能譜分析、熱重分析等得出，交河故城遺址土加固前后物相沒有發(fā)生變化，PS材料固化的主要物理作用包括：一方面填塞遺址土顆粒間的空隙；另一方面增加土顆粒間的膠結強度,形成較致密的網狀結構。此外，PS在遺址土顆粒之間的化學鍵合作用、PS與遺址土礦物的羥基發(fā)生縮合反應，把遺址土進一步結合起來，形成了整體的聯(lián)結[8]。遺址土風蝕試驗結果表明：風蝕主要發(fā)生在挾沙風條件下，未加固樣的風蝕態(tài)勢表現為隨挾沙風風速增大，風蝕強度增加；加固樣品的抗風蝕強度（抗風蝕強度指試樣抵抗風蝕的能力。在相同風蝕條件下，抗風蝕強度與風蝕量成反比例關系，即試樣風蝕量的減小，表征抗風蝕強度的提高）明顯高于未加固樣品，其抗風蝕強度是原樣的 2.9～38倍[9]。

以上研究結果初步揭示了 PS加固遺址土的加固機制和加固前后的風蝕特性，給顆粒元模擬提出了建模概念。而接下來的研究方向應該是針對大量不同賦存條件、不同特性遺址土分別展開專項研究，驗證并深入研究 PS加固機制及加固后各耐環(huán)境因素。將土力學研究領域中的顆粒元模擬應用到土遺址保護加固中，為 PS加固機制分析乃至現場加固工程提供理論依據。顆粒元的模擬應用將增加 PS加固機制研究的途徑，增加試驗的針對性和施工的預見性，減小試驗取樣和施工反復對遺址的不利影響。但顆粒元在遺址土中要得到有效運用，首要工作是進行顆粒元與遺址土室內試驗結果的擬合與分析，這是本文擬進行的首要研究內容，本文選取具有代表性的交河故城遺址土作為研究對象，進行顆粒元模擬研究工作，最終對加固前后遺址土進行風蝕模擬，研究加固前后的風蝕特征。

顆粒元方法最早由Cundall提出[10]，采用顆粒集合模擬巖土體在近 10 年受到廣泛關注。該模型不需要輸入巖土體的本構模型，只需要指定顆粒間的接觸剛度和接觸強度參數，顆粒間的計算只是采用簡單的牛頓運動定律，便能模擬出巖土體的非線性力學特征和微觀破壞機制，與連續(xù)網格的有限元計算相比具有獨特的優(yōu)勢。通過在顆粒表面增加壓力，能模擬風壓、液壓的作用。張曉平等[11－12]采用顆粒元模擬含了軟弱夾層土樣的變形破壞機制，并嘗試模擬滑坡過程[13－15]及抗滑樁加固[16]；吳順川[17]利用顆粒元模擬了壓力注漿；劉洋[18]利用顆粒元流-固耦合分析模型模擬了飽和砂土的液化現象；張剛[19]利用顆粒元模擬了管涌現象等。遺址土多為人工夯土、版筑土等，不同于自然環(huán)境的土體，顆粒元在土遺址 PS加固模擬中的運用，必將拓寬其應用范圍，也為土遺址加固機制研究及各耐風蝕環(huán)境影響分析提供一種新的思路和方法。

2 遺址土的物理力學性質

2.1 遺址土的物理特性

交河故城本體土的類型有4種：生土、夯土、垛泥和土坯。主要為粉土和粉質黏土。顆粒分析結果表明，粉粒含量普遍在80%左右，垛泥的粒度成分變化極大，級配不佳，而生土和夯土的級配良好。李最雄等[1]曾統(tǒng)計了故城遺址建筑材料的物理性質指標；2008年和法國等[20]又進行了大量的室內常規(guī)試驗。比較兩作者所給出的物性指標數據基本吻合，現僅給出文獻[20]的最新試驗資料，如表 1所示。

2.2 PS加固前后遺址土的力學性質

將原狀土用5 %濃度PS溶液噴灑3次，每次噴灑時間間隔為24 h，全部噴灑完畢后在室內養(yǎng)護至自然風干，測試其單軸抗壓和抗拉強度，選取典型的土樣 JH-029、JH-030、JH-011和 JH–012，其測試結果如表2?？梢钥闯?，經過PS材料加固后，生土和垛泥的強度得到了明顯提高，提高幅度120%～180%，生土的強度要高于垛泥[20]。

3 PS加固前后遺址土的力學性質顆粒元(PFC)模擬

本研究首先采用顆粒元對生土加固前后力學性質進行擬合，然后模擬加固前后生土的風蝕過程。

3.1 顆粒元模型簡介

遺址土模擬的基礎是顆粒元中的連接顆粒模型(bonded particle model，簡稱BPM)，因此，首先對BPM模型及其參數進行介紹。

本次模擬采用 Itasca公司開發(fā)的顆粒元程序PFC2D，該程序自1995年首次發(fā)布以來，在許多研究項目中得到了應用。PFC2D通過有限厚度的圓盤集合或單層（垂直于平面方向）球來描述土體或巖石，顆粒間通過黏結和摩擦連接，連接為一定剛度和顆粒間的有限區(qū)域。當指定的連接強度超過施加在連接上的應力（如拉力、剪力或由于顆粒旋轉而形成的彎矩），連接發(fā)生破裂形成破裂面，如裂紋。在這些過程中，PFC不需要任何本構關系來描述屈服后的響應和斷裂韌度來控制破裂行為，而只需要顆粒運動定律、簡單的顆粒變形公式和顆粒連接斷裂的準則，有關PFC的詳細情況參見文獻[21－23]。平行連接模型是PFC中連接模型(BPM)的一種，近似為兩個顆粒在接觸點處膠結作用，能承受拉力、剪力和彎矩，是模擬巖土材料的理想模型。確定平行連接模型的主要參數如下：

描述礦物顆粒參數有 { Rmin, Rmax/Rmin, ρ,Ec，(kn/ks), μ}；描述膠結作用的平行連接參數有

其中：Rmin和 Rmax/ Rmin描述顆粒半徑，在模型中滿足以Rmax和Rmin為上下限的均勻分布；ρ為顆粒的密度；Ec和分別為顆粒和平行連接的彈性模量；和分別為顆粒和平行連接的法向剛度與切向剛度的比率；為設置平行連接半徑R的半徑系數，為顆粒半徑；μ為顆粒間的摩擦系數；和分別為平行連接的抗拉強度和抗剪強度。在給定彈性模量和法向剛度與切向剛度的比率后，通過式（1）、（2）來計算顆?；蚱叫羞B接的法向和切向剛度。

對法向剛度來講，可以想象成如圖1所示的彈性梁。

圖1 顆粒間接觸和平行連接的等效連續(xù)Fig.1 Equivalent continuum material of parallel bond particle system

平行連接顆粒間相對運動會形成軸向力 T、切向力V和彎矩M。連接材料能承擔的最大法向應力和切向應力由式（3）、（4）表達。

A和I是平行連接截面（如圖1(b)）的面積和轉動慣量，T為正值表示抗拉。當或超過相應的連接強度連接就會破壞，在兩個顆粒間形成拉裂紋和剪裂紋。

PFC的顆粒參數和連接參數可以通過單軸壓縮試驗和巴西劈裂試驗來擬合。通過單軸壓縮試驗可得到徑向應變、軸向應變和峰值強度?？估瓘姸韧ㄟ^巴西劈裂計算得到。PFC試樣的生成和加載方法在文獻[21－23]中有詳細論述。由于PFC中引入了隨機模型，因而顆粒半徑和強度參數也都服從統(tǒng)計的隨機模型，雖然同樣的命令程序，但沒有兩次PFC的計算結果會完全相同。因此，在參數研究中，每種參數計算10次，求平均值。這樣就會增加很大的計算工作量，但如同實際的巖石力學試驗一樣，可避免由于個別樣品的差異而得到不相干的結論。

3.2 顆粒元模型

顆粒元模擬土體，顆粒粒徑不能與土體粒徑一一對應，但可以按照級配曲線進行近似粒徑比例模擬。實際上交河故城粉粒（粒徑0.005～0.075 mm）含量普遍在80%左右，因此，出于減小計算量，使得計算具有實用性，故不考慮顆粒級配，而擬采用粒徑0.005～0.075 mm范圍內均勻分布顆粒模型來模擬遺址土。實際顆粒元計算中發(fā)現，粒徑0.005～0.075 mm范圍內均勻分布顆粒，計算50 mm ×50 mm的原狀土壓縮試樣或 70.7 mm×70.7 mm的重塑土樣，仍然存在顆粒數巨大計算機無法執(zhí)行的情況。故考慮將土體粒徑擴大 20倍進行模擬，即在 PFC模型中，顆粒直徑在 0.1～1.5 mm 范圍內均勻分布。如圖2所示，模型尺寸50 mm×50 mm，顆粒直徑0.1～1.5 mm范圍均勻分布，模型包含4 177個顆粒。

PFC顆粒模型中，指定顆粒密度是指顆粒本身的密度。由于顆粒之間有孔隙，計算模型的密度比顆粒密度要小。因此，在擬合實際土樣密度中，先將實際土樣密度乘以計算模型體積（二維模型中為單位厚度），計算出總質量；然后用總質量除以模型顆粒的實際體積，能得到應賦顆粒密度值。這一段通過獨立編寫擴展Fish語言來實現，即只需輸入實際土樣的密度，便能計算出PFC顆粒密度值，并自動賦值。

3.3 生土PS加固前后抗壓、抗拉強度模擬

生土壓縮試樣計算模型尺寸為50 mm×50 mm，巴西劈裂試樣為φ50 mm的圓盤。生土取天然密度平均值，參照表1生土密度范圍為1.51～1.77 g/cm3，故近似取生土PFC模型密度為1.64 g/cm3。采用Fish計算得到對應顆粒密度為1.833 g/cm3。

通過反復試算逐步逼近生土的抗壓和抗拉強度參數，最后取表3所示參數作為生土加固前、后對應的顆粒元模型參數。對比表中的顆粒參數可以看出，PS加固后生土的強度參數比加固前要大。也就是說，PS加固土體的作用是通過增加顆粒間平行連接強度來實現的。如3.1節(jié)所述，出于統(tǒng)計考慮，每種試樣抗壓和抗拉試驗分別計算10次（見表4），取10次計算結果的平均值列于表 5。模擬得到加固前后抗壓和抗拉強度參數與試驗結果非常接近。

4 PS加固前后遺址土風蝕顆粒元(PFC)模擬

4.1 風蝕模型

遺址土風蝕試驗結果表明，風蝕主要發(fā)生在挾沙風條件下。目前對風沙運動的力學機制研究，在沙漠化過程研究中較多，且集中在風沙躍移運動發(fā)展過程。大量的研究者從離散的角度對其進行了數值模擬[24－28]，并對風沙的啟動機制和顆粒間的相互作用進行了相關研究，這些研究吻合于自然風沙躍移運動過程的基本特征。本文的研究內容主要集中在分析加固前、后遺址土的抗風蝕能力變化，因此，將風沙的運動簡化為一部分以一定初速度撞向土體表面的顆粒。撞擊的結果使得土體表面部分顆粒間連接斷開，從土體中剝落，形成風蝕（見圖3）。圖中左側為土體，右側豎線處為隨機產生不同粒徑的挾沙風顆粒，以一定的水平初速度（箭頭所示）撞向土體。

表3 模擬加固前生土顆粒參數Table 3 Particle parameters of raw soil before reinforcement

表4 生土加固前、后抗壓強度和抗拉強度模擬結果統(tǒng)計表Table 4 Numerical simulation of compressive strength and tensile strength before and after reinforcement of raw soil

表5 生土試驗結果和模擬結果對比表Table 5 Comparison between experimental and numerical outputs

圖3 挾沙風模擬示意圖Fig.3 Illustration of sand-driving wind flow simulation

4.2 風蝕模擬計算及分析

文獻[9]中采用顆粒速度為 6、10、15、20 m/s的挾沙風進行了風洞吹蝕試驗。發(fā)現20 m/s挾沙風條件下，加固前后土體的抗風蝕能力差別顯著。為了便于比較，本顆粒元模擬中，也選取20 m/s的顆粒速度模擬挾沙風。挾沙風顆粒從圖3所示右側豎線位置隨機產生，每個顆粒的初始位置都在該豎線上，但高低位置隨機分布。挾沙風顆粒直徑服從0.1～1.5 mm之間的均勻分布，與土體顆粒分布保持一致。挾沙風顆粒密度及彈性參數與生土相同：ρ=1.833 g/cm3，Ec=0.08 GPa，kn/ks=1，μ=0.5。計算過程設置每循環(huán)5 000個時步，產生一個新的挾沙風顆粒。因此，挾沙風顆粒數與循環(huán)步數成正比例的關系，可以用挾沙風顆粒數來代表吹蝕時間的長短。

圖4為生土加固前、后不同吹蝕時間（用不同挾沙風顆粒數表示）條件下土樣顆粒元模型的吹蝕圖片。加固前，在挾沙風顆粒為300和600個時，試樣右側面已經有較明顯的風蝕；挾沙風顆粒為900個時，風蝕變得非常顯著；挾沙風顆粒繼續(xù)增加到1 200個和1 500個時，風蝕面繼續(xù)加大、風蝕掏蝕深度加深。

加固后，在挾沙風顆粒數為300、600、900個時，試樣右側表面風蝕還非常輕微；到挾沙風顆粒達到1 200個以后，表面才出現較顯著的風蝕現象。

圖4 生土加固前、后不同吹蝕時間（即不同挾沙風顆粒數）條件下土樣顆粒元模型的吹蝕Fig.4 Snapshots of erosion under different times (i. e. different numbers of sand-driving wind particles)before and after reinforcement of raw soil

為了定量比較加固前后的風蝕量，即風蝕掉的顆粒面積，圖5列出生土加固前、后風蝕掉的顆粒面積與挾沙風顆粒數的對比曲線?？梢钥闯?，隨著挾沙風顆粒數的增加，即風蝕時間的增長，加固前土樣的風蝕量一直呈較快的增長趨勢，比加固后的顯著偏大。加固后的土樣隨挾沙風顆粒數的增加，增幅緩慢而且趨于穩(wěn)定?？梢娂庸毯笊镣翗涌癸L蝕能力得到了顯著的增強。與風洞試驗結果統(tǒng)計的規(guī)律一致[9]。即加固樣品的抗風蝕能力明顯高于未加固樣品，其抗風蝕強度是原樣的2.9～38倍[9]。

圖5 生土加固前、后風蝕量與挾沙風顆粒數對比曲線Fig.5 Erosion mass of soil vs. erosion time before and after reinforcement of raw soil

5 結 論

（1）顆粒元模擬通過增加顆粒間平行連接強度，實現對 PS加固的模擬。模擬得到土樣加固前后抗壓和抗拉強度參數與試驗結果接近。

（2）隨機生成挾沙風顆粒，以一定的速度撞向土體，模擬挾沙風的吹蝕作用。挾沙風顆粒數與循環(huán)步數成正比例，可以用挾沙風顆粒數來代表吹蝕時間的長短。而挾沙風顆粒的速度則代表挾沙風風速。

（3）風蝕模擬結果表明，在20 m/s的挾沙風吹蝕作用下，風蝕程度隨吹蝕時間的增加而增大，未加固土樣的風蝕程度增幅度遠大于加固土樣；同樣吹蝕時間條件下，加固土樣的抗風蝕強度明顯高于未加固土樣。PS加固后生土的抗風蝕能力得到顯著增強，與風洞試驗結果的統(tǒng)計規(guī)律一致。

（4）本研究擬合的顆粒元模型，可以應用于遺址土加固的進一步研究，如加固機制研究及耐風蝕、雨蝕、凍融等諸環(huán)境影響分析研究。

[1]李最雄. 絲綢之路古遺址保護[M]. 北京: 科學出版社,2003.

[2]李最雄, 王旭東, 張志軍, 等. 秦傭坑土遺址的加固試驗[J]. 敦煌研究, 1998, (4): 151－158.LI Zui-xiong, WANG Xu-dong, ZHANG Zhi-jun, et al.Consolidation tests of Qinyong site[J]. Dunhuang Research, 1998, (4): 151－158.

[3]李最雄, 王旭東, 郝利民. 室內土建筑遺址的加固試驗—半坡土建筑遺址的加固試驗[J]. 敦煌研究, 1998, (4):144－150.LI Zui-xiong, WANG Xu-dong, HAO Li-min. Research on the conservation of indoor earth structure sites—Experimentation of chemical consolidation on ancient earth—Structure Site of Ban Po and Qin Yong Keng[J]. Dunhuang Research, 1998, (4): 144－150.

[4]孫滿麗. 吐魯番交河故城保護加固研究[D]. 蘭州: 蘭州大學, 2006.

[5]王旭東, 張魯, 李最雄, 等. 銀川西夏 3號陵的現狀及保護加固研究[J]. 敦煌研究, 2002, (4): 64－72.WANG Xu-dong, ZHANG Lu, LI Zui-xiong, et al. The study of existing condition and consolidation project of NO.3 tomb of the Western Xia Mausoleums[J].Dunhuang Research, 2002, (4): 64－72.

[6]蘇伯民, 李最雄, 胡之德. PS與土遺址作用機理的初步探討[J]. 敦煌研究, 2000, 63(1): 3－35.SU Bo-min, LI Zui-xiong, HU Zhi-de. Study on consolidating mechanism between PS and earth relics[J].Dunhuang Research, 2000, 63(1): 3－31.

[7]王旭東等. 國家十一五科技支撐計劃課題: 大遺址保護關鍵技術研究與開發(fā)—土遺址保護關鍵技術研究報告[R]. 敦煌: 敦煌研究院, 2009.

[8]LI Li, WANG Si-jing, SHAO Ming-shen, et al. Impact of environment factors on consolidating earthen architecture sites with PS[C]//Proceedings of International Symposium Conservation of Ancient Sites &ISRM-Sponsored Regional Symposium. Beijing: Science Press, 2010.

[9]趙海英, 李最雄, 汪稔, 等. PS 材料加固土遺址風蝕試驗研究[J]. 巖土力學, 2008, 29(2): 392－396.ZHAO Hai-ying, LI Zui-xiong, WANG Ren, et al. Wind erosion experiment of ancient earthen site consolidated by PS material[J]. Rock and Soil Mechanics, 2008, 29(2):392－396.

[10]CUNDALL P A. A computer model for simulating progressive large scale movements in blocky rock system[C]//Proceedings of the Symposium of the International Society of Rock Mechanics. Nancy, France:[s. n.], 1971.

[11]張曉平, 吳順川, 張志增, 等. 含軟弱夾層土樣變形破壞過程細觀數值模擬及分析[J]. 巖土力學, 2008, 29(5):1200－1204.ZHANG Xiao-ping, WU Shun-chuan, ZHANG Zhi-zeng,et al. Numerical simulation and analysis of failure process of soil with weak intercalated layer[J]. Rock and Soil Mechanics, 2008, 29(5): 1200－1204.

[12]張曉平, 吳順川, 張兵, 等. 軟弱夾層幾何參數對試樣力學行為影響顆粒元模擬研究[J]. 工程地質學報, 2008,16(4): 539－545.ZHANG Xiao-ping, WU Shun-chuan, ZHANG Bing,et al. Particle flow code simulation of specimen mechanical behavior with different geometric parameters of weak seams[J]. Journal of Engineering Geology,2008, 16(4): 539－545.

[13]張曉平, 吳順川, 王思敬. 類土質路塹邊坡動態(tài)監(jiān)測及數值模擬分析[J]. 巖石力學與工程學報, 2008, 27(增刊2): 3431－3440.ZHANG Xiao-ping, WU Shun-chuan, WANG Si-jing.Dynamic monitoring and numerical analysis of soil-like cut slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(Supp.2): 3431－3440.

[14]張曉平. 類土質邊坡破壞機理分析及工程應用研究[D].北京: 北京科技大學, 2007.

[15]吳順川, 張曉平, 劉洋. 基于顆粒元模擬的含軟弱夾層類土質邊坡變形破壞過程分析[J]. 巖土力學, 2008,29(11): 16－21.WU Shun-chuan, ZHANG Xiao-ping, LIU Yang.Analysis of failure process of similar soil slope with weak intercalated layer based on particle flow simulation[J].Rock and Soil Mechanics, 2008, 29(11): 16－21.

[16]張曉平, 王思敬, 王幼明, 等. 二維離散元模擬抗滑樁的折算方法研究[J]. 巖土工程學報, 2010, 32(2): 271－278.ZHANG Xiao-ping, WANG Si-jing, WANG You-ming,et al. Conversion of anti-sliding piles into 2-dimensional discrete element simulation[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(2): 271－278.

[17]吳順川. 壓力注漿復合錨固樁地基處治理論研究及工程應用[D]. 北京: 北京科技大學, 2004.

[18]劉洋. 砂土液化破壞的細觀力學機制與數值模擬[D].上海: 同濟大學, 2006.

[19]張剛. 管涌現象細觀機理的模型試驗與顆粒流數值模擬研究[D]. 上海: 同濟大學, 2007.

[20]和法國, 諶文武, 張景科, 等. PS 材料加固交河故城土體試驗研究[J]. 敦煌研究, 2007, (5): 38－41.HE Fa-guo, CHEN Wen-wu, ZHANG Jing-ke, et al.Experimental research on PS reinforcing Jiaohe cultural heritage soil[J]. Dunhuang Research, 2007, (5): 38－41.

[21]POTYONDY D O, CUNDALL P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(8): 1329－1364.

[22]CHO N, MARTIN C D, SEGO D C, et al. A clumped particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(7): 997－1010.

[23]Itasca Consulting Group, Inc. PFC2D (Particle Flow Code in 2 Dimension)Version 3.1[M]. Minneapolis,Minnesota: Itasca Consulting Group, Inc., 2004.

[24]李萬清. 風沙躍移運動發(fā)展過程的離散動力學模擬[J].中國沙漠, 2006, 26(1): 47－53.LI Wan-qing. Discrete dynamics simulation on developing process of aeolian sand saltation[J]. Journal of Desert Research, 2006, 26(1): 47－53.

[25]HAFF P K, ANDERSON R S. Grain scale simulations of loose sedimentary beds: The example of grain-bed impacts in aeolian saltation[J]. Sedimentology, 1993,40(2): 175－198.

[26]亢力強, 郭烈錦. 風沙運動的 DPM 數值模擬[J]. 工程熱物理學報, 2006, 27(3): 441－444.KANG Li-qiang, GUO Lie-jin. Simulation of windblown sand movement by discrete particle model[J]. Journal of Engineering Thermophysics, 2006, 27(3): 441－444.

[27]亢力強, 郭烈錦. 風沙躍移中顆粒與多粒徑床面碰撞的數值模擬[J]. 工程熱物理學報, 2006, 27(1): 82－84.KANG Li-qiang, GUO Lie-jin. Numerical simulation of particles impacting multiple-grain size bed in saltation[J].Journal of Engineering Thermophysics, 2006, 27(1):82－84.

[28]亢力強, 郭烈錦. 風沙躍移中顆粒沖擊起動的數值模擬[J]. 自然科學進展, 2005, 15(2): 252－256.