大應變靜力觸探數(shù)值模擬及錐形因子影響因素分析

肖先波，李 波，王 婷，邱 敏

(1.湖州職業(yè)技術(shù)學院，浙江 湖州 313000；2.長江科學院，湖北 武漢 430010；3.中國地質(zhì)大學(武漢)工程學院，湖北 武漢 430074)

圓錐靜力觸探(CPT)作為主要的原位測試技術(shù)，已成功進行了大量工程應用，其實用性、簡便性及經(jīng)濟性得到了工程界的認可[1]，尤其適用于遠距離輸水工程、高速公路、鐵路這種線形分布的大型工程[2]。利用靜力觸探試驗可以直接獲取土體錐尖阻力、側(cè)壁摩阻力和孔隙水壓力等參數(shù)，由于側(cè)壁摩阻力受摩擦筒材料特性影響較大，而孔隙水壓力測量值經(jīng)常出現(xiàn)很大的誤差，所以對測量成果的應用以錐尖阻力為主[3]。貫入黏性土和其它細粒土時，對錐尖阻力的研究主要集中在錐形因子Nc(承載力系數(shù))上，即：

(1)

式中：qc——錐尖阻力；

su——黏性土不排水強度；

σv0——上覆土壓力。

由于土體大應變和非線性等原因，對錐尖阻力的嚴格分析很困難，通常只能做一些近似理論分析。目前，主要的近似理論方法有承載力理論、孔穴擴張理論、穩(wěn)態(tài)變形理論、應變路徑法和運動點位錯方法等[4]。

在發(fā)現(xiàn)彈塑性介質(zhì)中開挖深孔所需壓力與相同條件下擴張相同體積的孔所需壓力之間存在一定的比例關(guān)系后，BISHOP等[5]首次建立起孔穴擴張和靜力觸探之間的聯(lián)系。YU等[6]研究表明，由于孔穴擴張理論考慮了土體的剛度、壓縮性(或剪脹性)和貫入過程中水平應力的減少等因素的影響，孔穴擴張理論得出的預測結(jié)果較承載力理論更為精確。這就使得利用孔穴擴張理論研究靜力觸探機理在更多的學者當中流行起來[7～9]。利用孔穴擴張理論，采用理想彈塑性模型，可以得出簡單的錐形因子精確解，但不能考慮貫入速率、土體應力歷史的影響，不能模擬垂直穩(wěn)態(tài)貫入的連續(xù)性。事實上，錐形因子還會受到圓錐粗糙程度、原位應力狀態(tài)的影響。要考慮以上因素的影響，就需要采用有限元方法。有限元方法分析錐尖阻力主要有小變形和大變形兩種分析方法。傳統(tǒng)的小變形方法容易產(chǎn)生較大誤差，因為在探頭的貫入過程中，探桿附近的側(cè)向應力增大了，探桿周圍的應力變化將導致錐尖阻力比小應變情況下更大。除此之外，傳統(tǒng)的小應變有限元分析不能產(chǎn)生必要的殘余應力場，所以不能得到一個合適的極限錐尖阻力。一般來說，當應變超過10%時就不再滿足小變形理論的變形條件。特別是當土體承受較大荷載且發(fā)生較大的壓縮變形時，土體的變形具有明顯的非線性、大變形的特點。研究表明，在靜力觸探過程中，探頭周圍土體的平均應變常超過10%，有的剪切應變甚至達到40%。因此在利用有限元方法研究靜力觸探貫入過程時，經(jīng)典的彈塑性小變形理論會產(chǎn)生較大的誤差，引入大變形理論很有必要。

