地基加固模拟与仿真

Foundation Reinforcement Simulation and Modeling

摘要: 利用FLAC-3D软件，对深夯实技术在地基强化中的作用进行了详细的有限元数值模拟研究，该软件的基础概念和操作特性在此得到简要概述。尽管有关强夯强化地基的深度影响理论在众多文献中有所讨论，但由于现场条件的多样性，至今未建立一致的控制标准或共识。因此，本部分依据实际工况，对含有砾石的粘性土层地基进行强夯法的数值模拟，重点探讨了在这种特定地质条件下强夯强化的有效影响范围。

Abstract: Using FLAC-3D software, a detailed finite element numerical simulation study on the role of deep compaction technology in foundation strengthening is carried out. The basic concepts and operational characteristics of the software are briefly summarized here. Although the depth influence theory of dynamic compaction has been discussed in many literatures, there is no consistent control standard or consensus due to the diversity of field conditions. Therefore, according to the actual working conditions, the numerical simulation of dynamic compaction method for the cohesive soil foundation with gravel is carried out in this part, and the effective influence range of dynamic compaction under this specific geological condition is discussed.

1. 引言

数值分析领域涵盖了一系列方法，其中包括经典的欧拉法和创新的拉格朗日法。拉格朗日法独树一帜，因为它是一种显式计算技术，这与常规的隐式计算方法有所区别，特别是在与边界元法和传统的有限元法比较时[1]。在实践中，FLAC-3D是一个典型的数值模拟工具，本研究选择了一个特定的公路地基处理实验区作为研究对象，利用该软件结合实地试验数据，对含砾粘土的强夯过程进行了深入的数值模拟分析。

近年来，国内强夯技术发展迅速，应用范围更加广泛，主要集中在饱和软粘土地基的强夯技术研究上。一般认为强夯法特别适合于粗颗粒非饱和土、含水量不大的杂填土与湿陷性黄土，低饱和度粘性土与粉土可采用。对于饱和粘性土，可采用铺设砂垫层和打设排水板或袋装砂井的方法来加速孔压的消散。铺设砂垫层有利于强夯能量的有效吸收，均匀传播，使波的传播范围扩大，从而在较大土体范围内形成毛细管裂隙，促进内部水分的排出。但是对于我国特别是沿海地区广泛分布的淤泥、淤泥质土等软粘土地基，强夯法的加固效果和适用性目前尚存在争议。

2. FLAC3D软件特征

FLAC(快速拉格朗日有限差分分析程序)源于流体力学，诞生于1970年代中期的美国。自此，它已渗透到交通、水力、地质、环境科学以及土木建筑等众多工程领域，特别是在岩土工程中发挥了重要作用。作为一种显式有限差分法，FLAC无需构建整体刚度矩阵，因此节省了内存，特别适合处理大变形、非线性和不稳定问题。在应对大形变、非线性行为、结构不稳定性及大规模塌陷这类极具挑战的问题时，其性能远超同类有限元软件[2]。接下来，我将从应用特性和计算特性两个角度探讨其特点。

计算特征

本FLAC-3D这款有限差分程序在数值模拟领域展现出独特的优点，主要体现在：

(1) 它运用“混合离散方法”，能更精确且合乎逻辑地描绘物料的塑性破损和流动行为。

(2) 该软件倾向于使用显式差分技术来解决微分方程，这种方法的独特之处在于无需构建繁琐的刚度矩阵，显著减少了内存占用；它具有通用性，无论是线性还是非线性材料响应，都可通过应变增量的简单处理，轻松推算出应力增量和作用力，从而实时跟踪系统动态演变过程。针对大规模变形问题，它巧妙地采用分步策略简化求解，有效降低了计算的复杂度，极大地提高了效率[3]。

FLAC的网格构造原理与众不同，不同于常见的有限元网格设计。在FLAC框架内，网格和节点的布局依据一个独特的(I, J)坐标体系构建，其中I象征着水平的X轴方向，而J则指向垂直的Y轴。如图1中的数学网格示例清晰标注了(I, J)坐标，每个黑色标记点都对应着一个精确的(I, J)坐标对。实际上，如阴影区域所示的网格结构，都是围绕这个(I, J)坐标系展开的。

Figure 1. Differential grid and node diagram

图1. 差分网格和节点示意图

3. FLAC3D模型建立及参数选取

3.1. 本构模型

