基于離散元的南海軟黏土剪切變形模擬

程 升，朱超祁，單紅仙,2*，劉曉磊,2，賈永剛,2

(1.中國海洋大學(xué)山東省海洋環(huán)境地質(zhì)工程重點(diǎn)實(shí)驗(yàn)室,青島 266100;2.青島海洋科學(xué)與技術(shù)國家實(shí)驗(yàn)室海洋地質(zhì)過程與環(huán)境功能實(shí)驗(yàn)室,青島 266061)

隨著中國對南海深海油氣資源以及天然氣水合物的勘測及開采的進(jìn)行，對于南海的工程地質(zhì)性質(zhì)研究也逐漸增多。軟黏土的剪切變形特性是決定海洋構(gòu)筑物、管線能否安全存在的重要因素，因此亟需開展這方面的研究。朱超祁等[1-2]對南海北部陸坡淺表層的沉積物的物理力學(xué)性質(zhì)進(jìn)行了初步研究，分析得到陸架坡折帶淺表層沉積物以粉質(zhì)黏土和黏土質(zhì)粉砂為主，工程地質(zhì)特性較差。蔣明鏡等[3]研究了南海東北部陸坡軟黏土的微觀結(jié)構(gòu)和力學(xué)特性，土體呈開放式絮凝結(jié)構(gòu)，具有一定的結(jié)構(gòu)性；任玉賓等[4-5]針對南海西部深海沉積物的物理性質(zhì)、全流動循環(huán)軟化特性和微觀結(jié)構(gòu)進(jìn)行了研究，認(rèn)為生物硅顆粒的破碎會加劇土體軟化程度同；郭興森等[6]考慮海底低溫環(huán)境，利用南海北部軟黏土研究得到低溫環(huán)境下泥流的剪應(yīng)力與表觀黏度會顯著提高；馬永政等[7]采用高壓三軸試驗(yàn)儀模擬了海洋土取樣應(yīng)力釋放的過程，分析得到高壓卸載過程會使得土體體積膨脹?？傮w而言，對于南海沉積物的變形特性的研究較少，且對南海沉積物的細(xì)觀變形特性的研究尚未見報(bào)道。因此基于室內(nèi)試驗(yàn)獲取南海神狐海域原狀軟黏土的物理力學(xué)參數(shù)，結(jié)合離散元數(shù)值模擬來研究其細(xì)觀變形特性。

離散元數(shù)值模擬廣泛應(yīng)用于巖土體的工程力學(xué)性質(zhì)研究中[8-10]，采用的是顆粒流數(shù)值模擬法。Anandarajah[11]基于傾斜顆粒之間的雙層排斥力的研究，且考慮了顆粒之間存在物理化學(xué)力，將每個粒子劃分為多個互連的離散元素，作為整體分析的一部分模擬了粒子的彎曲；Zhang等[12]為了準(zhǔn)確地模擬和分析黏性土的行為，考慮毛細(xì)管和土壤顆粒間水的存在引起的動態(tài)黏性力，建立了顆粒間平行黏結(jié)的黏性土DEM力學(xué)模型。前人對黏土的研究主要集中在微觀機(jī)理性方面。近年來采用顆粒流法對于黏土的宏觀變形特性研究增多。李力等[13]驗(yàn)證了黏性沉積物顆粒流數(shù)值模型和數(shù)值模擬方法的正確性，顆粒間黏結(jié)強(qiáng)度決定沉積物宏觀剪切強(qiáng)度；寧孝梁[14]開展黏性土三軸試驗(yàn)的細(xì)觀模擬，得到了與室內(nèi)土工試驗(yàn)變化趨勢一致的應(yīng)力-應(yīng)變曲線，同時研究了不同細(xì)觀參數(shù)對黏性土材料宏觀力學(xué)響應(yīng)；邢煒杰等[15]采用顆粒流離散元法對黏性土的三軸剪切進(jìn)行了初步模擬，分析了細(xì)觀參數(shù)變化對應(yīng)力應(yīng)變曲線的影響。

