高壓條件下石英砂中的靜力觸探試驗及大變形模擬

孔壯壯,王棟,張民生,裴會敏

（中國海洋大學 山東省海洋環境地質工程重點實驗室，山東 青島 266100）

深基礎在陸上和海洋工程中的應用日益增多，例如，近海導管架平臺的樁基埋深可達150 m 以上，樁基底端的上覆有效應力水平約在1～2 MPa 之間[1]；中國南海第二輪天然氣水合物試采的儲層埋深為208 m，水平開采套管承受的上覆有效應力約為1.5 MPa[2]。工程實踐中常采用原位試驗確定土層的物理和力學性質，但大多針對中等或低應力水平土層，已有的原位試驗數據解譯方法也很少考慮高應力水平的影響。

靜力觸探（CPT）試驗是海洋工程勘察中最常用的原位試驗方法，可通過量測的錐尖阻力經驗性地預測砂層的相對密實度和內摩擦角[3-5]。為建立預測公式，需要進行一定數量的室內標定罐試驗：在標定罐中制備均勻砂樣，對砂樣施加預定的上覆壓力或圍壓，量測施加壓力后砂樣的相對密實度；將靜力觸探儀貫入砂樣中，記錄貫入阻力隨深度的變化；改變砂土的相對密實度和壓力水平，進行多個試驗。Lunne 等[6]總結了不同地區石英砂的標定罐試驗和現場測試數據，建立起錐尖阻力qc和砂土相對密實度Dr值及平均有效應力p′之間的關系，但該依托試驗的p′值不超過500 kPa。Ahmadi 等[7]開展了不同K0條件下的標定罐試驗，考慮不同K0狀態對錐尖阻力的影響，也提出過類似公式，但其研究中的p′值不超過200 kPa。目前公開報道的高應力條件下砂土標定罐試驗和現場試驗的數據極少。

標定罐試驗砂樣準備時間長、成本高，數值模擬可以作為有效輔助手段，以完善錐尖阻力與相對密實度之間的關系。然而，傳統的小變形或有限變形框架內的有限元方法很難模擬CPT，原因是探頭在土中的貫入導致周圍土體網格的嚴重扭曲，使得計算無法進行。近年來，多種大變形數值方法已被用于追蹤觸探儀在砂土中的貫入過程，如任意拉格朗日—歐拉方法（Arbitrary Lagrangian Eulerian method，簡稱ALE）[8]、耦合歐拉—拉格朗日方法[9]和物質點法[10]。ALE 方法結合了拉格朗日和歐拉步驟，允許網格獨立于材料移動，保持剩余的網格拓撲不變。與傳統的拉格朗日分析相同，滿足平衡方程、邊界條件、外部載荷和接觸條件，然后重新定位節點，并將所有變量從舊網格映射到新網格。因此，當發生大變形時，可以保持高質量的網格。Kouretzis等[11]采用ALE 法模擬砂土中靜力觸探試驗，并與離心機試驗結果相互驗證，建立了相對密實度和錐尖阻力之間的關系，但其試驗中p′值不超過340 kPa。

筆者采用高圍壓標定罐試驗和ALE 大變形有限元方法，進行砂土中靜力觸探貫入，建立高應力條件下錐尖阻力與石英砂相對密實度的關系。大變形模擬中采用修正摩爾—庫倫模型（MMC）本構模型描述砂土內摩擦角和剪脹角隨累積塑性剪應變的變化。

1 標定罐試驗

采用某種青島石英砂進行標定罐試驗，砂的顆粒級配曲線如圖1 所示，平均粒徑D50=0.173 mm，不均勻系數Cu=1.70，最大孔隙比emax=0.949，最小孔隙比emin=0.490。

圖1 青島某石英砂顆粒級配曲線Fig.1 Particle size distribution curve of a silica Qingdao sand

早期標定罐側壁是剛性的，當罐體直徑不夠大時，容易出現明顯的邊界效應。目前，標定罐都采用柔性雙壁式標定罐，在剛性外壁和柔性內壁之間施加水頭形成圍壓，從而有效消除剛性邊界對貫入阻力的增強效應。筆者使用的高壓標定罐（圖2）也屬于柔性雙壁式：試樣直徑為600 mm、高度為750 mm，試樣四周被柔性囊包裹，圍壓最大為5 MPa。采用砂雨法制備試樣，罐口上部安裝水平往復運動的砂斗，運動速度自動控制。砂斗底部開口寬度及其與砂面的相對高度為手動調節，試樣的相對密實度取決于撒砂高度、開口寬度和運行速度[12]。