配置模型的核心在于揭示土壤应力与形变之间的宏观行为法则，它在构造模拟平台时，倾向于通过抽象的度量，如压实度比、各向异性参数和孔隙结构，来概括颗粒材料的基本力学特性，而略去了微粒间精细的相互作用描述。在实际模拟运算中，依赖于高级数学工具和算法，以这些状态参数为基础，构建出土壤单元的应力–应变响应关系。

目前的配置理论框架广泛包括非线性弹性、弹塑性和粘弹性三个核心领域[4]。在常规动力学应力的数值解析中，常见的配置模型有弹性模型和莫尔–库仑理论，本研究的数值模拟探讨则特地选择了莫尔–库仑模型作为工作基础。

3.2. 冲击荷载形式

针对强夯施工过程中锤击与地基交互的冲击力特性，全球范围内的学者们从理论上和实际操作层面进行了深入透彻的研究。高有斌采用动力学的应力–速度关系原理，巧妙地将应力变化的时域描绘转化为速度响应的时域理解，他的工作着重于理论转化[5]。王桂尧则通过边界元技术，深入解析了动态应力与接触时间之间的复杂关系，形成了一套独特的关联图谱[6]。周世良的研究聚焦于高填土地层的强化夯击处理，他运用非线性有限元模拟手段，揭示了夯坑结构、土壤应力分布以及塑性区域演化的过程，为揭示这一工程实践的独特现象提供了重要参考[7] [8]。郭见扬的研究强调了强夯能量主要沿垂直方向传递，从而引发了土壤显著的垂直位移效应[9]。这些丰富的研究成果为冲击负荷作用机制的理解提供了坚实的科学基础。

3.3. 网格模型建立

在处理力学问题时，最常用的是弹性模型。这个模型运算效率高，可以揭示应力应变关系，并揭示应力集中区域，确定网格的精细度，只需要两个关键参数——体积模量和剪切模量G。基于已有的信息评估模型的有效性是很重要的，随着数据的积累，可能需要转向更复杂的摩尔库伦模型，它能更精确地模拟土壤的应力应变行为[2]。

在深度强化处理粘土层中嵌有砾石的区域时，利用了先进的FLAC3D仿真工具进行了理论探讨。研究的核心焦点是设定一个具体的力学模型：夯锤的直径设定为1.5米，其作用点精确地定位在土壤表层的中央。该模拟环境构建在一个三维空间内，规模设定为10米的长度、10米的宽度以及25米的高度。模型的边界条件极具现实性，顶部即地表被定义为自由运动边界，这意味着它不受任何限制。另一方面，模型底部和侧壁均设置了零位移约束，以模拟实际工程中的刚性支撑。为了实现精细的网格划分，采用了包含20,000个节点的网格结构，以确保计算结果的准确性和细致性。

3.4. 模型参数选取

地基强化过程采用强夯技术，其本质是通过重锤的自由下落对土壤施加瞬时冲击力。在FLAC3D的模拟场景中，自始至终关注的是从锤子与土壤接触的那一刻，每一级加载阶段维持0.071秒[10]。通过深入分析实验室中的强夯试验场地原土样本，成功提取并构建了一套详尽的地基土力学特性模型，以获取地基土的特性。这些特性包括初始弹性模量E₀、压缩模量E_s、剪切阻力G、体积响应K、土壤的固有密度$\rho$ 、泊松比$\nu$ 以及孔隙介质中的水压p。此外，通过整合先前的地质探测资料，精确获得了渗透率k和孔隙填充程度e₀等核心参数，具体数值见表1中的详尽数据记录。

Table 1. Subsoil model parameters

表1. 地基土模型参数

土体模型参数	土层深度(m)
	2	4	8	25
泊松比$\nu$	0.441	0.441	0.441	0.441
变形模量E0 (MPa)	2.511	3.414	4.202	5.106
体积模量Ev (MPa)	5.232	7.113	8.741	10.632
剪切模量G (MPa)	0.883	1.198	1.481	1.797
粘聚力C (kPa)	25.425	37.506	44.230	46.542
内摩擦角$\phi$ (˚)	15.129	17.440	21.222	23.533
密度$\rho$ (g/cm^3)	1.702	1.912	2.007	2.217
孔隙水压力p (kPa)	7.354	9.445	12.607	3.152
渗透系数k (10^−6 cm/s)	3.656	4.055	4.549	4.444
孔隙比e0	0.842	0.804	0.755	0.748

其中剪切模量G、体积模量K、变形模量E₀见下：

$G=\frac{E_0}{2\left(1+\mu \right)}$ (3.1)

$K=\frac{E_0}{3\left(1-2\mu \right)}$ (3.2)

