土壤干縮開裂三維模型構建及其參數敏感性分析

農 睿 ，萬愉快 ，孫伯顏 ，孫振源 ，朱 磊

（1.寧夏大學土木與水利工程學院, 銀川 750021；2.寧夏大學 寧夏回族自治區黃河水聯網數字治水重點實驗室，銀川 750021）

0 引 言

在干旱條件下，土壤干縮開裂是在自然和人造物體破壞過程中可以觀察到的最常見的物理現象之一[1]。水分蒸發產生的裂隙對土壤的力學性能和水力性能有顯著影響，受到了農業、地質、環境及工程等學科研究者的廣泛關注。裂隙網絡為水和溶質運輸提供了優先的流動路徑，會顯著增加土體的導水性[2]，從而增加地下水污染的可能性[3]。同時，裂隙破壞了土壤的完整性和穩定性，會導致土體強度受損、過度變形和壓縮性增加等一系列問題[4]，對土方工程[5]、堤壩[6]、路堤[7]、工程屏障[8]的穩定構成重大威脅。此外，干燥裂隙的存在使得土壤表面粗糙度增加，加速了土壤風化和侵蝕，推動了局部或區域尺度沙塵暴和荒漠化的形成[9-11]。因此，對土壤干縮開裂進行系統研究，對于防治地質和巖土工程問題具有重要意義。

為了模擬土壤干縮開裂的發展過程，學者們已經建立了大量的土壤干裂數值模型[12-15]，但目前還缺乏關于裂隙深度的全面建模方法。土壤干縮開裂的數值建模方法主要分為有限元法（finite element method, FEM）、離散元方法（discrete element method, DEM）和線彈性斷裂力學理論。VO等[16]采用水力-力學耦合和內聚裂隙單元的有限元程序，研究了黏土干縮裂隙的幾何形狀和發育過程，數值分析結果在水力擴散、裂隙萌生、收縮演化及干燥階段時序方面與試驗數據較為吻合。沈珠江等[17]基于FEM搭建了非飽和土簡化固結理論數值模型，數值分析結果表明簡化固結理論方法具有實用性。然而，由于需要預先定義少量裂紋或將模型邊界作為潛在裂紋，FEM模擬對于裂隙產生和發育的預測能力較弱[18]。

DEM將土壤視為離散元的組合，在模擬干燥開裂方面具有良好的潛力[19]。PERON等[20]使用DEM對土壤的干縮開裂進行仿真，重現了細粒土壤裂隙的萌生、擴展和幾何形狀。AMARASIRI等[21]采用通用離散元代碼（universal distinct element code, UDEC）對黏土的超固結現象進行建模，該模型可以捕捉到不同長度、高度和材料模具中包含的土壤裂隙演變特征。SIMA等[22]基于三維離散元方法搭建數值模型，分析了土基界面特性、試樣厚度、收縮參數和微觀力學參數等材料特性對黏土干縮開裂的影響。需要指出的是，傳統離散元模型將離散單元視為剛體，無法預測干縮過程中的水分蒸發、熱量和質量交換等多物理過程[19]。

在已有的研究中，KONRAD等[23]利用線彈性斷裂力學理論構建了平均裂縫深度和間距預測框架，簡要描述了裂紋擴展現象。VOGEL等[24]開發了一個基于Hookean彈簧網格的物理模型，以重現土壤裂縫模式的真實整體外觀，模擬水分蒸發引起的裂縫動態拓展。與Konrad模型相比，它能捕捉到自然界中可觀測到的裂縫網絡發展非線性動力學的顯著特征（即裂隙的特征形狀和分叉的特征角度）。一般情況下，彈簧網絡模型一個周期的仿真計算時間（即彈簧網絡對節點的計算力和節點位置的更新）與模型的復雜性（即節點和彈簧的數量）成正比增加，而有限元模型的計算時間（即計算模型的邊界條件和求解剛度矩陣定義的聯立方程所用時間）則與節點數量的立方成正比[25]。因此，彈簧網格模型在模擬土壤開裂過程上可以節省大量計算時間。然而，Vogel模型只能再現土壤表面的裂隙形態特征，不能模擬裂隙彈簧網格沿深度的破壞，故需要進行三維彈簧網格建模。

