高心墻堆石壩心墻水力劈裂的顆粒流模擬

楊 艷 ，周 偉，常曉林，花俊杰

（1. 武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072；2. 西北農林科技大學 水利與建筑工程學院，陜西 楊凌 712100）

1 引 言

在眾多壩型中，土石壩因其具有廣泛的適應性、可就地取材、施工方法簡單且抗震性能良好等優點而成為高壩的首選壩型，在壩工建設中得到了極其廣泛的應用。特別是隨著我國西南水電開發進程的加快，在建和擬建的200～300 m 級高土石壩的數量在逐漸增加。但在高土石壩的建設中，心墻堆石壩作為土石壩中常用的壩型之一，還存在著一些急需解決的重要問題，其中心墻的水力劈裂問題是土石壩設計和建設中倍受關注的焦點之一。

目前，已有很多學者針對心墻的水力劈裂問題進行了研究，并取得了不少成果。黃文熙[1]指出，心墻中任何一點處的孔隙水壓力如果使該點處的最小主應力的有效值降低至心墻土料的抗拉強度，心墻就會沿著這個最小主應力面產生水力劈裂。殷宗澤等[2]分析了常用的有效應力法和總應力法計算水力劈裂時結果產生較大差異的原因，認為應該用心墻外水壓力是否超過心墻上游面處土中的中主應力來判別水力劈裂發生的可能。沈珠江等[3]模擬了Teton 壩水力劈裂的離心模型試驗，并得出結論：深截水槽心墻未必會發生水力劈裂破壞，常用的總應力法過大地估計了水力劈裂發生的可能性。張丙印等[4]研制了一種新型的水力劈裂試驗裝置，采用糯扎渡高心墻堆石壩心墻混合土料進行了水力劈裂試驗，結果證實了土石壩心墻中可能存在的滲水弱面以及水庫在快速蓄水的過程中產生的弱面水壓楔劈效應是導致水力劈裂發生的重要條件。李全明 等[5]采用彌散裂縫理論描述水力劈裂裂縫的發展過程，建立了用于描述水力劈裂發生和擴展過程的數學模型及有限元計算模式。朱俊高等[6]研究認為，水庫蓄水初期是水力劈裂的危險期；完全均質的心墻內不會發生水力劈裂；裂縫或局部的缺陷及迅速蓄水的初期是土石壩心墻發生水力劈裂的兩個重要條件，水力劈裂發生的根本原因是局部高水力梯度的存在。陳五一等[7]指出，目前水力劈裂分析常用的有效應力法存在不足，并提出改進方法，定義了一種判定水力劈裂的安全系數。曹雪山等[8]提出了研究心墻水力劈裂問題的非飽和土固結簡化計算的有效應力分析方法，研究發現，提高心墻的滲透系數和心墻填筑土的初始飽和度及在初次蓄水時放慢蓄水速度等均可防止心墻水力劈裂的發生。畢慶 濤[9]提出了用非飽和心墻料做固結不排水試驗得出總應力法計算模型參數的改進方法以及一種新的水力劈裂判斷標準，用心墻內緊靠上游表面單元的組合應力(32σ -1σ )與心墻前庫水壓力進行比較來判定水力劈裂的發生。

綜上所述，可以看出現有的成果大部分是通過試驗和有限元數值模擬等手段，從宏觀的角度研究水力劈裂的發生、發展過程及其力學機制，至今尚未形成一致的觀點，仍存在許多爭議。正是由于從宏觀的角度難以完全揭示水力劈裂的破壞機制，因此，本文在考慮流-固耦合的前提下，采用PFC2D軟件的顆粒流方法，從細觀的角度初步探究了心墻水力劈裂發生的機制，并將該成果與現有的室內試驗成果[4]進行比較分析，證實了該方法的可行性，對今后進一步研究心墻的水力劈裂具有重要的參考價值。

2 顆粒流流-固耦合的基本原理

顆粒流流-固耦合的基本原理[10]基于2 個假設：①假想固體中流體的滲流路徑是由顆粒間接觸處的平行板通道組成，該通道稱為“管道”（如圖1 中黑色線段所示，圖中灰色圓形表示顆粒，白色線段表示顆粒間的接觸），其管徑的大小與接觸處顆粒間的法向距離成正比；②本文采用的PFC2D顆粒流程序沒有模擬流體的存在，僅僅通過流體壓力來體現流體的作用，于是假設計算模型中存在可以存儲壓力的單元，該單元在PFC 中用域（如圖1 中由白色線條圍成的一個閉合多邊形區域，黑色圓點表示域的中心）來定義，相鄰域之間通過管道連接，于是可以根據域間的壓力差實現流體流動的模擬，而流體的流動將會引起域內壓力的改變。

