基于渐近均匀化理论的黄土高原草本植物固土效果模拟

黄建坤，王学林，及金楠，陈丽华，张之伟

基于渐近均匀化理论的黄土高原草本植物固土效果模拟

黄建坤，王学林，及金楠，陈丽华，张之伟

（1. 山西吉县森林生态系统国家野外科学观测研究站，北京林业大学，吉县 042200；2. 北京林业大学水土保持学院水土保持国家林业局重点实验室，北京 100083；3. 北京林业大学水土保持学院，北京 100083）

黄土高原是中国水土流失和浅层滑坡灾害最严重的地区之一，而植被根系能有效抑制浅层滑坡。黄土高原草系分布广泛，根系在空间配置上具有鲜明的周期重复特征，计算大量根系的边坡固土效果非常复杂。该文利用渐进均匀化理论，构建“根系-土体”复合土体本构关系，通过与已有研究对比说明，该文所述方法计算三维本构关系具有较高精度。为说明均匀化方法在含根群边坡中的准确性和高效性，该文基于有限元软件ANSYS，进一步模拟了含草边坡和无草边坡的应力和应变分布规律。数值模型设置为3组，分别是无草边坡模型、含草边坡的分离式模型，以及基于该文均匀化理论的含草边坡模型。该文进一步采用观察等效塑性应变区是否贯通作为边坡失稳临界状态的标准，利用强度折减法计算了边坡安全系数。研究表明：1）渐进均匀化理论不仅能准确地构建根土复合体本构关系，而且极大地减少了模型计算工作量（减少了95.58%的单元数）。2）披碱草根系能改善边坡浅层土体的应力场，使得根系分布区内的剪应力更加均匀，使坡体浅层更趋于稳定。3）当边坡坡度较小时（坡角为30°），无草边坡安全系数较大（=4.28），根系对边坡的稳定性加持较小（平均仅提高2.92%）；当坡角增大到45°时，无草边坡安全系数下降为2.90，而含草边坡平均安全系数提升了13.45%，根系固土效果更加显著。

边坡稳定性；根系；固土；均匀化理论；有限元法

0 引 言

中国是世界上水土流失和滑坡灾害最严重的国家之一，山地面积占国土面积的2/3。黄土高原是中国水土流失最严重的地区，滑坡等重力侵蚀相当严重。传统的土木工程护坡方式不仅造价高昂，而且破坏了原有的生态环境。植被是防治水土流失最积极、有效的因素，特别是植物根系固持表层土壤、防止浅层滑坡的作用更是不可忽视。因此，植被措施是一种经济、可靠、环保的滑坡防治措施，具有不可替代的重要作用，在国内外越来越受到重视[1-9]。

通过根系形态分析、拉伸试验、剪切试验和数值模拟等方式，研究人员探讨了根系力学特性[3,10-12]、分布特点[13-14]以及根土相互作用[15]，推动了根系固土力学机理的揭示过程。新材料、新技术的引入，加快了根系固土的研究历程[16-20]。植被根系属于生物有机材料，具有不连续、非均匀、各向异性等特点，其本构关系要比普通工程材料复杂得多。目前，关于根系固土的研究多关注于单根或根束对土体的影响，而此方面的研究还不能对现实中广泛存在的根群固土效应做出全面而准确的评价。植物根系固土护坡能力的研究主要包括了对根强度的研究和根土复合体强度的研究，根与土之间呈弱耦合关系；自然界中根土复合体像钢筋混凝土一样，可以看成一个整体，即一种“新材料”，根土之间呈强耦合关系。在科学和工程方面都迫切需要深入研究根系与土体之间的力学关系，为边坡稳定性评价提供准确理论支持。