$E_0=\left(1-\frac{2{\mu }^2}{1-\mu }\right)E_s$ (3.3)

4. 数值模拟计算及结果分析

4.1. 数值模拟计算

本数值模拟重点探讨了夯锤与土壤接触后的土壤应力应变状态，未涉及夯锤下落的动态过程。仅在夯锤与土壤接触瞬间赋予夯锤一个起始速度。通过计算夯锤的质量与其下落距离的乘积，可以确定夯击能量。接触时的初始速度可根据自由落体运动的法则来估算。当夯击能量等于2100 kN·m时，依据自由落体运动的规则，初速度为v = 13.56 m/s。

夯锤接触面作用力F及接触面作用应力σ：

$F_t=mv$ (4.1)

可以推算出夯锤与接触面之间的力F以及应力σ的表现。

当时间t等于0.071秒时，计算得出力F的峰值为294,560牛顿，而应力σ的集中度则达到了9,379千牛顿每平方米。通过精细操控落距参数，可以调控并释放不同程度的夯击能量。在FLAC3D的高级数值模拟环境中，对这些应力参数的微调引发了土壤基础层的不同响应效应。下表2列出了在不同能量级别下强夯操作的相关参数，这些参数是基于上述的计算方法得出的。

Table 2. Dynamic compaction parameters corresponding to different dynamic compaction actions

表2. 不同强夯作用下对应的强夯参数

强夯能级(kN·m)	接触面速度(m/s)	接触面作用应力(kPa)	接触面作用力(kN)
1150	11.1494	7606	26,408
1500	12.8507	8867	28,967
2100	14.5092	10,036	31,518
2800	17.334	13,170	37,727

4.2. 模拟结果分析

(1) 解析

利用先进的FLAC3D仿真平台，细致探讨了Q1至Q4四个实验场地在实施强夯加固过程中的动态行为。通过深入分析不同冲击能量条件下的土壤模型运动，一个显著的发现是，土壤的最大响应点位于夯锤核心底部，其位移模式呈现出独特的“梨形”分布特征，这与现场实操中对土体动态应变研究的结果相吻合[11]。

另一方面，Kuhlemye等人在其研究中强调了有限元建模网格设计的关键性，他们建议网格间距应小于对应波长的八分之一，这是确保模拟精度和物理现象准确再现的重要参数。

网格大小应满足$\Delta x\le 1/6{\lambda }_{\min }$ ；

网格最小尺寸为$\Delta x\le 1/4\cdot 4{\lambda }_{\min }$ 。

网格单元尺寸基本可按下式选取

$\Delta x\le \left(1/6\text{-}1/12\right)\cdot 4{\lambda }_T$ (4.2)

上式中：

$\Delta x$ ——网格尺寸；

${\lambda }_T$ ——所需考虑截止频率对应的波长。

根据项目工程的地质特性及其计算简化的考量，本文设定网格基础尺寸遵循公式(4.1)，取波长的十分之一作为基准，最终选定每个网格单元的大小为1.0米。特别强调的是，为了确保有限元软件在预设时间内能输出精确且可靠的分析结果，时间步长$\Delta t$ 的选择至关重要。时间步幅过大可能导致计算精度减损，而过小则会导致计算过程显著拉长。基于对众多学者关于时间步长研究的整合，刘波在其工作中提出了一种优化策略，具体体现在方程(4.3)中。

$\begin{array}{l}\Delta t=\Delta x/V,\Delta x\le y,\Delta x/V\\ \Delta t=\Delta y/V,\Delta y\le \Delta x,\Delta x/V\end{array}$ (4.3)

4.3. 荷载输入及参数

根据项目施工计划中的详细设计，所采用的夯击设备重量达到了惊人的3000千牛顿，考虑到地球重力加速度为10米每平方秒，其锤头直径为2.6米，而落下的距离则设定为20米。基于过往项目的实践经验及初步试夯成果，预估每次夯击对土壤产生的下沉效应L大约为0.6米。在进行强夯作业时，锤头对地基土壤产生影响的时间可以通过公式(3.4)精确估算得出。然而，这个计算过程并不单纯依赖于单一参数，而是综合了多个因素的动态交互作用。

$\Delta t=\sqrt{\frac{2}{gh}}L$ (4.4)

同时查阅文献可知，通过利用动量定理和重力势能做功[12]，推导出了接触面之间的应力公式：

$P_{\max }=\frac{2W}{A}\left(1+\frac{h}{L}\right)$ (4.5)

式中：

Δt——接触时间(s)；

g——重力加速度(m/s²)；