研究裂隙沿深度產生與拓展的規律及機理對于研究土壤干縮開裂過程非常重要。但迄今為止，學界對土壤三維干縮開裂過程的研究較少，也鮮有考慮裂隙深度分布的土壤干縮裂隙數值模型，土壤表面下干燥裂隙的萌生和擴展過程尚不清楚。缺乏對土體內部開裂過程或三維開裂特征的認識，將阻礙對土體水力學行為作出準確評價。因此，為了模擬裂隙沿深度方向的形成與拓展，本文基于Vogel模型引入重力項，劃分三棱柱狀彈簧網格結構，并考慮土壤豎向彈性對裂隙深度的影響，建立三維土壤干縮開裂模型，通過Monte-Carlo模擬驗證模型參數值與裂隙模式之間的關系，并通過最小化裂隙面積、長度、歐拉數密度和深度頻率與實測值之間的差值來設定優化參數，最后分析模擬結果對縱向彈性系數的敏感性，以期為研究土壤彈性變化對裂隙深度的影響提供理論參考。

1 模型構建

1.1 三維模型重力項引入及三維彈簧網格劃分

本文基于線彈性斷裂理論[26]構建了三維土壤干裂模型，模擬因水分蒸發而導致土壤干燥開裂的動態。與Vogel模型不同的是，三維模型考慮了重力的作用，每個節點的重力被確定為重力加速度（g，固定為9.8 m/s2）和節點上的質量（mi）的乘積。模擬區域的土壤被劃分為三棱柱（原為三角形）組成的彈簧網格結構，如圖1所示。

如圖1所示，節點間連接著長度為d和彈簧系數為K的彈簧，任意2 個相鄰節點的應變εij，彈簧和重力所產生的合力Fi以及每個節點處的總能量Ei分別為[27]：

式中│xi-xj│代表網格中相鄰節點間的距離，x=(x,y,z)T。d是彈簧的松弛長度，將其從初始值為1（完全松弛狀態）逐漸減小來模擬土壤干縮開裂過程。K為彈簧的彈性系數，與土壤彈性有關。VOGEL等[24]將水平彈簧的K設置為1。當彈簧受力F大于摩擦力f時，節點才會移動，并且節點根據式（3）向能量最小位置移動，由式（2）可以求解節點新位置。此時判斷節點應變是否大于彈簧臨界應變，若滿足上述條件，則彈簧斷裂形成裂隙。通過彈簧的臨界應變引入土壤的異質性，臨界應變閾值是均值（））和方差（δ2）的高斯概率分布，在空間隨機分布。邊界條件指定為邊界節點不允許垂直邊界移動，這是線彈性斷裂理論模型的前提條件[25]。

在本文構建的三維模型中，一個公共節點連接著n個相鄰三棱柱，每個三棱柱結構的質量均勻分布在其6個節點上。且每個三棱柱結構的質量由土壤的密度（ρ）和結構的體積（V）求得，因此，節點i的質量可表示為

式中n表示由節點i連接的三棱柱結構，Vn是三棱柱n的體積。

1.2 土壤結構收縮特性

VOGEL等[24]模型中沒有涉及時間的概念，但將彈簧張力u=1-d解釋為時間或土壤含水率的函數。同時TANOUE等[25]觀察到，黏土塊的表面比內部干縮得更快。由于收縮速率的不同，膨脹應力導致土壤表面干縮開裂。為了描述這種現象，假設表層土壤（與空氣直接接觸）的彈簧結構處于不斷收縮的狀態。也就是說，表層土壤的彈簧松弛長度d1隨著時間的推移而減小，減小步長即土壤脫水時間的離散化設置為10-4，減小至最小松弛長度d1min（代表土壤開裂程度）。而內部土壤的收縮率保持在初始值，即d2=1。

1.3 模型求解及裂隙深度提取