LU等[10]將黏土視為理想彈塑性材料，服從Tresca屈服準則，采用RITSS (網(wǎng)格重劃分和插值技術(shù))大位移方法量化了剛度指數(shù)、原位應力狀態(tài)和錐土摩擦系數(shù)對錐形因子的影響，并與應變路徑法進行了對比。Tolooiyan等采用Drucker-Prager模型及任意拉格朗日-歐拉(ALE)技術(shù)模擬了砂土中CPT貫入，并將模擬結(jié)果與實測錐尖阻力展開對比，獲得了較好的一致性，證明了ALE技術(shù)的適用性。本文采用有限元商業(yè)軟件Abaqus模擬探頭在不排水黏土中的貫入過程，不分析孔壓場的變化，采用Mises屈服準則，相關(guān)聯(lián)流動法則模擬不排水黏土的彈塑性行為。采用 ALE網(wǎng)格劃分技術(shù)對探頭的貫入進行復雜的大應變分析，并將大應變分析結(jié)果與傳統(tǒng)的小應變分析結(jié)果進行了對比。研究了探頭貫入達到穩(wěn)定狀態(tài)時土體剛度指數(shù)、原位應力狀態(tài)和錐尖粗糙程度對錐形因子與塑性區(qū)半徑的影響，并由此建立錐形因子表達式，利用該表達式反演不排水強度，驗證該表達式的精確性。

1 ALE方法在Abaqus中的實現(xiàn)形式

1.1 有限元法方法的優(yōu)點

對于巖土貫入問題的求解，相比于其他方法，有限元方法有如下優(yōu)點[11]，這些優(yōu)點與孔穴擴張理論的優(yōu)點正好對應起來，便于對二者進行分析對比。

(1) 可以模擬土的剛度和壓縮性；

(2) 可以定義初始應力場；

(3) 可以采用不同的本構(gòu)模型；

(4) 不需要預先假設破壞模式；

(5) 可以準確模擬探頭的幾何形狀。

1.2 ALE方法的實現(xiàn)形式

盡管有限元方法具有許多優(yōu)點，在利用有限元方法模擬靜力觸探貫入過程中存在幾個特別需要解決的問題：(1) 錐尖附加的高應變區(qū)，存在大應變、大變形；(2) 高度扭曲的網(wǎng)格；(3) 邊界條件明顯改變。而解決以上問題，需要用到大變形分析理論。Abaqus提供了一系列自適應技術(shù)處理以上問題，其中最主要的是任意拉格朗日—歐拉自適應網(wǎng)格劃分方法(ALE)。在ALE描述中，計算網(wǎng)格可以在空間中以任意的形式運動，即可以獨立于物質(zhì)坐標系和空間坐標系運動。這樣通過規(guī)定合適的網(wǎng)格運動形式就可以準確地描述物體的移動界面，并維持單元的合理形狀，即便網(wǎng)格發(fā)生極大扭曲變形，ALE方法也能在分析過程中保證高質(zhì)量的網(wǎng)格。

ALE將純拉格朗日分析和純歐拉分析的優(yōu)點結(jié)合起來。一個完整的ALE分析包括兩個步驟：(1) 建立一個新的網(wǎng)格；(2) 將舊網(wǎng)格的解答與狀態(tài)變量傳輸?shù)叫戮W(wǎng)格上。ALE通過反復掃掠自適應域以創(chuàng)造一個平滑的網(wǎng)格。在網(wǎng)格掃掠過程中，分析域中節(jié)點根據(jù)相鄰節(jié)點和網(wǎng)格中心位置按照一定的網(wǎng)格光滑化方法重新定位。一般采用體積光滑方法減少扭曲，并保持原有網(wǎng)格的變化梯度。如圖1所示，M點的新位置由周邊4個單元中心點C的體積加權(quán)平均值來確定。加權(quán)量會試圖將M點朝遠離C1、指向C3的方向推進，從而減少4個單元的扭曲。

圖1 網(wǎng)格掃描過程中的一個節(jié)點的重新定位Fig.1 Relocation of a node during a mesh sweep

自適應網(wǎng)格ALE方法引入對流動量平衡方程與質(zhì)量守恒方程，這就導致了獨立的網(wǎng)格和物質(zhì)運動，Abaqus就是通過耦合物質(zhì)運動和網(wǎng)格運動來解答修正后的方程。方程網(wǎng)格建立以后，需要將解答和狀態(tài)變量從舊網(wǎng)格傳輸?shù)叫戮W(wǎng)格上。Abaqus提供兩種傳輸方法，一是Van Leer提出的二階方法，二是donor cell提出的一階方法。其中二階方法對準靜態(tài)問題、瞬時動力沖擊都比較適用。由于靜力觸探測試被當作準靜態(tài)問題而且在模擬探頭的貫入過程中需要進行頻繁自適應網(wǎng)格劃分，所以采用二階方法。其傳輸方式如圖2所示，主要可以分為4步。

