数值模拟报告(DOC)

数值模拟报告(DOC)
第一部分：数值模拟技术研究文献综述
浅析数值模拟技术
1．引言
近年来，随着我国大规模地进行“西部大开发”和“南水北调”等巨型工程，越来越多的岩土工程难题摆在我们面前，单纯依靠经验、解析法显然已不能有效指导工程问题的解决，迫切需要更强有力的分析手段来进行这些问题的研究和分析。

自R.W. Clough 上世纪60年代末首次将有限元引入某土石坝的稳定性分析以来，数值模拟技术在岩土工程领域取得了巨大的进步，并成功解决了许多重大工程问题。

特别是个人电脑的普及及计算性能的不断提高，使得分析人员在室内进行岩土工程数值模拟成为可能。

在这样的背景下，数值模拟特别是三维数值模拟技术逐渐成为当前中国岩土工程研究和设计的主流方法之一，也使得岩土工程数值模拟技术成为当今高校和科研院所岩土工程专业学生学习的一个热点。

采用大型通用软件对岩土工程进行数值模拟计算，在目前已成为项目科研、工程设计、风险评估等岩土类项目的必须，学习和掌握Ansys、FLAC3D、UDEC 等数值计算软件已成为学校、科研院所对工程从业人员的基本要求。

数值模拟方法主要有限元法、边界元法、加权余量法、半解析元法、刚体元法、非连续变形分析法、离散元法、无界元法和流形元法等，各种方法都有其对应的软件。

2.数值模拟的发展趋势
可以说, 继理论分析和科学试验之后, 数值模拟已成为科学技术发展的主要手段之一。

随着软件技术和计算机技术的发展, 目前国际上数值模拟软件发展呈现出以下一些趋势：
（1）. 由二维扩展为三维。

早期计算机的能力十分有限，受计算费用和计算机储存能力的限制，数值模拟程序大多是一维或二维的，只能计算垂直碰撞或球形爆炸等特定问题。

随着第三代、第四代计算机的出现, 才开始研制和发展更多的三维计算程序。

现在，计算程序一般都由二维扩展到了三维，如LS-DYNA2D和LS - DYNA3D、AUTODYN2D 和AUTO-DYN3D。

（2）.从单纯的结构力学计算发展到求解许多物理场问题。

数值模拟分析方法最早是从结构化矩阵分析发展而来，逐步推广到板、壳和实体等连续体固体力学分析，实践证明这是一种非常有效的数值模拟方法。

近年来数值模拟方法已发展到流体力学、温度场、电传导、磁场、渗流等求解计算，最近又发展到求解几个交叉学科的问题。

例如内爆炸时，空气冲击波使墙、板、柱产生变形，而墙、板、柱的变形又反过来影响到空气冲击波的传播，这就需要用固体力学和流体动力学的数值模拟结果交叉迭代求解。

（3）.由求解线性问题进展到分析非线性问题。

随着科学技术的发展，线性理论已经远远不能满足设计的要求。

诸如岩石、土壤、混凝土等，仅靠线性计算理论就不足以解决遇到的问题，只有采用非线性数值算法才能解决。

众所周知，非线性的数值模拟是很复杂的，它涉及到很多专门的数学问题和运算技巧，很难为一般工程技术人员所掌握。

为此，近年来国外一些公司花费了大量的人力和资金，开发了诸如LS- DYNA3D、ABAQUS和AU-TODYN等专长求解非
线性问题的有限元分析软件，并广泛应用于工程实践。

这些软件的共同特点是具有高效的非线性求解器以及丰富和实用的非线性材料库。

3．数值模拟的基本原理
一般而言，岩、土体处于三向受力状态，其破坏模式往往表现为压-剪破坏和拉伸破坏。

要分析和预测岩、土体在外力作用下的变形、破坏，就需要对其变形、破坏情况进行较为直观地再现。

岩土工程数值模拟正是从岩、土体的受力状态出发，来分析和预测岩、土体破坏情况的一种手段。

