基于有限元遺址土應力-應變關系的方法研究1)

岳建偉 黃軒嘉 趙麗敏 孔慶梅 陳 穎 王自法

*(河南大學土木建筑學院，河南開封 475004)

?(開封市不可移動文物修復與安全評價重點實驗室，河南開封 475004)

土遺址作為古代建筑物的遺留部分，受組成材料的限制及外部因素的影響，易造成失水收縮，力學性能降低，最終產生劣化現象[1]。這些裂隙不僅會對土遺址的工程特性產生影響[2]，還會對遺址的穩定性及長期保存構成嚴重威脅，這也是土遺址保護領域公認的工程“疑難雜癥”。

學者們對遺址土的宏觀力學性能進行了一定的研究[3]，而細觀尺度下的物理力學性能研究成果較少，遺址土的劣化機理缺乏系統的理論基礎支撐。復雜力學行為數值仿真研究方向的飛速發展，能夠有效避免物理試驗中人為及環境因素的影響，滿足工程領域研究遺址土基礎特性的目的，而建立遺址土組成材料的本構關系及數值模型是滿足工程應用的首要工作，也是影響數值模擬結果的重要因素。

巖土材料的力學性能極易受工程環境的影響，由此產生的裂縫顯著影響土壤的力學特性[4]。馮永等[5]基于室內試驗及有限元數值模擬研究了土的細觀多相變形機理及其宏、細觀變形特性之間的相互關系；郭夢圓等[6]通過對土進行單軸抗壓強度實驗，并對S–M 模型中各參數的意義進行修正，得到土在不同溫度及載荷加載等級下的應力-應變曲線和改進的S–M 蠕變顯式模型。徐安全等[7]通過壓縮試驗研究了壓縮模量的變化規律并利用MATLAB 三維擬合得出壓縮模量的表達式，更直觀表現壓縮模量的變化規律。

上述研究成果在一定程度上推動了宏、細觀土樣力學特性研究的發展，但鑒于土樣自身的復雜性及多變性[8]，建立模型時，均對顆粒間的膠結物進行了一定的簡化，同時也簡化了不同環境作用下土樣顆粒間膠結物的性能變化情況，這會使土樣和土樣各組成成分的應力變化規律存在偏差[9]。由于不能在既有土遺址本體上進行大量的取樣以及進行反復的試驗測試，本工作基于材料類比思想和廣義相似理論，對比分析遺址土細觀尺度和混凝土宏觀尺度組成的特點，獲得材料宏觀和細觀間的內在聯系，推導顆粒間膠結物的計算公式，依據遺址土掃描電子顯微鏡(scanning electron microscope， SEM) 圖像及少量的力學測試結果，以探索細觀尺度下遺址土顆粒間膠結物的力學性能的確定方法。

1 研究思路

細觀尺度下(10-6～10-9m) 的遺址土由土顆粒、顆粒間的膠結物、氣泡三部分組成，而宏觀尺度下(10-2m) 混凝土也由石子、石子間的砂漿、氣泡三部分組成。石子和土顆粒為彈性體，受力過程中假設為處于彈性狀態，相對于砂漿而言，土顆粒間膠結物的力學性能不清楚，本文以混凝土、水泥砂漿及遺址土的力學性能參數為依據，基于材料特性類比思想，根據相似理論及Hooke 定律推導出膠結材料力學性能參數[10-11]，為研究遺址土的力學性能和劣化奠定基礎。

試驗用土取自河南省開封市州橋遺址，按照《土工試驗方法標準》(GB/T 50123–2019) 進行土工試驗，土樣基本物理指標液限為32.8%，塑限為18.7%，塑性指數為14.1，孔隙率為41%，天然干密度為1.67 g/cm3，天然含水率為12.3%，不同含水率下遺址土和混凝土的應力- 應變曲線如圖1 和圖2所示。

圖1 遺址土應力-應變曲線

圖2 混凝土應力-應變曲線[12]

混凝土與遺址土兩者的力學性能滿足相似理論的第一定理，主要有如下相似特征：

(1)分析圖1 和圖2 發現，含水率升高對兩類材料的力學性能均有影響，并呈現出含水率越大，應力及應變值越小的情況。

