考慮摩擦系數和顆粒強度劣化效應的堆石體濕化細觀數值模擬

黃緒武，周 偉，馬 剛，陳 遠，常曉林

(1. 武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072；2. 武漢大學水工巖石力學教育部重點實驗室，武漢 430072)

堆石體作為一種重要的建筑材料，已經在堆石壩、道路路基等基礎設施建設中得到廣泛應用。工程監測資料表明，堆石壩上游堆石區在庫水位變動，尤其是水庫初次蓄水時會產生較大的沉降，當遇到暴雨時，堆石壩下游堆石區也能觀察到明顯的沉降變形[1,2]。

為了研究堆石體濕化變形的規律，國內學者作了大量的試驗研究，并基于試驗結果提出了一系列濕化模型[3-5]。然而，這些模型多是建立在濕化變形與圍壓和應力水平的統計關系上，并沒有深入到細觀顆粒層面。因此，導致濕化變形的細觀力學機理一直沒有得到充分的認識。

在堆石體濕化變形機理方面，國內學者進行了一些試驗研究。左永振等[6]對雙江口壩殼料進行大型三軸濕化試驗發現，固結排氣剪的黏聚力c值最大，濕化后的c值次之，飽和排水剪的c值最小。粗粒料的c值是顆粒之間相互錯動時表現出的咬合力，浸水濕化后，由于顆粒浸水軟化和水的大量存在使顆粒間的潤滑作用增強，導致其c值比風干態的c值低。魏松等[7]在大型三軸儀上完成濕化試驗并結合顆粒破碎檢測試驗發現，濕化變形與濕化顆粒破碎之間有很好的線性關系。然而，這些試驗研究只是通過堆石體的一些宏觀物理量的變化來分析濕化變形的原因，無法得到每種因素對于濕化變形的影響程度以及影響機理。

浸水對堆石體的影響主要有2個：①堆石體遇水之后，水在顆粒接觸處會形成一層水幕，起到潤滑的作用，使得顆粒之間的有效摩擦系數降低；②水會使得巖石變脆變軟，導致堆石體單顆粒的破碎強度降低[8]。基于此，學者們開展了豐富的試驗研究。Ulusay等[9]通過傾斜試驗得到了不同巖石在干燥和浸水時的摩擦系數，與干態時的摩擦系數相比，浸水時的摩擦系數會降低。降低比例如下：安山巖降低9.46%，石灰巖降低13.07%，花崗巖降低16.3%。尤明慶等[10]通過巴西圓盤劈裂試驗發現，巖石在飽和狀態下的強度比干燥狀態下的強度低，降低比例分別為：大理巖降低15.3%，粗砂巖降低11.2%，細砂巖降低22.3%。Mohamad等[11]對不同風化程度的花崗巖進行點荷載試驗發現，浸水后Ⅰ類風化花崗巖強度平均降低7.1%。Hashiba等[12]通過不同巖石單軸試驗發現，浸水濕化后，巖石的抗拉、抗壓強度都會降低。降低比例如下：安山巖抗拉強度降低34%，抗壓強度降低15%；花崗巖抗拉強度降低7.6%，抗壓強度降低3.2%。

除了物理試驗外，學者們也在數值模擬方面進行了一些研究，楊貴等[13]運用PFC3D模擬堆石體三軸濕化試驗，通過改變顆粒剪切模量和摩擦系數產生濕化變形，但其并未考慮顆粒破碎的影響。Silvani等[14]運用PFC2D模擬水庫蓄水對于堆石壩變形的影響，通過考慮水對于堆石體骨架的浮力以及對摩擦系數的影響建立了濕化機理概化模型。然而在建立濕化機理概化模型時，并沒有結合試驗資料，只是假定了一些濕化后的強度和摩擦系數值，然后分析濕化變形的規律。

本文考慮水對于堆石體顆粒的兩種主要作用：①降低顆粒之間的摩擦系數；②低顆粒破碎強度，并結合室內試驗資料建立了濕化機理概化模型，然后將其引入連續-離散耦合分析方法(FDEM)[15,16]，來研究不同圍壓下濕化變形的規律以及顆粒破碎的規律，并從細觀組構[17]層面研究了濕化前后顆粒之間接觸力的變化，最后通過對比試驗分析了細觀摩擦系數和顆粒強度對于堆石體濕化變形的影響。