圖2 中國海洋大學研制的高壓標定罐Fig.2 High-pressure calibration chamber developed in Ocean University of China

砂中靜力觸探儀的貫入基本處于排水條件，Kluger 等[13]進行了飽和砂和干砂中的CPT 比較試驗，認為當應力水平很低時，干砂樣的錐尖阻力會略高于飽和砂樣；但當平均有效應力p′增加到200 kPa 時，干砂樣和飽和砂樣的錐尖阻力無明顯差異。因此，采用干砂雨法制備不同相對密實度的試樣：為了制備兩種不同相對密實度的砂樣，砂樣下落高度分別采用0.5、1.3 m，每次撒砂的砂斗開口保持3 mm 不變，運行速度控制在3 mm/s。在標定罐內部每隔0.1 m 進行標記，撒砂厚度一旦達到標記線，即升高砂斗0.1 m，以保證砂樣的下落高度基本不變。施加有效圍壓σ′3分別為0.5、1.0、2.0 MPa，采用等向固結，即靜止土壓力系數K0=1。提前在罐底布設4 個體積為200 mm3的帶底環刀，用于試驗后確定砂樣的相對密實度。

施加圍壓后砂樣大致可分為兩組，相對密實度分別約為46%和60%。共完成6 個不同相對密實度和圍壓條件下的標定罐試驗，對應的條件見表1。采用的CPT 探頭直徑D=20 mm，錐角為60°，探頭貫入速度為20 mm/s。

表1 標定罐中砂樣相對密實度Table 1 Relative density of sand sample in calibration chamber

2 有限元模擬

2.1 本構模型及參數標定

砂土在剪切時表現出的應變硬化或軟化行為取決于砂土的相對密實度與應力水平。傳統摩爾—庫倫模型規定砂土內摩擦角和剪脹角為定值，不能描述砂土的軟化和剪脹性。這里采用Hu 等[14]提出的MMC 模型，通過引進內摩擦角和剪脹角隨累積塑性剪應變γ的變化來描述砂土的軟化和剪脹性，其變化規律如圖3 所示，圖中φi、φp和φcv分別為初始、峰值和臨界內摩擦角。因為摩爾—庫倫類模型無法描述剪縮，所以假定當γ<γ1時，剪脹角為0。此后剪脹角迅速線性增大，在γ=γ2時達到峰值剪脹角ψp；剪脹角保持峰值至γ=γ3，然后線性減小，直至在臨界應變γ4處減小到0。當砂樣表現為硬化和剪縮性時，MMC 模型實質上退化為傳統的摩爾—庫倫模型：內摩擦角保持為φcv，剪脹角始終為0。對于一般石英砂，按照Hu 等[14]和Zheng 等[15]的建議，取φi=φcv、γ1=1%、γ2=1.2%、γ3=5%和γ4=15%。

圖3 MMC 模型中內摩擦角和剪脹角的變化Fig.3 Variation of friction and dilation angles of modified Mohr-Coulomb model

MMC 模型包含的參數有：泊松比ν，取定值0.35；黏聚力c為0；剪切模量G；峰值內摩擦角φp；臨界內摩擦角φcv和峰值剪脹角ψp。通過彎曲元試驗推算G，利用三軸排水剪切試驗獲得φp、φcv和ψp。

2.1.1 彎曲元試驗 彎曲元試驗用于確定土體在微小應變（應變約為10-6量級）時的最大剪切模量Gmax。將彎曲元量測設備安裝在三軸儀上，在三軸儀中分別制備初始孔隙比為0.694 和0.764 的砂樣，逐級施加100～1 500 kPa 的有效圍壓。每一級的固結完成后，首先測定排水量，用于計算固結后的孔隙比e，然后進行彎曲元試驗，獲得微小應變水平的Gmax。Gmax是e和p′的函數[16]。

圖4 為彎曲元試驗得到的歸一化最大剪切模量Gmax/pa隨p′/pa的變化，試驗數據可擬合為

圖4 青島某石英砂的小應變剪切模量Fig.4 Small-strain shear modulus of a silica Qingdao sand

彎曲元試驗確定Gmax時的應變約為10-6量級，而MMC 模型中剪切模量G的應變范圍大致在10-2量級。Loukidis 等[17]提出一種由Gmax估計G的方法。

式中：T為退化因子。按照Papadimitriou 等[18]和Pei等[19]的建議，采用ALE 方法試算CPT 在標定罐中的貫入，根據試算反推，確定T=5.5。