(2) 圖3 為宏觀尺度混凝土材料示意圖(10-2m)，圖4 為細觀尺度遺址土材料示意圖(10-6～10-9m)，其中，ΔL為對應材料沿F方向的單位變形量。分析圖 3 和圖 4 發現，宏觀尺度下混凝土整體性能是粗骨料及水泥砂漿共同作用的結果[13]，在細觀尺度下遺址土整體性能是土顆粒及膠結物共同作用的結果[14]，兩類材料在組分劃分和構成規律、結構受力特點、力學特性等方面均較為相似。

圖3 宏觀尺度混凝土示意圖

圖4 細觀尺度遺址土示意圖

(3) 雖然兩類材料的力學性質在數值上存在差距，但是兩者的應力-應變歷程規律較為相似，均屬于準脆性材料，因為這類材料在載荷作用下，當自身應力、應變和彈性模量超過某一臨界值后，材料會出現明顯軟化現象[15]，在實際工程中，對于遺址土和混凝土兩類材料，當本體出現微裂縫后，其力學特性和工程質量會發生改變，為了使兩類材料的力學性能發展規律更為接近，因此本文主要針對材料發生破壞前的性能進行分析，不再考慮后期的材料軟化階段，也滿足實際工程需求。

(4)遺址土和混凝土結構在受力破壞后具有分形特性[16-17]。

2 模型參數的確定

基于遺址土細觀尺度和混凝土宏觀尺度組成的特點，通過材料類比思想和廣義相似理論，推導顆粒間膠結物的計算公式，依據遺址土SEM 圖像及少量的力學測試結果，通過下列方法探索細觀尺度下遺址土顆粒間膠結物的力學性能的確定方法。

混凝土宏觀尺度平衡方程與遺址土細觀尺度平衡方程可表示為

式中，σx與σy分別為材料在X和Y方向的正應力分量，τxy為材料在X和Y方向的剪切應力分量，X與Y分別為材料在X和Y方向的體力。

兩類材料的宏觀力學與自身的材料組成密切相關，表1 為各物理力學參量相似常數，其中，Cl1為混凝土與遺址土幾何長度l的相似常數，Cσ1為混凝土與遺址土應力σ的相似常數，Cε1為混凝土與遺址土應變ε的相似常數，Cγ1為混凝土與遺址土容重γ的相似常數?；趶V義相似理論將表1 中相關式代入式(1)，可得

表1 各物理力學參量相似常數

并且為了讓遺址土能反映出混凝土的應力狀態，要求式(1)與式(2)平衡方程必須一致，并且強度極限與應力的量綱一致，因此混凝土與遺址土強度相似模數為

同理，可得水泥砂漿與膠結材料的平衡方程為

式中，Cσ2為水泥砂漿與膠結物應力相似比，Cγ2為水泥砂漿與膠結物容重相似比，Cl2為水泥砂漿與膠結物幾何長度相似比，σ水泥砂漿為水泥砂漿材料應力，γ水泥砂漿為水泥砂漿材料容重，l水泥砂漿為水泥砂漿材料體長，σ膠結物為膠結物材料應力，γ膠結物為膠結物材料容重，l膠結物為膠結物材料體長，K2為水泥砂漿與膠結物相似模數。

基于相似理論第一定理可得K1=K2，則式(3)與式(4) 可聯立為

對式(5) 進行簡化變形可得遺址土顆粒間膠結物的應力關系式為

基于Hooke 定律可得材料中的應力與應變之間的關系為

式中，F為作用力，A為受力面積，ΔL為圖3 所示各相材料沿F方向的單位變形量值，L為各相材料沿F方向的總變形量值。

由于粗骨料比水泥砂漿的彈性模量大了近一個量級，導致混凝土的宏觀變形主要體現在粗骨料之間水泥砂漿接觸狀態的變化，可將混凝土材料沿F方向的單位變形量值ΔL混凝土視為由粗骨料及水泥砂漿材料兩部分的沿F方向的單位變形量值ΔL構成，即

聯立式(9)與式(10)并同時除以沿F方向的總變形量值L可得水泥砂漿與混凝土材料的應變比為

針對本工作的試驗對象具有不可重復性，不能在既有本體上進行大量的取樣，考慮借助于SEM 及CT 圖像等細觀結構圖像結果進行研究(圖5)。所以對式(11) 中沿F方向的總變形量值L取電鏡圖片單位寬度得到

圖5 混凝土圖像處理

式中，idem為混凝土中水泥砂漿占混凝土CT 圖片的單位寬度面積的倒數。當分段數為1，2，···，10 時，水泥砂漿與混凝土電鏡圖片單位寬度面積比分別為0.35，0.37，0.48，0.51，0.49，0.42，0.57，0.41，0.47，0.53?？芍猧dem= 0.46。

