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

孔壯壯,王棟,張民生,裴會(huì)敏

（中國(guó)海洋大學(xué) 山東省海洋環(huán)境地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室，山東 青島 266100）

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

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

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

筆者采用高圍壓標(biāo)定罐試驗(yàn)和ALE 大變形有限元方法，進(jìn)行砂土中靜力觸探貫入，建立高應(yīng)力條件下錐尖阻力與石英砂相對(duì)密實(shí)度的關(guān)系。大變形模擬中采用修正摩爾—庫(kù)倫模型（MMC）本構(gòu)模型描述砂土內(nèi)摩擦角和剪脹角隨累積塑性剪應(yīng)變的變化。

1 標(biāo)定罐試驗(yàn)

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

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

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

圖2 中國(guó)海洋大學(xué)研制的高壓標(biāo)定罐Fig.2 High-pressure calibration chamber developed in Ocean University of China

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

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

表1 標(biāo)定罐中砂樣相對(duì)密實(shí)度Table 1 Relative density of sand sample in calibration chamber

2 有限元模擬

2.1 本構(gòu)模型及參數(shù)標(biāo)定

砂土在剪切時(shí)表現(xiàn)出的應(yīng)變硬化或軟化行為取決于砂土的相對(duì)密實(shí)度與應(yīng)力水平。傳統(tǒng)摩爾—庫(kù)倫模型規(guī)定砂土內(nèi)摩擦角和剪脹角為定值，不能描述砂土的軟化和剪脹性。這里采用Hu 等[14]提出的MMC 模型，通過引進(jìn)內(nèi)摩擦角和剪脹角隨累積塑性剪應(yīng)變?chǔ)玫淖兓瘉?lái)描述砂土的軟化和剪脹性，其變化規(guī)律如圖3 所示，圖中φi、φp和φcv分別為初始、峰值和臨界內(nèi)摩擦角。因?yàn)槟枴獛?kù)倫類模型無(wú)法描述剪縮，所以假定當(dāng)γ<γ1時(shí)，剪脹角為0。此后剪脹角迅速線性增大，在γ=γ2時(shí)達(dá)到峰值剪脹角ψp；剪脹角保持峰值至γ=γ3，然后線性減小，直至在臨界應(yīng)變?chǔ)?處減小到0。當(dāng)砂樣表現(xiàn)為硬化和剪縮性時(shí)，MMC 模型實(shí)質(zhì)上退化為傳統(tǒng)的摩爾—庫(kù)倫模型：內(nèi)摩擦角保持為φcv，剪脹角始終為0。對(duì)于一般石英砂，按照Hu 等[14]和Zheng 等[15]的建議，取φi=φcv、γ1=1%、γ2=1.2%、γ3=5%和γ4=15%。

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

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

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

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

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

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

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

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

表2 三軸排水試驗(yàn)方案和結(jié)果Table 2 Scheme and results of drained triaxial tests

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

式中：dε1和dε3為主應(yīng)變?cè)隽?；三軸條件下k=2。

為了后續(xù)的大量變動(dòng)參數(shù)分析，需要建立砂土峰值內(nèi)摩擦角或峰值剪脹角與相對(duì)密實(shí)度、平均有效應(yīng)力之間的關(guān)系。參考經(jīng)典Bolton 公式[20]，擬合表2 的三軸試驗(yàn)結(jié)果，得到式（4）～式（6）。

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

2.2 有限元模型

采用通用有限元軟件ABAQUS 中的大變形模塊ALE 模擬CPT 試驗(yàn)。建立的ALE 模型如圖5 所示。參考Fan 等[21]和Mahmoodzadeh 等[22]的成果，在錐尖下人為設(shè)置直徑為D/20 的剛性小管，此剛性小管與CPT 一起移動(dòng)，從而防止最左邊的土體單元向?qū)ΨQ軸的內(nèi)部移動(dòng)，有利于保證土體網(wǎng)格調(diào)整的有效性，同時(shí)不會(huì)影響計(jì)算結(jié)果。

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

按照二維軸對(duì)稱計(jì)算，土體區(qū)域的寬和高分別是28D和42D。為了減少計(jì)算量，實(shí)施網(wǎng)格調(diào)整的ALE 區(qū)域?qū)?0D，根據(jù)前期試算，ALE 區(qū)域中單元大小設(shè)置為D/8，這樣能在保證計(jì)算精度的同時(shí)提高計(jì)算效率。ALE 區(qū)域需要保持節(jié)點(diǎn)位置調(diào)整后網(wǎng)格的質(zhì)量，以避免單元畸變，選用ABAQUS 內(nèi)置幾何加強(qiáng)形式的網(wǎng)格平滑方法[11]。

