土壤干縮開裂三維模型構(gòu)建及其參數(shù)敏感性分析

農(nóng) 睿 ，萬愉快 ，孫伯顏 ，孫振源 ，朱 磊

（1.寧夏大學(xué)土木與水利工程學(xué)院, 銀川 750021；2.寧夏大學(xué) 寧夏回族自治區(qū)黃河水聯(lián)網(wǎng)數(shù)字治水重點(diǎn)實(shí)驗(yàn)室，銀川 750021）

0 引 言

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

為了模擬土壤干縮開裂的發(fā)展過程，學(xué)者們已經(jīng)建立了大量的土壤干裂數(shù)值模型[12-15]，但目前還缺乏關(guān)于裂隙深度的全面建模方法。土壤干縮開裂的數(shù)值建模方法主要分為有限元法（finite element method, FEM）、離散元方法（discrete element method, DEM）和線彈性斷裂力學(xué)理論。VO等[16]采用水力-力學(xué)耦合和內(nèi)聚裂隙單元的有限元程序，研究了黏土干縮裂隙的幾何形狀和發(fā)育過程，數(shù)值分析結(jié)果在水力擴(kuò)散、裂隙萌生、收縮演化及干燥階段時(shí)序方面與試驗(yàn)數(shù)據(jù)較為吻合。沈珠江等[17]基于FEM搭建了非飽和土簡化固結(jié)理論數(shù)值模型，數(shù)值分析結(jié)果表明簡化固結(jié)理論方法具有實(shí)用性。然而，由于需要預(yù)先定義少量裂紋或?qū)⒛Ｐ瓦吔缱鳛闈撛诹鸭y，F(xiàn)EM模擬對(duì)于裂隙產(chǎn)生和發(fā)育的預(yù)測(cè)能力較弱[18]。

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

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

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

1 模型構(gòu)建

1.1 三維模型重力項(xiàng)引入及三維彈簧網(wǎng)格劃分

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

如圖1所示，節(jié)點(diǎn)間連接著長度為d和彈簧系數(shù)為K的彈簧，任意2 個(gè)相鄰節(jié)點(diǎn)的應(yīng)變?chǔ)舏j，彈簧和重力所產(chǎn)生的合力Fi以及每個(gè)節(jié)點(diǎn)處的總能量Ei分別為[27]：

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

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

式中n表示由節(jié)點(diǎn)i連接的三棱柱結(jié)構(gòu)，Vn是三棱柱n的體積。

1.2 土壤結(jié)構(gòu)收縮特性

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

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

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

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

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

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

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

2 案例分析

2.1 試驗(yàn)區(qū)概況

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

表2 試驗(yàn)區(qū)土壤物理性質(zhì)Table 2 Soil physical properties in the test area

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

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

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

2.2.1 土壤表面裂隙圖像處理

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

圖4 試驗(yàn)區(qū)S1圖像二值化Fig.4 Image binarization of experimental area S1

2.2.2 土壤裂隙網(wǎng)絡(luò)量化方法

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

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

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

2.3 三維土壤干裂模型應(yīng)用

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

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

②出入線明挖段為三線矩形斷面，結(jié)構(gòu)斷面寬度15～19m，托換梁需跨過明挖隧道，托換梁最大跨度達(dá)23.3m，工程施工難度大，實(shí)施風(fēng)險(xiǎn)高。

3 結(jié)果與分析

3.1 試驗(yàn)實(shí)測(cè)和計(jì)算結(jié)果

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

圖5 試驗(yàn)區(qū)S1～S4裂隙Minkowski密度函數(shù)Fig.5 Minkowski density function of crack for experimental area S1-S4

通過鋼針法測(cè)量了試驗(yàn)的裂隙深度，結(jié)果如圖6所示。

圖6 試驗(yàn)區(qū)S1～S4裂隙深度相對(duì)頻率統(tǒng)計(jì)Fig.6 Statistics of crack depth relative frequency for Experimental area S1-S4

由圖6可知，試驗(yàn)裂隙深度主要分布在土深1～6 cm。4次獨(dú)立試驗(yàn)S1～S4的裂隙深度相對(duì)頻率隨土壤深度均呈現(xiàn)先增大后減小的趨勢(shì)，但峰值及其所處深度不同，當(dāng)?shù)竭_(dá)各峰值后相對(duì)頻率值均迅速下降。

3.2 模型率定結(jié)果

基于如表3所示初始值，根據(jù)Monte-Carlo仿真結(jié)果，計(jì)算各輸入?yún)?shù)組合的評(píng)價(jià)指標(biāo)，確定模型的參數(shù)率定值。

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

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

表4 模型率定的評(píng)價(jià)指標(biāo)Table 4 Evaluation indicators of model calibration

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

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

3.3 模型驗(yàn)證結(jié)果