基于此，針對在南海神狐海域取得的大量原狀樣開展室內(nèi)試驗(yàn)，獲取典型軟黏土的綜合物理力學(xué)參數(shù)，這對南海工程地質(zhì)研究具有重要意義。通過建立雙軸壓縮模型研究軟黏土的剪切變形特性，是對進(jìn)一步通過數(shù)值模擬手段研究南海地質(zhì)災(zāi)害的初步探索工作。

1 研究區(qū)概況

神狐海域位于中國南海北部陸坡中段的神狐暗沙東南海域附近，屬于典型的海底峽谷區(qū)，涵蓋了2017年水合物試采海域[16]。該區(qū)域水深200～300 m位上陸坡和陸架的過渡區(qū)，地形平緩，水深超過500 m后地形變陡。水深等值線總體呈NE-SW分布，區(qū)域位置圖如圖1中(a)所示。圖1(b)中紅色圓點(diǎn)為過去五年取得的柱狀樣品站位點(diǎn)(115.0°～115.4°E，19.7°～20.1°N)，水深為300～2 000 m。開展固結(jié)不排水試驗(yàn)的柱狀樣長為5.25 m，水深1 532 m，取其中1.0～2.0 m進(jìn)行三軸試驗(yàn)及其他室內(nèi)土工試驗(yàn)。

圖1 研究區(qū)及取樣站位Fig.1 Research sea area and sampling stations

2 試驗(yàn)原理及方法

2.1 室內(nèi)土工試驗(yàn)

圖2 南海神狐海域原狀沉積物樣品及GDS三軸試驗(yàn)儀Fig.2 Undisturbed sediment sample from shenhu area,South China Sea and GDS triaxial instrument

南海神狐海域原狀樣品如圖2(a)所示，將樣品運(yùn)回實(shí)驗(yàn)室后存放在4 ℃環(huán)境下。樣品含水率超過100%，因此在制備樣品時要盡量避免擾動樣品。使用GDS應(yīng)力路徑靜三軸儀[圖2(b)]進(jìn)行固結(jié)不排水試驗(yàn)，有效圍壓分別為50、100、150 kPa。由于軟黏土的滲透系數(shù)較低，因此采用CO2飽和、反壓飽和聯(lián)合的方法，反壓均為100 kPa，試驗(yàn)過程中進(jìn)行B值檢測確保試驗(yàn)飽和度達(dá)到規(guī)范要求，剪切速率為0.06 mm/min。具體試驗(yàn)過程嚴(yán)格按照《土工試驗(yàn)方法標(biāo)準(zhǔn)》(GB/T 50123—1999)。通過開展系列的室內(nèi)土工試驗(yàn)，包括含水率、滲透率、顆粒級配、液塑限、比重、壓縮系數(shù)和固結(jié)不排水試驗(yàn)等，能夠系統(tǒng)描述南海軟黏土的物理力學(xué)指標(biāo)。

2.2 離散元數(shù)值模擬

采用顆粒流理論基于離散元方法，從細(xì)觀角度研究材料的變化來模擬介質(zhì)的受力及變形[13]。顆粒和墻體需要滿足牛頓第二定律和力-位移定律，通過每一時步更新顆粒、墻的位置，形成新的接觸和接觸力條件[17]。

黏土在微觀結(jié)構(gòu)上呈片狀、針狀等，粒徑在微米級，研究著重于從細(xì)觀角度研究軟黏土的變形特性，因此選擇顆粒粒徑0.2～0.5 mm。模型箱尺寸為39.1 mm×80 mm，與實(shí)際三軸試驗(yàn)土體的橫截面積相同，且在剪切階段保證體積不發(fā)生改變。顆粒之間為接觸黏結(jié)模型，該模型的兩顆粒之間能夠傳遞力，但是不能傳遞力矩。通過對比其他研究者的模型基本參數(shù)[18-20]，確定基本模型參數(shù)如表1所示。

表1 模型參數(shù)Table 1 Parameters of particle flow model

為了使得試樣剪切過程的可視化，通過設(shè)置一系列的測量圓來監(jiān)測試樣內(nèi)的孔隙率、應(yīng)變率、配位數(shù)、滑動接觸比等細(xì)觀參量的變化，試驗(yàn)建立的測量圓如圖3所示，一共有405個測量圓。這與Gu等[19]在雙軸壓縮試驗(yàn)中采用的1 200個測量圓，用于研究顆粒物質(zhì)剪切帶的形成演變過程是一致的。