1 考慮顆粒破碎的連續離散耦合分析方法

1.1 顆粒破碎模型

顆粒破碎采用內聚力模型進行模擬[18]。使用內聚力模型時，需要定義初始起裂準則和材料損傷演化準則。達到初始起裂準則前，材料的應力與相對位移滿足線性關系；達到初始起裂準則后，材料剛度逐漸劣化并伴隨能量耗散。建立數值模型時，在所有毗連的實體單元之間插入無厚度的界面單元(Cohesive Interface Element，CIE)，見圖1。實體單元定義為線彈性模型，界面單元采用內聚力模型。

圖1 典型顆粒模型Fig.1 A typical particle model

內聚力模型起裂準則有3種：①拉應力達到材料抗拉強度ft；②剪應力達到材料抗剪強度fs；③拉應力、剪應力滿足二次起裂準則的混合模式，其表達式如下：

(1)

式中：〈〉為Macaulay括號；ts表示界面單元所受的切向應力。

材料達到初始起裂準則后繼續加載時，采用Benzeggagh和 Kenane[19]提出的基于能量的復合損傷演化準則，即：

(2)

式中：Gc為復合斷裂能；Gcn為Ⅰ型斷裂能；Gcs為Ⅱ型斷裂能；GS為2個切向斷裂能之和；GT為總斷裂能；Gn，Gs，Gt分別為法向和2個切向的斷裂能；δn為法向相對位移；δs，δt為2個切向的剪切滑移量；G0為初始損傷時的彈性應變能；t0為達到損傷準則時的應力；δm表示法向張開和切向滑移的等效位移量；D表示損傷因子；δ0m，δfm分別為初始損傷和界面完全失效時的等效位移量；η為表征界面拉伸損傷與剪切損傷相互作用的因子，對于堆石體顆粒材料通常取η=2。

在損傷演化過程中，D從0單調增加到1，這一過程中界面單元剛度表示為：

(3)

式中：kcn，kcs分別表示損傷演化過程中的界面單元的法向、切向剛度；k0n，k0s分別表示界面單元未損傷時的法向、切向剛度。

當材料剛度降為0時，界面單元完全失效，界面單元處出現自由面，實體單元在此處變為接觸關系。

1.2 顆粒接觸模型

離散顆粒之間的接觸算法采用通用接觸罰函數法[20]。法向接觸模型的表達式如下：

(4)

式中：p為顆粒之間的法向接觸力；Kn為顆粒之間的法向接觸剛度；δn為顆粒之間的法向接觸位移。

當相互接觸的顆粒之間存在相對滑移的趨勢時，接觸面之間可以傳遞切向接觸力，即：

τ=Ksδs

(5)

式中：Ks為顆粒之間的切向接觸剛度；δs為顆粒之間切向相對位移。

當摩擦力小于某一臨界值時，接觸面處于黏結狀態，此時接觸面之間可以發生很小的彈性滑移；當摩擦力超過臨界值之后，接觸面開始出現相對滑移變形。采用庫倫摩擦定律計算臨界摩擦力：

τcr=μp

(6)

式中：μ為摩擦系數。

1.3 濕化機理概化模型

由于堆石體濕化變形較為復雜且影響因素眾多，所以目前還沒有成熟的濕化機理模型存在。Zhao等[8]進行堆石體濕化數值模擬時，引入了用來表征水對堆石體摩擦系數和強度劣化作用的劣化因子的概念，劣化因子分別用λ1和λ2表示，它們代表了堆石體濕化后的摩擦系數和強度相比于濕化前的降低比例。本文在其基礎上結合相關室內試驗資料，建立了濕化機理概化模型。λ1和λ2的選取參考文獻[9,11,12]，文獻[9]進行傾斜試驗研究發現花崗巖試塊浸水時的摩擦系數相對于干燥狀態降低了16.3%，文獻[11]對花崗巖進行點荷載試驗發現，浸水后花崗巖強度平均降低7.1%，文獻[12]通過巖石單軸試驗發現，浸水濕化后花崗巖抗拉強度降低7.6%。文獻[21]選取雙江口壩殼料(主要為花崗巖)進行了一系列室內大三軸濕化試驗研究，為了與其進行比較，以巖性一致為原則選取以上花崗巖的試驗結果，取λ1=16%、λ2=7.5%。由于堆石體是由粒徑不同的堆石顆粒組成，堆石顆粒對應的劣化因子存在差異性，所以以上取值雖然不能做到絕對準確，但作為一個算例來獲取堆石體濕化變形的規律認識還是可行的。堆石體遇水后，顆粒的抗拉強度fwett、抗剪強度fwets以及顆粒之間的摩擦系數μwet由劣化因子表示為：