黄土高原地区草本和灌木受地质和气候条件的影响，根系具有如下明显的形态特征：根系以竖向根系为主，且较为发达，在地表以下有较长的延伸；根系生长方向较为一致，排列比较均匀，侧根比较少，根系分布大致呈现周期性的特征。根土复合体可以看成长纤维加强的周期复合材料，形成强耦合关系。三轴压缩试验表明，基于均匀化理论构建的根土复合体模型具有计算高效、结果准确的特点[18]。与简单地将根土复合体视为各向同性材料不同，均匀化理论的思想是基于对根系-根周土“单胞”的分析，推导根土复合体的各向异性本构关系，能够将微观根系形态特征量与宏观边坡力学分析结合，建立不同尺度之间的联系，有利于揭示根系固土力学机制。

本文引入周期复合材料的渐进均匀化理论，针对黄土高原浅层滑坡问题，开发用于探讨草本植物根土复合体本构关系的三维力学模型，从周期复合材料的视角推进理解含林（草）边坡的固土机制，解决准确评价根群固土效果的困难，为含林（草）边坡稳定性分析提供理论支撑。

1 渐进均匀化理论

1.1 控制方程

根据山西省黄土高原优势固土护坡植物的生长特性，草本根系的长度远大于根系的直径，因此可以将根系和土壤看成周期性分布的长纤维根土复合材料（如摄于山西吉县森林生态系统国家野外科学观测研究站的图 1）。针对该地区植物根系生长的特点，可以忽略侧根（或须根）的影响，并假定根系在土壤中呈单向周期性分布，因而可以把植物根系和土壤简化成周期性分布的长纤维增强复合材料。因此，该文提出基于均匀化理论的根土强耦合固土机制，先从较为简单的直根系出发，不考虑复杂形态的根系等问题，根土复合体可以由单根-根周土“单胞”表示，如图2所示。

图1 具有周期性长纤维复合材料特征的根系分布

注：a为单胞尺寸，cm。

以具有周期性细观“单胞”结构Y的根土复合体为例，假设在空间坐标系某处存在此复合材料Ω，其全部边界为，并分别受体力和面力作用。根据虚位移原理，控制方法为[21]

式中0与细观尺度无关，描述了根土复合体宏观尺度的位移，1描述了细观尺度的位移，由方程（4）可得

引入[22]

式中,1,2,3代表微观坐标系的坐标，以示与和表示的宏观坐标系区别，根系-土体复合体等效密度可以表示为

1.2 等效刚度求解

将函数用插值函数近似表示为

式中是形函数矩阵；是节点的广义坐标；表示有限元系统中的自由度总数。因此式（6）可写成有限元的标准形式[23]

其中和分别为应变和弹性矩阵。有效刚度矩阵，即方程（8）可表示成

1.3 三维本构关系具体形式

二维平面应变问题的本构关系可以表示为[24]

而自然根系为有限长度，犹如微型锚杆，是典型三维空间问题。其本构关系具有典型的正交各向异性特征，形式为[24]

（15）

为了将二维应变问题升维为三维问题，根据方程[25]

最终可以确定三维根系-土体复合体的正交各向异性的本构关系。

1.4 方法验证

表1 材料参数

表2 材料参数均匀化结果对比

注：EE、E分别指、、方向弹性模量；υ、υ υ分别指、、方向泊松比。下同。

Note: EE, Eare elastic modulus in,,direction; υ υυarePossion’s ratioin、,plane, respectively. The same below.

2 实例数值建模

2.1 植被根系和研究区域选择

本文以披碱草根系为研究对象，播种当年，节根入土深度可达70 cm，第二年能达110 cm以上[18]。在灌溉条件下，虽然根深可达100 cm，但是50 cm土层以下，根系极少，约占总根量的4%[18]。为突出根系固土效果，本文根长统一取为50 cm，根半径取为0.11 cm，根系体积填充率为0.01，“单胞”尺寸取为=2 cm。由于“单胞”具有对称性，计算中仅取1/4“单胞”进行计算，如图3所示。根系和土体的材料参数主要通过文献[18]获得，如表1所示。

图3 “单胞”有限元模型

2.2 本构关系计算