本文在CodeLite平臺采用 C++語言編寫程序代碼進行計算求解，模型的模擬結果由平均臨界應變、方差δ2、松弛參數nit、摩擦力f、縱向彈性系數Kz共同控制，其物理性質如表1所示。

表1 三維模型參數Table 1 List of parameters of the three-dimensional Model

圖2展示了根據裂隙點坐標生成的三維可視化圖像。圖3展示了沿土壤深度分層的裂隙二值.tiff圖像。該算例設置400×400×40個節點（代表100 cm × 100 cm× 10 cm的土壤體積）的長方體彈簧網格結構，2節點間隔為0.25 cm，均值和方差分別為0.1和0.20，摩擦力f為8×10-4。選取開裂時間為18 h的裂隙圖像進行展示。如圖3所示，裂隙面積、長度形態土表至深處呈現逐漸減小的變化特征。STEWART等[28]假設裂隙垂直于土壤表面發育，截面形態呈三角形或楔形。作為參考，該模型假設裂隙垂直于土壤表面發育，與試驗中測量的裂隙豎向深度進行對比。讀取每層的裂隙分布信息，篩選土壤表層裂隙點（以1表示），再逐層豎直對比裂隙點的變化，若在下層土壤中該裂隙點轉變為非裂隙點（以0表示）則記錄該點深度，循環結束后統計模擬裂隙的深度頻率分布特征。

圖2 三維裂隙圖像可視化示例Fig.2 An example of visualization of 3D crack images

圖3 模擬裂隙圖像隨土壤剖面不同深度（h）的變化趨勢Fig.3 Variation trend of simulated crack image along different depth (h) of soil proile

2 案例分析

2.1 試驗區概況

試驗場地位于寧夏回族自治區吳忠市同心縣（105°54'E, 36°58'N），地處鄂爾多斯臺地南部黃土高原，屬北溫帶大陸性季風氣候區。干旱少雨，年平均降水272 mm，蒸發量大，年蒸發量2 325 mm。年平均氣溫8.6 ℃，多年平均日照3 024 h。土壤為粉質黏土，田間持水率質量分數為23.82%，0～50 cm土層平均容重為1.58 g/cm3，試驗區土壤的物理性質如表2所示。

表2 試驗區土壤物理性質Table 2 Soil physical properties in the test area

該試驗于2022年8月在同心縣試驗區進行。前期調研發現，田間土壤裂隙平均最大深度小于10 cm。結合前人對農田裂隙特征參數的研究經驗[29]，選取深度為10 cm，表面尺寸100 cm×100 cm的4個土壤長方體試驗區作為研究對象，編號為S1～S4。試驗前，對土體進行松動消除田間土壤的原有裂紋；隨后于土壤表面施加積水，保持恒定的注水速率（2 L/s）。利用時域反射技術（time domain reflectometry, TDR）測量土壤含水率，以確保土壤在0 ~ 10 cm完全飽和。飽和時，土體表面膨脹而無裂隙，至此開始記錄土壤收縮時間，使其自然脫水，期間日均溫為20.0～32.0 ℃，空氣濕度為43%～60%，試驗過程中天氣狀況良好，無降雨。

運用Canon EOS 5DII數碼相機（最高分辨率5 616×3 744像素）每隔2 h拍照記錄裂隙形成和發育狀況，試驗時間為08:00-18:00，連續觀測3 d。夜間采用防雨塑料鋪設試驗區，防止雨水降落和水分夜間蒸發。土壤開裂完全后，采用長為15 cm，直徑為1.5 mm，測量精度為1mm的鋼針測量裂隙豎向深度，將鋼針分別依次插至裂隙底部。沿裂隙長度方向布置測量點，2個測量點的距離間隔均為1 cm，每個測量點重復測量5次并取平均值，每個試驗區分別統計100個點。

2.2 基于圖片的裂隙量化方法

2.2.1 土壤表面裂隙圖像處理