圖3 測量圓分布Fig.3 Distribution of measure circles

通過室內(nèi)試驗(yàn)獲取軟黏土的物理力學(xué)參數(shù)，建立顆粒流數(shù)值模擬模型。顆粒流數(shù)值模擬的主要過程包括生成樣品、樣品固結(jié)和雙軸剪切。通過施加初始圍壓，模擬原始地應(yīng)力條件。樣品來自海底3～4 m，因此設(shè)置初始地應(yīng)力為30 kPa。通過伺服機(jī)制能夠使得圍壓達(dá)到預(yù)定值，且不平衡力小于10-5。模型設(shè)置的剪切速度為0.1 mm/s，而顆粒流模擬中上下墻移動1 mm需要運(yùn)算105～106時間步，因此能夠保證在準(zhǔn)靜態(tài)下進(jìn)行剪切試驗(yàn)，這個建模過程與Gu等[19]的研究一致。

3 試驗(yàn)結(jié)果

3.1 室內(nèi)試驗(yàn)

通過分析沉積物的粒度得到沉積物主要為細(xì)粒黏土類，顆粒粒徑為0.001～0.1 mm，粒徑小于0.01 mm的累積質(zhì)量分?jǐn)?shù)約為85%，顆粒粒徑級配如圖4所示。

圖4 深度100～110 cm處沉積物粒徑級配累計(jì)曲線Fig.4 Accumulation curve of sediment particle size grading at depths of 100～110 cm

通過開展室內(nèi)土工試驗(yàn)得到典型軟黏土的物理力學(xué)特性如表2所示。其中含水率、滲透系數(shù)、顆粒級配、比重、液塑限和壓縮系數(shù)是一個綜合統(tǒng)計(jì)分析結(jié)果，樣品數(shù)據(jù)涉及多個站點(diǎn)，指標(biāo)的獲取方法如表1所示。而雙軸壓縮試驗(yàn)選用的固結(jié)不排水試驗(yàn)結(jié)果選自2018年其中一個樣品的試驗(yàn)數(shù)據(jù)，具體指黏聚力(Ccu)、內(nèi)摩擦角(φ)指標(biāo)。通過開展綜合性的室內(nèi)試驗(yàn)得到，南海神狐沉積物具有的典型特征為高含水率、高液塑限、低密度、高壓縮性、低強(qiáng)度。這一結(jié)論與朱超祁等[1-2]、蔣明鏡等[3]和郭興森等[6]的試驗(yàn)結(jié)論相同，但是對這些參數(shù)進(jìn)行綜合性的統(tǒng)計(jì)分析，對其他學(xué)者在該區(qū)域開展工程地質(zhì)研究提供了必要的基礎(chǔ)參數(shù)。

表2 神狐海域軟黏土物理力學(xué)參數(shù)Table 2 Physical and mechanics properties of the sediment in the Shenhu area

通過對原狀樣進(jìn)行固結(jié)不排水試驗(yàn)，得到在不同圍壓下的應(yīng)力-應(yīng)變曲線，如圖4所示。

圖5為數(shù)值模擬結(jié)果與室內(nèi)三軸試驗(yàn)的結(jié)果對比圖。通過分析室內(nèi)試驗(yàn)數(shù)據(jù)可得，在50、100、150 kPa下的初始剪切模量分別為2.6、5.0、6.7 MPa，峰值強(qiáng)度分別為36、72.5、104 kPa。在軸向應(yīng)變達(dá)到3%時候，偏應(yīng)力均達(dá)到最大值，隨后都呈現(xiàn)應(yīng)變硬化型(圖5)。由于端部約束的存在，試樣出現(xiàn)微鼓脹現(xiàn)象。

圖5 數(shù)值模擬結(jié)果與室內(nèi)三軸試驗(yàn)結(jié)果對比Fig.5 Comparison of numerical simulation and indoor triaxial test results

3.2 數(shù)值模擬