2.1.2 三軸排水試驗 為了獲取砂土的φ和ψ，進行6 個三軸排水剪切試驗，其相對密實度約為46%和60%，施加有效圍壓為0.5、1.0、2.0 MPa，剪切速率為0.5 %/min，具體方案和試驗結果見表2。由于相對密實度約46% 的砂土在圍壓1.0、2.0 MPa 下表現為硬化和剪縮現象，故φp=φcv，ψp=0。

表2 三軸排水試驗方案和結果Table 2 Scheme and results of drained triaxial tests

根據三軸排水剪切試驗的應力—應變曲線峰值強度和臨界強度可以分別求得砂土的φp和φcv。φcv基本不依賴應力水平與初始孔隙比，在6 個試驗中大致為定值34.4°。剪脹角利用式（3）根據體變—應變曲線得出。

式中：dε1和dε3為主應變增量；三軸條件下k=2。

為了后續的大量變動參數分析，需要建立砂土峰值內摩擦角或峰值剪脹角與相對密實度、平均有效應力之間的關系。參考經典Bolton 公式[20]，擬合表2 的三軸試驗結果，得到式（4）～式（6）。

式中：IR為剪脹系數。由于φcv為定值，一旦確定Dr和p′后，即可根據式（4）～式（6）計算φp和ψp。

2.2 有限元模型

采用通用有限元軟件ABAQUS 中的大變形模塊ALE 模擬CPT 試驗。建立的ALE 模型如圖5 所示。參考Fan 等[21]和Mahmoodzadeh 等[22]的成果，在錐尖下人為設置直徑為D/20 的剛性小管，此剛性小管與CPT 一起移動，從而防止最左邊的土體單元向對稱軸的內部移動，有利于保證土體網格調整的有效性，同時不會影響計算結果。

圖5 有限元軸對稱模型Fig.5 Finite element axisymmetric model

按照二維軸對稱計算，土體區域的寬和高分別是28D和42D。為了減少計算量，實施網格調整的ALE 區域寬10D，根據前期試算，ALE 區域中單元大小設置為D/8，這樣能在保證計算精度的同時提高計算效率。ALE 區域需要保持節點位置調整后網格的質量，以避免單元畸變，選用ABAQUS 內置幾何加強形式的網格平滑方法[11]。

采用硬接觸算法模擬靜力觸探儀與土體之間的相互作用。假定二者交界面上的摩擦應力為0，即忽略貫入過程中探桿上的側摩阻力，這是因為砂土中探桿上的側摩阻力遠小于錐尖阻力。

設置CPT 貫入速度為20 mm/s 進行勻速貫入，貫入深度為0～20D，采用ALE 大變形方法避免網格扭曲，從而實現靜力觸探在砂土中的連續貫入。為了與標定罐試驗對比，對土體施加K0=1的壓力。

3 ALE 有限元模型的驗證

3.1 試驗和模擬結果對比

圖6 對比了有限元模擬和標定罐試驗獲得的錐尖阻力，其中H代表錐尖深度。從標定罐試驗結果曲線可以看到，每條錐尖阻力曲線都分為3 個階段：1）H=0～1.5D時，qc隨深度幾乎呈線性增加，增長斜率和圍壓與相對密實度呈正相關；2）H=1.5D～10D時，qc的增長斜率變緩并逐漸趨于穩定；3）H>10D時，qc略微減小或基本保持不變。qc略微減小的原因可能是在試樣上放置的頂蓋使得淺層試樣的密實度略高于深部。

圖6 有限元模擬結果和標定罐試驗結果對比Fig.6 Comparison of finite element simulation results and calibration chamber test results

高應力條件下錐尖阻力仍滿足隨相對密實度和圍壓增大而增大的現象。對比Dr約為46%的砂樣，圍壓1.0 MPa 或2.0 MPa 時的錐尖阻力大約為0.5 MPa 時的1.76 倍和2.66 倍；對比Dr約為60%的砂樣，圍壓1.0 MPa 或2.0 MPa 時的錐尖阻力大約為0.5 MPa 時的1.58 倍和2.66 倍。這表明對于不同相對密實度的砂樣，高應力條件下的錐尖阻力仍依賴圍壓水平。如果圍壓相同，Dr約為60%砂樣的錐尖阻力是Dr約為46%砂樣的1.17～1.3 倍。

整體來說，ALE 方法能夠較好地模擬高圍壓石英砂中的CPT 貫入。從圖6 中可以看出：

