[发明专利]一种基于边坡稳定性的山区公路线位优化方法

权利要求书

1.一种基于边坡稳定性的山区公路线位优化方法，其特征在于：所述方法具体过程为：

步骤1：确定边坡和山区公路横断面几何参数；

步骤2：基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的模型；

步骤3：确定边坡的物理力学参数；

步骤4：基于步骤2和步骤3，将网格化的模型导入FLAC 3D软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算，得到地应力平衡后的三维数值计算模型；

步骤5：在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型参数中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1×10^-5，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；

步骤6：基于最后一次计算收敛后模型的动能与弹性剪切应变能数值，计算“动剪比”系数；

步骤7：改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值；

步骤8：对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度；

所述步骤1中确定边坡和山区公路横断面几何参数；具体过程为：

几何参数包括边坡相对高度b、边坡坡度α、路堤填方坡度β、公路横断面宽度s、公路横断面嵌入边坡深度a，0≤a≤s；

其中公路横断面宽度s和公路横断面嵌入边坡深度a为整数；

基于边坡坡度α确定未作处理的边坡；

基于路堤填方坡度β确定公路路堤；

基于边坡相对高度b确定未作处理的边坡、公路路堤、公路路面、公路路堑；

沿未作处理的边坡段前端向前延长30m，形成第一平台段；

在第一平台段的最前端向下延长20m，形成稳定基岩段；

从稳定基岩段最底端向右延长L，形成数值计算模型的底部边界；

从底部边界最右端向上垂直延长至边坡最高的位置，形成数值计算模型的右侧边界，右侧边界最顶端与公路路堑最顶端形成了第二平台段；

与数值计算模型的底部边界对应的顶部从左到右依次为：第一平台段、未作处理的边坡、公路路堤、公路路面、公路路堑、第二平台段；

所述步骤2中基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的三维网格数值计算模型；具体过程如下：

步骤21、基于几何参数，在XZ平面内确定数值计算模型的平面控制点坐标及其表达式；

步骤22、基于XZ平面内控制点形成的平面，沿Y轴方向拓展1m，形成三维数值计算模型；

步骤23：进行三维数值计算模型的网格划分；

所述步骤21中基于几何参数，在XZ平面内确定数值计算模型的平面控制点坐标及其表达式；具体过程为：

在XZ平面内，将数值计算模型底部边界最左端设置为控制点O，数值计算模型底部边界最右端设置为控制点H，稳定基岩段最底端设置为控制点O，稳定基岩段最顶端设置为控制点A，第一平台段最左端设置为控制点A，第一平台段最右端设置为控制点B，未作处理的边坡最底端设置为控制点B，未作处理的边坡最顶端设置为控制点C，公路路堤最底端设置为控制点C，公路路堤最顶端设置为控制点D，公路路面最左端设置为控制点D，公路路面最右端设置为控制点E，公路路堑最底端设置为控制点E，公路路堑最顶端设置为控制点F，第二平台段最左端设置为控制点F，第二平台段最右端设置为控制点G，右侧边界最顶端设置为控制点G，右侧边界最底端设置为控制点H；

其中AB段为第一平台段、AO段为稳定基岩段、OH段为数值计算模型的底部边界L、BC段为未作处理的边坡、CD段为公路路堤、DE段为公路路面、EF段为公路路堑、FG段为第二平台段、GH段为数值计算模型的右侧边界；

以控制点O为原点，OA方向为Z轴，OH方向为X轴；

将控制点B在底部边界的投影点设置为控制点X1；

将控制点C在底部边界的投影点设置为控制点X2；

将控制点D在底部边界的投影点设置为控制点X3；

将控制点E在底部边界的投影点设置为控制点X4；

将控制点F在底部边界的投影点设置为控制点X5；

将控制点G在底部边界的投影点设置为控制点X6，控制点X6和控制点H为同一控制点；

将控制点A在右侧边界的投影点设置为控制点Z1；

将控制点B在右侧边界的投影点设置为控制点Z1；

将控制点C在右侧边界的投影点设置为控制点Z2；

将控制点D在右侧边界的投影点设置为控制点Z3；

将控制点E在右侧边界的投影点设置为控制点Z3；

将控制点F在右侧边界的投影点设置为控制点Z4，控制点Z4和控制点G为同一控制点；

由几何关系可得各坐标的表达式如下：

X1＝30，X6＝L；

Z1＝20，Z4＝20+b；

所述步骤22中在XZ平面内基于控制点形成的平面，沿Y轴方向拓展1m，形成三维数值计算模型；

所述步骤23中进行三维数值计算模型的网格划分，具体过程为：

将三维数值计算模型进行网格划分，划分为八个网格体；

O、A、B、X1构成网格体1；