目前,對于細(xì)觀參數(shù)的選取沒有統(tǒng)一的方法，一些學(xué)者提出了細(xì)觀參數(shù)的標(biāo)定方法[18-20]，但都缺乏通用性，因此多數(shù)學(xué)者還是采用試錯法來確定模型所需的細(xì)觀參數(shù)[13-14,19]。采用試錯法，即通過不斷地調(diào)整細(xì)觀參數(shù)使得雙軸模擬結(jié)果與室內(nèi)固結(jié)不排水試驗(yàn)相吻合，來確定適用南海典型軟黏土的細(xì)觀參數(shù)組合，具體如表3所示。

表3 接觸黏結(jié)模型參數(shù)Table 3 Parameter of linear contact bond model

顆粒間的接觸模量為6.5 MPa，初始剪切模量分別為2.4、4.3、5.7 MPa。顆粒間摩擦系數(shù)為0.3，顆粒間剛度比為1.25。

采用接觸黏結(jié)模型得到的結(jié)果與實(shí)際相比總體比較吻合(圖5)。在圍壓為50 kPa情況下，應(yīng)變在前5%時吻合度很好；軸向應(yīng)變?yōu)?%～7.5%時偏應(yīng)力呈現(xiàn)增大趨勢，最大誤差率在18.4%。可能是低圍壓且等體積條件下，顆粒之間的相對轉(zhuǎn)動使得偏應(yīng)力呈現(xiàn)增大。在圍壓為100 kPa下，模擬結(jié)果總體小于實(shí)際結(jié)果，在軸向應(yīng)變5.5%～6.5%，偏應(yīng)力出現(xiàn)一個低谷階段；圍壓為150 kPa下的擬合度相對最好，同時軸向應(yīng)變?yōu)?%～7%時偏應(yīng)力也呈現(xiàn)一個低谷。因此在軸向應(yīng)變達(dá)到5%時，達(dá)到峰值應(yīng)力階段，顆粒之間的相對轉(zhuǎn)動使得低圍壓下的應(yīng)力呈現(xiàn)一個增大趨勢，在圍壓100、150 kPa下，出現(xiàn)一個應(yīng)力減小的階段。綜合來看，選取的細(xì)觀參數(shù)組合能夠有效表征南海神狐海域軟黏土。結(jié)合數(shù)值模擬過程中測量圓監(jiān)測到的數(shù)據(jù)，對剪切過程中的變形特性進(jìn)行分析。

4 剪切變形分析

4.1 局部孔隙率演化分析

對于分析土體變形而言，孔隙率的變化尤為關(guān)鍵。模擬的是固結(jié)不排水試驗(yàn)，因此總的孔隙率變化為0。但可通過局部孔隙率變化，來分析剪切變形過程。

圖6 軸向應(yīng)變分別為0.0%、5.3%和10%時的孔隙率分布Fig.6 The distribution of porosity when the axial strain was 0.0%,5.3% and 10%,respectively

如圖6所示，在初始時刻，孔隙率分布并不均勻，這也與實(shí)際土體情況一致；在應(yīng)變達(dá)到5.3%時候，土體發(fā)生破壞，形成一條與X軸負(fù)方向?yàn)?9°的大孔隙率帶；應(yīng)變?yōu)?0%時，可形成“X”形的大孔隙率區(qū)帶，與X軸負(fù)方向的夾角為57°，隨著應(yīng)變增大，剪切帶傾角呈現(xiàn)增大趨勢。根據(jù)摩爾庫倫理論計(jì)算剪切帶傾角公式為π/4+φ/2，其中φ為內(nèi)摩擦角，由表2的內(nèi)摩擦角數(shù)據(jù)可以計(jì)算理論剪切帶傾角為60°，與試驗(yàn)結(jié)果比較接近。這與Gu等[19]的研究結(jié)果一致，即理論預(yù)測剪切帶傾角比實(shí)際大，但是總體上比較吻合。

4.2 顆粒接觸組構(gòu)分析

土的接觸組構(gòu)可量化土體細(xì)觀特性對宏觀力學(xué)響應(yīng)的影響。Rothenburg等[24]建議采用傅里葉函數(shù)分別表示接觸法向分布函數(shù)E(θ)、法向平均接觸力函數(shù)fn(θ)和切向接觸力函數(shù)ft(θ)，具體計(jì)算公式如式(1)所示：

(1)