為了深入研究表面裂隙模式的幾何形態，采用以下4個步驟處理圖像：1）將圖像中的非試驗區域裁切去除，并統一圖像尺寸為1 024×1 024像素；2）將裂隙彩色圖像轉化為灰度圖像；3）采用閾值分割法，將裂隙灰度圖像轉化為二值圖像，若像素點的值大于閾值則劃分為1，小于閾值劃分為0。4）通過去噪去除二值圖像中的雜質。以S1為例，處理前后裂隙圖像對比如圖4所示。

圖4 試驗區S1圖像二值化Fig.4 Image binarization of experimental area S1

2.2.2 土壤裂隙網絡量化方法

采用描述裂隙拓撲性質的Minkowski數（Mk） 量化土壤表層裂隙網絡，其包含3個基本幾何特征，即裂隙的面積（A，cm2）、長度（L，cm）和歐拉數（E）。對于二維表面裂隙，Mk（k=0，1，2）的計算式如下[30]：式中X表示黑色像素集合；δX表示結構X的邊界；r是邊界的曲率半徑；ds表示邊界單元；M0、M1、M2分別表示二維結構中裂隙的面積、裂隙的邊界長度和結構聯通性的拓撲度量。A(X)、L(X)、E(X)分別表示黑色像素集合X的面積、長度和歐拉數。歐拉數E定義為孤立對象的數量N（封閉凸）減去孔洞數H（封閉凹），即E=N-H，描述裂隙網絡的連通性。

閔可夫斯基密度（Minkowski densities）可以用來衡量不同尺寸圖像的裂隙結構特征，其計算式[30]如下：

式中Δ為圖像面積；m0、m1、m2分別表示Minkowski面積密度、長度密度和歐拉數密度。

2.3 三維土壤干裂模型應用

選擇S1試驗區率定模型參數，包括均值、方差δ2、松弛參數nit、摩擦力f和縱向彈性系數Kx，S2～S4試驗區進行模型驗證。模型率定的目的是通過調整輸入參數來捕捉現場記錄的裂縫模式。本研究采用Monte-Carlo法率定輸入參數，使用24個核（Intel Xeon Gold 6 226 CPU，運行頻率為2.9 GHz，內存為256 GB）的工作站進行運算，按照1.3節進行模型求解和裂隙深度提取，與基于圖片的裂隙量化方法計算的結果進行對比。

ZHU等[31]通過最小化試驗和模擬的Minkowski密度差異來校準模型時間。本研究采用ZHU的方法，通過試驗的4個時間點（4、14、24、34 h）對模型時間進行匹配，每個時間點對應模型中一定的彈簧松弛長度。通過式（8）對試驗和模擬裂隙圖像的Minkowski密度進行計算。在率定和驗證過程中，采用決定系數（R2）、均方根誤差（SRMSE）、一致性指標（SIA）和偏差（SBIAS）4個指標[32-33]對模擬結果Minkowski密度和裂隙深度頻率的準確性和一致性進行檢驗，其中R2和SIA值越接近1、SBIAS和SRMSE越小表示模擬精度越高。一般認為，R2大于0.5，SRMSE在20%以內時模型達到率定要求[34]。

②出入線明挖段為三線矩形斷面，結構斷面寬度15～19m，托換梁需跨過明挖隧道，托換梁最大跨度達23.3m，工程施工難度大，實施風險高。

3 結果與分析

3.1 試驗實測和計算結果

圖5展示了S1～S4試驗區裂隙Minkowski密度在4個時間點（4、14、24、34 h）的變化趨勢。在土壤干縮裂隙隨時間變化的過程中，由于裂隙不斷產生與發育，面積密度、長度密度與歐拉數密度均呈現增函數的變化，這符合自然條件下干縮裂隙的發育規律[3]。隨著裂隙繼續拓展，土壤含水率減少，裂隙模式逐漸固定，長度密度的增長趨勢隨時間變緩，而面積密度仍在持續增大。

圖5 試驗區S1～S4裂隙Minkowski密度函數Fig.5 Minkowski density function of crack for experimental area S1-S4

通過鋼針法測量了試驗的裂隙深度，結果如圖6所示。