圖1 管道和域的示意圖 Fig.1 Sketch of pipes and domains

PFC2D顆粒流程序在模擬固體和流體的耦合作用時，是通過下面的流體計算公式和壓力計算公式實現的。

2.1 流體計算公式

在PFC 中用管道概化流體的滲流路徑，該管道相當于一個平行板通道，其長度為L，管道孔徑為a，在垂直于平面的方向上取單位厚度；另外，域將流體計算區域離散后，計算域內的“流體”是通過相鄰域之間的壓力差產生流動。于是，由平行板均勻流的立方定律可以得到管道內的流量公式為

式中：k 為滲透系數；a 為管道的孔徑；P2-P1為兩個相鄰域的壓力差；L 為管道的長度。

對于具有一定滲透性的無裂縫材料而言，其PFC 模型中的滲流通道用管道概化后，a 的大小將影響模型的滲透性。假設當顆粒間的接觸力為0 時對應的管道孔徑為a0，并稱a0為殘余孔徑。當給模型施加實際的應力邊界條件后，接觸處法向可能的受力情況分別是受壓、受拉和不受力。

當法向接觸力為壓力時，a 將隨著接觸處法向壓力的增大而逐漸減小，此時a 與a0間的經驗關系如下：

式中：F0為a 減小到a0/2 時的法向壓力；F 為荷載作用下的法向接觸壓力。

當法向接觸力為拉力或者0 時：

式中：g 為兩顆粒表面間的法向距離；m 為兩顆粒表面間法向距離的縮放因子，其作用是調節a 的大小。

當接觸力為0 時，兩接觸顆粒表面間的法向距離g 為0，從式（3）中可以得出此時顆粒間接觸處的a =a0，這與前面的假設是一致的。

2.2 壓力計算公式

在一個時間步 tΔ 內，域中流體壓力的變化可以通過下式計算：

式中：Kf為流體的體積模數；Vd為域的表觀體積；∑q 為每個域從周圍管道中獲得的總流量，以流入為正；dVΔ 是考慮力的作用引起域的體積改變。

2.3 流-固耦合方式

在PFC2D中，流-固耦合的實現可以采用以下3種方式：（1）通過接觸的張開與閉合或者接觸力的改變實現管道孔徑的變化；（2）作用在顆粒上的力改變了域的體積，從而引起域內壓力的變化；（3）域內的壓力差導致顆粒上作用著滲透體積力。

以上3 種形式中，前兩點已經通過式（2）～（4）實現。針對第3 點，為了簡化施加在顆粒上的力，假設域內壓力施加在周圍顆粒上的作用力沿其周圍接觸顆粒間的連線均勻分布（如圖2 所示，圖中黑色圓點表示域的中心，白色線段表示相鄰接觸顆粒間的連線，黑色箭頭表示均勻分布的水壓力P）。

圖2 域內壓力分布 Fig.2 Distribution of pressure in domain

由此可知，作用在顆粒上的力矢量可以表示為

式中：ni為相鄰接觸顆粒間連線的法向單位矢量；s為顆粒圓心到接觸點的距離。

2.4 求解方法

PFC2D進行流-固耦合計算時采用顯示求解方法，對每個管道和域分別進行流量和壓力計算，并且在整個求解過程中兩者的計算是交替循環進行的。

然而，為了保證滲流計算的穩定性，滲流計算時步不應大于臨界時步，下面推導臨界時步的計算公式。假設某個域存在擾動壓力 pPΔ ，由于擾動流入域內的流量可以根據式（1）計算得出：

式中：R 為該域周圍顆粒的平均半徑；N 為連通該域的管道個數。

由式（4）可以計算出該流量在域內產生的壓力變化 rPΔ 為

要保證滲流計算穩定，那么擾動流量引起的壓力變化 rPΔ 必須小于擾動壓力 pPΔ ，當兩者相等時，由式（6）、（7）可求出臨界時間步長為

在計算中，整體的時間步長必須取所有局部時間步長中的最小值，同時需乘以一個小于1 的安全系數。