圖2 二階傳輸方法Fig.2 Second order advection

(1)根據(jù)當前單元和2個相鄰單元上積分點的變量值Φ(單元內(nèi)為恒定值)構(gòu)造一個二次函數(shù)。

(2)對二次函數(shù)在當前單元積分點處求導，得到一個線性分布模式計算式Φtrial。

(3)較小二次函數(shù)斜率，直到其最小值和最大值落到兩側(cè)單元原始值范圍之內(nèi)為止。

(4)所有單元線性分布模式確定以后，在每個新單元上對變量積分，處于新單元的體積之后即為映射支后的新變量值。

2 黏土中圓錐靜力貫入的數(shù)值模擬

2.1 基本參數(shù)

CPT的貫入可以簡化為軸對稱問題，選擇直徑為0.036 m、錐角為60°的常規(guī)CPT探頭作為貫入圓錐。圓錐貫入深度為20倍的圓錐直徑，即0.72 m。取錐尖附近1.5 m×1.5 m作為計算域以保證模型邊界在塑性變形區(qū)之外，模型左側(cè)0.2 m×1.5 m范圍內(nèi)使用ALE方法對網(wǎng)格進行重劃分。

本研究對不排水黏土采用理想彈塑性模型，不分析孔壓場的變化，屈服面采用Mises模型，流動法則采用相關(guān)聯(lián)的法則。土體彈性模量保持為2.98 MPa，泊松比為0.49(模擬不排水狀態(tài))，屈服強度cu=10 kPa，土體飽和重度取20 kN/m3，土體剛度指數(shù)Ir=G/su。

圓錐的貫入過程分兩步進行，第一個Dynamic/Explicit分析步確定圓錐貫入前的應力狀態(tài)，為了克服Dynamic/Explicit在計算初始加速度時因為不包含初始應力對解答造成波動，將對應的分析步時間取一個極小值，并固定所有節(jié)點的所有自由度。將第二個Dynamic/Explicit分析步時間總長設置為36 s，CPT的標準貫入速度為2 cm/s，36 s意味著圓錐貫入了0.72 m，即20倍的探頭直徑。

在定義荷載、邊界條件環(huán)節(jié)，通過定義模型所受豎向應力模擬不排水黏土所受上覆壓力，本次模擬中豎向應力取100 kPa。在貫入過程中，約束模型區(qū)域底部兩個方向的位移和左、右邊界的水平位移，模型初始網(wǎng)格和邊界條件如圖3所示。

圖3 模型邊界條件與受力狀態(tài)Fig.3 Model boundary conditions and loading

圓錐貫入過程中，將圓錐的外表面設定為主控面，模型左上角的面定義為從屬面，表明錐面可以穿透土體，而土體不能進入錐面。在Tangential/Behavior項定義接觸面切向特性，通過定義圓錐表面的摩擦系數(shù)分析貫入過程中探頭的粗糙程度對土體力學行為的影響。

ALE重劃分技術(shù)成功的效率很大程度取決于初始網(wǎng)格的形態(tài)，對于普通區(qū)域，采用四邊形單元形狀，結(jié)構(gòu)化網(wǎng)格劃分技術(shù)即可實現(xiàn)。對于ALE重劃分區(qū)，考慮到圓錐貫入的變形形態(tài)，還需要對網(wǎng)格控制做一些調(diào)整，增加網(wǎng)格密度，確保圓錐貫入相當深的距離后網(wǎng)格質(zhì)量仍能令人滿意。圓錐貫入前后模型網(wǎng)格如圖4所示。

2.2 土體剛度指數(shù)的影響

2.2.1穩(wěn)定狀態(tài)的辨識

圖4 圓錐貫入前后模型網(wǎng)格Fig.4 Model grid before and after cone penetration

圖5 錐形因子隨貫入深度變化散點圖Fig.5 Scatter diagram for the cone factor vs penetration depth