式(1)中:θ為顆粒間接觸法向的角度；f0為所有法向接觸力的平均值；θa、θn、θt分別為顆粒間接觸法向、法向接觸力和切向接觸力各向異性分布主方向的角度；θn、θt分別為顆粒法向接觸法向和切向各項(xiàng)異性的二階主方向；an、at分別為二階平均法向和切向力各向異性系數(shù)。

數(shù)值試驗(yàn)得到軸向應(yīng)變分別為0%、5.3%和10%時試樣的法向接觸力和切向接觸力分布圖，粗實(shí)線為式(1)中傅里葉函數(shù)擬合的曲線，如圖6～圖8所示。

如圖7所示，初始時接觸法向分布比較均勻，說明樣品生成比較均勻。在應(yīng)變5.3%對應(yīng)峰值階段和應(yīng)變?yōu)?0%試驗(yàn)結(jié)束階段，接觸法向的主軸方向發(fā)生較小的旋轉(zhuǎn)，形狀只有微小的改變。由于重力和軸向力的加載，這也與實(shí)際相吻合。

如圖8所示，初始狀態(tài)，法向接觸力分布比較均勻；應(yīng)變?yōu)?.3%時，主應(yīng)力角度為80°～90°，發(fā)生小角度偏轉(zhuǎn)；應(yīng)變?yōu)?0%時，主應(yīng)力角度為90°～100°，發(fā)生主應(yīng)力的旋轉(zhuǎn)。

如圖9所示，初始狀態(tài)，切向接觸力分布比較均勻；應(yīng)變?yōu)?.3%時，主應(yīng)力角度為-40°～40°，發(fā)生小角度偏轉(zhuǎn)；應(yīng)變?yōu)?0%時切向應(yīng)力狀態(tài)基本與5.3%時一致，未發(fā)生旋轉(zhuǎn)。

在應(yīng)變?yōu)?時，傅里葉函數(shù)擬合效果較好；隨著應(yīng)變增大，對于法向接觸力的擬合效果較好。該函數(shù)未能反應(yīng)出接觸法向在90°附近出現(xiàn)的一個峰值，對于切向接觸力的擬合未能反應(yīng)出切向力的主優(yōu)勢方向。

圖7 接觸法向分布Fig.7 Normal direction distribute of contact

圖9 切向接觸力分布Fig.9 Distribute of shear contact force

4.3 局部應(yīng)變率演化分析

室內(nèi)三軸試驗(yàn)測得的應(yīng)變?yōu)檩S向應(yīng)變，未能測得某一狀態(tài)下的應(yīng)變場。通過測量圓監(jiān)測到每一狀態(tài)下的總應(yīng)變場來分析剪切變形過程，如圖10所示。圖10中應(yīng)變?yōu)?0%的應(yīng)變場中標(biāo)注了3個測量圓的位置，分別為C213、C217、C220號測量圓。

10%應(yīng)變場內(nèi)的三個藍(lán)色原點(diǎn)為三個測量圓，從左到右分別為C213、C217、C220圖10 應(yīng)變分別為0.8%、5.3%和10%時的應(yīng)變場Fig.10 Strain field at strains of 0.8%,5.3%,and 10%,respectively

某一點(diǎn)的應(yīng)變狀態(tài)可以用一個二階的應(yīng)變張量表示，因此通過設(shè)置測量圓監(jiān)測到每一個位置的X方向和Y方向的應(yīng)變量來求取總的應(yīng)變場。應(yīng)變?yōu)?.8%時，開始出現(xiàn)連續(xù)應(yīng)變剪切帶；應(yīng)變?yōu)?.3%時出現(xiàn)一條與X軸負(fù)方向?yàn)?7°的連續(xù)應(yīng)變帶；應(yīng)變?yōu)?0%時，出現(xiàn)兩個“X”形連續(xù)應(yīng)變帶，以及很多次要的剪切帶，這些次要的剪切帶在室內(nèi)試驗(yàn)中是無法觀測或者測量到的，這與Gu等[19]對顆粒土模擬中發(fā)現(xiàn)在主要剪切帶附近還存在其他次要剪切帶是一致的。隨著應(yīng)變增加，從單一應(yīng)變帶演化為共軛應(yīng)變帶，這也與試驗(yàn)得到的剪切帶演化形成機(jī)制一致。