通过本文方法，根土复合体本构关系计算如表3所示，为了便于比较，本文同时给出了文献[25]所述简化方法的计算结果，使简化方法和本文提出的方法相互校核，可以验证结果控制在一定误差范围内，从而避免得到错误的结果。从表3中可以看出，本文方法和简化方法除了在υ计算中存在一定偏差之外，其他弹性参数误差均较小。平面内的泊松比υ对应力分布影响较小，这一点将在数值模拟中得到验证，因此误差在可接受范围。

表3 根土复合体本构参数

注：G、G、G分别表示平面剪切模量。

Note: G,G,Gare shear modulus inplane, respectively

2.3 有限元模型建立

在披碱草根系固土效果数值分析中，将边坡坡度设为30°，坡长5 m。本文所建的模型为3D模型，边坡剖面位于平面，方向上为一个“单胞”厚度=0.02 m，有限元模型如图4。边坡体的底面完全约束，坡面完全自由，两侧平面上的面采用向约束，厚度方向平面上采用对称边界[28]。

模型统一采用SOLID45单元，充分模拟边坡及根系的三维特性，模型假设所有界面理想连接，仅受自重应力影响。为对比含草边坡和无草边坡的应力分布规律，以及本文所述方法的准确性，数值模型设置3组，分别是无草边坡模型、含草边坡的分离式模型，以及基于均匀化理论的含草边坡模型。在含草边坡的分离式模型中，根系分布如图5所示。

注：O、x和y分别为全局坐标系原点和坐标轴；O′、x′和y′分别为局部坐标系原点和坐标轴。

图5 根系边坡分离式模型

虽然该文建立的是单排草根的直根系工况建模，该模型两侧设置为对称边界。对称边界可以模拟边坡在方向无限延伸。因此，该模型本质上模拟的是左右无限扩展的多行多排根系情况。该文分离式模型中，根系分布区已经含有250根草根，能体现含根边坡的真实性和典型性，也有利于说明均匀化理论的高效性。如果继续扩大根系规模，势必难以计算（收敛）。同时，调查发现，黄土高原灌草植被存在一定量的直根系，如图1。该地区的草本和灌木受土壤和气候条件影响，不仅具有发达的竖向根系，可从坡面延伸至深层土坡，总体呈现出生长方向一致、排列均匀的周期性特征。因此，直根系模型也常被用于数值模拟中[15,18]。

在基于均匀化理论的含草边坡模型中，图5中的矩形区域视为根土复合体，采用本文均匀化理论计算得到的参数，即表3中数据，等效密度H= 1 363.6 kg/m3。因此，均匀化参数具有方向性，在数值模拟中，该模型使用了局部坐标系以准确定义计算得到的均匀化参数。

3 数值模拟结果及分析

3.1 边坡应力分布

图6为边坡平面剪应力分布等值线图。从图中可以看出，无草边坡应力分布比较均匀，且与含草边坡模型结果不同；含草边坡分离式模型和含草边坡均匀化模型的结果几乎是一致的，3种工况都在坡脚产生了一定的应力集中现象。由于草根的存在，含草边坡应力分布较为分散。因为根系有较强的抗拉和抗剪刚度，边坡自重荷载产生的应力向根系转移，根系分担了较多的土体应力，并将其扩散到根周土中，因此，改变了自然边坡的应力场。根系明显对边坡应力场分布的影响集中在边坡浅层，特别是根系分布区，即图6中虚线矩形区，使得根系分布区内的剪应力更加均匀，提高了浅层边坡土体稳定性。对比图6b和6c可以看出，本文提出的均匀化理论计算结果和根土分离式模型计算结果具有高度的一致性，充分验证了本文均匀化理论计算结果的准确性，也说明了υ的误差对应力分布影响不显著。

图6 边坡xy平面内剪应力等值线图

图7为边坡von Mises应力分布等值线图。从图中可以看出，3种工况都在坡脚产生了一定的应力集中现象，根系的存在对von Mises应力场的影响比对平面剪应力场的影响小。根系的影响主要在根系分布区，使得根系分布区内的von Mises应力分布更加均匀。

图7 边坡von Mises应力等值线图

3.2 边坡安全系数分析