(7)

進行細觀數值模擬時，還需要作2個假設：①堆石體浸水過程是瞬時完成的，濕化對整個試樣的影響是均勻的、各向同性的；②不考慮水的浮力影響，相關研究表明，浮力的影響可以忽略不計[14]。

2 雙軸壓縮數值試驗

2.1 數值試樣

為了反應堆石顆粒形狀的多樣性，采用隨機均布Voronoi算法生成多邊形顆粒[22]。數值試樣的顆粒等效粒徑為4～20 mm(等效粒徑是指多邊形的等面積圓的直徑)。將松散的顆粒集合體等向壓縮至需要的尺寸和孔隙比，在壓縮制樣過程中，將重力加速度與顆粒之間的摩擦系數設為0，并使顆粒破碎不發生。最終生成的數值試驗試樣含有5 535個多邊形顆粒，劃分為43 221個實體單元、55 570個界面單元，孔隙比為0.145。數值試樣及顆粒級配曲線見圖2。

圖2 數值試樣及顆粒級配曲線Fig.2 Numerical specimen and its particle size distribution curve

數值試驗的步驟按照單線法濕化試驗的步驟進行：①水平和軸向方向等向加壓至初始圍壓；②保持水平向圍壓不變，繼續增大軸向壓力至指定的應力水平；③保持水平向、軸向的壓力不變，按照上述濕化機理概化模型調整細觀力學參數模擬濕化過程；④保持水平向壓力不變，繼續增加軸向壓力至試樣破壞。

2.2 細觀力學參數

數值試驗包含多個細觀力學參數，按照其用途可以劃分為2組：實體單元{ρ,E,υ,Kn,Ks,μ}，界面單元{k0n,Le,k0s,φ,c,ft,Gcn,Gcs}。其中ρ為材料密度；E，υ分別為彈性模量和泊松比；Le為實體單元網格的平均尺寸；其他參數意義同前。細觀摩擦系數的取值參照文獻[23,24]取為0.5，其他參數按照能反應堆石體力學特性的原則并結合已有的研究成果[20]進行取值，詳見表1。

表1 細觀力學參數取值Tab.1 Values of mesomechanical parameters

3 數值試驗結果

3.1 應力～應變關系

為了驗證模型以及細觀參數取值的合理性，首先進行了3組不同圍壓(0.8、1.6和2.4 MPa)下的雙軸剪切數值試驗，得到的偏應力～軸向應變，體積應變～軸向應變關系曲線見圖3。進行這3組數值試驗時，沒有引入濕化機理概化模型，所以得到的結果可視為室內試驗中的風干樣剪切試驗結果。從試驗結果來看，細觀數值模擬能再現堆石體應力～應變曲線的非線性、剪脹等特點，與室內試驗得到的規律相似。

圖3 不同圍壓下偏應力、體積應變與軸向應變的關系Fig.3 Relation curves of principal stress difference and volume strain with axial strain under different confining pressures

進行濕化數值試驗時，保持濕化過程中圍壓0.8、1.6和2.4 MPa不變，偏應力1.854、2.976和3.942 MPa不變，此時對應的濕化應力水平均為0.6。不同圍壓下濕化數值試驗得到的Δεa、Δεv見表2，得到的偏應力～軸向應變、體積應變～軸向應變關系曲線見圖4。其中Δεv>0表示剪縮體積應變，空心符號表示沒有引入濕化機理概化模型的干樣剪切數值試驗，實心符號表示引入濕化機理概化模型的濕化剪切數值試驗。

表2 不同圍壓下的濕化軸向應變和體積應變Tab.2 Axial strain and volume strain in the wetting process under different confining pressures

圖4 不同圍壓下濕化試驗偏應力、體積應變與軸向應變的關系曲線Fig.4 Relation curves of principal stress difference and volume strain with axial strain in the wetting tests under different confining pressures