3 顆粒流水力劈裂模擬

由于土石或堆石壩的心墻近似滿足平面應變的特性，在模擬時可以按照平面應變問題考慮，因此，本文采用PFC2D進行顆粒流水力劈裂的二維模擬。

3.1 顆粒流模型的建立

本文參考張丙印等[4]的室內試驗數據資料，將本次顆粒流模型的尺寸取為20 cm×20 cm，顆粒半徑的取值范圍為2.0～3.0 mm，半徑取值的分布形式采用均勻分布，孔隙率為0.15。根據以上參數，本次模擬一共生成顆粒1 665 個（見圖3，圖中灰色圓形表示顆粒，四周的黑色線段表示墻）。 另外，考慮到心墻土為黏性土，本次模擬在顆粒間設置了接觸黏結。但不足之處是，目前PFC 中的細觀和宏觀參數之間還沒有建立起對應的關系，然而黏結強度（包括法向和切向）的取值將直接影響數值試驗中水力劈裂的起始壓力，因此，在模擬過程中，為了使得數值計算結果符合室內試驗的規律，經多次調試，最終選取的法向和切向黏結強度n_bond 和s_bond 均為1.05 kN。當接觸力超過黏結強度時，相應的接觸黏結發生破壞，從而可以模擬裂縫的形成。

圖3 顆粒流模型 Fig.3 Particle flow model

為了保證計算模型的邊界條件符合心墻壩上游面實際的應力條件，本文采用PFC2D中的伺服控制機制控制邊界墻的速度，從而使得計算模型達到給定的應力邊界條件，如圖3 所示。

3.2 心墻水力劈裂發生的必備條件

目前很多學者對于產生水力劈裂的條件有著一致的看法，認為產生水力劈裂必須具備以下條件[11]：①心墻中存在初始裂縫或缺陷；②心墻材料具有較低的滲透性；③能夠產生楔劈效應；④水庫快速蓄水。本文為了能夠較真實地模擬水力劈裂的發生和發展過程，盡量滿足上述的4 個必備條件。

本文生成模型時，在上游側設置了長為0.08 m，厚為0.005 m 的水平軟弱帶，并通過在模型的上游面和水平軟弱帶區域內設置相應的庫水壓力來實現水庫的快速蓄水過程，如圖4 中深灰色區域所示（為了更加清晰地顯示上游水壓力和軟弱帶的布置，圖中只標出了顆粒圓心的位置并未顯示顆粒，淺灰色線條表示顆粒間的接觸，圓圈1 表示測量圓1，本文中設置該測量圓是用于測量模型內部的應力值）。如果在初設的上游庫水壓力作用下未發生水力劈裂，則逐漸增加庫水壓力，直到水力劈裂發生。本文在模擬過程中設置水平初始缺陷的主要原因是，與豎向初始缺陷相比土石或堆石壩在逐層碾壓的施工過程中較易產生水平缺陷。

圖4 數值模型 Fig.4 Numerical model

3.3 計算參數的選取

采用顆粒流模擬土的實際宏觀力學行為，如心墻的水力劈裂，在模擬中能否實現細觀參數與宏觀力學參數的統一將直接影響結果的正確與否。然而，這是一個十分復雜的過程，目前還沒有統一的方法。本文為了減少在選取參數時的盲目性，首先對流-固耦合中的流體計算參數逐一進行敏感性分析，結果見圖5。

圖5 中的5 組曲線圖均顯示的是計算模型中平均應力(1σ +2σ )/2 隨時步的變化過程，此處的1σ 和2σ 是指通過模型內部的測量圓1 得到的主應力。由于在水力劈裂研究中十分關注裂縫的繼續擴展，因此，本文只跟蹤記錄了水平軟弱帶之后模型內部平均應力的變化過程，將測量圓1 布置在如圖4 所示的位置。而本文之所以將心墻中的平均應力作為衡量指標，主要是考慮到目前很多學者認為，當心墻前的外水壓力值大于心墻內的某一主應力（1σ 、2σ或3σ ）時，將導致水力劈裂的發生。由此可見，心墻內的應力值是判斷水力劈裂發生的重要數據之一。另外，圖5 中只顯示了5 000 步之后平均應力的變化過程，這是因為在此之前是模型生成過程以及在給定的應力邊界條件下達到平衡所經歷的時步，而本文研究的重點是快速蓄水之后模型內部的應力變化及裂縫產生和發展的過程，因此，圖5 僅從施加水荷載之后開始記錄平均應力的變化過程。