通過試驗(yàn)區(qū)S2～S4的裂隙驗(yàn)證了模型參數(shù)的可靠性，在率定參數(shù)組合下，每組獨(dú)立試驗(yàn)通過調(diào)整“隨機(jī)種子數(shù)”生成100組土壤臨界應(yīng)變隨機(jī)場(chǎng)，得到100組模擬圖像。比較每組試驗(yàn)裂隙（Minkowski密度及裂隙深度頻率）與對(duì)應(yīng)100組模擬裂隙（Minkowski密度及裂隙深度頻率平均值）間的差異，如表5所示。結(jié)果顯示，S2～S4試驗(yàn)與模擬Minkowski密度和裂隙深度相對(duì)頻率的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，總體模擬結(jié)果較好，可見，所構(gòu)建模型的模擬結(jié)果較好，土壤表面收縮的假設(shè)與邊界條件的設(shè)置合理，模型可靠性較高，能夠應(yīng)用于實(shí)際模擬。由圖10可知，模擬裂隙深度為1～6 cm的占比為0.90，模擬裂隙深度為7～10 cm的占比為0.10，大部分裂隙的深度集中在1～6 cm之間。不同深度下裂隙試驗(yàn)與模擬結(jié)果的深度頻率平均值基本一致，擬合線變化趨勢(shì)和范圍基本相同，表明構(gòu)建的三維模型對(duì)于模擬土壤開裂結(jié)束后的裂隙深度頻率分布特征具有較高的一致性。

表5 模型驗(yàn)證結(jié)果評(píng)價(jià)指標(biāo)Table 5 Evaluation indexes of model validation results

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

3.4 敏感性分析

因?yàn)榭v向彈性系數(shù)Kz的變化代表不同性質(zhì)（縱向彈性）土壤下裂隙形態(tài)以及裂隙深度分布的變化，所以本文在率定參數(shù)不變的前提下，選取土壤開裂結(jié)束時(shí)間（26 h），選擇縱向彈性系數(shù)作為敏感性分析的參數(shù)，研究其變化對(duì)裂隙形態(tài)、裂隙深度頻率以及Minkowski密度的影響。

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

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

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

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

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

歐拉數(shù)反映的是裂隙網(wǎng)絡(luò)連通性的指標(biāo)，歐拉數(shù)密度函數(shù)的峰值隨著縱向彈性系數(shù)的增大而減小（圖13c），這表明孤立裂隙的數(shù)量是隨著縱向彈性系數(shù)的增大在減小。正如圖11所呈現(xiàn)的，縱向彈性系數(shù)越大，相同土層裂隙骨架的長度和寬度越小，孤立裂隙數(shù)越少。從而可以初步判斷在相同的土壤深度下，由于縱向彈性系數(shù)的減小使得土壤結(jié)構(gòu)彈性減小，增強(qiáng)了土體沿深度開裂的趨勢(shì)。

4 結(jié) 論

本文基于Vogel模型引入了重力項(xiàng)，劃分了三維彈簧網(wǎng)格，建立三維土壤干縮開裂模型，探究了縱向彈性系數(shù)對(duì)裂隙深度的影響規(guī)律。主要結(jié)論如下：

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

2）縱向彈性系數(shù)越小，Minkowski密度函數(shù)峰值越大，同一深度下對(duì)應(yīng)的土壤開裂程度越大，孤立裂隙數(shù)越多?？v向彈性系數(shù)的減小減弱了土壤彈性，使得下層兩節(jié)點(diǎn)位移增大，深層土壤網(wǎng)格更容易因達(dá)到彈簧臨界應(yīng)變值而開裂，總體上促進(jìn)了裂隙沿深度的發(fā)育。

本研究基于田間土壤試驗(yàn)區(qū)情況（土壤表面平整）建立了長方體彈簧網(wǎng)格結(jié)構(gòu)，但未考慮復(fù)雜結(jié)構(gòu)的情況，針對(duì)具有傾角的斜坡土壤可在后期構(gòu)建不同類型網(wǎng)格結(jié)構(gòu)進(jìn)行研究。模型中以彈簧收縮描述水分蒸發(fā)引起的土壤收縮，但未考慮土壤含水率的變化。而現(xiàn)實(shí)中含水率變化的程度和速率會(huì)顯著影響裂隙發(fā)育速率。在后續(xù)的研究中，可以嘗試在三維模型的基礎(chǔ)上耦合土壤的真實(shí)含水率，建立含水率與裂隙演化的關(guān)系，以提高模型的預(yù)測(cè)準(zhǔn)確性及可靠性。此外，該模型未考慮復(fù)雜發(fā)育（彎曲、傾斜等）裂隙的情況，在未來需推廣三維模型面對(duì)復(fù)雜裂隙情景下的應(yīng)用。