將宏觀試驗結果[19]中混凝土與水泥砂漿的應變比(圖6)與利用式(12)推導所得到的應變比進行對比，發現理論與試驗結果存在偏差。為修正混凝土與水泥砂漿應變關系存在的偏差，將宏觀試驗結果與理論推導結果的比值α定義為

圖6 混凝土與水泥砂漿宏觀試驗應變關系

式中，εs，p為混凝土宏觀試驗應變值，εs，m為水泥砂漿宏觀試驗應變值，εl，p為混凝土理論應變值，εl，m為水泥砂漿理論應變值。

在土力學中，一般認為土顆粒是剛性的，土的宏觀變形和破壞主要體現在土顆粒之間膠結物的接觸狀態的變化，則同理式(12)，膠結物與土樣的應變關系為

式中，idem為遺址土中膠結物占遺址土電鏡圖片的單位寬度面積的倒數，由遺址土圖像(圖7)可知，當分段數為1， 2，···， 10 時，膠結構與遺址土電鏡圖片單位寬度面積比分別為0.64， 0.67， 0.57，0.74，0.78，0.65，0.61，0.51，0.63，0.53，可得idem=0.63。

圖7 遺址土圖像處理

依據式(13) 及式(14) 可得遺址土顆粒間膠結物的應變關系式為

將遺址土基本物理指標：遺址土樣重度γ土=18.57 kN/m3，土顆粒重度γ顆粒=24.31 kN/m3，混凝土重度γ混凝土= 24 kN/m3，水泥砂漿重度γ水泥砂漿=20 kN/m3分別代入式(16)、式(6)及式(15)可求得遺址土顆粒間膠結物應力-應變曲線(圖8)。

圖8 遺址土與膠結物應力-應變曲線

3 細觀尺度仿真模型的構建

土遺址大多處于地上屬于無圍壓狀態，并且土遺址的劣化大多集中于土體外表面，加之準脆性材料(concrete damage plasticity， CDP) 本構的計算結果不存在宏細觀尺度上的差異[20]，所以本文基于SEM圖像，并結合遺址土與膠結物應力-應變曲線進行二維單軸壓縮試驗數值模擬，以驗證本文提出的遺址土顆粒間膠結物理論模型的有效性。

3.1 建模依據

圖9 為宏細觀尺度土顆粒結構分布，由圖9 可知500 倍電鏡圖片土體顆粒清晰，不規則裂縫較為明顯，在該尺度下可以觀測到直徑約5～25 μm 的土顆粒、孔隙的分布及無序排列的黏性顆粒；當提高電鏡倍數到2000～5000 倍時，顆粒形狀清晰，但顆粒數目較少，不具有普遍性；倍數提高到10 000～20 000倍時細觀結構中顆粒形狀模糊。綜上，本工作依據500 倍狀態下的電鏡圖建立遺址土的細觀數值模型。

圖9 宏細觀尺度土顆粒結構分布

為避免拍攝電鏡圖片的隨機性影響研究結果，在拍攝SEM 圖片時，在滿足篩分試驗級別的基礎上，對1 mm×1 mm 的面積域按從左至右及從下到上的原則進行取景，依據500X 下鏡頭拍攝的照片像素尺寸面積可得需拍攝電鏡圖片張數為4 張。

對圖像進行二值化處理的關鍵是閾值的選擇與確定，由于各種圖像上目標物性質的不同以及灰度變化的多樣性，使得傳統的二值化方法難以取得理想的處理效果。為了使處理后的圖像與真實土樣結構吻合，應用Image-Pro Plus 軟件對500X 放大倍數下的SEM 照片進行處理，為縮小單一閾值下計算結果隨機性帶來的誤差，通過調整閾值大小所代表的切平面高低，以壓汞試驗孔隙率結果為依據確定黑白色度的面積比，找到最能真實反應土體細觀結構的最佳閾值[21]，處理方式如圖10，圖中填充部分為土體顆粒，未填充部分為膠結物及孔隙，最終處理結果如圖11 所示。

圖10 不同閾值下的切片處理

圖11 圖像處理結果

將圖像軟件處理后的圖片導入到AutoCAD 軟件并依據以下規則對圖11 進行各相材料邊界劃分：