如圖11所示，分別在剪切帶外、剪切帶內(nèi)和剪切帶附近設(shè)置測量圓C213、C217和C220，用于定量分析局部應(yīng)變與軸向應(yīng)變的關(guān)系。在剪切帶外的C213區(qū)域，只發(fā)生了較小的應(yīng)變；在剪切帶內(nèi)(C217)初始應(yīng)變速率較大，在達(dá)到4%應(yīng)變后不再繼續(xù)發(fā)生變形；而在剪切帶邊緣區(qū)域，隨著應(yīng)變的增大，C220測得的應(yīng)變在一直增大，且增長速率基本不變。一般的室內(nèi)三軸試驗(yàn)得到的是整體的應(yīng)力——應(yīng)變曲線，而離散元能夠得到任意一點(diǎn)的應(yīng)變變化情況。Cheng等[21]將數(shù)字圖像測量三軸儀和離散元相結(jié)合起來，消除了端部約束的影響，這對于研究剪切帶形成機(jī)理及影響因素有重要作用。

圖11 C213、C217、C220號測量圓測得的總應(yīng)變曲線Fig.11 Total strain curve measured by measuring circle C213,C217,C220

4.4 剪切帶演化分析

應(yīng)變率在0.8%時，對應(yīng)從彈性變形轉(zhuǎn)變?yōu)閺椝苄宰冃危?.3%對應(yīng)土體達(dá)到臨界破壞狀態(tài)，偏應(yīng)力接近最大值，同時分析了應(yīng)變?yōu)?.3%和10%下的顆粒旋轉(zhuǎn)以及顆粒位移，從細(xì)觀角度解釋剪切帶的形成過程。

圖12 不同軸向應(yīng)變時的顆粒累計(jì)轉(zhuǎn)動量Fig.12 Total amount of rotation of the particles at different axial strain

圖13 不同軸向應(yīng)變時的顆粒位移Fig.13 The particles displacement at different axial strain

如圖12、圖13所示，將顆粒累計(jì)轉(zhuǎn)動量與顆粒累計(jì)位移圖對比分析可知，在較小應(yīng)變(0.8%)下，土體是從邊緣局部位置開始發(fā)生破壞，總體上顆粒是由上下向中間靠攏。原因是初始條件下，顆粒之間孔隙較多，在荷載作用下兩端土體向中間積聚，兩側(cè)土體位移相對較小。在應(yīng)變?yōu)?.3%時，顆粒累計(jì)轉(zhuǎn)動量與顆粒位移圖都顯示存在兩條比較明顯的旋轉(zhuǎn)帶，即可能形成剪切帶的區(qū)域。在應(yīng)變?yōu)?0%時候，形成X形的剪切帶，這與室內(nèi)試驗(yàn)以及Batiste等[25]的研究結(jié)果一致。

5 結(jié)論

基于南海神狐海域原狀沉積物開展室內(nèi)土工試驗(yàn)，獲取典型軟黏土的物理力學(xué)特性，基于此開展離散元模擬研究，分析軟黏土剪切變形特性。主要結(jié)論如下。

(1)通過試錯法標(biāo)定了適用于南海軟黏土的細(xì)觀參數(shù)組合：顆粒間接觸模量Emod=6.5 MPa、剛度比Kratio=1.25、摩擦系數(shù)fric=0.3、臨界法向阻尼比dp_nratio=0.5、抗剪強(qiáng)度cb_shears=1.5 kPa、抗拉強(qiáng)度cb_tens為0 kPa。

(2)在軸向應(yīng)變?yōu)?.8%時開始出現(xiàn)局部應(yīng)變集中現(xiàn)象，軸向應(yīng)變達(dá)到5.3%時出現(xiàn)一條與X軸反方向?yàn)?9°的連續(xù)應(yīng)變帶；隨應(yīng)變增加10%時，剪切帶傾角增加到57°。

(3)不同部位土體的變形不同。剪切帶外土體發(fā)生較小的應(yīng)變；剪切帶內(nèi)土體初始應(yīng)變率較大，在應(yīng)變達(dá)到4%后不再繼續(xù)發(fā)生變形；而剪切帶邊緣區(qū)域土體應(yīng)變持續(xù)增大，且增長速率基本不變。