離散元-有限差分跨尺度耦合下樁與土體穩定性細觀數值分析

高濤 陳云娟 閆甫 李艷龍 劉效智 王樂寧

摘要：為準確分析樁體宏觀變形與周圍土顆粒細觀力學行為，采用離散元-有限差分的跨尺度耦合進行樁和土體接觸過程中的穩定性分析，研究樁下沉過程中周圍土體的細觀變形、應力分布和樁體自身變化情況，通過FLAC3D建立樁和外部土體有限差分網格單元，對樁周圍側土體應用PFC3D離散元建立土顆粒微觀結構模型。研究結果表明：離散元與有限差分耦合方案能夠模擬樁周圍土體細觀力學行為;外部區域土體位移場呈包裹式C形對稱分布，應力場以豎直應力為主導，離散元土顆粒接觸力鏈網格以樹根狀向周圍遞減削弱分布;樁彈性單元體下沉時，應力場和位移場均有分層現象，應力場以壓應力為主，但存在局部拉應力即反彈現象。

關鍵詞：多尺度耦合;樁;位移場;應力場;接觸力鏈;數值模擬

中圖分類號：TP391.99;TU473.11

文獻標志碼：B

文章編號：1006-0871（2021）03-0024-07

DOI：10.13340/j.cae.2021.03.005

Abstract：To?accurately?analyze?the?macro?deformation?of?pile?and?the?micromechanical?behavior?of?surrounding?soil?particles，?discrete?element-finite?difference?multi-scale?coupling?is?used?to?study?the?stability?analysis?of?the?pile?and?soil?in?the?contact?process，?and?then?the?micro-deformation，?stress?distribution?and?changes?of?the?pile?itself?during?the?pile?sinking?process?is?studied.?The?finite?difference?mesh?of?the?pile?and?the?external?soil?is?established?by?FLAC3D.?The?soil?particle?microstructure?model?of?the?soil?around?the?pile?is?established?by?PFC3D?discrete?element.?The?results?show?that?the?coupling?scheme?of?discrete?element?and?finite?difference?can?be?used?to?simulate?the?micromechanical?behavior?of?the?soil?around?the?pile.?The?soil?displacement?field?in?the?outer?region?is?symmetrically?distributed?in?a?"wrapped?C"?shape，?and?the?stress?field?is?dominated?by?vertical?stress.?The?soil?particle?contact?force?chain?mesh?is?gradually?weakened?and?distributed?in?the?shape?of?tree?roots.?While?the?pile?elastic?element?body?sinks，?the?stress?field?and?displacement?field?are?both?stratified，?and?the?stress?field?is?mainly?compressive，?but?there?is?local?tensile?stress，?that?is?rebound?phenomenon.

Key?words：multi-scale?coupling;pile;displacement?field;stress?field;contact?force?chain;numerical?simulation

0?前?言

我國很多地域的地層土質狀況不穩定，由于雜填土成分復雜且持力層較深等因素，解決建筑施工地基層面問題常利用樁的徑長和自身強度高、穩定性好等優點，但施工過程主要依靠大量實踐衡量樁體的承載能力，樁下沉過程中底部土體的穩定性很難監測，樁土之間的載荷傳遞機理也不夠完善。[1-3]

為進一步了解樁土受力機理，已對其相互作用關系展開許多室內試驗，但大多數室內試驗采取比例控制模型，成本花費高、地質情況單一，造成試驗數據與現場監測差異較大。因此，國內外許多學者采用數值模擬技術進行此類研究，并取得一些進展。[4-6]周文龍等[7]應用有限差分數值軟件FLAC3D模擬樁體在復合地基情況下的自身沉降量和地表土體沉降，得到樁土在不同工況條件下應力比，為實際工程確定最優墊層厚度。邱瑞成等[8]對單樁進行參數敏感分析，改變數值模型中黏聚力和內摩擦角參數，結果認為黏聚力對計算結果影響較顯著，而內摩擦角參數對結果影響基本可忽略。季璇等[9]應用有限元MIDAS/GTS分析雙排樁支護過程中的變形情況，分析其水平位移和樁身彎矩，選取單排樁支護最佳優化方案。李新偉等[10]從樁對土層影響的角度出發，結合樁周圍土體應力場和位移場，分析施工過程中周圍土體的變形。劉文白等[11]應用流離散元PFC2D模擬樁上拔過程中的土體承載力，研究樁周圍土顆粒細觀力學變化，分析樁上拔過程中土體結構單元的變形情況，并與物理模擬試驗結果對比，認為顆粒流離散元軟件對工程模擬具有可靠性。