圖6 試驗區S1～S4裂隙深度相對頻率統計Fig.6 Statistics of crack depth relative frequency for Experimental area S1-S4

由圖6可知，試驗裂隙深度主要分布在土深1～6 cm。4次獨立試驗S1～S4的裂隙深度相對頻率隨土壤深度均呈現先增大后減小的趨勢，但峰值及其所處深度不同，當到達各峰值后相對頻率值均迅速下降。

3.2 模型率定結果

基于如表3所示初始值，根據Monte-Carlo仿真結果，計算各輸入參數組合的評價指標，確定模型的參數率定值。

表3 三維土壤干縮開裂模型參數率定值Table 3 Parameter calibration values of 3D soil shrinkage cracking model

對比不同時間土壤表面試驗與模擬裂隙圖像圖7，以及模擬裂隙圖像與試驗裂隙圖像的Minkowski密度匹配情況圖8，可以看出模擬與試驗裂隙的Minkowski密度發展規律基本一致，4個時間點（4、14、24、34 h）的面積、長度、歐拉數密度值基本相同。對比土壤開裂結束后模擬與試驗裂隙深度頻率分布圖9可知，裂隙深度（4～5 cm）的頻率占比最大，深度（8～10 cm）的頻率占比最小，試驗與模擬的裂隙深度相對頻率基本一致。量化率定結果如表4所示，模擬結果的R2均不小于0.944，SIA均大于等于0.986，SBIAS均小于等于0.103，SRMSE均小于等于0.052，可認為模型滿足率定要求。

表4 模型率定的評價指標Table 4 Evaluation indicators of model calibration

圖8 Minkowski密度試驗值與模擬值比較Fig.8 Comparison between experimental and simulated values of Minkowski density

圖9 裂隙深度試驗值與模擬值比較Fig.9 Comparison between experimetnal and simulated values of crack depth

3.3 模型驗證結果

通過試驗區S2～S4的裂隙驗證了模型參數的可靠性，在率定參數組合下，每組獨立試驗通過調整“隨機種子數”生成100組土壤臨界應變隨機場，得到100組模擬圖像。比較每組試驗裂隙（Minkowski密度及裂隙深度頻率）與對應100組模擬裂隙（Minkowski密度及裂隙深度頻率平均值）間的差異，如表5所示。結果顯示，S2～S4試驗與模擬Minkowski密度和裂隙深度相對頻率的R2在0.849～0.959之間，SIA為0.965～0.988，SBIAS為0.103～0.189，SRMSE為0.005～0.083。R2和SIA均大于0.8，SBIAS和SRMSE均小于0.19，總體模擬結果較好，可見，所構建模型的模擬結果較好，土壤表面收縮的假設與邊界條件的設置合理，模型可靠性較高，能夠應用于實際模擬。由圖10可知，模擬裂隙深度為1～6 cm的占比為0.90，模擬裂隙深度為7～10 cm的占比為0.10，大部分裂隙的深度集中在1～6 cm之間。不同深度下裂隙試驗與模擬結果的深度頻率平均值基本一致，擬合線變化趨勢和范圍基本相同，表明構建的三維模型對于模擬土壤開裂結束后的裂隙深度頻率分布特征具有較高的一致性。

表5 模型驗證結果評價指標Table 5 Evaluation indexes of model validation results

圖10 試驗與模擬裂隙深度頻率分布比較Fig.10 Comparison between simulated and measured values of crack depth frequency

3.4 敏感性分析

因為縱向彈性系數Kz的變化代表不同性質（縱向彈性）土壤下裂隙形態以及裂隙深度分布的變化，所以本文在率定參數不變的前提下，選取土壤開裂結束時間（26 h），選擇縱向彈性系數作為敏感性分析的參數，研究其變化對裂隙形態、裂隙深度頻率以及Minkowski密度的影響。

圖11 不同縱向彈性系數(KZ)下模擬裂隙沿土壤深度(h)的分層圖像Fig.11 The layered image of simulated crack along soil depth (h)under different longitudinal elastic coefficient (Kz)