圖5 平均應力-時步關系曲線 Fig.5 Curves of average stress-step

從圖5(a)、(b)可以看出，殘余孔徑a0和滲透系數k 取不同值對試樣平均應力有較大影響。當a0值在0.010 0～0.000 5 mm 之間變化時，試樣平均應力的最大差值達到0.16 MPa，是a0=0.010 0 mm 時的1.3 倍；當k 值在1.0×10-6～1.0×10-9mm/s 之間變化時，試樣平均應力的最大差值達到0.13 MPa，是k = 1.0×10-6mm/s 時平均應力值的0.79 倍。這主要是由于上述2 個參數直接反映了心墻的滲透性，而心墻發生水力劈裂的實質是在高水力梯度的作用下，裂縫形成并發展的過程。因此，在模擬過程中應選取適當的滲透參數，確保心墻材料的低滲透性，從而模擬高水力梯度的作用。而a =a0/2 時的法向壓力F0、流體體積模數kf以及顆粒表面間法向距離的縮放因子m 分別取不同值時，試樣平均應力的最大差值僅為0.04 MPa，由此可見，上述3 個參數在合理地取值范圍內變化時，對平均應力的結果影響較小。因此，參數選取時應重點調整a0和k 的取值。

此外，圖5 中5 組曲線的變化規律基本一致，大約在5 500 步之前，壓應力隨時間步的增加而逐漸減小，這是因為在水壓力作用下，初始的水平軟弱帶內形成了水平裂縫，在裂縫尖端附近區域的接觸黏結受拉（見圖6，圖中黑色線段表示顆粒間的接觸黏結承受拉力，線段越粗表示其對應的拉力值越大，灰色線段表示顆粒間的接觸），這部分拉力抵消了一部分豎向初始壓力，這種效應被稱為水壓楔劈效應[3]。而隨后壓應力又隨時間步的增長而逐漸增加，其原因是當接觸黏結破壞導致顆粒間產生裂縫后，這部分拉力被釋放，作用在顆粒上的水壓力使得顆粒向裂縫兩側移動，增加了裂縫兩側顆粒間的接觸壓力。但是當a0和k 取值較大時，由于模型的滲透性較大未能形成高水力梯度的作用，因此，沒有發生水力劈裂，之前產生的拉應力也就無法釋放。

圖6 裂縫尖端區接觸拉力分布 Fig.6 Distribution of contact tension around crack tip

本文在敏感性分析的基礎之上，經過反復調整參數取值，最終選取的流體計算參數值見表1。

表1 流體計算參數 Table 1 Computational parameters of fluid

3.4 數值模擬結果及分析

本文在模擬心墻水力劈裂過程時，首先取一初始水壓力作用在模型的上游表面和水平弱面內，如果在給定的水壓力值下運行2 000 步后未能在模型內部產生新裂縫（本文以接觸黏結發生破壞作為相應接觸處產生裂縫的判定標準），則在計算過程中逐漸增加水壓力值，直至在其內部產生連通上下游的貫穿性裂縫（如圖7 所示），并在數值模擬中以此作為判別水力劈裂發生的依據，此時對應的庫水壓力值即為劈裂水壓力Pf。

圖7 PFC2D 模擬結果 Fig.7 Simulation result of PFC2D

本文針對表2 中1σ 、3σ 的3 種組合方案分別進行了顆粒流模擬，并將PFC2D的模擬結果與室內試驗成果[3]做了比較（見表2 和圖8）。從圖可以看出，由數值模擬得到的結果，其規律與室內試驗的基本一致，即劈裂水壓力Pf隨著豎向應力的增大而增大，且兩者基本成線性關系。但由于在PFC2D的數值模擬過程中還存在著許多不確定的因素，如計算結果受顆粒的大小以及顆粒的組裝形式等因素的影響較大。另外，在數值模擬過程中還做了一些簡化和假定，從而導致計算結果均較室內成果偏大。