采用單一數值軟件對樁和樁周圍土體分析，樁和土體均為同一種計算原理，如有限元法、非連續性變形法、離散元法等，但實際上兩者結構具有很大差異性，樁自身強度大、變形小，而土體性質復雜且呈顆粒離散狀態，因此從不同角度研究樁和土體性質更為重要[12-13]。本文將FLAC3D和PFC3D進行耦合，將土體松散不穩定等特性定義為離散元體，樁體和周圍宏觀土層視為彈塑性材料，通過數值模擬分析樁體宏觀變形與周圍土顆粒細觀力學行為。

1?FLAC3D-PFC3D耦合理論基礎

1.1?耦合原理

在FLAC3D和PFC3D數值軟件相互耦合分析中，FLAC3D從宏觀上模擬連續區域內介質的力學行為，采用顯示拉格朗日算法和混合離散分區原理，運用動態運動方程進行求解計算，拉格朗日算法將大變形問題轉化為小變形本構關系，通過離散集成方法求解應力增量和不平衡力。PFC3D用于模擬細觀離散元顆粒的微觀力學行為，計算原理為顆粒間相互關系的離散元方法，顆粒間力的傳遞與顆粒變形通過二者接觸剛度（法向剛度和切向剛度）進行計算，接觸剛度由強度參數（摩擦系數、黏結和抗拉強度）控制[14]。兩者相互耦合發生在域的接觸邊界上，不同域之間的相互關系主要借助Socket?O/I接口進行數據傳輸與交換，FLAC3D-PFC3D耦合計算原理見圖1。

1.2?耦合方案分析

為探討FLAC3D和PFC3D的相互耦合作用，模擬連續區域內FLAC3D與離散介質PFC3D中球、簇的相互作用，在邊界區域創建獨立的耦合體域，耦合體域使球、簇等結構單元與實體單元相互作用，耦合體邊界采用WALL單元作為相互交換媒介。WALL單元是由點構成三角形面區域的組成體，三角形區域位置可由時間函數確定，以此作為耦合邊界上的媒介體。因此，在耦合分析中，WALL單元必須與實體單元表面協調一致。

初步確立耦合邊界介質后，通過媒介體傳遞和接受單元體力、速度等信息，形成區域通道內的耦合邏輯（通過墻面接觸力、力矩確定三角面頂點處的等效力系統）。確定耦合邏輯后定義相應的耦合工作，即取接觸力和力矩與壁面，確定等效的力系統在面上的頂點，這些力與剛度發送信息至媒介體區域，相關的網格節點和單元體接受媒介體傳遞的信息。數值模擬耦合中力的傳遞示意見圖2。

若沒有附加約束條件，式（14）和（15）構成的欠正定方程組可能存在無數組解集。控制方程并沒有提供明確的約束性質，因此欠正定方程組需得到一個特解，同時需要提供約束實體單元或結構單元的額外約束條件，通過該額外約束條件確立欠正定方程組的等效系統，等效系統需要的額外約束條件為

當點OP外推三角形區域的節點和結構單元觸發更新時，其網格點與節點能夠直接添加剛度屬性，即分別對FLAC3D于PFC3D設置實體單元的相關參數，通過三角形區域這一介質實現傳遞力與速度信息。

2?連續-離散耦合的樁與土體數值模擬

2.1?耦合模型建立

以FLAC3D軟件作為耦合計算方案的平臺，通過設置軟件Load?PFC打開耦合開關系統，建立的數值耦合模型主要包括樁和土體2個模塊，樁結構體采用FLAC3D實體單元，下部土體分為2個模塊。外部土體采用FLAC3D實體單元作為“基床”，該實體單元利用有限差分數值原理計算周圍大范圍土體，處理大變形，較真實反映材料實際的動態行為，有效模擬外部土體隨時間、空間變化產生的系統相關問題。與樁周圍接觸的土體采用PFC3D離散元結構體，離散元模型不僅不受變形限制，而且處理非連續性介質力學問題更精確，符合實際土顆粒物理力學性質。外部土體模型尺寸為15?m×10?m×8?m，土體單元采用FLAC3D中的Mohr-Coulomb計算模型，樁周圍局部土體區域范圍為5?m×5?m×5?m的?PFC3D顆粒（BALL），顆粒直徑分布區間為0.2～4.0?mm。為使樁的數值模型與施工現場更好吻合，樁局部直徑為0.6?m、高度為5?m，計算模型采用FLAC3D彈性體單元，整體系統數值計算模型見圖3。從多個不同單元體模型體現不同的物理性質關系，可有效模擬土體結構分離、變形等非連續特性，優化分析土體與樁的細觀力學行為。