本文采用强度折减法[5,29-30]对边坡安全系数进行计算。边坡失稳可以通过数值迭代计算是否收敛、特征部位的位移拐点和等效塑性应变区是否贯通进行判断。本文中采用观察等效塑性应变区是否贯通作为边坡失稳临界状态的标准。

无草模型和均匀化模型都使用规则的六面体单元，计算收敛性稳定可靠。但是在分离式模型中，根系直径仅为0.11 cm，划分网格时不可避免地出现三棱锥等低质量单元，极易造成收敛性提前失稳，无法得到准确的安全系数。因此，在本小节分析中，采用Link单元代替实体单元模拟根系状况，提高计算稳定性。根土复合体试验研究表明[16,26,31]，根系对提高土体的内摩擦角十分有限，主要是提高黏聚力，根据文献[18]，本节计算中，设定土体黏聚力=9.56 kPa，内摩擦角=24.716°，根土复合材料黏聚力提高为[31]：c=42.08 kPa，内摩擦角保持不变。

边坡最大等效塑性应变随折减系数的变化，如图8所示。从图8a中可以看出，最大等效塑性应变随折减系数增大而增大。无草（素土）模型在折减系数为4.28时，出现峰值拐点。而3种含草模型的峰值拐点对应的折减系数都比无根模型略大。含草边坡分离式模型由于收敛性问题，最大等效塑性应变曲线与其他含草模型差别较大，且不易得到拐点，因此，不适合用于安全系数求解。

从图8b中可以看出，当坡角从30°增大到45°时，无根边坡最大等效塑性应变随折减系数增大而迅速出现拐点，且最大等效塑性应变的数值也从0.47增至0.77。而含草模型最大等效塑性应变的拐点对应的折减系数都比无根模型大很多，且最大等效塑性应变的数值较低。再次说明根系的存在使得边坡应力分布更加均匀，滑坡的潜在风险降低。图8a和8b中，均匀化模型与含草边坡Link模型的应变曲线变化规律基本一致，说明了均匀化模型的准确性。

本文定义首次出现峰值拐点的折减系数为安全系数，根据图8中的实（空心）点标注，不同工况下边坡的安全系数如表5所示。边坡安全系数与边坡几何参数和材料参数有密切关系[32]。当坡角为30°时，由于坡角较小，边坡安全系数较高。无草边坡的安全系数为4.28，含草边坡均匀化模型与含草边坡Link模型的安全系数分别为4.40和4.41，两者相对误差为0.22%，后者安全系数比前者提高了2.92%，说明草根起到一定的固土作用。受制于草根根系较浅且根系直径较小，根系对边坡稳定性的影响有限。

图8 最大等效塑性应变随折减系数的变化规律

当坡角为45°时，无草边坡安全系数下降为2.90，含草边坡均匀化模型与含草边坡Link模型的安全系数分别为3.28和3.30，两者相对误差为0.69%，而平均安全系数提升了13.45%，说明在高角度边坡情况下，根系固土效果更加显著。由于边坡坡度较小时，无草边坡安全系数较大，根系对边坡的稳定性加持较小。当边坡坡度增加时，无草边坡安全系数降低，根系固土的作用更加明显。

表5 边坡安全系数

图9a-9c分别给出了坡角为30°时，无草边坡和两种含草边坡模型的等效塑性应变云图。从图9中可以看出，无草边坡的滑动面较浅，含草边坡的滑动面较深。无草边坡坡脚变形较小，坡顶变形较大；2种含草边坡的坡脚处都有较大的隆起，坡顶有较大的竖向和水平位移。塑性区由坡角发展至坡顶，并从坡脚贯通至坡顶。由于根系的存在，塑性区主要分布在根系分布区以下部分，且更加分散，因此提高了含草边坡的稳定性。含草边坡均匀化模型与含草边坡Link模型的等效塑性应变云图分布规律基本一致，说明了本文所述均匀化方法的准确性。

图9 等效塑性应变分布