傳統(tǒng)的小應變模擬很難辨別深層貫入過程中的極限承載力，常用方法是設定一個初始剛度，一旦增量剛度減小到一個特定的數(shù)值后即認為達到穩(wěn)定狀態(tài)，該方法會大大低估錐尖阻力。對于大應變ALE模擬，由于土體是均質(zhì)的，圓錐貫入一定深度后錐尖的擠土效應就會達到穩(wěn)定狀態(tài)，即錐尖阻力達到穩(wěn)定狀態(tài)。這里，錐尖阻力qc以錐尖所受豎向應力表示，通過改變不排水強度確定土體剛度指數(shù)Ir。圖5給出了錐形因子Nc(歸一化錐尖阻力)隨圓錐歸一化貫入深度h/D的關(guān)系。如圖5所示，剛開始歸一化錐尖阻力均隨著貫入深度的增加而增大，貫入一定深度后，傳統(tǒng)小應變歸一化錐尖阻力停止增長，大應變歸一化錐尖阻力趨于穩(wěn)定。在傳統(tǒng)小應變模擬分析中，由于網(wǎng)格發(fā)生過度扭曲，導致小應變曲線與大應變曲線發(fā)生分離，最終停止增長。盡管歸一化小應變錐尖阻力和大應變錐尖阻力都隨著剛度指數(shù)的增大而增大，二者的極限錐尖阻力還是有很大的差別。Ir=100，小應變歸一化錐尖阻力可取7.48，大應變歸一化錐尖阻力穩(wěn)定在8.6左右，大應變歸一化錐尖阻力比小應變高15%；Ir=300，小應變歸一化錐尖阻力可取8.25，大應變歸一化錐尖阻力穩(wěn)定在10.8左右，大應變歸一化錐尖阻力比小應變高31%；Ir=500，小應變歸一化錐尖阻力可取8.43，大應變歸一化錐尖阻力穩(wěn)定在11.8左右，大應變歸一化錐尖阻力比小應變高40%。

模擬結(jié)果表明，傳統(tǒng)小應變分析不適合計算圓錐貫入過程中的極限承載，大應變極限承載力隨著貫入深度的增大而增大，在貫入深度超過12D以后，3條大應變曲線都達到穩(wěn)定狀態(tài)，超過12D以后承載力的增加可以忽略不計。

2.2.2土體剛度指數(shù)的影響

以往的研究表明，剛度指數(shù)在50～500之間時，剛度指數(shù)對錐形因子影響較大，本次研究中土體剛度指數(shù)Ir分別取100，300，500。繪制錐形因子Nc隨剛度指數(shù)Ir變化的對數(shù)圖(圖6)。由圖6可知，可以將3個點的連線擬合成1條直線，直線的對數(shù)方程為：

圖6 錐形因子隨土體剛度指數(shù)變化曲線Fig.6 Variation curves for the cone factor vs rigidity index

Nc=1.988lnIr-0.554

(2)

這樣就得到了原位應力狀態(tài)各向同性均質(zhì)土體中、錐尖完全光滑條件下錐形因子的關(guān)系式，為了驗證該關(guān)系式的正確性，將其他學者獲得的錐形因子擬合直線繪制于圖6內(nèi)，獲得了較好的一致性，該擬合直線與Abu-Farsakh[12]等人提出擬合直線變化規(guī)律一致，且處于中間值。

2.2.3塑性區(qū)的變化

圓錐貫入土體時，最初只有錐尖周圍出現(xiàn)塑性區(qū)，隨著貫入深度的增加，塑性區(qū)向探桿側(cè)壁擴展。塑性區(qū)大小受到很多因素的影響，如錐頭的粗糙程度、土體剛度和初始應力各向異性等。YU等[13]研究表明，錐頭和探桿的粗糙程度對塑性區(qū)的發(fā)展微乎其微，而原位應力各向異性對塑性區(qū)影響較小，當土體剛度指數(shù)較大的，各向同性原位應力會使塑性區(qū)域豎向范圍最大化、徑向范圍最小化。塑性區(qū)的擴展主要取決于土體剛度，剛度指數(shù)的增加會使得塑性區(qū)發(fā)生巨大變化，本文主要對此展開研究。

