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

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

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

(1)

式中：qc——錐尖阻力；

su——黏性土不排水強(qiáng)度；

σv0——上覆土壓力。

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

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

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

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

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

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

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

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

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

(4) 不需要預(yù)先假設(shè)破壞模式；

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

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

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

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

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

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

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

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

(2)對(duì)二次函數(shù)在當(dāng)前單元積分點(diǎn)處求導(dǎo)，得到一個(gè)線性分布模式計(jì)算式Φtrial。

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

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

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

2.1 基本參數(shù)

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

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

圓錐的貫入過程分兩步進(jìn)行，第一個(gè)Dynamic/Explicit分析步確定圓錐貫入前的應(yīng)力狀態(tài)，為了克服Dynamic/Explicit在計(jì)算初始加速度時(shí)因?yàn)椴话跏紤?yīng)力對(duì)解答造成波動(dòng)，將對(duì)應(yīng)的分析步時(shí)間取一個(gè)極小值，并固定所有節(jié)點(diǎn)的所有自由度。將第二個(gè)Dynamic/Explicit分析步時(shí)間總長設(shè)置為36 s，CPT的標(biāo)準(zhǔn)貫入速度為2 cm/s，36 s意味著圓錐貫入了0.72 m，即20倍的探頭直徑。

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

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

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

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

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

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

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

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

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

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

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

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

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

Nc=1.988lnIr-0.554

(2)

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

2.2.3塑性區(qū)的變化

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

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

圖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 圓錐與探桿粗糙程度的影響

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

(3)

對(duì)于給定的Mises屈服強(qiáng)度cu，su的表達(dá)式為：

(4)

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

(5)

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

大量研究表明[14,16～17]，探桿的摩擦力對(duì)錐尖阻力影響甚微，因?yàn)閷?duì)于Mises材料而言，探桿周圍剪應(yīng)力的增加不可能使得錐尖附近的平均應(yīng)力增加。因此，本文只研究圓錐的摩擦力對(duì)錐形因子的影響。取αc=0，1/4，1/2，3/4，1等5種情況分析圓錐摩擦系數(shù)對(duì)錐形因子的影響，這里αc=0表示圓錐完全光滑，αc=1表示圓錐完全粗糙。對(duì)于均值各項(xiàng)同性，剛度指數(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é)果表明，試驗(yàn)前必須測(cè)量圓錐的摩擦系數(shù)，因?yàn)槭欠窨紤]錐面的摩擦系數(shù)會(huì)使得模擬結(jié)果差別高達(dá)25%。

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

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

(7)

式中：σv0——原位豎向應(yīng)力；

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

σh0——原位水平應(yīng)力。

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

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

圖10 錐形因子隨土體原位應(yīng)力狀態(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ù)、初始應(yīng)力狀態(tài)參數(shù)的關(guān)系式為：

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

(8)

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

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

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

3 結(jié)論

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

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

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

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

[1] 李小雪,彭祖武,滕偉福.CPTU 試驗(yàn)在土體動(dòng)力學(xué)參數(shù)計(jì)算中的應(yīng)用研究[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] 陳維家,汪火旺.靜力觸探中土層界面效應(yīng)試驗(yàn)研究[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ù)理論與工程應(yīng)用[M].北京:科學(xué)出版社,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.