圖12為裂隙發育結束后不同縱向彈簧系數Kz對裂隙深度頻率分布的影響。結果表明，不同縱向彈性系數下，裂隙深度頻率峰值不同，但變化形態基本一致。顯然，隨著縱向彈性系數的增大，在深度（h=1～4 cm）下，裂隙深度相對頻率增大，這表明深度1～4 cm的裂隙占比越大。反之，在深度（h=5～10 cm）下，裂隙深度相對頻率減小，這表明深度5～10 cm的裂隙占比越小。從整體上看，縱向彈性系數越小，裂隙沿深度方向發育的趨勢越明顯，深裂縫（h=5～10 cm）的占比越大。

圖12 不同縱向彈性系數KZ下裂隙深度頻率分布Fig.12 Frequency distribution of crack depth under different longitudinal elastic coefficient (Kz)

由Minkowski 密度函數定量描述不同縱向彈性系數下裂隙圖像的變化（圖13a和圖13b），其面積、長度和歐拉數密度具有相同的變化規律，在同一深度，面積密度和長度密度隨縱向彈性系數的增大而減小，表示裂隙數目隨縱向彈性系數的增大而減小。其原因可能為土壤表面收縮相同的情況下，結構產生的總能量相同，縱向彈性系數的增大使得土壤結構彈性增強，導致下層節點位移減小，裂隙骨架的長度與寬度減小。

圖13 不同縱向彈性系數（KZ）下的Minkowski函數Fig.13 Minkowski functions under different longitudinal elastic coefficient (KZ)

歐拉數反映的是裂隙網絡連通性的指標，歐拉數密度函數的峰值隨著縱向彈性系數的增大而減小（圖13c），這表明孤立裂隙的數量是隨著縱向彈性系數的增大在減小。正如圖11所呈現的，縱向彈性系數越大，相同土層裂隙骨架的長度和寬度越小，孤立裂隙數越少。從而可以初步判斷在相同的土壤深度下，由于縱向彈性系數的減小使得土壤結構彈性減小，增強了土體沿深度開裂的趨勢。

4 結 論

本文基于Vogel模型引入了重力項，劃分了三維彈簧網格，建立三維土壤干縮開裂模型，探究了縱向彈性系數對裂隙深度的影響規律。主要結論如下：

1）本文所建立的三維模型從土壤裂隙的表面動態形態特征和深度分布2個方面描述三維土壤裂隙。不僅較好地再現農田土壤表面干裂網絡的動態形態特征，并且能夠模擬土壤水分蒸發結束后的裂隙深度。應用Minkowski函數描述的農田土壤表面裂隙面積、長度、歐拉數密度隨時間呈現增函數變化，符合自然裂隙的發育規律。面積、長度、歐拉數密度及裂隙深度頻率的決定系數在0.849～0.959之間，一致性指標大于0.940，偏差在0.103～0.190之間，均方根誤差在0.005～0.083之間，模擬結果較好。

2）縱向彈性系數越小，Minkowski密度函數峰值越大，同一深度下對應的土壤開裂程度越大，孤立裂隙數越多。縱向彈性系數的減小減弱了土壤彈性，使得下層兩節點位移增大，深層土壤網格更容易因達到彈簧臨界應變值而開裂，總體上促進了裂隙沿深度的發育。

本研究基于田間土壤試驗區情況（土壤表面平整）建立了長方體彈簧網格結構，但未考慮復雜結構的情況，針對具有傾角的斜坡土壤可在后期構建不同類型網格結構進行研究。模型中以彈簧收縮描述水分蒸發引起的土壤收縮，但未考慮土壤含水率的變化。而現實中含水率變化的程度和速率會顯著影響裂隙發育速率。在后續的研究中，可以嘗試在三維模型的基礎上耦合土壤的真實含水率，建立含水率與裂隙演化的關系，以提高模型的預測準確性及可靠性。此外，該模型未考慮復雜發育（彎曲、傾斜等）裂隙的情況，在未來需推廣三維模型面對復雜裂隙情景下的應用。