模擬錐尖完全粗糙、探桿完全光滑，土體剛度指數(shù)Ir分別為100，300，500條件下，圓錐在原位應力各向同性均質(zhì)土體中貫入達到穩(wěn)定狀態(tài)時的塑性區(qū)(圖7)。如圖7所示，網(wǎng)格的形態(tài)都比較好，體現(xiàn)了ALE方法的效果?？梢杂^察到塑性區(qū)隨著剛度指數(shù)的增大而增大，最大塑性應變主要集中在錐尖肩部與探桿一側(cè)一定范圍內(nèi)，因為錐尖肩部有明顯轉(zhuǎn)折，容易出現(xiàn)應力集中現(xiàn)象，所以塑性應變最大。探桿一側(cè)的塑性區(qū)幾乎與探桿平行，且塑性應變由近及遠逐漸減弱，該結(jié)果與TEH等[14]、YU 等[15]模擬結(jié)果一致?？拷F尖的塑性區(qū)在徑向的擴張類似于一個半球形。圖8顯示了歸一化塑性區(qū)半徑隨土體剛度指數(shù)變化曲線，其中柱形孔與球形孔擴張曲線由VESIC計算得來。對比可知，LIYANAPATHIRANA[11]的模擬曲線、LU的模擬曲線位于柱形孔擴張曲線和球形孔擴張曲線之間。LU的模擬曲線和本文的模擬曲線變化規(guī)律與球形孔擴張一致，三者的大小比較接近，表明塑性區(qū)的徑向擴張更接近于球形孔擴。

圖7 穩(wěn)定狀態(tài)下塑性區(qū)大小(Ir分別為100，300，500)Fig.7 Size of the plastic zone in the steady state (Ir=100,300,500)

圖8 歸一化塑性區(qū)半徑隨土體剛度指數(shù)變化曲線Fig.8 Variation curves for the normalized plastic zone radius vs rigidity index

2.3 圓錐與探桿粗糙程度的影響

研究圓錐和探桿粗糙度度對錐形因子的影響時，由于顯式算法可以應對復雜的接觸條件，所以不需要設置特別的接觸單元。圓錐—土體接觸面和探桿—土體接觸面的摩擦狀態(tài)特性是一樣的，即假設土體附著于錐體上，直到接觸面上的剪應力達到預設的土體剪切強度土體才能發(fā)生自由滑移。Mises材料最大剪應力表達式如下：

(3)

對于給定的Mises屈服強度cu，su的表達式為：

(4)

接觸面的摩擦力與錐體的粗糙程度成正比，因此，為了模擬錐體不同的粗糙程度，需要改變接觸面最大允許剪應力，接觸面上最大允許剪應的分數(shù)形式(摩擦系數(shù))如下：

(5)

式中：τint——錐體表面與土體接觸面上的允許剪應力。

大量研究表明[14,16～17]，探桿的摩擦力對錐尖阻力影響甚微，因為對于Mises材料而言，探桿周圍剪應力的增加不可能使得錐尖附近的平均應力增加。因此，本文只研究圓錐的摩擦力對錐形因子的影響。取αc=0，1/4，1/2，3/4，1等5種情況分析圓錐摩擦系數(shù)對錐形因子的影響，這里αc=0表示圓錐完全光滑，αc=1表示圓錐完全粗糙。對于均值各項同性，剛度指數(shù)分別為100，300，500的土體，錐形因子隨圓錐摩擦系數(shù)變化的模擬結(jié)果如圖9所示。錐形因子隨圓錐摩擦系數(shù)的增大而增大，近似一種線性關(guān)系。通過最佳線性擬合，取擬合直線斜率的平均值，得到錐形因子與剛度指數(shù)、摩擦系數(shù)的關(guān)系式：

圖9 錐形因子隨圓錐摩擦系數(shù)變化圖Fig.9 Variation diagram for the cone factor vs cone friction coefficient

(6)

模擬結(jié)果表明，試驗前必須測量圓錐的摩擦系數(shù)，因為是否考慮錐面的摩擦系數(shù)會使得模擬結(jié)果差別高達25%。

2.4 原位應力狀態(tài)的影響

土體的原位應力狀態(tài)對錐形因子也會有很大的影響，TEH等[13]定義原位應力狀態(tài)參數(shù)Δ為：

(7)

式中：σv0——原位豎向應力；