采用硬接觸算法模擬靜力觸探儀與土體之間的相互作用。假定二者交界面上的摩擦應(yīng)力為0，即忽略貫入過程中探桿上的側(cè)摩阻力，這是因?yàn)樯巴林刑綏U上的側(cè)摩阻力遠(yuǎn)小于錐尖阻力。

設(shè)置CPT 貫入速度為20 mm/s 進(jìn)行勻速貫入，貫入深度為0～20D，采用ALE 大變形方法避免網(wǎng)格扭曲，從而實(shí)現(xiàn)靜力觸探在砂土中的連續(xù)貫入。為了與標(biāo)定罐試驗(yàn)對(duì)比，對(duì)土體施加K0=1的壓力。

3 ALE 有限元模型的驗(yàn)證

3.1 試驗(yàn)和模擬結(jié)果對(duì)比

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

圖6 有限元模擬結(jié)果和標(biāo)定罐試驗(yàn)結(jié)果對(duì)比Fig.6 Comparison of finite element simulation results and calibration chamber test results

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

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

1）對(duì)于不同相對(duì)密實(shí)度的砂樣，當(dāng)貫入深度H>10D時(shí)，有限元計(jì)算得到的qc基本不再變化，存在一定程度的噪音波動(dòng)，但在可接受范圍內(nèi)。

2）對(duì)于Dr約為46%的砂樣，試驗(yàn)與有限元結(jié)果基本吻合。

3）對(duì)于Dr約為60%的砂樣，H<10D時(shí)的試驗(yàn)結(jié)果比有限元結(jié)果約高15%，但二者給出的H>10D時(shí)的穩(wěn)態(tài)錐尖阻力很接近。

3.2 剪切模量敏感性分析

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

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

4 相對(duì)密實(shí)度和錐尖阻力的關(guān)系

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

圖8 給出了0.8、1.3 MPa 圍壓下不同相對(duì)密實(shí)度的錐尖阻力曲線。對(duì)比圍壓為0.8、1.3 MPa 的砂樣，相對(duì)密實(shí)度80%和60%時(shí)的qc大約是40%時(shí)的1.22 倍和1.56 倍左右。當(dāng)Dr=40%、60%和80% 時(shí)，圍 壓1.3 MPa 時(shí)，qc約 為0.8 MPa 時(shí)的1.35 倍。圖9 給出了CPT 貫入過程中的土體位移矢量圖，大變形條件下土體位移場(chǎng)可近似認(rèn)為是破壞場(chǎng)?？梢钥闯觯F面附近土體的位移大致沿錐面外法線方向，CPT 探頭下的土體處于三軸壓狀態(tài)，被徑向擠出。

圖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]總結(jié)的相對(duì)密實(shí)度和錐尖阻力之間的關(guān)系為

式中：參數(shù)C1、C2與土體性質(zhì)有關(guān)；Qtn為歸一化錐尖阻 力；m值 取0.6。Ahmadi 等[7]考慮不同K0狀態(tài)對(duì)錐尖阻力的影響，建議m值取0.5。

取貫入深度H=10D作為穩(wěn)定時(shí)的錐尖阻力進(jìn)行分析。總結(jié)18 個(gè)工況的ALE 模擬結(jié)果發(fā)現(xiàn)，在p′=0.5～2.0 MPa 的高應(yīng)力水平下，仍可以按照式（7）由錐尖阻力預(yù)測(cè)相對(duì)密實(shí)度，但式（7）的參數(shù)擬合為：C1=0.58、C2=1.91、m=0.64。即

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

圖10 砂土相對(duì)密實(shí)度預(yù)測(cè)值和實(shí)際值對(duì)比Fig.10 Comparison of the predicted and measured values of relative density of sand

5 結(jié)論

開展了高應(yīng)力條件下的靜力觸探標(biāo)定罐試驗(yàn)和大變形有限元模擬。大變形模擬中采用的本構(gòu)模型考慮了砂土內(nèi)摩擦角和剪脹角隨累積塑性剪應(yīng)變的變化。標(biāo)定罐試驗(yàn)和大變形模擬能夠很好地相互驗(yàn)證，主要結(jié)論如下：

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

2）靜力觸探的錐尖貫入到10倍的直徑時(shí)，錐尖阻力達(dá)到穩(wěn)定。錐尖阻力隨相對(duì)密實(shí)度和圍壓的增大而增大。

3）結(jié)合有效應(yīng)力水平在0.5～2.0 MPa 范圍內(nèi)的標(biāo)定罐試驗(yàn)和大量大變形模擬，提出了由錐尖阻力預(yù)測(cè)相對(duì)密實(shí)度的經(jīng)驗(yàn)關(guān)系式。