h——落距(m)；

W——夯锤质量(kN)；

L——单击沉降量(m)；

A——夯锤底面积(m²)。

将上述参数代入公式(4.4)、(4.5)求得在模型中接触时间为0.06 s，接触面应力为4.56 MPa，运用显示动力学模块将荷载输入。模型中各土层的物理力学参数如表3所示，结果如表4所示。

从表3和图表4的数据清晰显示，ABAQUS的有限元分析在模拟本项目中表现出色。模拟得出的单次沉降和累积沉降值与现场实测数据之间的差异微小，误差范围可接受，这验证了构建的模型与工程实际情况的高度吻合。此模型的建立为进一步探究强夯处理对加固效果的研究奠定了坚实的基础。

但随着打击能量增加，土体在夯击影响区内的水平和垂直位移出现无规则变化。

Table 3. Model soil parameters table

表3. 模型土体参数表

土类	厚度(m)	重度(kN/m^3)	压缩模量(MPa)	泊松比	摩擦角(˚)	膨胀角(0)	内黏聚力(kPa)
杂填土	8	18.3	6.2	0.34	14.3	5	20
粉质粘土	2	19.2	7.1	0.32	16.4		25
红粘土	12	18.8	8.7	0.3	17.3		35

Table 4. Table of vertical displacement of soil under dynamic compaction

表4. 强夯下土体竖向位移表

夯击次数	工程实际单击沉降量(cm)	工程实际累计沉降量(cm)	数值模拟单击沉降量(cm)	数值模拟累计沉降量(cm)
1	36	84	73	65
2	54	131	52	131
3	51	176	51	157
4	40	187	40	192
5	25	221	23	212
6	13	234	12	223

当能量从1150 kN提升至1500 kN时，位移区域逐渐扩大。当能量进一步提升时，土壤在横向的移动明显加剧，这种情况与实地重击施工中的观察结果相吻合——能量增加往往造成夯实区域附近产生裂隙和更深的沉降。这验证了模拟研究与实际工程现象的相符性。同时，这也暗示对于含砾黏土地基而言，并非重击能量越大效果越好，存在一个最为适宜的强夯能量等级。

(2) 压力分布图解析

观察图2，可以看出在各种打击能量影响下土壤模型的主要应力：随着深度增加，最大主应力逐渐降低，最高点位于锤子与土壤接触点之下。随着强夯能量级别的提高，夯坑侧壁开始出现张应力，并且随着能量的增强这一趋势持续。当打击能量从1150 kN提升到1500 kN，最大主应力的影响区域有所扩大，但进一步增加能量时，这一影响范围反而缩小。可能的成因可能涉及到土壤表面结构的深刻扰动，使得土壤与水分的相互作用显著加剧，并产生了显著的耦合效应，这在某种程度上分散了原本应被强化的能量，从而对强夯的强化效果产生了负面影响。土壤中最大主应力的非线性演变进一步突显了土壤强夯加固过程的动态不规则性。利用FLAC3D专业软件，详细记录并分析了不同冲击力条件下强夯所能达到的加固深度，详尽的数据已在表5中展示。

Figure 2. Stress cloud map of dynamic compaction Q2 zone

图2. 强夯Q2区应力云图

Table 5. FLAC3D numerical simulation of effective reinforcement depth

表5. FLAC3D数值模拟有效加固深

夯击能(kN·m)	1150	1500	2100	2800
有效加固深度(m)	5.0	5.5	4.3	4.1

5. 有效加固深度分析

在强夯地基处置领域，确定有效的加固区域是至关重要的设计考量，尽管目前对于这一范围的理解尚未达成共识。本文界定的强化作用域是指土壤在经历强夯处理后，其应力表现出明显转型的区域界限。通常来说，相对于实际的强化深度，有效影响范围的地理覆盖更为广泛[7]。

对于强化深度的有效度测定，工程实施中采用的策略多变[13]：(1) 依据强夯前后场地的常规贯入试验中锤击次数的变化作为指标；(2) 通过对土壤室内压缩模量的前后对比进行分析；(3) 依赖地质雷达探测技术的数据评估；(4) 通过考察土体液化程度的消除情况来判断；(5) 借助静力触探实验中的贯入阻力读数；(6) 结合Menard公式进行理论上的预测计算。值得注意的是，每项工程对强化深度的评估方法可能会有所差异，这主要取决于土壤的独特力学特性。

5.1. 有效加固深度计算

L.Menard其表达式为：

$H=\alpha \sqrt{\frac{Mh}{10}}$ (5.1)