主要教學內(nèi)容：通過參加工業(yè)工程專業(yè)的相關(guān)競賽或完成校級、省級大學生創(chuàng)新項目，從而提升學生對所學知識和能力的深度和廣度。教學地點：實驗室和企業(yè)實習基地等；教學課時：3～12個月。

σh0——原位水平應力。

原位應力狀態(tài)參數(shù)反映土體的原位應力各向異性，其變化范圍在-1～1之間。當Δ=1時，原位豎向應力最大；當Δ=-1時，原位水平應力最大；當Δ=0時，原位豎向應力等于原位水平應力，即原位應力各向同性。

本次模擬中，對于特定的Ir，保持σv0=100 kPa不變，通過改變σh0輸入值大小來實現(xiàn)Δ的變化。繪制不同土體剛度指數(shù)和圓錐摩擦系數(shù)組合情況下，錐形因子隨原位應力狀態(tài)參數(shù)變化的散點圖(圖10)。

圖10 錐形因子隨土體原位應力狀態(tài)參數(shù)變化圖Fig.10 Variation diagram for the cone factor vs in situ stress state parameter

如圖10所示，6條擬合曲線的斜率都比較接近，分別為-1.85，-1.7，-1.94，-1.99，-1.85，-1.7，其最大差值為0.29，表明擬合曲線的斜率受土體剛度指數(shù)和圓錐摩擦系數(shù)的影響較小，取曲線的平均斜率1.84，得到錐形因子與剛度指數(shù)、摩擦系數(shù)、初始應力狀態(tài)參數(shù)的關(guān)系式為：

Nc=1.988lnIr-0.554+2.47αc-1.84Δ

(8)

該公式適用范圍為土體剛度指數(shù)Ir=0～500，圓錐摩擦系數(shù)αc=0～1，原位應力狀態(tài)參數(shù)Δ=-1～1。為了驗證本公式的精確性，取Ir=200，αc=0.5，Δ=0.5模擬圓錐在穩(wěn)態(tài)貫入時的錐尖阻力，利用飽和黏土不排水強度理論公式 (1) 和錐形因子關(guān)系式(8)反演不排水抗剪強度su，將反演的su值與數(shù)值模擬設定的su值展開對比。結(jié)果顯示，su最大差值為3.2%，標準差為2.1%。

最后，將由數(shù)值模擬關(guān)系式與孔擴張理論表達式得到的錐形因子繪制于表1內(nèi)。其中，LADANYI等[8]采用的是球形孔擴張方法，BALIGH[9]和YU[18]采用的是柱形孔擴張方法，S表示圓錐完全光滑，R表示圓錐完全粗糙。對比可知，利用數(shù)值模擬方法得到的錐形因子比較接近，體現(xiàn)了大應變有限元法的優(yōu)勢；對比孔穴擴張方法和大應變有限元法，發(fā)現(xiàn)無論錐頭完全光滑還是完全粗糙，利用大應變有限元方法得到的錐形因子處于LADANYI等、BALIGH和YU值之間，更接近于LADANYI 等的計算結(jié)果。表明圓錐貫入過程中，土體變形處于球形孔擴張和柱形孔擴張之間，更接近球形孔擴張，這與圖8顯示的結(jié)果一致。

表1 不同計算方法下的錐形因子Table 1 Cone factor from different calculation methods

3 結(jié)論

(1)利用傳統(tǒng)小應變分析無法得到圓錐貫入過程中的極限承載力，大應變ALE分析中，極限承載力隨著貫入深度的增大而增大，并在到達一定深度后達到穩(wěn)定狀態(tài)。

(2)圓錐貫入達到穩(wěn)定狀態(tài)后，塑性區(qū)隨著剛度指數(shù)的增大而增大,錐尖周圍土體塑性區(qū)的徑向擴張?zhí)幱谥慰讛U張和球形孔擴張之間，更接近于球形孔擴張。

(3)錐形因子隨土體剛度指數(shù)、錐尖粗糙程度的增大而增大，隨土體原位應力狀態(tài)參數(shù)的增大而減小。采用大應變ALE模擬方法獲得的錐形因子表達式能夠量化土體剛度指數(shù)、原位應力狀態(tài)和錐尖粗糙程度的影響，具有較高的精確度。