2.2?模型參數標定

選取濟南開源路地區粉質黏土力學參數，按照現場試驗數據設置模型參數，在樁周圍局部土體單元設置線性接觸模型，對2種不同土體之間連接處設置WALL單元，作為FLAC3D有限差分網格與PFC3D離散元顆粒耦合的邊界。對于離散元計算模型的細觀參數，可通過一系列三軸數值試驗繪制不同圍壓強度下的莫爾強度包絡曲線，并進行多次標定計算，得到離散元模型細觀參數，樁周圍土體離散元模型材料參數見表1。FLAC3D區域內的材料參數根據現場測定得到，樁周圍土體的材料參數見表2，樁體的變形模量為4.86?MPa，泊松比為0.15。

3?數值模擬結果分析

采用離散元（PFC3D）與有限差分（FLAC3D）耦合模擬樁與土體的接觸變形關系，主要通過模擬樁接觸土體顆粒，分析樁身和周圍土顆粒的應力和變形情況。2種不同的單元體受力形態變化見圖4，可驗證PFC3D與FLAC3D相互耦合的可行性。PFC3D顆粒流單元的應力邊界傳遞至FLAC3D實體網格中，同一工況加載狀態下其應力云擴散性效果較為一致，PFC3D和FLAC3D的應力云梯度基本相同，均無網格單元體或球顆粒單元因應力傳遞問題導致應力點集中現象，二者在邊界處具備耦合相互傳遞性能，說明邊界連接體單元WALL對顆粒程序和有限差分網格的融接性好，同時也證明Socket?O/I接口進行計算數據傳輸與交換的可行性。

3.1?樁周圍外部有限差分模型土體變形特性

數值模擬加載后樁底周圍土體的水平位移場見圖5a。在樁加載下沉過程中，外部土體受到樁周圍土體傳遞應力產生一定的位移，土體接觸面附近水平方向產生4個較為顯著的聚集性位移場，最大位移量為3.3?mm。產生該現象主要原因是樁在下沉過程沖擊土體，造成樁周圍土體被迅速密實擠壓，產生沖擊力側壓現象，使水平位移云圖呈包裹式C形狀分布;坑底側土體在樁下沉擠壓作用下呈兩側包裹圓位移，但從縱向切片圖來看，底部較深處的一側土體位移并非連續的，該情況與現場土體位移監測較為符合。樁體下沉受載非一點集中力，在樁的豎直方向集中力與土體接觸過程中，產生中心線偏離引起樁身彎矩力效應，使下部土體位移不是等距變化。

數值模擬加載后樁底周圍土體豎直方向的位移場見圖5b。位移場基本呈現法向散狀分布，土體沉降量自中心向四周減緩，樁底中心處土體壓縮位移量約為四周的2倍左右，局部最大位移達19?mm左右，但從豎直方向位移切片圖看，局部最大位移區域范圍較小，說明越靠近樁底端部，土體壓縮位移越大，豎直位移分布云圖也可以更好地表明對土體建模劃分的科學性，將土體壓縮量較大區域用離散元顆粒表示更能貼合實際施工狀態。

樁周圍外部土體應力場變化云圖見圖6。

在外部上表層土體應力場中，水平應力場拉應力和壓應力對稱性分布，拉應力區域面積遠大于壓應力區域，整體以拉應力產生的土體變形為主導，但最大壓應力數值比拉應力約大0.8?kPa;豎直應力場中表面土體均為拉應力且無壓應力存在，主要原因是在樁豎向受載中，底部土體受壓引起上部土體產生被動受拉效應，而在樁底部的土體中，水平應力場呈現“兩壓一拉”現象，且壓應力區域遠遠超過拉應力區域范圍。這說明同等時間段內樁底中心處土體較周圍土體先壓縮固結，周圍土層固結量小，但中心豎向載荷不斷增加，同一深度應力場發生改變。

3.2?樁周圍離散元模型土體變形特性

