一种基于边坡稳定性的山区公路线位优化方法

1.本发明属于山区公路设计领域，具体涉及基于边坡稳定性的山区公路线位优化方法。


背景技术：

2.在山区公路设计选线阶段就应当充分考虑如滑坡等地质灾害的影响，尤其为了解决山区复杂艰险地区的道路线路安全问题，需对山区公路通过边坡路段的线位优化开展研究。
3.目前山区公路通过边坡路段的线位优化方法研究不多，现有研究又主要集中在公路线形优化或山区公路边坡治理等方面，并未充分考虑边坡路段对山区公路线位优化的影响，所以开展基于边坡稳定性的山区公路线位优化方法的研究有较好的理论意义与工程价值。


技术实现要素：

4.本发明目的是为了解决目前山区公路优化方法主要集中在公路线形优化或山区公路边坡治理等方面，并未充分考虑边坡路段对山区公路线位优化的影响，致使山区公路边坡路段的滑坡风险高，使用寿命低，维修频率高，易发生交通事故，且缺少一种综合化的山区公路线位优化方法的问题，从而提出一种基于边坡稳定性的山区公路线位优化方法。
5.一种基于边坡稳定性的山区公路线位优化方法具体过程为：
6.步骤1：确定边坡和山区公路横断面几何参数；
7.步骤2：基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的模型；
8.步骤3：确定边坡的物理力学参数；
9.步骤4：基于步骤2和步骤3，将网格化的模型导入flac 3d软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算，得到地应力平衡后的三维数值计算模型；
10.步骤5：在flac 3d软件中采用强度折减法，对地应力平衡后的三维数值计算模型参数中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1
×
10-5
，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
11.步骤6：基于最后一次计算收敛后模型的动能与弹性剪切应变能数值，计算“动剪比”系数；
12.步骤7：改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值；
13.步骤8：对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度。
14.本发明的有益效果为：
15.本发明一种基于边坡稳定性的山区公路线位优化方法，首先，确定边坡和山区公路横断面几何参数；基于几何参数，确定边坡与公路横断面的数值计算模型中的控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的模型；基于“莫尔—库伦”准则，确定边坡的物理力学参数；
16.然后，将网格化的模型导入flac 3d软件，对网格化的模型赋值、设置边界条件并施加重力场以进行弹性计算；在flac 3d软件中采用强度折减法，对作归零平衡处理后模型参数中的粘聚力和内摩擦角进行强度折减，设置收敛精度条件1
×
10-5
，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；基于最后一次计算收敛后模型的动能与弹性剪切应变能数值，计算“动剪比”系数；
17.最后，改变公路横断面嵌入边坡深度的数值，重复上述的过程，直至取遍公路横断面嵌入边坡深度的所有数值；对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度。
18.解决了目前山区公路优化方法主要集中在公路线形优化或山区公路边坡治理等方面，并未充分考虑边坡路段对山区公路线位优化的影响，致使山区公路边坡路段的滑坡风险高，使用寿命低，维修频率高，易发生交通事故，且缺少一种综合化的山区公路线位优化方法的问题；
19.本发明充分考虑边坡路段对山区公路线位优化的影响，提高了山区公路使用寿命，降低了维修频率，提高了通行效率，降低了交通事故发生率，且能根据山区地质情况和公路线形参数动态优化山区公路线位。
附图说明
20.图1为本发明一种基于边坡稳定性的山区线位优化方法的流程图；
21.图2为本发明边坡与公路横断面二维计算模型示意图；
22.图3为本发明边坡与公路横断面计算模型示意图；
23.图4为本发明xz平面的网格划分数量示意图；
24.图5为本发明三个复合滑动面示意图。
具体实施方式
25.具体实施方式一：本实施方式一种基于边坡稳定性的山区公路线位优化方法具体过程为：
26.步骤1：确定边坡和山区公路横断面几何参数；
27.步骤2：基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的模型；
28.步骤3：基于“莫尔—库伦”准则，确定边坡的物理力学参数；
29.步骤4：基于步骤2和步骤3，将网格化的模型导入flac 3d软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算，得到地应力平衡后的三维数值计算模型；
30.步骤5：在flac 3d软件中采用强度折减法，对地应力平衡后的三维数值计算模型参数中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1
×
10-5
，对每次折减
计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
31.步骤6：基于最后一次计算收敛后模型的动能与弹性剪切应变能数值，计算“动剪比”系数；
32.步骤7：改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值；
33.步骤8：对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度。
34.具体实施方式二：本实施方式与具体实施方式一不同的是，所述步骤1中确定边坡和山区公路横断面几何参数；具体过程为：
35.几何参数主要包括边坡相对高度b、边坡坡度α、路堤填方坡度β、公路横断面宽度 s(取整)、公路横断面嵌入边坡深度a，0≤a≤s(取整)；
36.其中公路横断面宽度s和公路横断面嵌入边坡深度a为整数(向下取整，比如3.4取3，3.9取3)；
37.基于边坡坡度α确定未作处理的边坡；
38.基于路堤填方坡度β确定公路路堤；
39.基于边坡相对高度b确定未作处理的边坡、公路路堤、公路路面、公路路堑；
40.沿未作处理的边坡段前端向前延长30m，形成第一平台段；
41.在第一平台段的最前端向下延长20m，形成稳定基岩段；
42.接着从稳定基岩段最底端向右延长l，形成数值计算模型的底部边界；
43.最后从底部边界最右端向上垂直延长至边坡最高的位置，形成数值计算模型的右侧边界，右侧边界最顶端与公路路堑最顶端形成了第二平台段；
44.与数值计算模型的底部边界对应的顶部从左到右依次为：第一平台段、未作处理的边坡、公路路堤、公路路面、公路路堑、第二平台段。
45.几何参数及单位如图2和表1所示。
46.表1边坡与公路横断面几何参数表
47.参数符号参数含义单位α边坡坡度度(
°
)β路堤填方坡度度(
°
)b边坡相对高度米(m)s公路横断面宽度(取整)米(m)a公路横断面嵌入边坡内侧深度(取整)米(m)l计算模型长度米(m)
48.其它步骤及参数与具体实施方式一相同。
49.具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤2中基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的三维网格数值计算模型；具体过程如下：
50.步骤21、基于几何参数，在xz平面内确定数值计算模型的平面控制点坐标及其表达式；
51.步骤22、基于xz平面内控制点形成的平面，沿y轴方向拓展1m，形成三维数值计算模型；
52.步骤23：进行三维数值计算模型的网格划分。
53.其它步骤及参数与具体实施方式一至二之一相同。
54.具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤21中基于几何参数，在xz平面内确定数值计算模型的平面控制点坐标及其表达式，如图3和表2所示：
55.在xz平面内，将数值计算模型底部边界最左端设置为控制点o，数值计算模型底部边界最右端设置为控制点h，稳定基岩段最底端设置为控制点o，稳定基岩段最顶端设置为控制点a，第一平台段最左端设置为控制点a，第一平台段最右端设置为控制点b，未作处理的边坡最底端设置为控制点b，未作处理的边坡最顶端设置为控制点c，公路路堤最底端设置为控制点c，公路路堤最顶端设置为控制点d，公路路面最左端设置为控制点 d，公路路面最右端设置为控制点e，公路路堑最底端设置为控制点e，公路路堑最顶端设置为控制点f，第二平台段最左端设置为控制点f，第二平台段最右端设置为控制点g，右侧边界最顶端设置为控制点g，右侧边界最底端设置为控制点h；
56.其中ab段为第一平台段、ao段为稳定基岩段、oh段为数值计算模型的底部边界 l、bc段为未作处理的边坡、cd段为公路路堤、de段为公路路面、ef段为公路路堑、 fg段为第二平台段、gh段为数值计算模型的右侧边界；
57.以控制点o为原点，oa方向为z轴，oh方向为x轴；
58.将控制点b在底部边界的投影点(控制点b向底部边界做垂线，在底部边界的交点) 设置为控制点x1；
59.将控制点c在底部边界的投影点设置为控制点x2；
60.将控制点d在底部边界的投影点设置为控制点x3；
61.将控制点e在底部边界的投影点设置为控制点x4；
62.将控制点f在底部边界的投影点设置为控制点x5；
63.将控制点g在底部边界的投影点设置为控制点x6，控制点x6和控制点h为同一控制点；
64.将控制点a在右侧边界的投影点(控制点a向右侧边界做垂线，在右侧边界的交点) 设置为控制点z1；
65.将控制点b在右侧边界的投影点设置为控制点z1；
66.将控制点c在右侧边界的投影点设置为控制点z2；
67.将控制点d在右侧边界的投影点设置为控制点z3；
68.将控制点e在右侧边界的投影点设置为控制点z3；
69.将控制点f在右侧边界的投影点设置为控制点z4，控制点z4和控制点g为同一控制点；
70.控制点o的x轴坐标为0，z轴坐标为0；
71.控制点a的x轴坐标为0，z轴坐标为z1；
72.控制点b的x轴坐标为x1，z轴坐标为z1；
73.控制点c的x轴坐标为x2，z轴坐标为z2；
74.控制点d的x轴坐标为x3，z轴坐标为z3；
75.控制点e的x轴坐标为x4，z轴坐标为z3；
76.控制点f的x轴坐标为x5，z轴坐标为z4；
77.控制点g的x轴坐标为x6，z轴坐标为z4；
78.控制点h的x轴坐标为x6，z轴坐标为0；
79.表2数值计算模型控制点坐标表达式
80.控制点x坐标z坐标o00a0z1bx1z1cx2z2dx3z3ex4z3fx5z4gx6z4hx60
81.由几何关系可得各坐标的表达式如下：
[0082][0083][0084][0085]
其它步骤及参数与具体实施方式一至三之一相同。
[0086]
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述步骤22中在xz平面内基于控制点形成的平面，沿y轴方向拓展1m，形成三维数值计算模型，在此过程中拓展方向应为朝着y轴正方向。
[0087]
其它步骤及参数与具体实施方式一至四之一相同。
[0088]
具体实施方式六：本实施方式与具体实施方式一至五之一不同的是，所述步骤23中进行三维数值计算模型的网格划分，具体过程为：
[0089]
网格划分单元应为八面体，每边尺寸不宜超过0.5m，xz平面的网格划分数量如图4 所示：
[0090]
由几何关系和坐标可计算各边划分的网格数量，int()表示取整数：
[0091]
将三维数值计算模型进行网格划分，划分为八个网格体；
[0092]
o、a、b、x1构成网格体1；
[0093]
线段az1与线段ex4的交点为δ，x1、b、δ、x4构成网格体2；
[0094]
x4、δ、z1、x6构成网格体3；
[0095]
线段cz2与线段ex4的交点为φ，b、c、φ、δ构成网格体4；
[0096]
δ、φ、z2、z1构成网格体5；
[0097]
c、d、e、φ构成网格体6；
[0098]
φ、e、z3、z2构成网格体7；
[0099]
e、f、g、z3构成网格体8；
[0100]
网格体1的ab段划分的网格数量为w1，网格体1的oa段划分的网格数量为w4，网格体1的ox1段划分的网格数量为w1，网格体1的x1b段划分的网格数量为w4；
[0101]
网格体2的x1x4段划分的网格数量为w2，网格体2的bδ段划分的网格数量为w2，网格体2的x1b段划分的网格数量为w4，网格体2的x4δ段划分的网格数量为w4；
[0102]
网格体3的x4x6段划分的网格数量为w3，网格体3的x4δ段划分的网格数量为w4，网格体3的δz1段划分的网格数量为w3，网格体3的z1x6段划分的网格数量为w4，
[0103]
网格体4的bδ段划分的网格数量为w2，网格体4的bc段划分的网格数量为w5，网格体4的cφ段划分的网格数量为w2，网格体4的φδ段划分的网格数量为w5；
[0104]
网格体5的φδ段划分的网格数量为w5，网格体5的φz2段划分的网格数量为w3，网格体5的z2z1段划分的网格数量为w5，网格体5的z1δ段划分的网格数量为w3；
[0105]
网格体6的cd段划分的网格数量为w6，网格体6的de段划分的网格数量为w2，网格体6的eφ段划分的网格数量为w6，网格体6的φc段划分的网格数量为w2；
[0106]
网格体7的eφ段划分的网格数量为w6，网格体7的ez3段划分的网格数量为w3，网格体7的z3z2段划分的网格数量为w6，网格体7的z2φ段划分的网格数量为w3；
[0107]
网格体8的ef段划分的网格数量为w7，网格体8的fg段划分的网格数量为w3，网格体8的gz3段划分的网格数量为w7，网格体8的z3e段划分的网格数量为w3；
[0108]
其中
[0109][0110][0111][0112]
沿y轴方向的网格划分数量wy＝2，网格尺寸为0.5m(y轴方向拓展1m，划分为2，每个网格y轴方向尺寸为0.5m)，int()表示取整数。
[0113]
其它步骤及参数与具体实施方式一至五之一相同。
[0114]
具体实施方式七：本实施方式与具体实施方式一至六之一不同的是，所述步骤3中基于“莫尔—库伦”准则，确定边坡的物理力学参数，如表3所示：
[0115]
边坡的物理力学参数包括：弹性模量e1、泊松比υ、土体密度ρ、内摩擦角粘聚力c。
[0116]
表3边坡的物理力学参数
[0117][0118][0119]
其它步骤及参数与具体实施方式一至六之一相同。
[0120]
具体实施方式八：本实施方式与具体实施方式一至七之一不同的是，所述步骤4中基于步骤2和步骤3，将网格化的模型导入flac 3d软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算，得到地应力平衡后的三维数值计算模型；具体过程为：
[0121]
步骤41、基于步骤3中的弹性模量e1和泊松比υ计算体积模量k1和切变模量g1，将网格化的模型导入flac 3d软件中，并将体积模量k1、切变模量g1、土体密度ρ、内摩擦角和粘聚力c的参数值赋予flac 3d软件中的网格模型；
[0122]
其中体积模量k1和剪切模量g1计算公式如下：
[0123][0124][0125]
式中：k1为体积模量，单位为pa；g1为切变模量，单位为pa；e1为弹性模量，单位为pa；v为泊松比；
[0126]
步骤42、针对三维网格数值计算模型设置边界条件，对模型的左侧边界(ao段及拓展形成的面，即稳定基岩段)和右侧边界(gh段及拓展形成的面)进行左右约束，即边界上的网格单元沿x轴方向不能移动，对底部边界(oh段及拓展形成的面)进行全约束，即底部边界的网格单元完全固定；对三维网格数值计算模型前后两个面进行拓展方向约束，即前后两个面沿y轴方向不能移动；
[0127]
步骤43、对模型施加重力场，取重力加速度g＝-9.8m/s2(沿z轴负方向)，计算收敛精度设置为1
×
10-5
，即表示计算过程中当网格体系的最大不平衡力与典型内力的比率小于1
×
10-5
时计算结束(所谓体系最大不平衡力，是指每一个计算时步中，外力通过网格节点传递分配到体系各节点时，所有节点的外力与内力之差中的最大值；所谓典型内力，则是指计算模型所有网格点力的平均值)，然后将已经赋值且约束完成的三维网格数值计算模型基于“莫尔—库伦”准则进行弹性求解(“莫尔—库伦”准则是flac 3d软件中自带的，此处弹性求解即为地应力平衡)，得到速度、位移和塑性区等，对速度、位移和塑性区作归零平衡处理，得到地应力平衡后的三维数值计算模型。
[0128]
其它步骤及参数与具体实施方式一至七之一相同。
[0129]
具体实施方式九：本实施方式与具体实施方式一至八之一不同的是，所述步骤5中在 flac 3d软件中采用强度折减法，对地应力平衡后的三维数值计算模型中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1
×
10-5
，对每次折减计算收敛后的三维网
格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；具体过程为：
[0130]
步骤51、在flac 3d软件中采用强度折减法，对地应力平衡后的三维数值计算模型中的粘聚力和内摩擦角进行强度折减(强度折减后内摩擦角与粘聚力c改变)，将未改变的体积模量k1、切变模量g1、土体密度ρ的参数值和改变后内摩擦角和粘聚力c的参数值赋予flac 3d软件中的网格模型(三维数值计算模型)，再次基于“莫尔-库伦”准则进行弹性求解，设置收敛精度为1
×
10-5
，即表示计算过程中当网格体系的最大不平衡力与典型内力的比率小于等于1
×
10-5
时计算结束(所谓体系最大不平衡力，是指每一个计算时步中，外力通过网格节点传递分配到体系各节点时，所有节点的外力与内力之差中的最大值；所谓典型内力，则是指计算模型所有网格点力的平均值)，保存计算模型与结果，之后再对内摩擦角与粘聚力c继续折减，重复步骤51步骤中的内容；
[0131]
强度折减法的原理是同步等比例折减边坡土体的内摩擦角与粘聚力c，方法如下：
[0132][0133][0134]
式中：fi为强度折减系数，且fi＞1；
[0135]ci
、分别为折减后的土体粘聚力(单位：pa)与内摩擦角(单位：
°
)；
[0136]
c、分别为折减前的土体粘聚力(单位：pa)与内摩擦角(单位：
°
)；
[0137]
步骤52、强度折减过程中，观察数值计算收敛精度的变化，直至数值计算不收敛为止，即flac 3d软件中模型网格体系最大不平衡力与典型内力的比率一直大于1
×
10-5
， (flac 3d软件中模型网格体系最大不平衡力与典型内力的比率一直大于1
×
10-5
，表示收敛)，即数值计算精度长时间无法收敛至1
×
10-5
以内，完成数值计算；
[0138]
步骤53、数值计算完成后，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
[0139]
所述步骤53中计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能，计算方法如下：
[0140]
基于最后一次折减计算收敛后的模型，提取第j个网格单元的质点速度(j为网格单元编号)，计算第j个网格单元的动能：
[0141][0142]
式中：为第j个网格单元的动能，单位为j；ρj为第j个网格单元的密度，等于土体密度ρ，即ρj＝ρ，单位为kg/m3；为第j个网格单元的质点速度，单位为m/s；
[0143]
逐一提取所有网格单元的质点速度并计算所有网格单元的动能，累计所有网格单元的动能，求得整个三维网格数值计算模型的动能，计算公式如下：
[0144][0145]
式中：uk为三维网格数值计算模型的动能，单位为j；n为三维网格数值计算模型的网格单元数量；为第j个网格单元的动能，单位为j；ρj为第j个网格单元的密度，等于土体密度ρ，即ρj＝ρ，单位为kg/m3；为第j个网格单元的质点速度，单位为m/s；
[0146]
所述步骤53中计算并提取最后一次折减计算收敛后三维网格数值计算模型的弹性剪切应变能，其定义和计算方法如下：
[0147]
针对岩土材料，若从能量角度解释其屈服行为，可在三个主应力构成的莫尔库伦圆中，将其代表性体积单元(模型中的网格单元)的内部看作存在三个复合滑动面(cmp)，即图5中和的作用面，同时假设岩土材料屈服时三个复合滑动面的弹性剪切应变能之和为常数。
[0148]
基于能量准则假设，三条莫尔库伦圆的切线在横轴上相交于三轴等拉强度应力点a，假定是土体的内摩擦角和可由几何关系分别求得；然后，基于几何关系得出三个复合滑动面上法向应力和剪切应力，接着可求得岩石材料屈服时三个复合滑动面的弹性剪切应变能p
13
、p
12
和p
23
，最终得到弹性剪切应变能之和p。
[0149]
基于网格单元的第一主应力(最大主应力)和第二主应力(中间主应力)确定第三半圆，基于网格单元的第一主应力(最大主应力)和第三主应力(最小主应力)确定第一半圆，基于网格单元的第二主应力(中间主应力)和第三主应力(最小主应力)确定第二半圆；
[0150]
假定是土体的内摩擦角第一半圆的切线在横轴上相交于点a，点a与第一半圆的切线在横轴上的交点夹角为
[0151]
令第二半圆的切线在横轴上的交点位于a点，确定a与第二半圆的切线在横轴上的交点夹角为
[0152]
令第三半圆的切线在横轴上的交点位于a点，确定a与第三半圆的切线在横轴上的交点夹角为
[0153]
切线与半圆相交的点为复合滑动面；
[0154]
步骤531、基于最后一次折减计算收敛后的模型，提取第j个网格单元的三个主应力σ
1-j
、σ
2-j
、σ
3-j
，计算第j个网格单元的lode角正切值：
[0155][0156]
式中：tanθj为第j个网格单元的lode角正切值；σ
1-j
、σ
2-j
、σ
3-j
分别为第j个网格单元的第一主应力(最大主应力)、第二主应力(中间主应力)和第三主应力(最小主应力)，单位为pa；第一主应力(最大主应力)、第二主应力(中间主应力)和第三主应力(最小主应力)是上面自动保存时每个网格就有的，“莫尔—库伦”准则得到的；
[0157]
步骤532、基于三个复合滑动面准则，计算第j个网格单元的三个复合滑动面的几何角度和的正弦值：
[0158][0159][0160][0161]
式中：分别为第j个网格单元三个复合滑动面的几何角度的正弦值；tanθj为第j个网格单元的lode角正切值；为第j个网格单元的内摩擦角，等于土体的内摩擦角即
[0162]
步骤533、计算第j个网格单元三个复合滑动面上的法向应力和剪切应力：
[0163][0164][0165][0166][0167][0168][0169]
式中：σ
12-j
、τ
12-j
分别为第j个网格单元第一个复合滑动面的法向应力和剪切应力，单位为pa；σ
23-j
、τ
23-j
分别为第j个网格单元第二个复合滑动面的法向应力和剪切应力，单位为pa；σ
13-j
、τ
13-j
分别为第j个网格单元第三个复合滑动面的法向应力和剪切应力，单位为pa；
[0170]
步骤534、分别计算第j个网格单元三个复合滑动面的弹性剪切应变能及其之和：
[0171][0172][0173][0174]
pj＝p
12-j
+p
23-j
+p
13-j
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0175]
式中：pj为第j个网格单元的总弹性剪切应变能，单位为j；p
13-j
、p
12-j
和p
23-j
分别为第j个网格单元中三个复合滑动面的弹性剪切应变能，单位为j；g
1-j
为第j个网格单元的切变模量，单位为pa，等于土体的切变模量，即g
1-j
＝g1；
[0176]
步骤535、逐一计算所有网格单元的总弹性剪切应变能并累计求和，得到三维网格数值计算模型总的弹性剪切应变能，计算公式如下：
[0177][0178]
式中：p为模型总的弹性剪切应变能，单位为j；n为模型的网格单元数量；pj为第 j个网格单元的总弹性剪切应变能，单位为j。
[0179]
其它步骤及参数与具体实施方式一至九之一相同。
[0180]
具体实施方式十：本实施方式与具体实施方式一至十之一不同的是，所述步骤6中基于最后一次计算收敛后模型的动能与弹性剪切应变能，计算“动剪比”系数up，其计算公式如下：
[0181][0182]
式中：up为“动剪比”系数；uk为整个模型的动能，单位为j；p为模型总的弹性剪切应变能，单位为j；
[0183]
所述步骤7中改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值，0≤a≤s；
[0184]
公路横断面嵌入边坡深度a取整数(向下取整数，比如5.2取5；5.9取5)，公路横断面s取整数(向下取整数，比如5.2取5；5.9取5)，a的取值为0，1，2，
·····
，s，合计s+1种取值方式。
[0185]
特别指出，公路横断面嵌入边坡深度a和公路横断面s的数值都取整数，所以a的取值为0，1，2，
·····
，s(合计s+1种取值方式)。
[0186]
所述步骤8中对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度；具体过程为：
[0187]
根据步骤7中公路横断面嵌入边坡深度a的s+1种取值，求得s+1个“动剪比”系数，以“动剪比”系数小于0.2为标准(up＜0.2)确定合理的嵌入深度。
[0188]
其它步骤及参数与具体实施方式一至九之一相同。
[0189]
采用以下实施例验证本发明的有益效果：
[0190]
实施例一：
[0191]
下面结合附图以及具体实施例子对本发明的技术方案做出详细说明。如图1所示，为本发明的一种基于边坡稳定性的山区线位优化方法的流程图，包括以下步骤：
[0192]
步骤1：确定边坡和山区公路横断面几何参数；如表4和图2所示
[0193]
表4边坡与公路横断面几何参数表
[0194]
参数符号参数含义取值α边坡坡度45
°
β路堤填方坡度56.3
°
b边坡相对高度60ms公路横断面宽度(取整)10ma公路横断面嵌入边坡内侧深度(取整)5ml计算模型长度140m
[0195]
步骤2：基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，如表5和图3所示，建立三维计算模型，确定网格划分数量，得到网格化的模型；
[0196][0197][0198][0199]
表5数值计算模型控制点坐标
[0200][0201][0202]
基于控制点形成的平面，沿y轴方向拓展1m，形成三维数值计算模型，在此过程中拓展方向应为朝着y轴正方向。
[0203]
网格划分单元应为八面体，每边尺寸不宜超过0.5m，xz平面的网格划分如图4所示：
[0204]
由几何关系和坐标可计算各边划分的网格数量，int()表示取整数：
[0205][0206]
[0207][0208]
沿y轴方向的网格划分数量wy＝2，网格尺寸为0.5m。
[0209]
步骤3：基于“莫尔—库伦”准则，确定边坡的物理力学参数，如表6所示：
[0210]
表6边坡的物理力学参数取值
[0211][0212]
步骤4：基于步骤2和步骤3，将网格化的模型导入flac 3d软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算；
[0213]
步骤41、基于步骤3中的弹性模量e1和泊松比v换算成体积模量k1和切变模量g1，将网格化的模型导入flac 3d软件中，并将体积模量k1、切变模量g1、土体密度ρ、内摩擦角和粘聚力c的参数值赋予flac 3d软件中的网格模型，其中体积模量k1和剪切模量g1计算如下：
[0214][0215][0216]
式中：k1—体积模量；g1—切变模量；e1—弹性模量；v—泊松比。
[0217]
步骤42、针对三维网格数值计算模型设置边界条件，对模型的左侧边界(ao段及拓展形成的面，即稳定基岩段)和右侧边界(gh段及拓展形成的面)进行左右约束，即边界上的网格单元沿x轴方向不能移动，对底部边界(oh段及拓展形成的面)进行全约束，即底部边界的网格单元完全固定；对三维网格数值计算模型前后两个面进行拓展方向约束，即前后两个面沿y轴方向不能移动；
[0218]
步骤43、对模型施加重力场，取重力加速度g＝-9.8m/s2(沿z轴负方向)，计算收敛精度设置为1
×
10-5
，然后将已经赋值且约束完成的三维网格数值计算模型基于“莫尔—库伦”准则进行弹性求解(“莫尔—库伦”准则是flac 3d软件中自带的)，得到速度、位移和塑性区，对速度、位移和塑性区作归零平衡处理，得到地应力平衡后的三维数值计算模型。
[0219]
步骤5：在flac 3d软件中采用强度折减法，对地应力平衡后的三维数值计算模型参数中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1
×
10-5
，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
[0220]
步骤51、在flac 3d软件中采用强度折减法，对地应力平衡后的三维数值计算模型参数中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1
×
10-5
，当flac 3d 软件中体系最大不平衡力与典型力的比值小于等于1
×
10-5
时表示收敛，对每次折减计算收
敛后的三维网格数值计算模型进行自动保存。强度折减法的原理是同步等比例折减边坡土体的内摩擦角与粘聚力c，方法如下：
[0221][0222][0223]
式中：fi为强度折减系数，且fi＞1；
[0224]ci
、分别为折减后的土体粘聚力(单位：pa)与内摩擦角(单位：
°
)；
[0225]
c、分别为折减前的土体粘聚力(单位：pa)与内摩擦角(单位：
°
)。
[0226]
步骤52、强度折减过程中，观察数值计算收敛精度的变化，直至数值计算不收敛为止，即flac 3d软件中体系最大不平衡力与典型力的比值一直大于等于1
×
10-5
，完成数值计算。
[0227]
步骤53、数值计算完成后，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
[0228]
根据数值计算结果，最后一次折减计算收敛时的折减系数为1.321，计算和提取得到的动能和弹性剪切应变能分别：
[0229]
动能：uk＝7.03923
×
10-7
j；弹性剪切应变能：p＝4.59857
×
107j。
[0230]
步骤6：基于最后一次计算收敛后模型的动能与弹性剪切应变能数值，计算“动剪比”系数：
[0231][0232]
式中：up—“动剪比”系数；
[0233]
uk—整个模型的动能(单位：j)；
[0234]
p—模型总的弹性剪切应变能(单位：j)。
[0235]
步骤7：改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值(a＝1m，2m，3m，4m，5m，6m，7m，8m，9m 和10m)。
[0236]
步骤8：对比分析公路横断面嵌入边坡不同深度的“动剪比”系数，如表7所示，以“动剪比”系数小于0.2为标准(up＜0.2)确定合理的嵌入深度。
[0237]
表7公路横断面嵌入边坡深度与“动剪比”系数
[0238][0239][0240]
根据表7所示，以“动剪比”系数小于0.2为标准(up＜0.2)确定合理的嵌入深度，即合理的嵌入深度a＝3m(up＝0.129)，a＝4m(up＝0.108)，a＝5m(up＝0.153)和a＝6m (up＝
0.129)。
[0241]
本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。