其基本原理是以典型试样的物理试验（室内试验或现场试验）获得的强度来表征整个地质体的岩、土体强度，以边界条件替代地质体周围所受的约束条件，借由本构关系表达岩、土体在外力作用下的应力-应变特性，最终了解、预测岩、土体变形破坏情况。

它具有鲜明的时代特征，以计算机为实现平台，是信息化时代的产物。

通过与其它方法（如人工智能、人工生命科学、随机模拟、模糊数学、灰色理论以及分形理论等）交叉共生、相互耦合嫁接，以获得更广阔的发展空间。

从广义上来说，岩、土体的室内试验和原位试验也是一种模拟手段，本文称之为物理模拟。

之所以如此称谓，是因为它们也是为较真实地近似再现岩、土体在其所赋存的环境中所处的受力状态所采用的一种手段。

从这个意义上来说，它与数值模拟的基本原理是相同的，因此，可以将数值模拟称为虚拟实验室模拟。

所不同的是，数值模拟除可以进行常规尺寸模型的模拟外，还可以进行宏观和细观两个层面尺寸模型的模拟，而其输入的参数则需通过物理模拟来提供。

因此，数值模拟是与物理模拟并行发展、相互补充和相互验证的试验系统。

相较于其它方法，数值模拟具有可重复和操作性强，费用低廉，不受模型尺寸控制，可视化程度高的优点，能有效延伸和扩展分析人员的认知范围，为分析人员洞悉岩、土体内部的破坏机理提供了强有力的可视化手段。

当然作为一种分析方法，它也有自身的缺点，主要是易受制于岩、土体结构的描述和模型概化的准确性及合理性；受制于岩、土体物理试验模拟结果的准确性；受制于岩、土体本构关系与实际岩、土体力学响应特性拟合程度的高低。

4.数值分析方法中存在的问题
到目前为止，研究计算工程的文章很多，但真正用于实际工程的数值分析方法(例如有限元法等)却较少。

部分原因在于有较多不成功应用的实例。

为什么会有这种情况，原因是多方面的，下面列出几条仅供参考：
(1)对岩土工程数值分析方法缺乏系统的知识和深入的理解，出现问题时不知道在什么情况下属于理论问题或数学模型问题；在什么情况下是属于计算方法问题或本构模型问题；在什么情况下是参数的确定问题或计算本身的问题等。

(2)各种本构模型固有的局限性。

具有多相性土的物理力学性质太复杂，难以准确地用数学模型和本构模型描述。

例如邓肯一张模型不能反映剪胀性，不能反映压缩与剪切的交叉影响；模型只能考虑硬化，不能反映软化；模型不能反映各向异性。

剑桥模型也仅能考虑硬化而不能反映软化，不能反映土的剪切膨胀和各向异性，不能用于超固结土等。

(3)现有的试验手段和设备不能提供适当、合理和精确的参数。

靠少数样本点所获得的参数难以准确地描述整个空间场地的物理力学性能；土的参数因土样扰动难以高质量的获取，其精度很差。

有些模型要求较多的参数，但这些参数用常规的试验手段和设备难以获取等。

岩土工程中如何应用精确的数学模型和本构模型是一个值得注意的问题。

在一般结构分析中，因材料的力学性质简单、均匀，不确定性较小，一般采用较精确的数学模型会得到较精确的分析结果。

但就土这种材料而言，因其不确定性非常大，其情况发生了很大的变化。

众所周知，场地土性及其参数勘察结果的精度和准确性是很差的，由此导致既使采用了很精确的数学模型，但因输入参数的精度不能与之相匹配，其计算结果同样会很差。

采用精确的数学模型还会给人造成一种错觉，让人觉得其计算结果也一定会更好、更可靠。

这样可能使人们忽略了精确的数学公式也照样会有出错的可能性。

只有当输入参数的质量和精度很高，并能与数学模型的精度相匹配时，才有可能得到较为准确的计算结果。

5.结语
20世纪60年代以后，由于电子计算机的飞速发展使岩土工程数值分析方法得到不断发展和完善，并用于岩土工程实践。