土顆粒的水平位移場和豎直位移場見圖7。由此可以看出，位移量之間的傳遞具有較好的銜接性，土顆粒位移基本與外側土體位移趨勢相同。離散元土顆粒水平位移量集中在土體表面兩側，最大位移量局部相差2?mm左右，豎直方向位移呈現中間大兩邊小的形態，顆粒位移遞減量約35?mm。

為進一步研究土體位移過程中的細觀變化，提取離散元模型的土體速度矢量云圖，見圖8。

從2個方向的速度矢量場可知，水平方向的速度矢量運動軌跡基本與豎直位移場相同，豎直方向的顆粒運動軌跡與其水平位移場相同。速度矢量場這一現場可解釋在上部樁載作用下，樁載壓縮土體導致水平方向的土體顆粒向中間靠攏，而中心處土體在樁直接壓縮下先固結，因此周圍土顆粒向下位移產生拉應力。這種運動軌跡符合豎直方向位移場中間大兩邊小的規律。

3.3?離散元模型土體細觀變形特性

通過離散元顆粒模型構建樁周圍土體，能夠更好地從細觀角度分析顆粒間受力形態的變化。應用FLAC3D與PFC3D耦合技術可得到接觸模型顆粒細觀接觸力，見圖9。顆粒與顆粒間通過接觸力鏈傳遞應力和位移，接觸力鏈網格約35?000個，其接觸力鏈相互連接形成交叉網格呈樹根狀分布，在樁體與顆粒接觸軸向處，這些接觸力鏈網格分布密集并向周圍散射。

樁底部豎向載荷的傳遞，使樁底側顆粒出現強度高但數量少的接觸力鏈網格。隨著深度的增加，高強接觸力鏈網格轉化為數目更多的低強度力鏈網絡。在高強度接觸力鏈網格周圍，水平方向的接觸力基本以拉應力為主，但壓應力出現最高峰值82.41?N;豎直方向的接觸拉應力和壓應力數量相同，但豎向接觸力鏈應力值遠大于水平方向，接觸應力最大值為1.97?kPa。這說明樁體下沉中豎向應力起主導作用。

接觸力鏈網格分布與應力傳遞形態可揭示樁體下沉中土體的細觀特征，樁頂載荷通過樁身彈性體傳遞至樁底側土顆粒，傳遞的應力以豎向應力為主導作用，形成樹根狀接觸力鏈網格，其中心處土體固結速度大于周圍土體壓縮量，使周圍產生水平方向的接觸力鏈網格。這些水平方向接觸網格相對強度較低，從而對豎向接觸力鏈網格起到輔助支撐作用。若這些土顆粒無法對豎向接觸力鏈網格起約束作用，則會導致樁體下沉過程中出現力鏈坍塌錯落現象，甚至引起樁體宏觀坍塌破壞;若接觸力鏈對豎向力鏈約束程度過強，則會導致水平應力集中，使土體表面出現隆起，導致原有土層破壞。

3.4?樁有限差分網格單元變形分析

應用FLAC3D切片功能，將樁有限差分網格單元沿縱向作切片處理，得到樁單元不同方向的位移場和應力場，見圖10。樁身上半部分水平位移場與應力場保持相對均勻，主要原因是樁單元模型采用彈性體模型，該模型自身強度高、變形較為一致;下半部分樁體位移場與應力場不同，主要原因是土體產生的水平方向的束縛應力和樁自身摩擦阻力，使樁下半側水平位移場呈相對擠壓的位移，豎直位移場具有明顯的層次階段變化，且越靠近樁底端部，位移場分布區域越聚集于角落。

因此，針對樁身各個部位網格單元進行位移變化監測，結果見圖11。由此發現，越靠近樁底端部，水平位移越小，樁的豎直方向位移變化越大，說明彈性體樁單元符合胡克變形準則，彈性體變形基本以軸向為主導作用，水平位移量變化主要是由于土顆粒間的接觸力鏈約束引起的。

圖10b說明樁身應力場變化主要集中于樁底側單元區域，水平應力場呈軸對稱分布形態。水平方向土顆粒間受接觸力鏈約束作用，豎直方向的應力場上部為拉應力作用，但拉應力數值較下部壓應力小（壓應力最大值為238?kPa，是拉應力近20倍），樁身上半部的拉應力主要是自身在壓縮受力過程中發生的應力回彈現象，因此上部出現較小的均勻拉應力分布。