模拟强夯加固深度H分别为5.0 m、5.5 m、4.3 m、4.1 m，计算结果如表6所示：

Table 6. Strengthening depth and correction factor corresponding to different dynamic compaction levels

表6. 不同强夯能级对应的加固深度及修正系数

夯击能(KN·m)	1150	1500	2100	2800
加固深度(m)	4.1~5.5m
L. Menard计算公式	$H=\alpha \sqrt{\frac{Mh}{10}}$
修正系数	0.41	0.52	0.32	0.38

L. Menard提出的校正因子范围广泛，从0.2到0.95，而此次计算中获取的校正系数浮动在0.37至0.58的区间内，恰好落在Menard方法的有效适用范围内。这证实了运用此公式来验证数值模拟精确性的合理性。

5.2. 有效加固深度验证

针对四种不同的强夯能量级别，即1150千牛顿米、1500千牛顿米、2100千牛顿米及2800千牛顿米的数值模拟结果显示，预估的土壤强化深度依次为5.0米、5.4米、4.3米和4.1米。一项全面的集成研究探讨了这些高强度强夯处理前后对地基性能的影响，涵盖了标准贯入试验的数据对比、瑞利表面波测试的评估、孔隙水压力变化的监测，以及土壤物理特性的演变动态。研究发现，区域Q1 (1150 kN·m)、Q2 (1500 kN·m)、Q3 (2100 kN·m)和Q4 (2800 kN·m)的实际有效加固深度分别为4.8米、5.6米、4.6米和4.4米。由此可见，数值模拟的数据与实地测量的结果显示出良好的一致性。

5.3. 有效加固总结

(1) 运用FLAC3D建模和参数设定，结合室内土工试验和地质勘查数据，分析了在不同能量级别下，强夯处理对含砾粘土基础的动力响应和变形情况。对比分析结果显示，随着能量级别的提升，土壤的位移和应力响应呈现出非线性特征，过高的强夯能量可能不仅未能增强土壤的稳定性，反而加剧了表面土体的潜在损伤，促使土体与孔隙水的相互作用更为显著，反而消耗了部分能量，削弱了加固效果。此外，高强度的强夯操作还可能引发周边区域的裂缝扩散和地面下沉，对表层土的压实过程产生负面影响。因此，对于砾质粘土基础而言，理想的强夯能量并非一味追求高能，而是存在一个最优的能量控制区间。

(2) 数值模拟深入分析了不同强夯能量水平下，如何修正并适应砾质粘土基础的Menard有效加固深度预测模型，以实现更精确的工程应用。

(3) 精确的数值模拟结果显示，当强夯能量设定为1150千牛顿米、1500千牛顿米、2100千牛顿米和2800千牛顿米时，对应的有效加固深度分别为5.0米、5.5米、4.3米和4.1米，这些数据与实际强夯工程中的观察结果高度契合，证实了模拟方法的有效性。

6. 结论

运用FLAC3D软件，基于现场测试和实验室土工实验的数据，对砾石粘土质地基的强夯处理进行了深入的数值模拟研究，以评估不同强夯能量水平下的强化效果。在1150 kN·m、1500 kN·m、2100 kN·m和2800 kN·m的强夯能量下，模拟得出的有效加固深度分别为5.0 m、5.5 m、4.3 m和4.1 m，这些结果与实际的强夯加固深度相符。

然而，本文在使用FLAC3D模拟强夯地基加固时，对孔隙水与土壤之间的交互作用的探讨尚显不足，这可能导致模拟结果与真实情况存在一定的偏差。因此，建议未来的研究应更全面地考虑各岩土层中孔隙水和土体的耦合效应，以提高模拟的精确度并使结果更加贴近实际工况。

尽管在强夯技术的理论研究方面已取得一些进步，但现有的知识仍远远不足以为实际的工程应用提供充分的理论支持。目前，设计计算方法尚不成熟且不全面，主要依赖于经验积累和实验验证。

此外，对于强夯作业中观察到的一些具体现象，我们还缺乏精确的理论解释。鉴于此，未来对强夯技术的研究应关注以下几个关键点：

(1) 针对饱和软粘土，深入理解强夯强化的机制以及强化后土壤强度随时间演变的规律；

(2) 研究强夯处理前后的饱和软粘土在微观结构上的变化；

(3) 优化施工参数，通过对比不同土壤类型，探索单击夯击能量、夯击次数、间歇时间等与土壤物理力学性质之间的最佳匹配关系；

(4) 开发和完善强夯加固效果的监测和评估方法。