线段AZ1与线段EX4的交点为δ，X1、B、δ、X4构成网格体2；

X4、δ、Z1、X6构成网格体3；

线段CZ2与线段EX4的交点为φ，B、C、φ、δ构成网格体4；

δ、φ、Z2、Z1构成网格体5；

C、D、E、φ构成网格体6；

φ、E、Z3、Z2构成网格体7；

E、F、G、Z3构成网格体8；

网格体1的AB段划分的网格数量为w₁，网格体1的OA段划分的网格数量为w₄，网格体1的OX1段划分的网格数量为w₁，网格体1的X1B段划分的网格数量为w₄；

网格体2的X1X4段划分的网格数量为w₂，网格体2的Bδ段划分的网格数量为w₂，网格体2的X1B段划分的网格数量为w₄，网格体2的X4δ段划分的网格数量为w₄；

网格体3的X4X6段划分的网格数量为w₃，网格体3的X4δ段划分的网格数量为w₄，网格体3的δZ1段划分的网格数量为w₃，网格体3的Z1X6段划分的网格数量为w₄，

网格体4的Bδ段划分的网格数量为w₂，网格体4的BC段划分的网格数量为w₅，网格体4的Cφ段划分的网格数量为w₂，网格体4的φδ段划分的网格数量为w₅；

网格体5的φδ段划分的网格数量为w₅，网格体5的φZ2段划分的网格数量为w₃，网格体5的Z2Z1段划分的网格数量为w₅，网格体5的Z1δ段划分的网格数量为w₃；

网格体6的CD段划分的网格数量为w₆，网格体6的DE段划分的网格数量为w₂，网格体6的Eφ段划分的网格数量为w₆，网格体6的φC段划分的网格数量为w₂；

网格体7的Eφ段划分的网格数量为w₆，网格体7的EZ3段划分的网格数量为w₃，网格体7的Z3Z2段划分的网格数量为w₆，网格体7的Z2φ段划分的网格数量为w₃；

网格体8的EF段划分的网格数量为w₇，网格体8的FG段划分的网格数量为w₃，网格体8的GZ3段划分的网格数量为w₇，网格体8的Z3E段划分的网格数量为w₃；

其中

w₁＝60，w₃＝2L-60-w₂，w₄＝40，

沿Y轴方向的网格划分数量w_y＝2，网格尺寸为0.5m，int()表示取整数；

所述步骤3中确定边坡的物理力学参数；具体过程为：

边坡的物理力学参数包括：弹性模量E₁、泊松比υ、土体密度ρ、内摩擦角粘聚力；

所述步骤4中基于步骤2和步骤3，将网格化的模型导入FLAC 3D软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算，得到地应力平衡后的三维数值计算模型；具体过程为：

步骤41、基于步骤3中的弹性模量E₁和泊松比υ计算体积模量K₁和切变模量G₁，将网格化的模型导入FLAC 3D软件中，并将体积模量K₁、切变模量G₁、土体密度ρ、内摩擦角和粘聚力c的参数值赋予FLAC 3D软件中的网格模型；

其中体积模量K₁和剪切模量G₁计算公式如下：

式中：K₁为体积模量，单位为Pa；G₁为切变模量，单位为Pa；E₁为弹性模量，单位为Pa；v为泊松比；

步骤42、针对三维网格数值计算模型设置边界条件，对模型的左侧边界和右侧边界进行左右约束，即边界上的网格单元沿X轴方向不能移动，对底部边界进行全约束，即底部边界的网格单元完全固定；对三维网格数值计算模型前后两个面进行拓展方向约束，即前后两个面沿Y轴方向不能移动；

步骤43、对模型施加重力场，取重力加速度g＝-9.8m/s²，计算收敛精度设置为1×10-⁵，即表示计算过程中当网格体系的最大不平衡力与典型内力的比率小于1×10^-5时计算结束，然后将已经赋值且约束完成的三维网格数值计算模型基于“莫尔—库伦”准则进行弹性求解，得到速度、位移和塑性区，对速度、位移和塑性区作归零平衡处理，得到地应力平衡后的三维数值计算模型；

所述步骤5中在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1×10-⁵，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；具体过程为：

步骤51、在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型中的粘聚力和内摩擦角进行强度折减，将未改变的体积模量K₁、切变模量G₁、土体密度ρ的参数值和改变后内摩擦角和粘聚力c的参数值赋予FLAC 3D软件中的网格模型，再次基于“莫尔-库伦”准则进行弹性求解，设置收敛精度为1×10-⁵，即表示计算过程中当网格体系的最大不平衡力与典型内力的比率小于等于1×10-⁵时计算结束，保存计算模型与结果，之后再对内摩擦角与粘聚力c继续折减，重复步骤51步骤中的内容；