由表2可知，數值試驗得到濕化軸向應變隨圍壓的增加而略有增加，在數值上較室內試驗結果稍偏小。數值試驗得到的濕化軸向應變與室內試驗結果在數值量級上接近，表明水對于堆石體顆粒之間的摩擦系數和顆粒強度的劣化作用是引起濕化軸向應變的主要原因。

濕化體積應變隨圍壓的增加明顯增加，與室內試驗得到的規律吻合。低圍壓下濕化體積變形表現為剪脹變形而高圍壓下表現為剪縮變形，這是因為低圍壓下濕化過程中的顆粒破碎數量較少，對應的剪縮變形小，而由于濕化后顆粒之間摩擦系數的降低引起的顆粒之間的滑移和重新排列較多，對應的剪脹變形大，所以低圍壓下濕化體積應變表現為剪脹變形，高圍壓下則剛好相反。

由圖4(a)可見，濕化后的峰值偏應力都較干樣峰值偏應力低，降低的幅度從低圍壓到高圍壓依次增大。2.4 MPa圍壓下濕化后的偏應力～軸向應變曲線都落在干樣的偏應力～軸向應變曲線下，0.8和1.6 MPa圍壓下濕化后的偏應力～軸向應變曲線除了峰后一小段外，其他也都落在干樣的偏應力～軸向應變曲線下，與室內試驗得到的規律大致吻合。由圖4(b)可見，濕化后試樣的剪縮性大于干樣的剪縮性，與室內試驗得到的規律吻合。

3.2 濕化過程中的顆粒破碎

圖5所示為不同圍壓下濕化過程中的顆粒破碎數量和濕化軸向應變的柱狀圖。由圖5可見，隨著圍壓的增加，顆粒破碎數量明顯增加，與文獻[7]得到的規律吻合。濕化軸向應變雖然也隨圍壓的增加而增加，但不如顆粒破碎數量顯著，1.6 MPa圍壓下的顆粒破碎數量約為0.8 MPa圍壓下的4.56倍，但濕化軸向應變只增加了0.012%。

圖6所示為濕化過程中試樣的顆粒破碎分布情況，其中每一個點代表一處顆粒破碎。由圖6可見，不同圍壓下濕化過程中的顆粒破碎分布呈現出局部集中的現象，即在一些局部區域內顆粒破碎數量明顯較周圍的多，如圖中橢圓標示區所示，與相關文獻研究發現的局部危機(Local Crisis)[14]現象吻合。這可以從以下2個方面解釋：①由于濕化后顆粒之間的摩擦系數降低，一些原本處于黏結狀態的接觸變成滑移狀態，顆粒之間發生重新排列，在這一動態調整過程中顆粒之間相互擠壓碰撞，使得接觸力突然增加，造成局部區域內顆粒破碎集中；②由于濕化后顆粒破碎強度降低，在周圍顆粒較大的外力作用下，局部區域內的顆粒不足以繼續承擔如此大的外荷載，出現較多的顆粒破碎，使顆粒之間的接觸點數和接觸面積增加，應力集中得到釋放。圖7所示為2.4 MPa圍壓下濕化前后試樣某一局部放大圖。由圖7可見，顆粒破碎現象可以用本文的方法顯示直觀地表達。

圖5 不同圍壓下濕化過程中的顆粒破碎數量和軸向應變Fig.5 Counts of broken CIEs with axial strain in the wetting process under different confining pressures

圖6 不同圍壓下濕化過程中的顆粒破碎分布Fig.6 Distribution of broken CIEs in the wetting process under different confining pressures

圖7 濕化前后的顆粒破碎局部放大圖Fig.7 Magnification of broken particle before and at the end of wetting process

3.3 顆粒間接觸力分布

(8)

圖8、圖9為2.4 MPa圍壓下濕化前后顆粒間法向和切向接觸力平均值的極坐標分布圖。由圖8、圖9可見，濕化后法向接觸力平均值較濕化前法向接觸力平均值略有增加，濕化后切向接觸力平均值較濕化前切向接觸力平均值明顯降低。這可以歸結為以下2個原因：①由于濕化后顆粒之間的摩擦系數降低，使得切向接觸力平均值明顯降低，此時外荷載主要由法向力鏈承擔，所以法向接觸力平均值有增大的趨勢；②濕化過程中存在一定量的顆粒破碎，導致顆粒之間的接觸點數和接觸面積增加，進而使得顆粒間的法向和切向接觸力平均值有減小的趨勢。由于前一個因素的影響較后一個因素的影響大，所以濕化后法向接觸力平均值較濕化前法向接觸力平均值略有增加，而在兩個因素的共同作用下，濕化后切向接觸力平均值較濕化前切向接觸力平均值明顯降低。