隨著靜力觸探技術(shù)在海洋工程中的應用，研究對象由普通黏土轉(zhuǎn)為軟黏土和海相黏土。此時，影響錐形因子的主要因素不再是土體剛度指數(shù)、原位應力狀態(tài)和錐尖粗糙程度，而是土體靈敏度和應變速率等。如何將更多的因素考慮進來，并以一定的科學方法處理這些因素，是未來研究的重點。

[1] 李小雪,彭祖武,滕偉福.CPTU 試驗在土體動力學參數(shù)計算中的應用研究[J].水文地質(zhì)工程地質(zhì),2013,40(3): 84-88.[LI X X,PENG Z W,TNEG W F. Analyses of calculating soil kinetic parameters with CPTU testing[J].Hydrogeology & Engineering Geology, 2013,40(3):84-88. (in Chinese)]

[2] 陳維家,汪火旺.靜力觸探中土層界面效應試驗研究[J]. 水文地質(zhì)工程地質(zhì),2003,30(2):25-27.[CHEN W J, WANG H W. Research of soil-layer interface effect in cone penetrafion test [J]. Hydrogeology & Engineering Geology, 2003,30(2):25-27. (in Chinese)]

[3] TOLOOIYAN A, GAVIN K. Modelling the cone penetration test in sand using cavity expansion and arbitrary Lagrangian Eulerian finite element methods [J]. Computers and Geotechnics, 2011, 38(4):482-490.

[4] 劉松玉,蔡國軍,童立元,等.現(xiàn)代多功能CPTU:技術(shù)理論與工程應用[M].北京:科學出版社,2013:31.[LIU S Y, CAI G J, TONG L Y,etal. Modern multifunctional CPTU: technology theory and engineering application [M]. Beijing: Science Press, 2013:31. (in Chinese)]

[5] BISHOP R F, HILL R, MOTT N F. The theory of indentation and hardness tests [J]. Proceedings of the Physical Society, 1945, 57(3).

[6] YU H S, MITCHELI J K. Analysis of cone resistance: review of methods[J]. Journal of Geotechnical and Geoenvironmental Engineering, ASCE, 1998, 124(2):140-149.

[7] VESIC A S. Design of pile foundations[R]. National Cooperative Highway Research Program, Synthesis of Highway Practice 42, Transportation Research Board, National Research Council, 1977, Washington, D.C.

[8] LADANYI B, JOHNSTON G H. Behavior of circular footings and plate anchors embedded in permafrost[J]. Canadian Geotechnical Journal, 1974, 11:531-553.

[9] BALIGH M M. Theory of deep static cone penetration resistance[R]. Report No. R75-56, Department of Civil and Environmental Engineering, 1975, MIT, USA.

[10] LU Q, RANDOLPH M F, HU Y,etal. A numerical study of cone penetration in clay [J]. Geotechnique, 2004, 54(4):257-67.

[11] LIYANAPATHIRANA D S. Arbitrary Lagrangian Eulerian based finite element analysis of cone penetration in soft clay[J]. Computers and Geotechnics, 2009,36: 851-860.

[12] Abu-Farsakh M, Tumay M, Voyiadjis G. Numerical parametric study of piezocone penetration test in clays[J]. Int J Geomech, 2003, 3(2):170-181.

[13] YU H S, HERRMANN L R, BOULANGER R W. Analysis of steady state cone penetration in clay[J]. J Geotech Geoenviron Eng, ASCE 2000, 126(7):594-605.

[14] TEH C I, HOULSBY G T. An analytical study of the cone penetration test in clay [J]. Geotechnique, 1991,41(1):17-34.

[15] YU H S. Cavity expansion methods in geomechanics[M]. Rotterdam: Balkema, 2000.

[16] HERRMAAN L R. Finite element analysis of contact problems [J]. J Engng Mech, ASCE 1978, 104:1043-1059.

[17] LU Q, HU Y, RANDOLPH M F. FE analysis for T-bar and ball penetration in cohesive soil [C]//Proc 10th Int Offshore and Polar Engineering Conf. ISOPE 2000, 2:617-623.

[18] David Durban,Norman A Fleck. Singular plastic fields in steady penetration of a rigid cone [J]. J Appl Mech, 1993, 59(4):706-710.