步骤52、强度折减过程中，观察数值计算收敛精度的变化，直至数值计算不收敛为止，即FLAC 3D软件中模型网格体系最大不平衡力与典型内力的比率一直大于1×10-⁵，完成数值计算；

步骤53、数值计算完成后，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；

所述步骤53中计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能，计算方法如下：

基于最后一次折减计算收敛后的模型，提取第j个网格单元的质点速度计算第j个网格单元的动能：

式中：为第j个网格单元的动能，单位为J；ρ_j为第j个网格单元的密度，等于土体密度ρ，即ρ_j＝ρ，单位为kg/m³；为第j个网格单元的质点速度，单位为m/s；

逐一提取所有网格单元的质点速度并计算所有网格单元的动能，累计所有网格单元的动能，求得整个三维网格数值计算模型的动能，计算公式如下：

式中：U^k为三维网格数值计算模型的动能，单位为J；n为三维网格数值计算模型的网格单元数量；为第j个网格单元的动能，单位为J；ρ_j为第j个网格单元的密度，等于土体密度ρ，即ρ_j＝ρ，单位为kg/m³；为第j个网格单元的质点速度，单位为m/s；

所述步骤53中计算并提取最后一次折减计算收敛后三维网格数值计算模型的弹性剪切应变能，计算方法如下：

基于网格单元的第一主应力和第二主应力确定第三半圆，基于网格单元的第一主应力和第三主应力确定第一半圆，基于网格单元的第二主应力和第三主应力确定第二半圆；

假定是土体的内摩擦角第一半圆的切线在横轴上相交于点A，点A与第一半圆的切线在横轴上的交点夹角为

令第二半圆的切线在横轴上的交点位于A点，确定A与第二半圆的切线在横轴上的交点夹角为

令第三半圆的切线在横轴上的交点位于A点，确定A与第三半圆的切线在横轴上的交点夹角为

切线与半圆相交的点为复合滑动面；

步骤531、基于最后一次折减计算收敛后的模型，提取第j个网格单元的三个主应力σ_1-j、σ_2-j、σ_3-j，计算第j个网格单元的Lode角正切值：

式中：tanθ_j为第j个网格单元的Lode角正切值；σ_1-j、σ_2-j、σ_3-j分别为第j个网格单元的第一主应力、第二主应力和第三主应力，单位为Pa；

步骤532、计算第j个网格单元的三个复合滑动面的几何角度和的正弦值：

式中：分别为第j个网格单元三个复合滑动面的几何角度和的正弦值；tanθ_j为第j个网格单元的Lode角正切值；为第j个网格单元的内摩擦角，等于土体的内摩擦角即

步骤533、计算第j个网格单元三个复合滑动面上的法向应力和剪切应力：

式中：σ_12-j、τ_12-j分别为第j个网格单元第一个复合滑动面的法向应力和剪切应力，单位为Pa；σ_23-j、τ_23-j分别为第j个网格单元第二个复合滑动面的法向应力和剪切应力，单位为Pa；σ_13-j、τ_13-j分别为第j个网格单元第三个复合滑动面的法向应力和剪切应力，单位为Pa；

步骤534、分别计算第j个网格单元三个复合滑动面的弹性剪切应变能及其之和：

P_j＝P_12-j+P_23-j+P_13-j       (20)

式中：P_j为第j个网格单元的总弹性剪切应变能，单位为J；P_13-j、P_12-j和P_23-j分别为第j个网格单元中三个复合滑动面的弹性剪切应变能，单位为J；G_1-j为第j个网格单元的切变模量，单位为Pa，等于土体的切变模量，即G_1-j＝G₁；

步骤535、逐一计算所有网格单元的总弹性剪切应变能并累计求和，得到三维网格数值计算模型总的弹性剪切应变能，计算公式如下：

式中：P为模型总的弹性剪切应变能，单位为J；n为模型的网格单元数量；P_j为第j个网格单元的总弹性剪切应变能，单位为J；

所述步骤6中基于最后一次计算收敛后模型的动能与弹性剪切应变能，计算“动剪比”系数up，其计算公式如下：

式中：up为“动剪比”系数；U^k为整个模型的动能，单位为J；P为模型总的弹性剪切应变能，单位为J；

所述步骤7中改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值，0≤a≤s；具体过程为：

公路横断面嵌入边坡深度a取整数，公路横断面s取整数，a的取值为0，1，2，·····，s，合计s+1种取值方式；

所述步骤8中对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度；具体过程为：

根据步骤7中公路横断面嵌入边坡深度a的s+1种取值，求得s+1个“动剪比”系数，以“动剪比”系数小于0.2为标准确定合理的嵌入深度。