圖8 濕化前后的法向接觸力平均值各向異性分布極坐標圖Fig.8 Polar representation of anisotropy of average normal contact force before and at the end of the wetting process

圖9 濕化前后的切向接觸力平均值各向異性分布極坐標圖Fig.9 Polar representation of anisotropy of average tangential contact force before and at the end of the wetting process

3.4 影響濕化變形的主要因素

為了進一步研究影響濕化變形的主要因素，本文設置了2組對比試驗：①假定堆石體濕化過程只影響顆粒之間的摩擦系數，即只考慮濕化后μwet變為 ，此概化模型為W1；②假定堆石體濕化過程只影響顆粒強度，即只考慮濕化后ft變為fwett，fs變為fwets，此概化模型為W2。本文只進行了2.4 MPa圍壓下的對比試驗，試驗結果見圖10。圖10中D表示干樣剪切試驗，W表示濕化剪切試驗，W1和W2分別為2組對比試驗。

W、W1和W2對應的濕化軸向應變分別為0.401%、0.338%和0.236%，W1對應的濕化軸向應變明顯大于W2對應的濕化軸向應變且與W對應的濕化軸向應變更接近，說明顆粒間摩擦系數的降低是引起濕化軸向應變的主要因素，顆粒強度的降低其次。由圖10(a)可見，W1對應的峰值偏應力與W對應的峰值偏應力更接近，表明顆粒間摩擦系數對濕化后峰值偏應力的影響更大。

W、W1和W2對應的濕化體積應變分別為0.046%、0.035%和0.055%，W2對應的濕化體積應變最大，這是因為W2建模時只考慮了濕化對于顆粒強度的影響，濕化過程中的顆粒破碎數量較多，破碎生成的小顆粒填充到大顆粒之間的孔隙中，引起較大的剪縮變形。由圖10(b)可見，當軸向應變小于5%時，W1對應的體積應變～軸向應變曲線與W對應的體積應變-軸向應變曲線接近，W2對應的剪縮變形較大；當軸向應變大于5%時，W2對應的體積應變～軸向應變曲線與W對應的體積應變-軸向應變曲線接近，W1對應的剪脹變形較大，這表明濕化后的體積應變-軸向應變曲線前段主要受顆粒間摩擦系數降低的影響，后段則主要受顆粒強度降低的影響。

圖10 不同濕化模型對應的偏應力、體積應變與軸向應變的關系Fig.10 Relation curves of principal stress difference and volume strain with axial strain of different wetting models

4 結 論

(1)數值試驗得到濕化軸向應變隨圍壓的增加而略有增加，數值在量級上與室內試驗結果接近，表明水對于堆石體顆粒之間的摩擦系數和顆粒強度的劣化作用是引起濕化軸向應變的主要因素。濕化體積應變隨圍壓的增加明顯增加，與室內試驗得到的規律吻合，但是數值明顯較室內試驗結果低。低圍壓下濕化體積變形表現為剪脹變形而高圍壓下濕化體積變形表現為剪縮變形。

(2)隨著圍壓的增加，顆粒破碎數量明顯增加，濕化軸向應變雖然也有所增加，但不如顆粒破碎數量增加顯著。不同圍壓下濕化過程中的顆粒破碎分布呈現出局部集中的現象。

(3)濕化后顆粒之間的法向接觸力平均值較濕化前法向接觸力平均值略有增加，濕化后切向接觸力平均值較濕化前切向接觸力平均值明顯降低。

(4)顆粒間摩擦系數的降低是引起濕化軸向應變的主要因素，顆粒強度的降低其次，且顆粒間的摩擦系數的降低對濕化后峰值偏應力的影響更大。只考慮濕化對于顆粒強度的影響時，濕化過程中剪縮體積應變最大。濕化后的體積應變～軸向應變曲線前段主要受顆粒間摩擦系數降低的影響，后段則主要受顆粒強度降低的影響。