(1)將灰色部分視為土顆粒，紅色部分視為膠結物作為遺址土數值模擬建模寬泛規則，并忽略距離土顆粒形心較近的紅色散點。

(2)在滿足計算精度及效率的前提上，粒間孔隙中點位置為各個土體顆粒的形心連線的垂線交點[22]，形狀為圓形[23]，孔隙大小范圍依據文獻[24] 進行設置。若存在交點處于顆粒邊緣或顆粒內部，則不進行孔隙位置定義，孔隙位置具體處理方法如圖12所示。

圖12 孔隙位置處理方法

在建模時不考慮土體顆粒表面微孔隙的吸水作用，主要有以下兩個原因：由于土體顆粒表面孔隙存在孔隙內氣壓高，顆粒間氣壓低的現象，并根據DLVO (Darjaguin Landau Verwey Overbeek) 理論，顆粒與水分子之間存在電荷層的相互作用，顆粒周圍的水分子受細觀力學作用僅在表面形成吸附水[25]。

3.2 模型建立

本文根據500X 電鏡圖尺寸大小，通過ABAQUS 軟件構建長為0.644 mm 的剛性板和0.644 mm× 0.431 mm 的二維幾何模型，以剛性板模擬試驗中的壓力設備，將幾何體設置為實體，設置平面厚度為0.020 mm 并選取二維實體單元賦予截面屬性，剛性板的彈性模量為206 GPa，泊松比為0.3。

數值模型由土顆粒、孔隙和膠結物三相組成，將土體顆粒定義為彈性剛體[26]，顆粒間膠結材料采用混凝土塑性損傷(concrete damaged plasticity，CDP)材料本構進行屬性定義。該數值仿真模型中各相材料參數依據具體工況數值進行標定，土顆粒：彈性模量為2.0 GPa，泊松比為0.3；膠結物：彈性模量為21 MPa，泊松比為0.3。統一硬化模型參數：膨脹角Ψ = 30°，流動勢偏移量?= 0.1，雙軸極限抗壓強度與單軸受壓極限強度之比fb0/fc0= 1.16，拉伸子午面上與壓縮子午面上的第二應力不變量之比K=0.666 7，黏性系數μ黏=0.000 5。

本文采用位移加載法，對上壓板頂部分別施加0.006 mm，0.009 mm 和0.012 mm 的豎直位移并采用平滑分析步進行加載。為確保土體與上壓板能夠真實地模擬接觸面之間的滑移以及破壞等現象，在定義相互作用時將剛性板下邊緣定義為接觸主面，土樣模型上邊緣設定為接觸從面，在法向行為中選為硬接觸，在切向行為中摩擦公式選為Penalty (罰函數)，摩擦系數μ摩取0.492[27]。有限元模型建立之后，將模型下邊界的X方向、Y方向及轉角設置為固定端約束，上壓板的X方向及轉角設置為固定端約束。模型網格單元類型采用四節點雙線性平面應變四邊形單元，單元庫選擇動力學、單元族選擇平面應變，網格全局尺寸為0.005。為防止給定載荷步過大，導致迭代發散，收斂困難，采用顯示動力學進行計算。圖13 為建模處理流程，圖14(b)(d)(f)(h) 分別為含水率10%，12%，14%，16%土樣的單軸壓縮仿真模型結果。

圖13 建模流程

圖14 單軸壓縮試驗仿真模型

4 單軸壓縮試驗仿真驗證及結果分析

由于篇幅有限以及工作的相似性，本文僅對模型1 進行不同階段下軸向位移云圖的分析。圖15 為模型1 在不同豎直壓縮位移下的軸向位移云圖，圖15 中圖例“+” 為正方向，“-” 為反方向，圖15 及圖16 中參數量綱為mm。由圖15 可得，土體位移變化比較有規律，整體位移呈現出分層依次變化，其中底部位移最小，而上部位移較大，因模型底部邊界是固支約束，所以底部軸向載荷位移必為0。試塊在單軸抗壓過程中，其內部位移分布不均衡，且這種位移的不均勻性，隨載荷的增加而加劇，在壓縮位移持續作用下，局部位移會發生突變。圖16 為模型1～模型4 在豎直壓縮位移12 μm 下的軸向位移云圖，通過對圖16 中4 個模型軸向位移云圖的分析發現，造成上述現象的原因是土體顆粒分布具有不均勻性，當有較大壓縮位移時，土體顆粒與膠結物之間會發生搓動，土體內部產生較大位移，試樣內部膨脹并造成內部應力過大，宏觀表現為靠近中部的土體最先發生脫落，仿真模型結果與宏觀試驗規律表現一致。