虽然在工程实际使用中数值分析方法存在一些问题，但只要认清问题的实质，并采取措施去解决它，相信随着岩土工程数值分析方法的不断发展及其工程经验的不断积累，在工程实践中将会得到越来越多的应用，它必将成为岩土工程分析中的有力工具。

参考文献
[1] 张森，言志信，段建. 边坡开挖数值模拟及其稳定性评价研究[J]. 西部探矿工程. (3).
[2] 汪军，刘海波. 边坡稳定性的有限元数值模拟建模[J]. 华北科技学院学报.
(0).
[3] 陈印东，刘叔灼. 基于强度折减法的边坡稳定性分析[J]. 科学技术与工程.
(0).
[4] 王浩. 类土质路堑高边坡典型失稳机制与加固工程对策的数值模拟研究[D]. 铁道部科学研究院, 2004.
[5] 张超，杨春和. 有限差分强度折减法求解边坡稳定性[J]. 土木工程与管理学报. (4).
[6] 郑颖人，赵尚毅，宋雅坤. 有限元强度折减法研究进展[J]. 后勤工程学院学报. (0).
[7] 邹宝祥，李明，唐伟华. 某大桥边坡稳定性FLAC3D数值模拟分析[J]. 山西建筑. (3).
[8] 郭辉. 山西晋城土质垂直高边坡稳定性计算及数值模拟研究[D]. 西安科技大学, 2011.
[9] 郭志柳，陈建东，吴鹏. 填土物理力学性质对路堤边坡稳定性影响的数值模拟[J]. 江西理工大学学报. (9).
第二部分：数值模拟技术FLAC 3D
上机报告 FLAC 3D
数值模拟上机题 计算模型分别如图1、2、3所示，边坡倾角
分别为30°、45°、60°，岩土体参数为：
密度ρ＝2500 kg/m 3, 弹性模量E ＝1×108 Pa ，泊松比μ＝0.3，
抗拉强度σt ＝0.8×106 Pa ，内聚力C ＝4.2×104 Pa ，摩擦角φ＝17°
试用FLAC 3D 软件建立单位厚度的计算模
型，并进行网格剖分，参数赋值，设定合理的边界条件，利用FLAC 3D 软件分别计算不同坡角情况下边坡的稳定性，并进行结果分析。