4?結?論

通過離散元-有限差分尺度耦合對樁與土體的穩定性進行細觀數值分析，得到以下結論。

（1）應用離散元和有限差分原理進行耦合計算，可實現PFC3D與FLAC3D數值模型耦合計算，其數值模型的位移場、應力場等變形傳遞在不同屬性單元中具有良好的融合性，為跨尺度耦合數值模擬提供參考。

（2）樁周圍外部土體位移場呈包裹式C形對稱性分布，豎直應力場起主導作用，且約為水平方向應力7倍;內部離散元土顆粒位移場與外部土體基本一致，土顆粒的速度矢量場可解釋土顆粒在固結作用下的運動變化，其細觀接觸力鏈網格分布呈樹根狀并向周圍遞減，接觸應力密集處為樁底部土體附近。

（3）對樁彈性體單元作切片處理可知，樁下沉應力場和位移場有分層現象，應力表現以壓應力為主但有較小的拉應力反彈現象，水平位移隨樁身增加而增加，而豎直位移則相反。

參考文獻：

[1]?劉漢龍，?趙明華.?地基處理研究進展[J].?土木工程學報，?2016，?49（1）：96-115.

[2]?STAHL?M，?KONIETZKY?H.?Discrete?element?simulation?of?ballast?and?gravel?under?special?consideration?of?grain-shape，?grain-size?and?relative?density[J].?Granular?Matter，?2011，?13（4）：417-428.?DOI：10.1007/s10035-010-0239-y.

[3]?朱怡，?魏祥，?趙登堅.?樁體貫入大變形計算分析[J].?建筑結構，?2018，?48（S2）：838-841.?DOI：10.19701/j.jzjg.2018.S2.171.

[4]?寇海磊，?張明義，?張吉坤.?層狀粘性土及砂土地基中靜力壓樁連續貫入的數值模擬[J].?工程力學，?2012，?29（12）：175-181.

[5]?HENKE?S.?Influence?of?pile?installation?on?adjacent?structures[J].?International?Journal?for?Numerical?and?Analytical?Methods?in?Geomechanics，?2010，?34（11）：1191-1210.?DOI：10.1002/nag.859.

[6]?GAO?G，?MEGUID?M?A.?Effect?of?particle?shape?on?the?response?of?geogrid-reinforced?systems：insights?from?3D?discrete?element?analysis[J].?Geotextiles?and?Geomembranes，?2018，?46（6）：685-698.?DOI：10.1016/j.geotexmem.2018.07.001.

[7]?周文龍，?劉玉柱，?王小鵬.?應用FLAC3D分析確定碎石樁復合地基的褥墊層厚度[J].?施工技術，?2014，?43（S1）：76-78.

[8]?邱瑞成，?艾健森.?基于FLAC3D的單樁靜載模擬樁土接觸面參數敏感性研究[J].?路基工程，?2019（2）：164-169.?DOI：10.13379/j.issn.1003-8825.2019.02.34.

[9]?季璇，?高全臣，?楊卓，?等.?H型雙排樁支護結構三維數值模擬分析與工程應用[J].?施工技術，?2016，?45（7）：70-73.?DOI：10.7672/sgjs2016070070.

[10]?李新偉，?曾昌.?樁基礎施工對地層影響的三維數值模擬[J].?四川建筑，?2015，?35（1）：95-97.?DOI：10.3969/j.issn.1007-8983.2015.01.036.

[11]?劉文白，?周健.?上拔荷載作用下樁的顆粒流數值模擬[J].?巖土工程學報，?2004，?26（4）：516-521.?DOI：10.3321/j.issn：1000-4548.2004.04.018.

[12]?譚鑫，?馮龍健，?趙明華，?等.?剛性基礎下筋箍碎石樁復合地基樁土應力比計算模型[J].?中國公路學報，?2019，?32（9）：42-50.?DOI：10.19721/j.cnki.1001-7372.2019.09.004.

[13]?傅金陽，?謝佳偉，?房雅楠，?等.?EPB盾構開挖面穩定性的PFC-FLAC耦合分析[J].?華中科技大學學報（自然科學版），?2019，?47（5）：116-121.?DOI：10.13245/j.hust.190522.

[14]?石崇，?張強，?王盛年.顆粒流（PFC5.0）數值模擬技術及應用[M].?北京：中國建筑工業出版社，?2018.

（編輯?武曉英）