加筋土技术是一种将土体和土工合成材料复合，用以改善土工结构力学性能的土工加固技术。该技术通过在土体中成层埋设高模量的土工格栅，利用拉筋与土之间的摩擦作用改善土体的变形条件并提高土体的工程性能，从而增强挡墙的稳定性。这种方法不仅能够使填方挡墙坡度增大，也能够减少土方用量，节约工程占地，并提高经济效益[1-2]。
尽管加筋土结构已应用数十年，但其理论研究滞后于实际工程经验。对于自锚式加筋土挡墙的设计，通常采用极限平衡法计算挡墙的稳定性[3]。然而，这种方法未考虑筋材与土体之间的接触以及土体内部的应力应变，无法准确评估挡墙的稳定性并获得各参数结果。为了更合理地分析加筋土挡墙力学效应和稳定性，国内外许多学者采用了数值模拟的方法对返包式加筋土挡墙进行研究。王多垠等[4]对E-B模型中7个参数进行调整，分析不同参数对筋带拉力和挡墙位移的影响程度；王宁等[5]采用有限元法（finite element method，FEM）研究了筋材长度、刚度和竖向间距对加筋土挡墙墙背土压力的影响；陈晓飞等[6]采用FEM，考虑了分级填土过程，对面板位移、格栅拉力的分布及抗滑安全系数进行了分析；张磊等[7]建立了加筋与不加筋的三维模型，从三维特性方面分析加筋挡墙的变形和应力状态；任鑫[8]采用有限差分软件对返包式加筋土挡墙的水平位移、剪应力分布规律及稳定性进行分析；曹磊等[9]应用有限元软件，分析了加筋层数、筋材位置、筋材长度、筋材模量对挡墙水平变形的影响；李永亮等[10]通过不同的网格划分技术探究不同的网格形式对挡墙稳定性的影响；李元松等[11]综述了边坡稳定性评价方法的研究现状并指出存在的问题，分析与归纳各类评价方法的优缺点及其适用条件。Amirhosseini等[12]使用FEM开发了一系列三维模型，以研究土工格栅在受限和非受限条件下的响应，提出了一种配置复杂的地理网格建模方法。
数值模拟方法在加筋土工程领域发挥着越来越重要的作用。在建立有限元模型时，选择合适的网格划分形式是保证模型运算结果的准确性的关键。因此，在自锚式加筋土挡墙工程中，依托有限元软件ABAQUS建立二维模型，并对加筋土挡墙模型的网格划分进行研究，分析不同网格尺寸对模型计算精度的影响，以确保建立的模型运算结果更加准确。
1 有限元计算模型的建立
为建立有限元模型，作出如下假设：（1）加筋土挡墙的回填土处于平面应变的应力状态；（2）墙后填土表面水平，土工格栅水平铺设；（3）在正常使用荷载作用下（无筋材蠕变），每层筋材的拉力沿水平方向连续分布；（4）挡墙仅受自重作用（不考虑顶部荷载），回填土的变形主要为塑性变形，弹性变形很小，可以忽略[13]。
1.1 几何模型与边界条件
本文以湖北省宜昌市某加筋土挡墙为研究对象，根据地勘资料，该加筋土挡墙的地基岩土体以中风化灰岩为主。根据设计资料和地勘资料，对比参考同类工程及相关文献数据选取模型材料计算参数，如表1所示。
表1 材料物理力学参数
Table 1 Physical and mechanical parameters of
each material
[材料 重度[γ] /
（kN/m3） 弹性模
量E / MPa 泊松比μ 黏聚力c / kPa 内摩擦
角φ / （°） 中风化灰岩 24 15 000 0.21 300 35 回填土 18 50 0.33 25 35 土工格栅 — 26 000 0.33 — — ]
挡墙采取台阶式断面，分3级放坡，每级高8 m，总高24 m，坡长为50 m，每级挡墙坡度为1∶0.5，马道长度1 m，返包段2 m。地基长度75 m，高度24 m。土工格栅在3级挡墙的长度依次为32、28和16 m。挡墙几何模型如图1所示。
将几何模型输入有限元软件中并施加边界条件和荷载，如图2所示。本研究只考虑在土体施加重力荷载的情况。地基底部设置固定约束，挡墙背部以及地基土的左右两侧设置水平约束，其他边界设置为自由边界。
1.2 本构模型
岩土的应力应变关系复杂，研究人员提出了多种本构模型，在岩土工程有限元分析中，选用适合的本构模型才能建立出更精确的有限元模型。常用的本构模型有Mohr-Coulomb、Duncan-Chang、剑桥模型和Drucker-Prager等[14]。本研究在进行数值分析时，对回填土和地基土采用Drucker-Prager理想弹塑性本构模型进行模拟。Drucker-Prager本构模型可以有效避免Mohr-Coulomb本构模型在棱角处屈服函数的不可微或奇异的不足[15]，同时考虑了中间主应力的影响，方便塑性增量的计算。
Drucker-Prager本构模型屈服准则表达式为
[FI1,J2=αI1+J2-k=0] （1）
[α=2sinφ33-sinφ] （2）
[k=6ccosφ33-sinφ] （3）
式（1）~式（3）中：[I1]为应力张量第一不变量，[J2]是应力偏张量第二不变量；[φ]为内摩擦角；[c]为黏聚力。同时，为简化计算和分析，采用了非相关联流动法则，采取极端条件即膨胀角等于0，避免材料膨胀或收缩。
土工格栅是一种高模量的抗拉材料，在土体中只承受拉力，可以近似认为其应力应变关系为一条直线。采用线弹性本构关系来模拟土工格栅，在计算时采用Truss单元（即桁架单元，也是杆单元）进行模拟。
1.3 筋带与土体接触
目前土工格栅与回填土接触问题数值模拟主要有如下几种方式：
（1）将加筋土中的土工格栅和回填土处理成不同的材料，分别用不同的本构关系进行有限元模拟计算，不考虑加筋土中的土工格栅和回填土接触面上的相互作用，将格栅和回填土分开处理。
（2）在土体和土工格栅之间设置接触单元，将加筋土体视为由土体、土工格栅及两者间的接触单元共同构成的整体，接触面的单元应力大小随土工格栅和回填土之间产生的相对位移而改变[16]。
（3）将筋材的作用等效为沿水平方向作用于土体的附加应力，并只用加筋土中的土体进行计算，在有限元计算中只出现土单元[17-18]。
本研究用ABAQUS对加筋土挡墙进行有限元建模分析时，在土体和土工格栅之间设置接触单元，通过ABAQUS中的“嵌入约束”（embedded）将回填土与土工格栅相连接实现筋土之间的位移协同。该方法本质上是一种位移协调约束，不考虑格栅与回填土之间的相对位移[19]。
2 有限元模型计算精度的影响因素
在解决实际工程问题时，网格划分极大程度地影响计算结果的精度。如何建立好的网格主要包含3个因素：合适的单元类型、良好的单元形状和适当的网格尺寸。在网格划分时，不同的单元形状、单元尺寸、单元类型等因素都会影响最终的计算结果，所以在建模过程中，网格质量直接关系到有限元分析计算能否顺利、快速地完成，也关系到能否得到高精度的分析结果。
本文研究的重点是如何在建立加筋土挡墙有限元模型时设置最适合的网格尺寸，从而在满足精度要求的同时提高工作效率。减小网格的尺寸，计算精度通常会提高，但是网格划分时间和模型计算时间也会成倍增长，甚至会因为计算误差累积导致有限元分析不收敛。所以在分析实际工程中，应选择合适的网格尺寸，提高计算精度的同时也避免产生过大的工作量。
3 有限元计算结果对比分析
对挡墙网格和土工格栅网格分别采用CPE4平面应变单元和T2D2桁架单元划分。控制不同的网格尺寸，对加筋土挡墙有限元模型计算结果进行分析，通过对挡墙安全系数、挡墙应力和筋带应力、挡墙水平位移等参数的研究，得到不同的网格尺寸对加筋土挡墙有限元模型计算结果的影响。
3.1 网格尺寸对挡墙安全系数的影响
首先选择网格尺寸为0.50 m（约等于坡高的1/50）进行试算，得出安全系数为1.61，而极限平衡法计算得到的安全系数为1.52，有限元计算结果误差较大。因此，考虑通过减小网格尺寸，即将网格尺寸从0.40 m（即坡高的1/60）逐渐减小至0.15 m（即坡高的1/160），重新进行有限元计算。不同网格尺寸的有限元模型如图3所示。
<G:\武汉工程大学\2025\第2期\段博文-3.tif>[（b）][（a）]
图3 不同网格尺寸下有限元模型：（a） 0.40 m，（b） 0.15 m
Fig. 3 Finite element model under different mesh sizes：
（a） 0.40 m，（b） 0.15 m
利用强度折减法计算的挡墙安全系数如图4所示。当网格尺寸较大时，减小网格尺寸可以增加计算精度，而计算时间未显著增长，表明挡墙安全系数随着网格尺寸的减小而不断减小，越来越接近理论值。然而，当网格尺寸达到0.25 m（约为坡高的1/100）时，划分的网格数成倍增长，导致计算时间大幅增加（增幅为117.1%），而计算精度的提升相对有限，0.15 m网格尺寸相比于0.25 m网格尺寸的相对误差仅减少了0.64%。减小网格尺寸并不能提高安全系数的计算精度，同时还会使软件计算速率降低。因此，在计算挡墙安全系数时，需要选择合适的网格尺寸平衡计算结果的精度和计算时间，以提高工作效率。
3.2 网格尺寸对挡墙最大主应力和xy方向剪应力的影响
如图5所示，随着网格尺寸逐渐减小，挡墙主应力的变化并不明显，在205 kPa上下浮动，网格尺寸达到一定值后，最大主应力值将会接近真实值。这说明网格尺寸对挡墙主应力的影响较低，即使使用尺寸较大的网格，同样能得到较准确的主应力值。因此，在建立有限元模型计算挡墙应力时，为了提高效率，可以优先选择使用网格尺寸较大的单元，避免不必要的网格细化，以减少计算时间和所需资源。同时，在选择使用网格尺寸较小的单元时，也需要注意单元的形状、大小等因素，以确保计算结果的准确性。
如图5所示，当减小网格尺寸时，受限于FEM的局限性，挡墙剪应力可能出现应力集中现象。当网格尺寸细化至0.20 m（即坡高的1/120）时，挡墙剪应力出现应力奇异现象，剪应力大幅增长，增幅超过12.5 kPa。但随着网格尺寸的减小，挡墙剪应力持续增长，并未呈现明显的收敛趋势。这是因为更小的网格能够捕捉到细节上的应力变化，特别是在应力高度集中的区域。剪应力最大值多出现于挡墙台阶处，由于台阶处存在尖角，当网格过于细化时，在台阶边界处引起了数值的不稳定性，导致在局部区域产生不现实的高应力值。因此，在建立挡墙有限元模型时，必须做出综合考虑，以满足精度的要求。在易发生应力集中的区域选择适当的网格尺寸，即可有效避免应力集中现象的发生，使挡墙剪应力不再异常增长。
3.3 网格尺寸对格栅应力的影响
土工格栅应力图如图6所示，初步观察结果显示，格栅最大应力出现在第二级挡墙第一层中后部，最大值为235.9 kPa。
<G:\武汉工程大学\2025\第2期\段博文-6.tif>[+2.359e +05
+2.162e+05
+1.966e+05
+1.769e+05
+1.573e+05
+1.376e+05
+1.180e+05
+9.830e+04
+7.864e+04
+5.899e+04
+3.933e+04
+1.967e+04
+1.847e+01
][Mises应力/Pa]
图6 土工格栅应力图
Fig. 6 Geogrid stress diagram
如图7所示，网格尺寸对格栅应力分布具有显著影响：当格栅的网格尺寸较大时，网格尺寸对格栅应力的改变影响较小。当网格尺寸减小到0.25 m时，格栅应力突然增大，这是因为随着网格尺寸的减小，有限元计算能够捕捉到细节上的应力变化，提高了计算结果的精度。在较大网格模型中，这些高应力区域可能被平均化或不被充分解析，从而导致应力低估。尽管在较大尺寸的网格划分下可呈现较均匀的应力，但在实际工程分析中，对关键受力区域进行网格细化至关重要，以更真实地模拟格栅的力学行为并提升计算效率，因此建议在计算中采用尺寸较小的网格。综上，对土工格栅的网格划分应当注重对受力较大区域的精细计算，并且要结合相应的工程实际问题综合考虑。
<G:\武汉工程大学\2025\第2期\段博文-7.tif>[0.15 0.20 0.25 0.30 0.35 0.40
网格尺寸 / m][237.0
236.5
236.0
235.5235.0
234.5
234.0233.5
][格栅应力 / kPa][15.0
14.8
14.6
14.4
14.2
14.0][边坡位移 / mn][格栅应力
边坡位移]
图7 不同网格尺寸下土工格栅应力与边坡水平位移点线图
Fig. 7 Plot of geogrid stress versus horizontal slope
displacement points for different mesh sizes
3.4 网格尺寸对水平位移的影响
不同网格尺寸下挡墙水平位移如图8所示。若网格尺寸较大，则减小网格尺寸会导致位移明显改变。但是随着网格尺寸的不断减小，水平位移的增量逐渐减小，且网格尺寸由0.25 m减小至0.15 m，水平位移变化仅为0.08 mm（图7）。通过逐步减小网格尺寸，有限元计算结果对挡墙水平位移的计算会越来越准确和稳定。但当网格尺寸细化到一定程度时，其对水平位移的影响逐渐减少。因此，在分析挡墙水平位移过程中，可以选择适当的网格尺寸进行有限元计算，以提高计算效率。当数据精度提升不足1%，可以停止减小网格尺寸，以提高计算效率。
4 结 论
通过对加筋土挡墙建立有限元分析，可以得出以下结论：
（1）在建立挡墙有限元模型时，应考虑网格的经济性，选择合理网格尺寸，在满足计算精度的同时尽量缩短计算时间，可以选择网格尺寸为坡高的1/50左右进行初步试算，将试算结果与监测结果进行比对，再改变网格尺寸。在对挡墙安全系数进行计算时，减小网格尺寸可以有效提高计算精度，但同时也会使计算时间成倍增长。因此在选择网格尺寸时需谨慎权衡，可以逐级减小，当耗费的时间超过精度提升所能带来的收益时（如精度提升不足1%时），就无需再减小网格尺寸。
（2）网格尺寸对挡墙水平位移和挡墙主应力的影响较小，主应力在网格尺寸改变下变化幅度不超过1.4%，水平位移变化值在1 mm内。对挡墙剪应力的影响较大，整体影响幅度超过15%，这是因为在减小网格尺寸的过程中，会出现应力集中效应，使剪应力不断增大且无法收敛。因此，在水平位移和挡墙应力满足计算精度的情况下，无需进一步细化网格以避免时间浪费。应在确保计算精度与控制计算成本之间寻求平衡，从而确定合理的网格尺寸。
（3）网格尺寸的改变对于格栅应力的影响较为显著。在研究格栅受力时，应根据工程需求，选择较小网格进行运算，以便达到研究格栅受力及位移的目的。