附 换算公式：
1 kN/m 3= 100 kg/m 3 剪切弹性模量：)
1(2μ+=E G =3.846e7 体积弹性模量：)
21(3μ-=E K =8.33e7 40
25.36
6040
10030°
图1 倾角为30°的边坡(单位：m) 计算命令流如下：
new
gen zone brick p0 0 0 0 p1 100 0 0 p2 0 1 0 p3 0 0 40 size 50 1 10
gen zone brick p0 40 0 40 p1 100 0 40 p2 40 1 40 p3 74.64 0 60 p4 100 1 40&
p5 74.64 1 60 p6 100 0 60 p7 100 1 60 size 30 1 10
fix x range x -0.1 0.1
fix x range x 99.9 100.1
fix y
fix z range z 0.1 -0.1
model elas
prop density 2500 bulk 3e9 shear 1e9
set gravity 0 0 -10
solve
ini xdisp 0 ydisp 0 zdisp 0
ini xvel 0 yvel 0 zvel 0
model mohr
prop density 2500 bulk 8.33e7 shear 3.846e7 c 42000 fric 17 ten 800000
solve fos file slope30.sav associated
计算结果如下：
图1-a，网格剖分图
图1-b，速度矢量图
图1-c，速度等值线图图1-d，位移等值线图
最终计算边坡的稳定性系数为：Fs＝1.49 分析：
30°边坡稳定性系数采用的是FLAC3D内置的强度折减法求解，稳定性系数 1.49>1，从稳定性系数系数可以判断该边坡处于安全状态。

坡面最大速度为0.25mm/s，随着深度的增加，竖
向应力逐渐增大。

坡肩处出现下沉，最大值达到
3.54m 。

100406040
40
45°
图2 倾角为45°的边坡(单位：m) 计算命令流如下：
new
gen zone bri p0 0 0 0 p1 100 0 0 p2 0 2 0 p3 0 0 40 size 50 1 10
gen zone bri p0 40 0 40 p1 100 0 40 p2 40 2 40 p3 60 0 60 p4 100 2 40 p5 60 2 60 p6 100 0 60& p7 100 1 60 size 30 1 10
model mohr
prop density 2500.0 bulk 8.3e7 shear 3.8e7 c 42000.0 tens 0.8E6 &
friction 17 dilation 20
fix x y z range z -0.1 0.1
fix x range x 99.9 100.1
fix x range x -0.1 0.1
fix y
set gravity = 10.0
plot add axe red
plot con dis
;定义循环终止条件
def calfos
ait1=0.02
k11=0.6
k12=2.0
loop while (k12-k11)>ait1
fs=(k12+k11)/2.0
refric=(atan((tan(17*pi/180))/fs))*180/pi recoh=42000/fs
;折减实现过程
command
ini sxx 0.0 syy 0.0 szz 0.0 sxy 0.0 sxz 0.0 syz 0.0
ini xvel 0.0 yvel 0.0 zvel 0.0
ini xdis 0.0 ydis 0.0 zdis 0.0
pro fric refric coh recoh
set mech ratio 1e-5
solve step 5000
print fs
end_command
aa=mech_ratio
if aa<1e-5 then
k11=fs
else
k12=fs
end_if
end_loop
end
calfos
SA VE results.sav
solve fos file slop45.sav associated 计算结果如下：
图2-a，网格剖分图
图2-b，速度矢量图
图2-c，速度等值线图图2-d，位移等值线图
最终计算边坡的稳定性系数为：Fs＝1.15625 分析：
45°边坡稳定性系数采用的是FLAC3D自编的强度折减法求解，稳定性系数 1.16>1，从稳定性系数系数可以判断该边坡处于安全状态。

坡面最大速度为0.64mm/s，随着深度的增加，竖向应力逐渐增大。

坡肩处出现下沉，最大值达到
5,34m 。

100
40
60
48.45
40
60°
图3 倾角为60°的边坡(单位：m) 计算命令流如下： new
gen zone brick p0 0 0 0 p1 100 0 0 p2 0 1 0 p3 0 0 40 size 50 1 10
gen zone brick p0 40 0 40 p1 100 0 40 p2 40 1 40 p3 51.55 0 60 p4 100 1 40 p5 51.55 1 60& p6 100 0 60 p7 100 1 60 size 30 1 10 attach face range z 39.9 40.1 fix x range x -0.1 0.1 fix x range x 99.9 100.1 fix y
fix z range z 0.1 -0.1
model elas
prop density 2500 bulk 3e9 shear 1e9
set gravity 0 0 -10
solve
ini xdisp 0 ydisp 0 zdisp 0
ini xvel 0 yvel 0 zvel 0
model mohr
prop density 2500 bulk 8.33e7 shear 3.846e7 c 42000 fric 17 ten 800000
solve fos file slope60.sav associated
计算结果如下：
图3-a，网格剖分图
图3-b，速度矢量图
图3-c，速度等值线图图3-d，位移等值线图
最终计算边坡的稳定性系数为：Fs＝0.94 分析：
60°边坡稳定性系数采用的是FLAC3D内置的强度折减法求解，稳定性系数0.94<1，从稳定性系数系数可以判断该边坡处于不安全状态。

坡面最大速度为0.576mm/s，随着深度的增加，竖向应力逐渐增大。

坡肩处出现下沉，最大值达到4.14m，故应采取措施以保证边坡安全。

开挖后坡面
原始坡面
345°
45°
100
40
60
48.45
40
60°
图4 边坡开挖算例分析
计算命令流如下： new
gen zone brick p0 0 0 0 p1 100 0 0 p2 0 1 0 p3 0 0 40 size 100 1 10
gen zone brick p0 40 0 40 p1 100 0 40 p2 40 1 40 p3 50 0 50 p4 100 1 40 p5 50 1 50 p6 100 0 50 p7 100 1 50 size 60 1 5
gen zone brick p0 53 0 50 p1 100 0 50 p2 53 1 50 p3 63 0 60 p4 100 1 50 p5 63 1 60 p6 100 0 60 p7 100 1 60 size 47 1 5 attach face range z 39.9 40.1 attach face range z 49.9 50.1 fix x range x -0.1 0.1
fix x range x 99.9 100.1
fix y
fix z range z 0.1 -0.1
model elas
prop density 2500 bulk 3e9 shear 1e9
set gravity 0 0 -10
solve
ini xdisp 0 ydisp 0 zdisp 0
ini xvel 0 yvel 0 zvel 0
model mohr
prop density 2500 bulk 8.33e7 shear 3.846e7 c 42000 fric 17 ten 800000
solve fos file slope60kaiwa.sav associated
计算结果如下：
图4-a，网格剖分图
图4-b，速度矢量图
图4-c，速度等值线图图4-d，位移等值线图
最终计算边坡的稳定性系数为：Fs＝1.26 分析：
60°开挖后边坡稳定性系数采用的是FLAC3D 内置的强度折减法求解，稳定性系数1.26>1，从稳定性系数系数可以判断该边坡处于安全状态。

坡面最大速度为0.51mm/s，随着深度的增加，竖向应力逐渐增大。

坡肩处出现下沉，最大值达到3.24m。

第三部分：数值模拟技术研究应用实例分析
基于FLAC3D在不同土体参数条件下的边坡稳定性研究
摘要：边坡稳定性的影响因素很多，最直接的是边坡岩土体的性质、地下水、边坡坡脚等。

对于人工填土的路堤边坡而言，主要是填土的物理参数。

本文介绍了数值模拟法的工程应用的历史及现状，数值模拟在边坡中的应用。

然后利用FLAC3D软件，采用控制变量法，在保持其他土体物理参数不变的情况下，分别依次改变弹性模量、泊松比、密度、凝聚力、内摩擦角，计算边坡的稳定性，观察边坡稳定性系数变化，以及各点位移量的变化，从而得出土体物理参数变化对路堤边坡稳定性的影响，为路堤边坡填土土类提供依据和建议。

关键词：边坡、物理参数、FLAC3D
1 研究背景
21世纪以来，中国处于快速发展的阶段，国内基础建设蓬勃发展。

随着水利工程、公路等基础设施建设的大力开展，尤其是我国西部大开发战略的实施，大量公路建设深入西部山区，严重地破坏了局部区域内地质环境的平衡，大量的工程活动对地质环境的改变日益加剧，导致了大量地质灾害的发生，所以边坡的稳定性研究显得越来越重要。

边坡稳定性评价一直是边坡工程的一项主要内容，也是边坡工程设计和施工的基础。

边坡稳定性计算理论和判别方法可靠与否，关系到工程的安全问题，一旦边坡失稳，不仅会给国家带来巨大的经济损失，而且会危及人民生命财产安全。

2 国内外研究状态
目前边坡稳定性分析方法研究主要集中在与计算机技术、岩土力学、数学模型结合上，也产生和发展了一些新方法，其中三维稳定性分析是研究热点之一，以有限元法为代表的数值分析法以及各种不确定性分析方法发展迅速，而传统的极限平衡法主要以改进为主。

根据分析认为边坡稳定性分析研究在以后需要解决以下主要问题：继续完善和发展现有理论和方法，扬长避短；建立具有普遍意义的边坡失稳机理和稳定性评价方法；统一评价标准，增加各方法之间的对比性；建立多因素的综合评价方法；建立反映边坡各个时段稳定状态的全过程评价方法；建立符合边坡稳定分析的理论体系或组合理论体系；重视人类活动动态，加强人类活动与边坡稳定的相互作用研究。

3 FLAC3D模拟计算
3.1基本模型
以下面的简单的模型为基础，计算不同物理参数下的稳定性系数，来进行本次研究，模型及形体参数见图1。

基本物理参数见表1。

20
40
20
40
80
路堤
地基
图1 基本模型及形体参数
土体物理参数参数
弹性模量E /M Pa 泊松比μ 密度ρ/（g ·cm -3） 凝聚力C /kP
a
内
摩
擦ψ /（°
）
角 抗
拉
强
度σt
/M Pa
剪切弹性模量/MPa
)1(2μ+=
E G 体积弹性模量/MP a
)
21(3μ-=
E
K 10
0.3
2.5 42
17
0.8
38.5
83.3
在FLAC 3D 中建立模型，并计算边坡稳定性 ;创建几何模型： new
gen zone brick p0 0 0 0 p1 80 0 0 p2 0 2 0 p3 0 0 20 size 40 1 10 gen zone brick &
p0 40 0 20 p1 80 0 20 p2 40 2 20 p3 60 0 40 &
p4 80 2 20 p5 60 2 40 p6 80 0 40 p7 80 2 40 &
size 20 1 10
;赋予材料模型属性
model elas
set gravity 0,0,-10
property bulk=8e9 shear=5e9 density=2500
;施加边界约束
fix x range x -0.1 0.1
fix x range x 79.9 80.1
fix y
fix z range z -0.1 0.1
;试算
solve
;设置重力场
ini xdis=0 ydis=0 zdis=0
ini xvel=0 yvel=0 zvel=0
;赋予材料模型属性
model mohr
property bulk=0.833e8 shear=0.385e8 friction=17
property cohesion=0.42e5 tension=0.8e6
;fos计算
solve fos file slopefos.sav associated
图2 基本模型的稳定性系数和剪应变增量图
3.2 填土物理力学性质对路堤稳定性的影响分析
3.2.1 填土弹性模量对路堤稳定性的影响
在表1 基本参数的情况下，改变路堤填土弹性模量E的大小，考虑弹性模量为10MPa、15 MPa、20 MPa、25MPa、30 MPa 等5种取值，采用FLAC3D软件分析路堤边坡安全系数Fs随弹性模量E的变化情况及其弹性模量E对路堤沉降的影响情况，计算得到的填土不同弹性模量E下路堤边坡的值Fs见下表2。

弹性模量
10 15 20 25 30
E/MPa
稳定性系
1.09 1.09 1.09 1.09 1.09
数Fs
从表2 可以看出安全系数为定值Fs=1.10，即路堤填土的弹性模量E 对路堤
的稳定性几乎没有影响.
图3 路堤表面在不同弹性模量E下的z方向位移情况
图3为路堤表面在不同弹性模量E下的z方向位移（即路堤表面沉降）情况，从中可以看出：
（1）同一弹性模量E下，路堤表面各点的沉降量由路堤中心向边缘逐渐减小，并不成线性变化。

（2）路堤表面距路堤边缘10 m内沉降量随填土弹性模量E的增加而大幅增大。

10 m 范围外，路堤表面沉降量缓慢增加，距路堤边缘位置13m 位置沉降量最大，即此处为路堤滑动面贯通位置。

图4 路堤表面各点在不同弹性模量E下的z方向位移情况
从图4当中可以得出结论：随着弹性模量的增大，各点沉降量逐渐增大，且增加速率也逐步增大。

3.2.2 填土泊松比对路堤稳定性的影响
在表1基本参数的情况下，改变路堤填土泊松比μ的大小，考虑泊松比为0.15、0.20、0.25、0.30、0.35等5种取值，采用FLAC3D软件分析路堤边坡随泊松比μ的变化情况，计算得到填土不同泊松比μ下路堤边坡的Fs值见表3。

表3 填土不同泊松比μ下路堤边坡的值Fs
泊松比0.15 0.20 0.25 0.30 0.35
稳定性系
1.09 1.09 1.09 1.09 1.09
数Fs
从表3可以看出安全系数为定值Fs=1.09，即路堤填土的泊松比μ对路堤边坡的稳定性几乎没有影响.
图5 为μ=0.3 下的路堤边坡剪应变增量及其位移图。

从中可以看出，边坡滑动位置贯穿于坡底和左半幅路基表面1/4 处的位置，边坡的最大水平位移在坡底位置。

图5 μ=0.3下的路堤边坡剪应变增量及其位移图
图6 和图7 为在不同泊松比μ下路堤表面的z方向位移（即路堤表面沉降）情况，从中可以看出：
（1）同一泊松比下，路堤表面各点的沉降量由路堤中心向边缘逐渐减小，并不成线性变化，在距路堤中心10 m内，z方向变化较平缓，几乎为水平直线；
10 m以外，沉降量变化较快且数值较大，最大沉降位置在距边坡边缘约5 m 位置，可以认为该位置为路堤边坡表面最危险破坏滑动位置.
（2）路堤表面的沉降量随填土泊松比μ的增加而减小，距路堤中心10 m 内的具有线性关系，距路堤中心10 m 外沉降量变化较大. 当泊松比μ＝0.2时，沉降量发生突变，其大小比相邻泊松比的情况下较小，总体趋势还是沉降量随泊松比增多而减小，且线性较陡。

图6 路堤表面在不同泊松比μ下的z方向位移情况
图7 路堤表面各点在不同泊松比μ下的z方向位移情况
3.2.3填土凝聚力对路堤稳定性的影响
在表1基本参数的情况下，改变路堤填土凝聚力c的大小，考虑凝聚力为
10kPa、20 kPa、30 kPa、42 kPa、50 kPa等5种取值，采用FLAC3D软件分析路堤边坡随凝聚力c 的变化情况，绘制c-Fs关系曲线及其拟合公式如图8. 图9 为c=42 时路堤边坡的安全系数、剪应变增量云图和速度矢量图。

图8 c-Fs关系曲线及其拟合公式
图9 c=42 时路堤边坡的剪应变增量云图、安全系数和速度矢量图
从图8和图9可以看出，随着粘聚力的增大，边坡的安全系数呈线性增大。

分析原因主要是粘聚力的增大，边坡的抗剪强度也增大，使得边坡的抗滑能力也增加，从而使边坡趋于稳定，即安全系数增大。

图9 还可以明显看出塑性贯通区域，即路堤边坡的潜在滑动面，速度矢量图表明滑动面外侧区域各网格点的速度明显大于其他区域，说明这一区域已经出现滑动，即该区域发生了破坏。

3.2.4填土内摩擦角对路堤稳定性的影响
在表1 基本参数的情况下，改变路堤填土内摩擦角ψ的大小，考虑填土内摩擦角为13、15、17、19、21等5种取值，采用FLAC3D软件分析路堤边坡安全系数Fs随填土内摩擦角ψ的变化情况，绘制ψ-Fs关系曲线及其拟合公式如图10。

图10 ψ-Fs关系曲线及其拟合公式
从图10可以看出，随着粘聚力的增大，边坡的安全系数呈线性增大。

对比图9和图10，尤其是其中的拟合公式，由图10的拟合公式中斜率大于图9中的拟合公式中的斜率，故知内摩擦角对稳定性系数的影响大于凝聚力对稳定性系数的影响。

4 结论
经过对影响路堤边坡稳定性的填土物理力学性质进行数值模拟分析，可以得出以下结论：
（1）路堤填土弹性模量和泊松比对路堤边坡的稳定性没有影响，安全系数为定值. 但对路堤的沉降有一定程度的影响，沉降量随填土弹性模量、泊松比的增加而成非线性减小。

（2）路堤边坡安全系数与填土粘聚力、内摩擦角分别成线性关系，安全系数随着粘聚力、内摩擦角的增大而成线性增大，边坡趋于稳定，通过分析c-Fs 和ψ-Fs的关系，得出摩擦角比粘聚力对路堤边坡的稳定性影响较大。

（3）路堤边坡的安全系数随边坡的坡率、填土高度的增加而减小。

考虑路堤填料的物理力学性质对路堤边坡稳定性的影响，对公路建设过程中路堤填料的选择具有一定的指导作用。