2024, Vol. 56
2. 中国科学院大学, 北京 100049;
3. 中国科学院大学南京学院, 南京 211135
土壤由固体颗粒、有机质和孔隙(水和气)组成。农田土壤压实是由于施加载重、振动或压力导致土壤颗粒重新排列，从而使土壤孔隙度降低和容重增加的过程[1-2]，主要是由于农业机械进行田间作业时，对土壤表面施加的外力作用造成的[3]。在土力学领域[4]，土的压缩(compression)是指土壤中的孔隙体积缩小；土的压实(compaction)是在动态荷载作用下，土壤颗粒克服粒间阻力而重新排列，土壤孔隙减小、密度增加，进而在短时间内得到土体新的结构强度。从土壤物理学(soil physics)和土力学(soil mechanics)的交叉领域出发，土壤压实可以看作压缩与压实的结合，因为土壤压实过程中孔隙度降低总是伴随着土体的体积减小，故本文不再区分压缩和压实。
近年来我国农村劳动力大规模转移到工业和服务业，为了满足提高农业生产力和保障粮食安全的需求，农业机械化得到了快速发展[5-6]。但是随着农机具田间作业次数增多、大中型机械普遍应用，导致了严重的土壤压实问题：机械压实破坏土壤结构，导致土壤孔隙减小，限制水分入渗和根系生长，降低作物产量潜力，从而影响农业生产的可持续发展[7-9]。20世纪50年代以来，国内外科研工作者围绕农田土壤压实问题开展了大量研究，在土壤压实行为和压实参数测量方面取得了许多成果，但这些研究往往局限于土力学或农业土壤学单一方面，而缺乏系统性的交叉研究和总结[6, 9-14]。
本文结合土力学和农业土壤机械压实研究，从土壤压缩行为及其特性、土壤形变机理和压实过程模拟3个方面对近年来土壤压实方面的研究进展进行系统性的梳理与评述。
1 土壤压缩行为土壤是由固、液、气组成的三相系统，农机具作业过程中对其表面施加应力，土壤应力通过颗粒接触传递的方式从表层土壤传递到深层，造成土壤变形。压应力传递过程中需要克服土壤自身强度才能破坏土壤结构，土壤受到外力作用后的压缩行为以土壤压缩特性为指标。室内实验常用不允许土样侧向膨胀的单轴压缩试验绘制土壤压缩曲线，并测定预固结压力、压缩指数和回弹指数等压缩特性参数[15]。
1.1 压缩曲线与压缩特性单轴压缩试验中，将原状土或扰动土制备成一定规格的土样，置于固结仪中施加不同的压力pi，测定相对稳定时的压缩变形量，计算出相应的孔隙比ei，由此可绘制出e-lgp曲线(即土壤压缩曲线)。典型的压缩曲线可分为弹性和塑性两部分：土壤承受的压力位于弹性曲线(Swelling Line, SL)段时，土壤可在外力撤去后自动恢复一部分，即发生弹性形变；而当土壤承受的压力超出弹性曲线，进入塑性曲线(Virgin Compression Line, VCL)段后，即使撤去外力土壤也不能或很少恢复，即发生塑性形变[16]，如图 1所示。
|
图 1 实测压缩曲线(e-lgp)示例(A)以及回弹指数Cs、压缩指数Cc、预固结压力值Pc示意图(B) Fig. 1 Example of a measured compression curve (A); swelling index (Cs), pre-compression stress (Pc) and compression index (Cc) (B) |
SL的平均斜率模量称为回弹指数(swelling index, Cs)，表示土壤被压实后恢复至压缩前状态的能力。Cs越大，土壤被压缩后膨胀回弹的能力越强[17]。VCL接近直线，其斜率模量称为压缩指数(compression index, Cc)，是判断土壤压实敏感性的指标。压缩指数越大，土的压缩性就越高，即压实敏感性越强，抗压实能力越弱。一般认为，Cc<0.2时为低压缩性土，0.2≤Cc≤0.4为中压缩性土，Cc > 0.4时为高压缩性土。压缩系数(compression coefficient，a)是孔隙减小量与竖向有效压应力增量的比值，即e-p曲线上任一点的切线斜率模量。压缩指数和压缩系数虽然都是反映土壤压缩性的指标，但a随着初始压力及后续压力变化范围的大小而改变，Cc在e-lg p塑性曲线段是常量，更方便实际使用[4]。
SL与VCL的交点对应的土壤承受压力，称为预固结压力(pre-compression stress，Pc)，代表土壤的承压能力极限。对于某一具体土壤而言，当p<Pc时，预计变形是弹性的；当p > Pc时，即使撤去外力也会导致永久变形(压实)，因此可以通过控制施加的压力不超过预固结压力来避免土壤压实[18]。但是也有研究表明，即使p<Pc，土壤变形也不是完全弹性的，仍存在残留变形[19]。
因此，压缩曲线可用来表达土壤压缩行为，量化土壤承载能力、压实敏感性和回弹能力，比较不同土壤发生板结的难易程度，还能预估变形后土壤结构的可恢复性。预固结压力、压缩指数和回弹指数这3个压缩特性指标还可作为土壤压实过程模型的输入参数。
1.2 预固结压力的确定方法测定土壤压缩性的方法和计算公式不同，得出的压缩特性指标(特别是预固结压力)可能不完全相同[20-21]。Keller等[17]发现，不同曲线表现出的土壤力学行为不完全相同：半对数图会将e-p图中的直线强制变为e–lgp图的曲线，造成数学假象。但e–lgp图可以较明显地将低应力范围的SL与高应力范围的VCL区分开来，使用起来更加方便。
Casagrande[22]提出了根据经验绘制e–lgp曲线并计算Pc的方法：确定最大曲率点并绘制水平线和切线夹角的平分线，VCL的反向延长线与平分线相交处的应力即为Pc。Casagrande法是国际上普遍采用的，也是我国《土工试验方法标准》[23]推荐的方法，但该方法通过人工绘图的方法得到试验压缩曲线，最大曲率点需主观判定，影响了Pc计算结果的准确性。此外，Casagrande法需采用原状土进行固结试验，在扰动土或破坏相对较大的原状土上存在局限性。
在对土壤压缩曲线的研究过程中，多种拟合方法被提出并应用。Baumgartl和Köck[24]使用van Genuchten方程对土壤压缩曲线进行建模，并通过试验证明了其可行性，但是与Casagrande法相比，该方法得到的压缩指数值偏高[21]。Gregory等[20]用四次多项式、Logistic和Gompertz方程拟合12种土壤的0 ~ 100 kPa和0 ~ 200 kPa单轴压缩试验数据，结果发现3种方法的拟合效果都比较好，但是相对而言，Gompertz方程得到的土壤压缩特性精度更高。Keller等[17]采用Gompertz方程拟合曲线研究了土壤颗粒和团聚体尺寸大小对一系列瑞典土壤压缩特性的影响，An等[10]、霍连飞等[14]也采用Gompertz方程分别研究了黑土和水稻土的压缩特性，曲线拟合效果较好。
估计预固结压力的不同方法其计算结果有所不同，如图 2所示，需要使用者仔细选择。Casagrande法：在最大曲率点绘制水平线和切线夹角的平分线，VCL的反向延长线与平分线相交处的应力[22]；初始截距法V–I(Virgin–Initial intercept)：原始压缩曲线VCL与零应变水平线的交点对应的应力[25]；最大曲率法MC(Maximum Curvature)：数学方法计算出拟合曲线的最大曲率，得到对应点的应力[20]。
|
图 2 三种预固结压力的计算方法示意图 Fig. 2 Three calculation methods of precompression stress (Pc) |
Cavalieri等[21]认为Casagrande法会高估预固结压力的值，而V–I法会造成低估，这与Gregory等[20]的结论相同。Gregory不仅用四次多项式、Logistic和Gompertz方程这3种方法拟合了12种土壤的压缩曲线，还用Casagrande法、V–I法和MC法分别估算了Pc，结果表明，Casagrande法的预测值较高，V–I法预测值较低，MC法位于两者之间，并且基于Gompertz方程的MC法预测值与已知应力的绝对偏差最小，预测准确性和稳定性较高。因此，Gompertz方程被认为是拟合压缩曲线和计算压缩特性的更优模型，由该函数计算最大曲率点从而得到的预固结压力精度最高。
1.3 压缩特性的影响因素土壤并非理想的弹塑性体，压缩特性是土壤具有的弹性、黏性、塑性的综合反映。弹性变形来自土壤结构、颗粒和孔隙水/薄膜水的弹性，封闭气体的压缩和溶解；塑性变形来自土壤颗粒相对位移及结构破坏、孔隙水和气体排出[26]。因此，影响压缩特性的土壤自身因素包括：质地、有机质含量、团聚体尺寸、黏粒含量及含水量等物质组成，以及容重、孔隙度、硬度等土壤性状。
多数研究[10, 27-29]认为，预固结压力Pc与容重正相关，与含水量、有机质和黏粒含量负相关。但也有研究[30]认为，Pc与土壤容重、含水量和黏粒含量显著相关，不受有机质含量影响。
对于压缩指数Cc的研究目前尚无定论。Reichert等[27]的研究表明，容重增大导致压缩指数降低，有机质和水分含量增加会导致压缩指数增加。林琳等[31]对配制不同有机质含量的重塑土进行压缩特性分析，结果表明压缩指数随有机质含量增加而增大，且高含水量时二者极显著正相关。An等[10]对中国东北的黑土和棕壤比较表明，压缩指数与容重呈负相关，但只有黑土的压缩指数与含水量呈正相关，棕壤的压缩指数与含水量的相关性并不显著。Imhoff等[30]认为压缩指数、土壤容重和黏粒含量间的关系具有统计学意义，但不受有机质和含水量影响。Arthur等[32]则发现土壤预固结压力值和压缩指数均随着有机质含量的增加而增大。
Gregory等[33]发现，对于草地土壤与耕地土壤，有机质含量和回弹指数Cs之间有强烈的正相关关系，并且在草地土壤中，较大的土壤颗粒可能会增加弹性恢复力。林琳等[31]的研究中，回弹指数随着有机质含量的增加表现出下降的趋势，但未达到显著水平。
综上所述，不同学者对于不同土壤的预固结压力Pc、压缩指数Cc和回弹指数Cs等压缩特性与有机质含量、含水量、容重和黏粒含量关系的研究结果不尽相同。有研究者[26]认为，这些影响因素相互影响、相互制约，表现出了极强的交互作用，因此某单一因素的影响程度可能会被其他因素掩盖，如林琳等[31]发现含水量不同，有机质与黏粒对黑土压缩–回弹特性的影响不同。因此对土壤压缩特性及其影响因素的关系、影响因素之间的交互作用和影响机理还需要进一步的研究。
2 土壤压实形变机理土壤是由固相、气相和液相组成的三相体，在外力作用下导致的土壤压实涉及到颗粒的重新排列和气体、液体的排出，其机理存在多种理论假说[34-37]。
Proctor[34]在1933年最先提出了毛管润滑理论，他指出土壤中的水在压实过程主要发挥毛细作用和润滑作用，认为含水量增加使颗粒表面结合水膜增厚，发挥润滑作用减小颗粒间摩擦，土体易于压实；而当含水量继续增加到占据空气孔隙时，压实过程中自由水不易排出，压实效果变差。一定压实能量下，土体干密度随含水量增加的变化曲线表现为单峰曲线，峰值点对应的纵坐标值为最大干密度，对应的横坐标值为最优含水量。
Hilf[35]在1956年将孔隙水压力理论应用于非饱和土，发现干密度随含水量的变化规律与孔隙比的变化相关，并提出当含水量位于最优含水量湿侧时，被封闭在土中的气体无法排出，承担了一部分压实能量，使土体干密度变小。
1958年，Lambe[36]根据土的结构和土粒的表面物理化学效应分析了压实作用对土体微观结构的影响，提出了基于土体物理化学特性的理论。该理论认为当含水量位于最优含水量干侧时，围绕黏粒的双电层不能充分发展，减小了颗粒间的溶质排斥作用，土体形成絮凝结构，粒间阻力大，不易压实；含水量继续增加，双电层扩大，颗粒间阻力减小，干密度增加；但当含水量位于最优含水量湿侧时，双电层进一步扩大，土体结构分散，孔隙水的增加使单位体积土体中颗粒相对减少，干密度降低。
Olson[37]在毛管润滑和表面物理化学等理论的基础上，提出了有效应力原理理论。该理论认为当含水量较低时，颗粒间黏结力强，颗粒排列方向不规则，土体即使在外力作用下也不易被破坏；当含水量增加时，土体中与外界连通的气体能够排出，土粒表面水膜逐渐变厚，发挥润滑作用，干密度增大；含水量超过最优含水量后，孔隙中水的体积增加，封闭气体不能排出，孔隙压力增大，抵消了部分外力作用，压实效果变差。
经过多年的理论发展，目前认为压实变形过程可以结合水膜润滑及电化学性质等理论来解释[4]。在土体中含水量较低时，由于土壤颗粒表面的结合水膜较薄，颗粒间距较小，颗粒间电作用力以引力占优势，土壤颗粒间的摩擦力、黏结力都较大，所以土壤颗粒的相对位移阻力大且难以克服，压实效果就差。随着土壤中含水量增加，结合水膜增厚，土壤颗粒间距也逐渐增加，水在颗粒间起到润滑作用，内摩擦力相对减小，外力比较容易克服颗粒间引力而使土壤颗粒相互位移，而趋于密实，压实效果较好。当土中含水量继续增加时，虽然可使颗粒间引力减小，但压实时孔隙中过多的水分不易立即排出，阻碍了土壤颗粒的靠拢，同时无法排出的气体以封闭气泡的形式存在于土体内部，很大一部分外力由孔隙气体承担并转化为孔隙压力，颗粒间所受的力减小，压实仅能导致土壤颗粒更高程度的定向排列，而土体几乎不发生体积变化，所以压实效果反而下降。
土壤压实形变机理十分复杂，因为土壤并非理想的弹塑性体。土壤压实及其压缩特性是土壤固有的弹性、黏性、塑性的综合反映，利用单一的土壤压实指标往往不能很好表征土壤的压缩特性，而且常用的预固结压力、压缩指数、压缩系数等难以应用于不同土壤类型间的比较。因此，需要充分运用数学模型等新方法，结合土力学和土壤物理学的研究成果，开展土壤压实过程模拟研究，深入理解压实形变的机理，并定量探究不同因素对压实的交互影响，有助于提出适用于不同土壤类型间可比较的压实指标。
3 土壤压实过程模拟在土壤压实研究中，模型模拟是评价和量化土壤压实的一种有效手段，既可以整合现有知识，量化土壤压实所涉及到的各个过程，还可用于预测特定土壤对压实的敏感性并评估其压实风险。
考虑到农业生产的实际需求，土壤压实过程模型所需的输入参数或变量主要是农机具的重量、轮胎尺寸与充气压力和土壤条件等。通过模型获得的输出主要包括：①农机具田间作业后的土壤压实程度及其最大值；②农机具施加载荷造成的土壤下陷深度及其破坏程度；③压实引起的土壤物理参数变化[38]。据此，研究者们提出了一系列的模拟模型，用于分析轮胎作用在土壤上的应力和应变特性。
根据其构建原理，土壤压实过程模型可分为机理模型和经验模型两类[39]。经验模型往往需要通过大量实验数据得出经验性的关系，通常仅适用于特定场景。机理模型则通过农机具车轮与土壤接触面的应力分析对压实情况进行模拟和预测，重点关注：①接触面积及其载荷分布；②应力在土壤中的传播；③应力引起的土壤体积变化。目前，所有机理模型均假设土壤是均质、各向同性且连续的介质。根据描述应力在土壤中传播的方式，机理模型可以分为近似分析模型(Pseudo-analytical Model)和有限元模型(the Finite Element Model, FEM)[40]。近年还有研究者将土壤看作由离散颗粒(土壤颗粒)、液体(土壤水分)和气体(土壤空气)组成的复杂的颗粒团聚体，将广泛应用于岩石力学的离散单元法(Distinct Element Method, DEM)引入了土壤力学领域[41]，但目前应用于土壤压实模拟的研究较少[42]。
3.1 近似分析模型现有的分析模型均基于Boussinesq方程。该方程由Boussinesq[43]在1885年提出，之后Fröhlich[44]在1934年引入了集中系数(concentration factor)，Söhne[45]在1953年将其应用于农业土壤研究。Boussinesq方程假设土壤是均质、各向同性、半无限、理想的弹性介质，后来Fröhlich考虑到实际土壤介质是弹塑性体甚至是黏弹塑性体，引入了集中系数以更精确地描述应力在土壤中的分布。
3.1.1 SOCOMO模型van den Akker等 [46]提出了SOCOMO模型，用于评估土壤失效状态。该模型综合了轮胎参数(数量、直径、宽度及充气压力)和土壤参数，根据Söhne方程计算土壤应力和底土承载力，通过比较两者大小来判断土壤是否失效，可为防止下层土壤压实和选择轮式农机具提供依据。但是在土体为干燥开裂的土壤时，SOCOMO模型评估效果不佳，因为干裂土不能简化为均匀介质[47]。
3.1.2 Compsoil模型O'Sullivan等[48]设计了Compsoil模型来预测土壤压实。首先根据轮胎参数和土壤条件，计算轮胎–土壤接触面积并转换为等效的圆形面积；然后用Söhne方程获得土壤垂直应力；最后在考虑土壤容重和压缩指数及含水率对其影响的基础上，计算出土壤受压后产生的形变即车辙深度(rut depth)。Compsoil模型为Excel表格形式，计算核心为VBA宏，便于非专业人士如学生和农民的使用。但是Söhne方程是针对均匀半无限介质建立的，容易低估松散表土覆盖在致密底土上的应力，因此该模型更适用于比较不同农机具的压实，不适用于预测特定土壤的实际容重。有研究者[49]使用Compsoil模型模拟田间轮耕试验，与实测数据对比表明，该模型可用于模拟田间干容重和车辙深度随土壤含水量增加的趋势，但是与实测值相比干容重和车辙深度均有实质性的差异。
3.1.3 SoilFlex模型Keller[50]在2005年提出了SoilFlex模型，是能够同时预测应力传播和体积密度变化的二维模型。该模型与Compsoil模型结构相同，首先描述轮胎–土壤接触面的应力分布(SoilFlex模型增加了剪切应力计算)，然后计算了应力在土体中的传播(Söhne方程)，最后根据应力计算土壤变形。但Compsoil模型将接触面视作圆形，SoilFlex模型用超椭圆描述。该模型同样是用Excel表格形式实现，并包含了多个VBA，用户可以根据需求选择相应的宏和输入输出参数，例如比较不同农机具的压实影响时，计算垂直应力即可。
SoilFlex模型较完整地解析了土壤压实过程，并提供了相应压实程度评价指标的计算方法，用户可以通过输入载荷、轮胎充气压力等参数，得到土体竖向应力分布、容重和车辙深度等，预测农机具压实土壤状况，帮助选择适宜的田间作业农机具。Lozano等[51]使用SoilFlex对巴西甘蔗收获期原状土样进行建模，结果表明表层0 ~ 20 cm土壤的板结风险较高，运输车辆应在含水率低于16% 时候通行以降低压实风险。在预测车辙深度时，SoilFlex模型[52]考虑了剪应变和侧向应变(之前模型未考虑的)的影响，但是土体变形的计算结果仍强烈依赖于土力学特性如压缩指数，而实际中往往难以直接获取这些有代表性的土力学特性。此外，SoilFlex模型的Excel表格包含70个子表，虽然增加了输出结果的选择性，但是使用过程较为复杂。
3.1.4 FRIDA模型Schjønning等[53]提出了FRIDA模型，该模型是对描述轮胎–土壤接触区域应力的Keller方程的细化和组合，可更详细地量化接触面积和整个接触区域的应力分布。接触区域同样用超椭圆描述，但是将应力方程标准化为单位接触面积和单位应力，使应力分布系数α、β不受轮胎形状影响，提高了模型应用范围，但是FRIDA模型在不同轮胎参数与土壤条件下的适应性尚不明确，仍需进一步优化参数和模型。王宪良等[54]验证了该模型在我国华北地区一年两熟壤土条件下的适应性，模拟了轮胎–土壤接触特性，通过预测土壤应力揭示轮胎气压和轴载对土壤压实的影响，为拖拉机轮胎气压及配套农机具选择提供了依据。
3.1.5 Terranimo模型Lassen等[55]在2013年基于Web开发了Terranimo模型，该模型需要农机具、场地特性、土壤质地和含水率等输入参数，输出轮胎–土壤接触面积、接触面应力和剖面应力分布等值线图。其中农机具参数可从数据库中选择，土壤质地和含水率数据可以由用户手动提供，也可以从系统提供的国家土壤和天气数据库中选择。网页模型使用更加方便，而且该模型可以模拟不同轮胎类型、轮胎充气压力和轮重的组合影响，因此可作为风险评估工具供农户和研究者使用。但是该模型输出参数中不包括土壤干容重和车辙深度等参数，无法直观评价土壤功能，仅能评估土壤应力是否会超过土壤机械强度上限；此外，关联的国家数据库仅包括丹麦、挪威、芬兰、比利时、英国和法国，对于其他国家的用户适应性较差。
3.1.6 soilphysicsda Silva和de Lima[56]使用R语言开发了soilphysics包，经过不断修改[57]，参考了SoilFlex模型计算土壤应力和容重、Terranimo模型评估压实风险，可用于模拟农田交通引起的土壤压实。其中，stressTrafic()函数计算接触应力、接触区域应力分布，以及应力在土壤剖面中的传播；soilDeformation()函数计算土壤容重的变化量，其中所需的土力学特性参数如压缩指数，可以用soilStrength()和compressive_ properties()函数计算得到。
soilphysics简化了计算农机具压实后轮胎–土壤接触面的应力分布和土体变形的步骤，并且它是开源的R语言包，用户可以方便地获得源代码，并根据自己需求调整代码获得更多输出。但是在soilDeformation()函数使用的压缩参数均是按照平均正应力条件计算的，这一点与Compsoil模型相似，而没有参照SoilFlex模型考虑剪应力。剪应力和侧向应变已经被证实对轮胎–土壤相互作用有显著影响[58]，因此在某些情况下该模型与现场实测值的吻合度不高，需要进一步完善。soilphysics的目的是简化模拟土壤压实的计算过程，方便非专家用户如学生和农户使用，如果更新模型细节以提高精度，可能会增加模型复杂度，反而降低其实用性。
近似分析模型软件的形式从公式到Excel表格、到Web网页再到R包，上边界条件中轮胎–土壤接触面形状的假设从圆形到超椭圆，应力传播的计算从只考虑平均正应力到考虑剪应力，但模型核心仍然是改自Boussinesq方程和Söhne方程，因此分析模型虽然计算相对简单，但是存在几个无法忽视的缺点：①集中系数需要估算；②难以计算多层土壤；③很少考虑剪应力和侧向应变对土壤变形的影响。因此未来的压实分析模型需要进一步考虑以下问题：模型上边界条件(即轮胎–土壤接触面的大小和形状，接触应力的大小和分布)及土壤强度对边界条件的影响；土体应力的计算需要考虑表土和心土强度差异的解析方法，以替代目前的单层分析方法；土壤压实评价指标需要进行更多研究，因为室内试验测得的应力–应变行为与原位传感器测量结果有较大差异，压缩特性如预固结压力在评估田间土壤压实程度时存在缺陷，压实评价过程需要更多地考虑剪切特性的影响等。
3.2 有限元模型(FEM)计算载荷产生的应力在土壤中的传播时，有限元模型将土壤离散为细小的单元，建立单元的方程再整合为描述土壤特性的近似方程，可以实现三维应力分析，也可以通过假设接触面形状轴对称来简化为二维问题。FEM所需参数较多但是对应力–应变行为的描述更加准确。
早期的土壤压实有限元模型主要是以单向、小变形和线弹性变形为主，通过杨氏弹性模量E和泊松比μ描述应力和应变状态。但是在压实过程中，耕作的土壤表现为非线性和塑性变形，这两个弹性参数不断变化，因此研究人员引入随应力水平变化的正切杨氏弹性模量Et和泊松比μ[59]。
研究者应用FEM时，多通过仿真软件建立土壤模型，进而模拟在外部载荷作用下土壤内部的动态过程。如Cui等[60]利用PLAXIS程序建立轮胎–土壤界面有限元模型，模拟垂直应力分布。Cui在土壤表面使用具有不同弯曲刚度的圆柱体模拟车轮的作用，圆柱的抗弯刚度将应力分布形状从均匀变为不均匀(砂土为抛物线形，黏土为U形)，可改进基于FEM模型的土壤压实预测效果。刘姝和李永奎[61]进行三轴压缩试验测定土壤内摩擦角φ、黏聚力c、弹性模量E以及应力应变曲线，利用Solidworks建立车轮和土壤模型，将实测土壤特性参数输入到模型中，利用插件Solidworks Simulation对车轮压缩土壤进行有限元分析，得到了车轮与土壤相互作用的应力、应变和位移分布图。
有限元模拟的缺点是使用者需要具有一定的弹塑性力学和土力学理论背景，一般使用者难以根据三轴试验参数运用软件准确建立土壤和轮胎模型，得到真实度高的模型和精确的结论。而且普通的有限元模型局限于连续介质理论，将土壤看作连续体进行分析，仅能对最终状态为静止的土壤进行整体分析。
3.3 离散单元法DEM一些学者在研究土壤行为时，将土壤看作由离散颗粒(固体颗粒)、液体(土壤水分)和气体(土壤空气)组成的复杂的颗粒团聚体，从土壤内部颗粒细观形态出发研究土壤复杂的动态行为。DEM是Cundall[62]在1971年提出的一种分析离散物质的不连续数值模拟方法，将介质看作由一系列离散的独立运动的单元所组成，单元运动受经典运动方程控制，整个介质的变形和演化由各单元的运动和相互位置来描述。土壤本身为离散结构，且耕作是土壤颗粒的动态变形和运动过程，因此采用DEM对耕作土壤进行分析模拟更符合实际。
DEM被提出后，1986年引入我国[63]，在采矿工程、岩土工程及水利水电工程等领域发展迅速。近年来，DEM在优化农机作业部件[64]、揭示农机触土部件与土壤颗粒的互作机理[65]等方面表现出独特优势。DEM也被应用于土壤压实相关研究中，赵振家等[66]通过离散元法模拟贯入试验和车轮土壤压实过程，分析车轮压实对土壤硬度的影响，经试验验证仿真结果可信，且仿真模拟可直观地分析压实土壤内部的受力分布和土壤颗粒间接触力情况。石林榕等[67]基于“快速填充法”(将有限元网格及离散元颗粒群建模结合)建立农田土壤单轴压缩模型，将仿真结果与试验结果对比，压缩曲线基本一致。王宪良[68]利用ECM模型(the Edinburgh Elasto-Plastic Cohesion Model)对轮胎–土壤相互作用进行模拟，研究轮胎运动过程中土壤应力传递规律，并基于仿真结果建立以轮胎轴载和压实次数为自变量的土壤应力传递模型，发现土壤压实次数对同一位置的土壤应力增大趋势呈对数函数增长。王宪良等[69]还基于敏感性分析和代理模型理论，提出了一种通用的离散元土壤模型的参数标定优化方法，使用优化后的参数建立的轮胎–土壤仿真模型预测的接触应力误差在5% 左右。
DEM克服了连续介质理论需要复杂边界条件的缺点，充分发挥了离散物质细观分析方法的优点，突破了将土壤看作连续体的限制，在分析土壤破裂变形的动态行为方面表现出了极大的优势，还能直观展现土壤受压后接触面的应力变化。但DEM应用于研究农机与土壤互作关系时，其重点和难点在于如何选择正确的土壤颗粒接触模型及合适的接触模型参数标定方法，这要求使用者对土壤及其颗粒的力学行为有深入的了解，并准确确定离散元仿真参数(主要包括材料参数和接触参数)。而且，DEM是为了分析离散颗粒的集合而建立的，需要使用者设定土壤颗粒大小和数量(通常以万为单位)，因此即使仅考虑有限颗粒数量的模型也需要大量的计算求解时间，最终的模型结果对实践应用往往仅具有参考价值。
综上所述，建立土壤压实模型是定量分析轮胎–土壤接触面的应力–应变关系和理解土壤压实过程的有效方法。基于Söhne方程的近似分析模型，可通过轮胎参数和土壤条件模拟预测二维应力传播和土壤变形过程。该模型所需参数相对较少，使用简单，是目前应用最广泛的模型，但是其未考虑剪应力和侧向应变对土壤变形的影响。有限元模型可实现三维分析、动态模拟，但是建立精确模型需要输入大量三轴试验参数，具体使用过程较为复杂。基于离散单元法的数值模型可更准确地分析土壤破裂变形的行为，揭示微观尺度农机与土壤的互作关系，但是同样需要输入大量的仿真参数，最终的模拟结果受限于颗粒大小和数量与实际之间的差别，往往仅有一定的参考意义。
4 总结与展望土壤压实影响土壤物理、化学和生物过程，与土壤肥力、作物生长和耕地功能直接相关，而在农业机械化快速发展的今天，农田土壤压实问题广泛存在，限制了土壤生产功能和农业可持续发展。国内外农业科研工作者已经关注到土壤压实问题并开展了许多研究，但是目前的研究以压实对土壤物理特性和作物产量的影响为主，尚无法系统全面地阐明农业土壤的压实过程及其机理。未来的研究工作需要重点围绕以下几个方面：
1) 开展细观颗粒/孔隙尺度的土壤压实过程和机理研究。土壤并非理想的弹塑性体，而是由固体颗粒、液体和气体组成的复杂的颗粒团聚体，土壤压实变形过程和宏观尺度上表现出的压缩特性(如预固结压力、压缩指数、回弹指数等)实质上是其弹性、黏性、塑性形变的综合反映。阐明外界作用力的大小、速率和方向等因素对颗粒排列组合、孔隙分布和结构特征等的影响机制，需借助新的观测、模拟技术从颗粒或孔隙尺度上开展细观和微观探究。
2) 改进土壤压缩曲线和压缩特性的数学描述和计算方法。土壤压缩特性可量化土体在外部荷载作用下的压缩和恢复行为，常利用室内单轴试验得到的土壤压缩曲线来估算。但目前尚无统一的、具有明确物理意义的土壤压缩曲线函数关系，以及可完全排除人为因素的由压缩曲线确定压缩特性(特别是预固结压力)的方法。此外，压缩特性同时受到外力作用方式和土壤自身性状(有机质含量、黏粒含量、容重、含水量等)等的控制，研究具有通用性的、可根据土壤性状估测压缩特性的土壤传递函数，有助于对不同土壤条件和受力场景下的压实特性进行定量分析。
3) 完善土壤压实过程的仿真模型和风险评估模型。当前农业土壤压实过程研究中，非侵入性和原位传感器监测方法受成本高昂、技术复杂等的限制，尚无法广泛应用，利用模拟模型来量化和评估压实风险具有重要意义。现有模型面临的主要问题是实现易用性和准确性之间的平衡。例如，近似分析模型虽然仅考虑了垂直应力的影响，但输入参数少、应用便捷；有限元和离散元方法可全面、直观地模拟压实形变过程，但模型参数众多、专业性强、使用门槛高。如何兼顾全过程仿真模拟和实际场景应用的需求，建立适于区域尺度农田压实风险评估的压实评价模型，是该领域下一步的研究重点。
| [1] |
Gürsoy S. Soil compaction due to increased machinery intensity in agricultural production: Its main causes, effects and management[M]//Ahmad F, Sultan M. Technology in Agriculture. IntechOpen, 2021.
(  0) |
| [2] |
Soil Science Society of America. Glossary of soil science terms[M].
American Society of Agronomy, Madison, 2008
( 0) |
| [3] |
Nawaz M F, Bourrié G, Trolard F. Soil compaction impact and modelling. A review[J]. Agronomy for Sustainable Development, 2013, 33(2): 291-309 DOI:10.1007/s13593-011-0071-8 ( 0) |
| [4] |
刘松玉. 土力学(5版)[M].
中国建筑工业出版社, 北京, 2020
( 0) |
| [5] |
Schjønning P, van den Akker J J H, Keller T, et al. Driver-pressure-state-impact-response (DPSIR) analysis and risk assessment for soil compaction—A European perspective[M]//Sparks D L. Advances in Agronomy. Amsterdam: Elsevier, 2015: 183–237.
( 0) |
| [6] |
周艳丽, 卢秉福. 农田机械压实对土壤物理特性的影响[J]. 中国农机化学报, 2018, 39(9): 66-70 ( 0) |
| [7] |
McKenzie R H. Agricultural soil compaction: Causes and management[EB/OL]//Alberta. (2010) [2022-12-27]. https://open.alberta.ca/publications/5333887.
( 0) |
| [8] |
严磊, 张中彬, 丁英志, 等. 覆盖作物根系对砂姜黑土压实的响应[J]. 土壤学报, 2021, 58(1): 140-150 ( 0) |
| [9] |
乔金友, 霍东旭, 张险峰, 等. 中型拖拉机压实对测试截面土壤坚实度和含水率的影响[J]. 东北农业大学学报, 2021, 52(6): 87-96 ( 0) |
| [10] |
An J, Zhang Y L, Yu N. Quantifying the effect of soil physical properties on the compressive characteristics of two arable soils using uniaxial compression tests[J]. Soil and Tillage Research, 2015, 145: 216-223 DOI:10.1016/j.still.2014.09.002 ( 0) |
| [11] |
李春林, 丁启朔, 陈青春. 水稻土的先期固结压力测定与分析[J]. 农业工程学报, 2010, 26(8): 141-144 ( 0) |
| [12] |
张军昌, 师帅兵, 朱瑞祥, 等. 轮式作业机械对农田土壤压实的模拟试验研究[J]. 农机化研究, 2012, 34(1): 161-164 ( 0) |
| [13] |
丁肇, 李耀明, 唐忠. 轮式和履带式车辆行走对农田土壤的压实作用分析[J]. 农业工程学报, 2020, 36(5): 10-18 ( 0) |
| [14] |
霍连飞, Belal E A, 丁启朔, 等. 采用Gompertz函数的水稻土压缩特性研究[J]. 农业工程学报, 2020, 36(12): 136-142 ( 0) |
| [15] |
Horn R, Lebert M. Soil compactability and compressibility[M]//Soane B D, van Ouwerkerk C. Developments in agricultural engineering. Amsterdam: Elsevier, 1994: 45–69.
( 0) |
| [16] |
Keller T, Arvidsson J. Compressive properties of some Swedish and Danish structured agricultural soils measured in uniaxial compression tests[J]. European Journal of Soil Science, 2007, 58(6): 1373-1381 DOI:10.1111/j.1365-2389.2007.00944.x ( 0) |
| [17] |
Keller T, Lamandé M, Schjønning P, et al. Analysis of soil compression curves from uniaxial confined compression tests[J]. Geoderma, 2011, 163(1/2): 13-23 ( 0) |
| [18] |
Lebert M, Horn R. A method to predict the mechanical strength of agricultural soils[J]. Soil and Tillage Research, 1991, 19(2/3): 275-286 ( 0) |
| [19] |
Keller T, Arvidsson J, Schjønning P, et al. In situ subsoil stress-strain behavior in relation to soil precompression stress[J]. Soil Science, 2012, 177(8): 490-497 DOI:10.1097/SS.0b013e318262554e ( 0) |
| [20] |
Gregory A S, Whalley W R, Watts C W, et al. Calculation of the compression index and precompression stress from soil compression test data[J]. Soil and Tillage Research, 2006, 89(1): 45-57 DOI:10.1016/j.still.2005.06.012 ( 0) |
| [21] |
Cavalieri K M V, Arvidsson J, da Silva A P, et al. Determination of precompression stress from uniaxial compression tests[J]. Soil and Tillage Research, 2008, 98(1): 17-26 DOI:10.1016/j.still.2007.09.020 ( 0) |
| [22] |
Casagrande A. The determination of pre-consolidation load and its practical significance[C]. 1st International Conference on Soil Mechanics and Foundation Engineering, Cambridge, MA: Harvard University, 1936(3): 60–64.
( 0) |
| [23] |
水利部水利水电规划设计总院, 南京水利科学研究院. 土工试验方法标准[M].
中国计划出版社, 北京, 2019
( 0) |
| [24] |
Baumgartl T, Köck B. Modeling volume change and mechanical properties with hydraulic models[J]. Soil Science Society of America Journal, 2004, 68(1): 57-65 DOI:10.2136/sssaj2004.5700 ( 0) |
| [25] |
Arvidsson J, Keller T. Soil precompression stress[J]. Soil and Tillage Research, 2004, 77(1): 85-95 DOI:10.1016/j.still.2004.01.003 ( 0) |
| [26] |
肖质秋, 虞娜, 安晶, 等. 土壤压实及有机质对其影响的研究进展[J]. 土壤通报, 2019, 50(5): 1253-1260 ( 0) |
| [27] |
Reichert J M, Mentges M I, Rodrigues M F, et al. Compressibility and elasticity of subtropical no-till soils varying in granulometry organic matter, bulk density and moisture[J]. Catena, 2018, 165: 345-357 DOI:10.1016/j.catena.2018.02.014 ( 0) |
| [28] |
Keller T, Lamandé M, Naderi-Boldaji M, et al. Soil compaction due to agricultural field traffic: An overview of current knowledge and techniques for compaction quantification and mapping[M]//Saljnikov E, Mueller L, Lavrishchev A, et al. Advances in understanding soil degradation. Cham: Springer, 2022: 287-312.
( 0) |
| [29] |
孙奥博, 安晶, 虞娜, 等. 铁矿尾砂配施有机物料对褐土压缩及回弹特性的影响[J]. 土壤, 2022, 54(6): 1300-1306 DOI:10.13758/j.cnki.tr.2022.06.026 ( 0) |
| [30] |
Imhoff S, Da Silva A P, Fallow D. Susceptibility to compaction, load support capacity, and soil compressibility of hapludox[J]. Soil Science Society of America Journal, 2004, 68(1): 17-24 DOI:10.2136/sssaj2004.1700 ( 0) |
| [31] |
林琳, 韩少杰, 王恩姮. 有机质与黏粒含量对黑土压缩-回弹特性的影响[J]. 土壤学报, 2016, 53(5): 1138-1147 ( 0) |
| [32] |
Arthur E, Schjønning P, Moldrup P, et al. Soil resistance and resilience to mechanical stresses for three differently managed sandy loam soils[J]. Geoderma, 2012, 173/174: 50-60 DOI:10.1016/j.geoderma.2012.01.007 ( 0) |
| [33] |
Gregory A S, Watts C W, Griffiths B S, et al. The effect of long-term soil management on the physical and biological resilience of a range of arable and grassland soils in England[J]. Geoderma, 2009, 153(1/2): 172-185 ( 0) |
| [34] |
Proctor R. Design and construction of rolled earth dams[J]. Engineering News Record, 1933, 31(8): 286-289 ( 0) |
| [35] |
Hilf J W. Compacting earth dams with heavy tamping rollers[J]. Journal of the Soil Mechanics and Foundations Division, 1957, 83(2): 1-28 ( 0) |
| [36] |
Lambe T W. The structure of compacted clays[J]. Journal of the Soil Mechanics and Foundations Division, 1958, 84(2): 1-34 ( 0) |
| [37] |
Olson R E. Settlement of embankments on soft clays: (the thirty-first terzaghi lecture)[J]. Journal of Geotechnical and Geoenvironmental Engineering, 1998, 124(8): 659-669 DOI:10.1061/(ASCE)1090-0241(1998)124:8(659) ( 0) |
| [38] |
Raper R L, Johnson C E, Bailey A C. Coupling normal and shearing stresses to use in finite element analysis of soil compaction[J]. Transactions of the ASAE, 1994, 37(5): 1417-1422 DOI:10.13031/2013.28222 ( 0) |
| [39] |
O'Sullivan M F, Simota C. Modelling the environmental impacts of soil compaction: A review[J]. Soil and Tillage Research, 1995, 35(1/2): 69-84 ( 0) |
| [40] |
任利东, 王丽, 林琳, 等. 农田土壤机械压实研究进展与展望[J]. 土壤学报, 2023, 60(3): 610-626 ( 0) |
| [41] |
张锐, 李建桥, 李因武. 离散单元法在土壤机械特性动态仿真中的应用进展[J]. 农业工程学报, 2003, 19(1): 16-19 ( 0) |
| [42] |
曾智伟, 马旭, 曹秀龙, 等. 离散元法在农业工程研究中的应用现状和展望[J]. 农业机械学报, 2021, 52(4): 1-20 ( 0) |
| [43] |
Boussinesq J. Application des potentiels : à l'etude de l'equilibre et du mouvement des solides élastiques[C]//Application of potentials in the study of the equilibrium and the movement of elastic bodies. Paris: Gauthier-Villars, 1885: 30.
( 0) |
| [44] |
Fröhlich O K. Druckverteilung im baugrunde: Mit besonderer berücksichtigung der plastischen erscheinungen[M].
Springer, Vienna, 1934
( 0) |
| [45] |
Söhne W. Druckverteilung im ackerboden und verformbarkeit des ackerbodens[J]. Kolloid-Zeitschrift, 1953, 131(2): 89-96 DOI:10.1007/BF01512930 ( 0) |
| [46] |
van den Akker J J H, Stuiver H J. A sensitive method to measure and visualize deformation and compaction of the subsoil with a photographed point grid[J]. Soil and Tillage Research, 1989, 14(3): 209-217 DOI:10.1016/0167-1987(89)90009-3 ( 0) |
| [47] |
van den Akker J J H. SOCOMO: A soil compaction model to calculate soil stresses and the subsoil carrying capacity[J]. Soil and Tillage Research, 2004, 79(1): 113-127 DOI:10.1016/j.still.2004.03.021 ( 0) |
| [48] |
O'Sullivan M F, Henshall J K, Dickson J W. A simplified method for estimating soil compaction[J]. Soil and Tillage Research, 1999, 49(4): 325-335 DOI:10.1016/S0167-1987(98)00187-1 ( 0) |
| [49] |
Défossez P, Richard G, Boizard H, et al. Modeling change in soil compaction due to agricultural traffic as function of soil water content[J]. Geoderma, 2003, 116(1/2): 89-105 ( 0) |
| [50] |
Keller T. A model for the prediction of the contact area and the distribution of vertical stress below agricultural tyres from readily available tyre parameters[J]. Biosystems Engineering, 2005, 92(1): 85-96 DOI:10.1016/j.biosystemseng.2005.05.012 ( 0) |
| [51] |
Lozano N, Rolim M M, Oliveira V S, et al. Evaluation of soil compaction by modeling field vehicle traffic with SoilFlex during sugarcane harvest[J]. Soil and Tillage Research, 2013, 129: 61-68 DOI:10.1016/j.still.2013.01.010 ( 0) |
| [52] |
Keller T, Défossez P, Weisskopf P, et al. SoilFlex: A model for prediction of soil stresses and soil compaction due to agricultural field traffic including a synthesis of analytical approaches[J]. Soil and Tillage Research, 2007, 93(2): 391-411 DOI:10.1016/j.still.2006.05.012 ( 0) |
| [53] |
Schjønning P, Lamandé M, Tøgersen F A, et al. Modelling effects of tyre inflation pressure on the stress distribution near the soil-tyre interface[J]. Biosystems Engineering, 2008, 99(1): 119-133 DOI:10.1016/j.biosystemseng.2007.08.005 ( 0) |
| [54] |
王宪良, 王庆杰, 李洪文, 等. 基于FRIDA模型的轮胎–土壤接触特性研究[J]. 农业机械学报, 2016, 47(9): 121-127 ( 0) |
| [55] |
Lassen P, Lamandé M, Stettler M, et al. Terranimo ®-A soil compaction model with internationally compatible input options[C]//Sustainable Agriculture through ICT Innovation. Turin, Italy, 2013: 24-27.
( 0) |
| [56] |
da Silva A R, de Lima R P. Aoilphysics: An R package to determine soil preconsolidation pressure[J]. Computers & Geosciences, 2015, 84: 54-60 ( 0) |
| [57] |
de Lima R P, da Silva A R, da Silva Á P. Soilphysics: An R package for simulation of soil compaction induced by agricultural field traffic[J]. Soil and Tillage Research, 2021, 206: 104824 DOI:10.1016/j.still.2020.104824 ( 0) |
| [58] |
Stettler M, Keller T, Weisskopf P, et al. Terranimo®-a web-based tool for evaluating soil compaction[J]. Landtechnik, 2014, 69(3): 132-138 ( 0) |
| [59] |
Raper R L, Erbach D C. Prediction of soil stresses using the finite element method[J]. Transactions of the ASABE, 1990, 33(3): 725-730 DOI:10.13031/2013.31392 ( 0) |
| [60] |
Cui K, Défossez P, Richard G. A new approach for modelling vertical stress distribution at the soil/tyre interface to predict the compaction of cultivated soils by using the PLAXIS code[J]. Soil and Tillage Research, 2007, 95(1/2): 277-287 ( 0) |
| [61] |
刘姝, 李永奎. 基于Solidworks车轮-土壤相互作用的有限元分析[J]. 农机化研究, 2012, 34(11): 42-45 ( 0) |
| [62] |
Cundall P A. A computer model for simulating progressive, large-scale movements in blocky rock systems[C]//Proceedings of the Symposium of the International Society of Rock Mechanics: Vol. 1. Nancy, France, 1971.
( 0) |
| [63] |
王泳嘉. 离散单元法——一种适用于节理岩石力学分析的数值方法[C]//第一届全国岩石力学数值计算及模型试验讨论会论文集. 成都: 西南交通大学出版社, 1988: 42–47.
( 0) |
| [64] |
李艳洁, 吴腾, 林剑辉, 等. 基于离散元法的贯入圆锥对沙土颗粒运动特性分析[J]. 农业工程学报, 2012, 28(24): 55–61, 359 ( 0) |
| [65] |
方会敏, 姬长英, Chandio F A, 等. 基于离散元法的旋耕过程土壤运动行为分析[J]. 农业机械学报, 2016, 47(3): 22-28 ( 0) |
| [66] |
赵振家, 邹猛, 薛龙, 等. 压实对土壤应力分布的影响仿真分析[J]. 农业机械学报, 2012, 43(S1): 311–313, 338 ( 0) |
| [67] |
石林榕, 吴建民, 赵武云, 等. 基于离散单元法农田土壤单轴压缩模型的建立及参数校核[J]. 中国农业大学学报, 2015, 20(4): 174-182 ( 0) |
| [68] |
王宪良. 农机作业土壤压实评价及组合铲松土技术研究[D]. 北京: 中国农业大学, 2018.
( 0) |
| [69] |
王宪良, 钟晓康, 耿元乐, 等. 基于离散元非线性弹塑性接触模型的免耕土壤参数标定[J]. 农业工程学报, 2021, 37(23): 100-107 ( 0) |
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. University of Chinese Academy of Sciences, Nanjing 211135, China