为达到必要的模拟精度，分离实体式模型含有155 513个三维实体单元，而本文提出的均匀化理论模型仅仅使用了6 873个单元，减少了95.58%的单元数。虽然两者在应力分布上可以取得类似精度，但是均匀化模型在计算速度上得到极大提高。当计算安全系数时，分离实体式模型由于异形单元的存在，极易造成收敛性失稳。而含草边坡Link模型在边坡工况不过分复杂的情况下可以解决此问题，但是依然使用了16 397个实体单元和250个Link单元，且建模过程复杂。本文仅列举了简单的直根系工况，并且建模时考虑到计算效率，在方向上采用了对称边界。如果按照自然界草根系多行多排方式直接建立模型，分离实体式和Link模型都将面临巨大困难，而本文方法将受干扰极小，优势将会进一步体现。因此，边坡根系含量越多，本文方法的优势将越明显。

4 结 论

本文提出了一种基于渐进均匀化理论的含草边坡本构关系的计算方法，通过对比研究和数值分析，验证了本文所述方法的准确性和高效性。本文选取黄土高原草本植物披碱草为对象，进一步研究了植物根系固土护坡的能力，得到以下结论：

1）本文所述方法研究含大量根系的边坡具有明显的效率优势，极大地减少了模型计算工作量（减少了95.58%的单元数），为含根群的边坡稳定性研究开拓了新思路。

2）草根根系的存在，改善了坡体浅层土体的应力状态，提高了边坡安全系数，保持了边坡的水土稳定性。

3）当边坡坡度较小时（坡角为30°），无草边坡安全系数较大（=4.28），根系对边坡的稳定性加持较小（平均仅提高2.92%）；当坡角增大到45°时，无草边坡安全系数下降为2.90，而含草边坡平均安全系数提升了13.45%，根系固土效果更加显著。

本文基于理想的竖向根系，下一步研究工作将以植物单株为“单胞”，考虑单株根系分布形态，形成真实三维“单胞”模型，提高计算精度和方法通用性。

[1] Moresi F V, Maesano M, Matteucci G, et al. Root biomechanical traits in a montane Mediterranean forest watershed: Variations with species diversity and soil depth[J]. Forests, 2019, 10(341). doi:10.3390/f10040341.

[2] 吴美苏，周成，王林，等. 根系和裂隙对土体水力和力学特性影响数值模拟[J]. 岩土力学，2019，40(S1)：519-526，534. Wu Meisu, Zhou Cheng, Wang Lin, et al. Numerical simulation of the influence of roots and fissures on hydraulic and mechanical characteristics of the soil[J]. Rock and Soil Mechanics, 2019, 40(S1): 519-526,534. (in Chinese with English abstract)