圖15 模型1 在不同豎直壓縮位移下的軸向位移云圖

圖16 模型1 ～模型4 在豎直壓縮位移12 μm 下的軸向位移云圖

圖17 為模型1 在不同豎直壓縮位移下的軸向應力云圖，圖17 中“+” 為拉應力，“-” 為壓應力，圖17 及圖18 中參數量綱為Pa。由圖17 可知，試件壓應力較大部分主要分布在中間部分，試件兩側壓應力較小，整體應力呈現向外側擴展的趨勢，這也是試樣中部發生較大位移的原因。此外，模型兩側的應力變化并不均勻，云圖分布特征呈現出豎條遞進關系，此外試塊中部及兩側應力數值正負號不同，表示中間的應力最大，因為與邊界區域相比，土壤是封閉的，在壓縮過程中土體會發生橫向膨脹的現象。圖18 為模型1～模型4 在豎直壓縮位移12 μm 下的軸向應力云圖，通過歸納圖18 中四個模型的軸向應力云圖，發現在受壓過程中，內部位移的不均勻性是導致試塊內部應力復雜的主要原因，其在宏觀表現為試件內部的破壞程度不同。

圖17 模型1 在不同豎直壓縮位移下的軸向應力云圖

圖18 模型1 ～模型4 在豎直壓縮位移12 μm 下的軸向應力云圖

圖19 為模型1 在不同豎直壓縮位移下的損傷云圖，通過對圖19 分析可得，試件損傷開始萌生并起始于試件底部兩端，此時土體內部粘聚力讓試塊變形可恢復，試塊表面無裂紋產生；隨著壓縮載荷的不斷增大，土體內部基體產生少量微裂紋，產生不可恢復變形；隨著損傷微裂紋進一步擴展、聚集，微裂紋不再是單一地某一條發展，而是沿著土體顆粒邊緣全斷面同時刻一起擴展。當應力達到一定值時，某些區域損傷程度開始加大，并且逐漸貫通，貫穿面外圍出現主裂縫，最后形成幾條貫通于整個試件的裂紋。圖20 為模型1～模型4 在豎直壓縮位移12 μm下的損傷云圖，從圖19 和圖20 中可以看出，孔隙對土體的破壞有一定的影響，即孔隙周圍容易出現應力集中現象導致損傷發生，而且孔洞對損傷路徑的發展有影響，并影響試件最終的破壞面。

圖19 模型1 在不同豎直壓縮位移下的損傷云圖

圖20 模型1 ～模型4 在豎直壓縮位移12 μm 下的損傷云圖

圖21 為宏觀試驗與仿真結果，圖中圖例括號內數字為土樣含水率。本文通過對比理論模型與單軸壓縮數值模型結果(圖21) 發現數值模擬與試驗應力-應變曲線結果吻合度高，初步驗證了所建立模型的可行性及參數選取的合理性。所求得的結果已滿足實際工程的需求，為下一步分析環境對土遺址劣化作用和土遺址安全評價奠定了基礎。

圖21 宏觀試驗與仿真結果

5 結論

本文基于材料類比思想、廣義相似理論推導了顆粒間膠結物的計算公式，探索了細觀尺度下遺址土顆粒間膠結物的力學性能的確定方法；基于遺址土500X 電鏡圖和圖像重建法，構建了顆粒、膠結物和孔隙三部分組成的細觀模型，采用位移加載方法獲得了不同含水率遺址土的應力-應變關系，并與宏觀試驗結果進行對比，從細觀角度對膠結顆粒力學特性進行了分析。主要結論如下：

(2)根據遺址土的膠結特性，將土體劃分為土顆粒、膠結物及孔隙組成的三相結構。通過SEM 圖像重建法確定了顆粒的形狀、大小及分布狀況并建立了反映真實土體的細觀有限元模型。

(3) 通過對細觀有限元模型進行單軸壓縮試驗，發現顆粒與膠結物接觸位置與孔隙周圍首先出現損傷，其主要原因是該區域存在應力集中現象；隨著載荷增大，損傷區域向外擴展延伸，最終導致試件破壞。通過對比發現數值分析結果與宏觀試驗應力-應變曲線結果具有很好的一致性，表明本工作關于遺址土的研究思路和研究方法是可行的，也為取樣困難的室內試驗提供了新方法。