當豎向應力為0.2 MPa 時，本文采用顆粒流數值模擬得到的劈裂水壓力為0.21 MPa，比室內試驗的結果大0.03 MPa，偏大部分占室內結果的16.7%；當豎向應力采用0.3 MPa 時，本文得到的劈裂水壓力為0.29 MPa，比室內試驗的結果大0.02 MPa，偏大部分占室內結果的6.9%；在第3 種方案中豎向應力取為0.4 MPa，此時本文計算得到的劈裂水壓力為0.38 MPa，與室內試驗的結果相比大0.03 MPa，偏大部分占室內結果的8.6%。經過上述的比較可以看出，采用PFC2D進行水力劈裂的模擬，其結果雖與室內試驗的成果有一定出入，但偏大的范圍僅在6.9%～16.7%之間，這說明采用顆粒流方法研究心墻的水力劈裂是一種可行的方法。

表2 PFC2D 數值計算結果 Table 2 Numerical results of PFC2D

圖8 數值結果與室內試驗成果 Fig.8 Numerical and laboratory results

產生水力劈裂后的模型見圖7，從該圖中可以清晰地看到模型在高水力梯度作用下，由初始水平軟弱帶向下游發展形成的水力劈裂面，見圖中近似水平的白色連通區域（其內部用黑色的線段表示顆粒間產生的裂縫）。此外，圖中還顯示了模型內部速度場的分布，從局部放大圖中可以看出，裂縫周圍顆粒的速度方向（黑色箭頭所示）基本垂直于裂縫的延伸方向，從而初步證明了心墻水力劈裂破壞應屬于顆粒間的張拉破壞，其破壞的主要力學原因是由于心墻中的張拉應力超過了土體的抗拉強度。

然而，目前關于水力劈裂發生的力學原因，主要有2 種觀點，即拉裂破壞和剪切破壞[11]。本文在模擬水力劈裂的發展過程中，除了根據裂縫周圍顆粒的速度方向定性地證明了水力劈裂破壞應屬于顆粒間的張拉破壞外，還通過PFC2D程序分別統計并記錄了法向和切向接觸黏結破壞數目的變化過程，圖9 顯示的是初始豎向應力為0.2 MPa 時的結果，另外2 種方案的結果規律與之類似。圖中初始的21條法向黏結破壞裂縫是由于在水平軟弱帶瞬時施加的水壓力產生的。當運行到6 622 步時，模型內部形成了貫穿性裂縫，總的法向黏結破壞裂縫數目達到57 條。從該圖可以看出，隨著水力劈裂過程的發展，顆粒之間的接觸黏結逐漸破壞，其中法向黏結破壞的數目隨著時步在逐漸增加，而切向黏結破壞的個數卻始終保持0，再一次證明水力劈裂現象的發生主要是由心墻內部法向拉裂破壞導致。

此外，從圖9 中還可以看出，形成的裂縫數目與時步之間并不呈線性關系。其原因是，只有在裂縫尖端處積聚的能量足夠使裂縫進一步擴展時，才能產生新的裂縫，因此，圖中線段斜率較小的時步段對應的是裂縫尖端處能量積聚以及內部調整的時段，一旦能量足夠大，就會形成大量新裂縫。

圖9 法向裂縫數目-時步關系曲線 Fig.9 Curves of number of normal crack-step

在本次數值模擬中，從裂縫的發展方向（見圖7）可以看出，水力劈裂面近似垂直于1σ 的作用方向，這表明在滿足高水力梯度作用下，產生的水力楔劈效應降低了裂縫尖端區附近的1σ ，從圖8 和表2 中可以看出，對于3 組不同的（1σ ，3σ ）組合，當1σ 小于或接近相同高程處心墻上游面的外水壓力時，就會發生水力劈裂。

4 結 論

（1）本文運用PFC2D模擬了心墻水力劈裂的發生和發展過程，并將模型中產生連通上下游的貫穿性裂縫作為判別水力劈裂發生的依據。文中將數值模擬的結果和室內試驗的成果做了對比，兩者的規律基本一致，即隨著豎向應力的增加所需的劈裂水壓力值大致呈線性增加趨勢。

（2）通過數值模擬的結果可以看出，水力劈裂發生的主要原因是心墻在高水力梯度作用下形成的水楔效應導致了心墻內裂縫尖端區的張拉破壞，并且破壞面是近似垂直于1σ 的作用方向，這表明水楔效應降低了心墻內原有的1σ ，當該值小于或接近心墻的外水壓力時將會發生水力劈裂。

（3）由于模型在發生水力劈裂破壞的過程中法向黏結破壞的數目隨著時步在不斷增加，切向黏結破壞卻始終未發生，進一步證明了水力劈裂的破壞形式屬于法向張拉破壞。