1）對于不同相對密實度的砂樣，當貫入深度H>10D時，有限元計算得到的qc基本不再變化，存在一定程度的噪音波動，但在可接受范圍內。

2）對于Dr約為46%的砂樣，試驗與有限元結果基本吻合。

3）對于Dr約為60%的砂樣，H<10D時的試驗結果比有限元結果約高15%，但二者給出的H>10D時的穩態錐尖阻力很接近。

3.2 剪切模量敏感性分析

為了進一步探討獲得的剪切模量G對錐尖阻力的影響，保持峰值摩擦角和剪脹角不變，人為變化G值。ALE 得到的錐尖阻力和剪切模量的變化曲線如圖7 所示。隨著G的增加，錐尖阻力逐漸增大，G=200 MPa 時對應的錐尖阻力是G=50 MPa 時的2.7倍，因此，數值模擬中不能直接采用利用彎曲元獲得的低應變條件下的剪切模量，而應該按照式（2）計算。

圖7 錐尖阻力隨剪切模量變化曲線Fig.7 Changing curve of cone tip resistance with shear modulus

4 相對密實度和錐尖阻力的關系

為了建立高應力下砂土錐尖阻力和相對密實度之間的定量關系，進行變動參數分析：圍壓分別為0.5、0.8、1.0、1.3、1.6、2.0 MPa；選擇典型相對密實度Dr=40%、60% 和80%。利用式（2）計算G，利用式（5）和式（6）計算φp和ψp。

圖8 給出了0.8、1.3 MPa 圍壓下不同相對密實度的錐尖阻力曲線。對比圍壓為0.8、1.3 MPa 的砂樣，相對密實度80%和60%時的qc大約是40%時的1.22 倍和1.56 倍左右。當Dr=40%、60%和80% 時，圍 壓1.3 MPa 時，qc約 為0.8 MPa 時的1.35 倍。圖9 給出了CPT 貫入過程中的土體位移矢量圖，大變形條件下土體位移場可近似認為是破壞場。可以看出，錐面附近土體的位移大致沿錐面外法線方向，CPT 探頭下的土體處于三軸壓狀態，被徑向擠出。

圖8 不同圍壓下錐尖阻力曲線Fig.8 Cone tip resistance curves under different confining pressures

圖9 CPT 貫入過程土體位移矢量圖Fig.9 Vector diagram of soil displacement during CPT penetration

Lunne 等[6]總結的相對密實度和錐尖阻力之間的關系為

式中：參數C1、C2與土體性質有關；Qtn為歸一化錐尖阻 力；m值 取0.6。Ahmadi 等[7]考慮不同K0狀態對錐尖阻力的影響，建議m值取0.5。

取貫入深度H=10D作為穩定時的錐尖阻力進行分析。總結18 個工況的ALE 模擬結果發現，在p′=0.5～2.0 MPa 的高應力水平下，仍可以按照式（7）由錐尖阻力預測相對密實度，但式（7）的參數擬合為：C1=0.58、C2=1.91、m=0.64。即

按照式（9），由6 個標定罐試驗得到的穩定錐尖阻力預測Dr，所得結果如圖10 所示。從圖10 中可以看出，式（9）能較好地預測高應力條件下（p′=0.5～2.0 MPa）砂土的Dr，其誤差在±10%以內；圖10 同時包含了Lunne 等[6]、Pournaghiazar 等[23]和Bolton 等[24]的標定罐和離心機試驗結果，其應力水平為p′=25～500 kPa，相對密實度Dr=20%～81%。可以看出，用于低應力水平時，對高密實度砂土的預測偏差相對較大，但其離散程度也能控制在±30%以內。

圖10 砂土相對密實度預測值和實際值對比Fig.10 Comparison of the predicted and measured values of relative density of sand

5 結論

開展了高應力條件下的靜力觸探標定罐試驗和大變形有限元模擬。大變形模擬中采用的本構模型考慮了砂土內摩擦角和剪脹角隨累積塑性剪應變的變化。標定罐試驗和大變形模擬能夠很好地相互驗證，主要結論如下：

1）土體剪切模量取值嚴重影響錐尖阻力，建議采用式（2）修正彎曲元試驗結果。

2）靜力觸探的錐尖貫入到10倍的直徑時，錐尖阻力達到穩定。錐尖阻力隨相對密實度和圍壓的增大而增大。

3）結合有效應力水平在0.5～2.0 MPa 范圍內的標定罐試驗和大量大變形模擬，提出了由錐尖阻力預測相對密實度的經驗關系式。