[3] 郑明新，黄钢，彭晶. 不同生长期多花木兰根系抗拉拔特性及其根系边坡的稳定性[J]. 农业工程学报，2018，34(20)：175-182.Zheng Mingxin, Huang Gang, Peng Jing. Tensile-pullout properties of roots ofin different growth stages and stability of slope with its root[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(20): 175-182. (in Chinese with English abstract)

[4] 及金楠，田佳，瞿文斌. 基于连续断裂过程的根系黏聚力Wu氏模型修正系数的确定[J]. 林业科学，2017，53(11)：170-178.Ji Jinnan, Tian Jia, Qu Wenbin. Determination of correction coefficients of Wu's model of root cohesion based on successive fracture process[J]. Scientia Silvae Sinicae, 2017, 53(11): 170-178. (in Chinese with English abstract)

[5] Fan C C, Tsai M H. Spatial distribution of plant root forces in root-permeated soils subject to shear[J]. Soil and Tillage Research, 2016, 156: 1-15.

[6] Leung F T, Yan W, Hau B C, et al. Root systems of native shrubs and trees in Hong Kong and their effects on enhancing slope stability[J]. Catena, 2015, 125: 102-110.

[7] Mao Z, Bourrier F, Stokes A, et al. Three-dimensional modelling of slope stability in heterogeneous montane forest ecosystems[J]. Ecological Modelling, 2014, 273: 11-22.

[8] 余芹芹，胡夏嵩，李国荣，等. 寒旱环境灌木植物根–土复合体强度模型试验研究[J]. 岩石力学与工程学报，2013，32(5)：1020-1031.Yu Qinqin, Hu Xiasong, Li Guorong, et al. Research on strength model test of shrub root-soil composite system in cold and arid environments[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(5): 1020-1031. (in Chinese with English abstract)

[9] 栗岳洲，付江涛，余冬梅，等. 寒旱环境盐生植物根系固土护坡力学效应及其最优含根量探讨[J]. 岩石力学与工程学报，2015，34(7)：1370-1383.Li Yuezhou, Fu Jiangtao, Yu Dongmei, et al.Mechanical effects of halophytes roots and optimal root content for slope protection in cold and arid environment[J].Chinese Journal of Rock Mechanics and Engineering, 2015, 34(7): 1370-1383. (in Chinese with English abstract)

[10] Zhang C B, Chen L H, Jiang J. Why fine tree roots are stronger than thicker roots: The role of cellulose and lignin in relation to slope stability[J]. Geomorphology, 2014, 206: 196-202.

[11] 陈丽华，及金楠，冀晓东，等. 林木根系基本力学性质[M]. 北京：科学出版社，2012.

[12] 肖本林，罗寿龙，陈军，等. 根系生态护坡的有限元分析[J]. 岩土力学，2011，32(6)：1881-1885.Xiao Benlin, Luo Shoulong, Chen Jun, et al. Finite element analysis of eco-protection slope through roots[J]. Rock and Soil Mechanics, 2011, 32(6): 1881-1885. (in Chinese with English abstract)

[13] 盖小刚，陈丽华，蒋坤云，等. 4种乔木根系不同埋根方式根-土复合体的抗剪特性[J]. 林业科学，2014，50(9)：105-111.Gai Xiaogang, Chen Lihua, Jiang Kunyun, et al. Shear characteristic research on root-soil composite in four kinds of roots of trees and different root buried ways[J]. Scientia Silvae Sinicae, 2014, 50(9): 105-111. (in Chinese with English abstract)

[14] 朱锦奇，王云琦，王玉杰，等. 基于试验与模型的根系增强抗剪强度分析[J]. 岩土力学，2014，35(2)：449-458.Zhu Jinqi, Wang Yunqi, Wang Yujie, et al. Analysis of root system enhancing shear strength based on experiment and model[J]. Rock and Soil Mechanics, 2014, 35(2): 449-458. (in Chinese with English abstract)

[15] 田佳，曹兵，及金楠，等. 花棒根-土复合体直剪试验的有限元数值模拟与验证[J]. 农业工程学报，2015，31(16)：152-158.Tian Jia, Cao Bing, Ji Jinnan, et al. Numerical simulation and validation test of direct shear test for root-soil composite ofusing finite element method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(16): 152-158. (in Chinese with English abstract)

[16] 孔纲强，文磊，刘汉龙，等. 植物根系分布形态及含根复合土强度特性试验[J]. 岩土力学，2019，40(10)：1-8.Kong Gangqiang, Wen Lei, Liu Hanlong, et al. Strength properties on root compound soil and morphological observation of plant root[J]. Rock and Soil Mechanics, 2019, 40(10): 1-8. (in Chinese with English abstract)

[17] 朱锦奇，王云琦，王玉杰，等. 基于植物生长过程的根系固土机制及Wu模型参数优化[J]. 林业科学，2018，54(4)： 49-57.Zhu Jinqi, Wang Yunqi, Wang Yujie, et al. Analyses on root reinforcement mechanism based on plant growth process and parameters optimization of Wu model[J]. Scientia Silvae Sinicae, 2018, 54(4): 49-57. (in Chinese with English abstract)

[18] 杨璞. 根土复合体极限载荷的数值计算方法和实验研究[D]. 北京：清华大学，2008.Yang Pu. Numerical and Experimental Research of Limit Loads of Soil-root Composites[D]. Beijing: Tsinghua University, 2008.

[19] Ji J, Mao Z, Qu W, et al. Energy-based fibre bundle model algorithms to predict soil reinforcement by roots[J]. Plant and Soil, 2020, 446: 307-329.

[20] Zhu J, Wang Y, Wang Y, et al. How does root biodegradation after plant felling change root reinforcement to soil?[J]. Plant and Soil, 2020, 446: 211-227.

[21] Hassani B, Hinton E. A review of homogenization and topology optimization I: Homogenization theory for media with periodic structure[J]. Computers & Structures, 1998, 69(6): 707-717.

[22] 董其伍，王哲，刘敏珊. 渐进均匀化理论研究复合材料有效力学性能[J]. 材料科学与工程学报，2008，26(1)：72-75.Dong Qiwu, Wang Zhe, Liu Minshan. Asymptotic homogenization theory for composites: Effective mechanics properties[J]. Journal of Materials Science & Engineering, 2008, 26(1): 72-75. (in Chinese with English abstract)

[23] Hassani B, Hinton E. A review of homogenization and topology opimization II-analytical and numerical solution of homogenization equations[J]. Computers and Structures, 1998, 69: 719-738.

[24] 陈惠发，萨里普 AF. 弹性与塑性力学[M]. 余天庆，王勋文，刘再华译. 北京：中国建筑工业出版社，2004.

[25] Omine K, Ohno S, Deformation analysis of composite ground by homogenization method[C] // Proceedings of the International Conference on Soil Mechanics and Foundation Engineering-International Society for Soil Mechanics and Foundation Engineering, AA BALKEMA, 1997, 719-722.

[26] 扈萍，宋修广，吴登高. 高速公路边坡植草护坡的根固效应试验研究[J]. 岩土力学，2008，29(2)：442-444.Hu Ping, Song Xiuguang, Wu Denggao. Experimental research on reinforcement mechanism of expressway slope protection with greensward[J]. Rock and Soil Mechanics, 2008, 29(2): 442-444. (in Chinese with English abstract)

[27] Borovkov A I, Sabadash V O, Finite element multiscale homogenization and sequential heterogenization of composite structures[C] // ANSYS conference, Pittsburgh. USA, 2002.

[28] 及金楠，张志强，郭军庭，等. 黄土高原刺槐和侧柏根系固坡的有限元数值模拟[J]. 农业工程学报，2014，30(19)：146-154.Ji Jinnan, Zhang Zhiqiang, Guo Junting, et al. Finite element numerical simulation of Black Locust () and Arborvitae () roots on slope stability on Loess Plateau of China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(19): 146-154. (in Chinese with English abstract)

[29] 江胜华，汪时机，李伟清，等. 基于位移变化率和强度折减有限元的边坡失稳判定方法[J]. 农业工程学报，2017，33(15)：155-161.Jiang Shenghua, Wang Shiji, Li Weiqing, et al. Slope instability evaluation method using finite element method of strength reduction and displacement rate[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(15): 155-161. (in Chinese with English abstract)

[30] 涂义亮，刘新荣，钟祖良，等. 三类边坡失稳判据的统一性[J]. 岩土力学，2018，39(1)：173-180，190.Tu Yiliang, Liu Xinrong, Zhong Zuliang, et al. The unity of three types of slope failure criteria[J]. Rock and Soil Mechanics, 2018, 39(1):173-180,190. (in Chinese with English abstract)

[31] 侍倩. 植被对斜坡土体土力学参数影响的试验研究[J]. 岩土力学，2005，26(12)：2027-2030.Shi Qian. Test research on influence of vegetation on mechanical parameters of soils[J]. Rock and Soil Mechanics, 2005, 26(12): 2027-2030. (in Chinese with English abstract)

[32] 田佳，及金楠，钟琦，等. 贺兰山云杉林根土复合体提高边坡稳定性分析[J]. 农业工程学报，2017，33(20)：144-152. Tian Jia, Ji Jinnan, Zhong Qi, et al. Analysis on improvement of slope stability in root-soil composite of Picea crassifolia forest in Helan Mountain[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(20): 144-152.

Numerical simulation of root reinforcement for herbs in Loess Plateau based on asymptotic homogenization theory

Huang Jiankun, Wang Xuelin, Ji Jinnan, Chen Lihua, Zhang Zhiwei

(1.,,,042200,; 2.,,,100083,; 3.,,100083,)

Loess Plateau is the most severe soil and water loss area in the world, as well subjected to shallow-landslide disaster in China. Currently, grass system has been widely distributed in the Loess Plateau, while the root system shows obviously periodic characteristics in spatial distribution. In order to accurately evaluate the influence of root groups on root reinforcement, and the coupling relationship between roots and soils, a constitutive relation of the “root-soil” composite was constructed via the “unit cell” of the root and soil based on the asymptotic homogenization theory. The “root-soil” composite in the nature similar to the reinforced concrete, can be regarded as a “new composite material”, with a strong coupling relationship between roots and soil. In the deduction, some numerical methods including the perturbation method, periodic boundary conditions, subsection integral method and divergence theorem, were used to derive the expression of the equivalent stiffness matrix, and the equivalent density of the root-soil composite. The detailed solution to these functions was also given based on the finite element method. The two-dimensional elastic parameters of the root-soil composites were calculated by MATLAB program. The three-dimensional equivalent elastic parameters of the root-soil composite were eventually obtained, where the two-dimensional plane strain problem can be extended to three-dimensional one by additional equations. The calculation accuracy and stability of the present method are better than those of the simplified method, particularly on calculating the three-dimensional constitutive relationship. There was a certain deviation (up to 29.4%) in the calculation of equivalent Poisson's ratioυ, whereas, the calculation errors of other equivalent parameters are less than 7.1%. To illustrate the accuracy and efficiency of the homogenization method in root reinforced slopes, the influence of the’s roots on the stress and strain field of the slope was analyzed based on the finite element software ANSYS. Three types of numerical models were constructed, including the slope model without grass, root reinforced slope with separated root elements, and root reinforced slope based on the present homogenization theory. The slope safety factor was calculated using the strength reduction method considering the penetrated equivalent plastic strain zone or not. The results show that: (1) The asymptotic homogenization theory can accurately construct the constitutive relation of “root-soil” composites, while reduce the calculation work (the element number reduced up to 95.58%). The error ofυhas little effect on the stress distribution of the simulated slopes with homogenized materials. (2)root system can modify the stress field of the shallow slope, indicating more uniform of the shear stress in the root distribution zone. Therefore, the slope stability can be improved. (3) If the slope angle is small (30°), the safety factor of the slope without grass is large (= 4.28). The root system has a small effect on the slope stability (an average increase is only 2.92%). If the slope angle increases up to 45°, the safety factor of the slope without grass reduces to= 2.90, while the averaged safety factor of slopes with grass increases by 13.45%, indicating the dominated reinforcement effect of root system on slopes. These findings can open up a new way to set “root-soil” periodic composites for root reinforced slopes.

slope stability; roots; reinforcement; homogenization theory; finite element method

黄建坤，王学林，及金楠，等. 基于渐近均匀化理论的黄土高原草本植物固土效果模拟[J]. 农业工程学报，2020，36(9)：168-176.doi：10.11975/j.issn.1002-6819.2020.09.019 http://www.tcsae.org

Huang Jiankun, Wang Xuelin, Ji Jinnan, et al. Numerical simulation of root reinforcement for herbs in Loess Plateau based on asymptotic homogenization theory[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(9): 168-176. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.09.019 http://www.tcsae.org

2019-11-6

2020-03-09

国家自然科学基金（31700637）；中央高校基本科研业务费专项资金（2015ZCQ-SB-01）

黄建坤，副教授，博士，主要从事防灾减灾、水土保持等方面的研究工作。Email：jiankunhuang@bjfu.edu.cn

10.11975/j.issn.1002-6819.2020.09.019

TU 441

A

1002-6819(2020)-09-0168-09