（4）本文的研究成果還存在不足之處，尚需進一步改進。如：顆粒的大小及組裝形式對計算結果有較大影響；計算參數的選取還沒有較嚴格的理論依據；流-固耦合中還沒有考慮到應力對滲透特性的影響；研究成果未能深入揭示水力劈裂發生的機制等。但本文的研究成果為今后進一步研究心墻的水力劈裂提供了新思路。

[1] 黃文熙. 對土石壩科研工作的幾點看法[J]. 水利水電技術, 1982, (4): 23－27.

[2] 殷宗澤, 朱俊高, 袁俊平, 等. 心墻堆石壩的水力劈裂分析[J]. 水利學報, 2006, 37(11): 1348－1352. YIN Zong-ze, ZHU Jun-gao, YUAN Jun-ping, et al. Hydraulic fracture analysis of rock-fill dam with core wall[J]. Journal of Hydraulic Engineering, 2006, 37(11): 1348－1352.

[3] 沈珠江, 易進棟, 左元明. 土壩水力劈裂的離心模型試驗及其分析[J]. 水利學報, 1994, (9): 67－77. SHEN Zhu-jiang, YI Jin-dong, ZUO Yuan-ming. Centrifuge model test of hydraulic fracture of earth dam and its analysis[J]. Journal of Hydraulic Engineering, 1994, (9): 67－77.

[4] 張丙印, 李娜, 李全明, 等. 土石壩水力劈裂發生機理及模型試驗研究[J]. 巖土工程學報, 2005, 27(11): 1277－1281. ZHANG Bing-yin, LI Na, LI Quan-ming, et al. Mechanism analysis and model test of hydraulic fracturing in embankment dams[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(11): 1277－1281.

[5] 李全明, 張丙印, 于玉貞, 等. 土石壩水力劈裂發生過程的有限元數值模擬[J]. 巖土工程學報, 2007, 29(2): 212－217. LI Quan-ming, ZHANG Bing-yin, YU Yu-zhen, et al. Numerical simulation of the process of hydraulic fracturing in earth and rockfill dams[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(2): 212－217.

[6] 朱俊高, 王俊杰, 張輝. 土石心墻水力劈裂機制研究[J]. 巖土力學, 2007, 28(3): 487－492. ZHU Jun-gao, WANG Jun-jie, ZHANG Hui. Study on mechanism of hydraulic fracturing in core of earth-rockfill dam[J]. Rock and Soil Mechanics, 2007, 28(3): 487－492.

[7] 陳五一, 趙顏輝. 土石心墻水力劈裂計算方法研究[J]. 巖石力學與工程學報, 2008, 27(7): 1380－1385. CHEN Wu-yi, ZHAO Yan-hui. Study of calculation method of hydraulic fracturing for core of earth-rockfill dam[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(7): 1380－1385.

[8] 曹雪山, 殷宗澤. 土石壩心墻水力劈裂的非飽和土固結方法研究[J]. 巖土工程學報, 2009, 31(12): 1851－1856. CAO Xue-shan, YIN Zong-ze. Consolidation method of unsaturated soils for hydraulic fracturing of core walls of rock-fill dams[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(12): 1851－1856.

[9] 畢慶濤. 心墻水力劈裂的總應力法分析[J]. 水力發電, 2010, 36(3): 51－54. BI Qing-tao. Hydraulic fracturing analysis of core wall by total stress method[J]. Water Power, 2010, 36(3): 51－54.

[10] NG A K L, SMALL J C. A case study of hydraulic fracturing using finite element methods[J]. Canada Geotechnical Journal, 1999, 36(5): 861－875.

[11] Itasca Consulting Group Inc. PFC2D User’s Manual(version 3.1)[M]. Minneapolis: Itasca Consulting Group Inc., 2004.

[12] Itasca Consulting Group Inc. PFC3D verification problems and example applications[M]. Minneapolis: Itasca Consulting Group Inc., 2005.

[13] 王俊杰, 朱俊高, 張輝. 關于土石壩心墻水力劈裂研究的一些思考[J]. 巖石力學與工程學報, 2005, 24(增刊2): 5664－5668. WANG Jun-jie, ZHU Jun-gao, ZHANG Hui. Some ideas of study on hydraulic fracturing of core of earth-rockfill dam[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(Supp.2): 5664－5668.