□

[1] 李全明, 于玉貞, 張丙印,等. 黃河公伯峽面板堆石壩三維濕化變形分析[J]. 水力發電學報, 2005,24(3):24-29.

[2] Neves E M D, Pinto A A V, Naylor D J, et al. A back-analysis of Beliche Dam[J]. Géotechnique, 1997,51(4):377-381.

[3] 程展林, 左永振, 丁紅順,等. 堆石料濕化特性試驗研究[J]. 巖土工程學報, 2010,(2):243-247.

[4] 彭 凱, 朱俊高, 王觀琪. 堆石料濕化變形三軸試驗研究[J]. 中南大學學報(自然科學版), 2010,41(5):1 953-1 960.

[5] 沈廣軍, 殷宗澤. 粗粒料浸水變形分析方法的改進[J]. 巖石力學與工程學報, 2009,28(12):2 437-2 444.

[6] 左永振, 程展林, 姜景山,等. 粗粒料濕化變形后的抗剪強度分析[J]. 巖土力學, 2008,29(S1):563-566.

[7] 魏 松, 朱俊高. 粗粒料三軸濕化顆粒破碎試驗研究[J]. 巖石力學與工程學報, 2006,25(6):1 252-1 258.

[8] Zhao Z, Song E X. Particle mechanics modeling of creep behavior of rockfill materials under dry and wet conditions[J]. Computers & Geotechnics, 2015,68:137-146.

[9] Ulusay R, Karakul H. Assessment of basic friction angles of various rock types from Turkey under dry, wet and submerged conditions and some considerations on tilt testing[J]. Bulletin of Engineering Geology & the Environment, 2015,28:1-17.

[10] 尤明慶, 陳向雷, 蘇承東. 干燥及飽水巖石圓盤和圓環的巴西劈裂強度[J]. 巖石力學與工程學報, 2011,30(3):464-472.

[11] Mohamad E T, Latifi N, Arefnia A, et al. Effects of moisture content on the strength of tropically weathered granite from Malaysia[J]. Bulletin of Engineering Geology and the Environment, 2016,75(1):1-22.

[12] Hashiba K, Fukui K. Effect of water on the deformation and failure of rock in uniaxial tension[J]. Rock Mechanics and Rock Engineering, 2015,48(5):1-11.

[13] 楊 貴, 劉漢龍, 朱俊高. 粗粒料濕化變形數值模擬研究[J]. 防災減災工程學報, 2012,(5):535-538.

[14] Silvani C, Bonelli S, Philippe P, et al. Buoyancy and local friction effects on rockfill settlements: a discrete modelling[J]. Computers & Mathematics with Applications, 2008, 55(2):208-217.

[15] 馬 剛, 周 偉, 常曉林,等. 考慮顆粒破碎的堆石體三維隨機多面體細觀數值模擬[J]. 巖石力學與工程學報, 2011,30(8):1 671-1 682.

[16] 馬 剛, 周 偉, 常曉林,等. 顆粒劣化效應的堆石料流變細觀數值模擬[J]. 巖土力學, 2012,(S1):257-264.

[17] Azéma E, Radja? F. Stress-strain behavior and geometrical properties of packings of elongated particles[J]. Physical Review E, 2010, 81:703-708.

[18] Nguyen O, Repetto E A, Ortiz M, et al. A cohesive model of fatigue crack growth[J]. International Journal of Fracture, 2001,110:351-369.

[19] Benzeggagh M L, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus[J]. Composites Science & Technology, 1996,56(4):439-449.

[20] Ma G, Zhou W, Chang X L. Modeling the particle breakage of rockfill materials with the cohesive crack model[J]. Computers & Geotechnics, 2014,61(61):132-143.

[21] 左永振. 粗粒料的蠕變和濕化試驗研究[D]. 武漢：長江科學院, 2008:40-52.

[22] Mollon G, Zhao J. Fourier-Voronoi-based generation of realistic samples for discrete modelling of granular materials[J]. Granular Matter, 2012,14(5):621-638.

[23] 日本土質工學會. 粗粒料的現場壓實[M]. 郭熙靈，文 丹，譯. 北京：中國水利水電出版社，1998.

[24] 李廣信. 高等土力學[M]. 北京: 清華大學出版社, 2004.