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

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

Info

Publication number
CN115168953B
CN115168953B CN202210811745.0A CN202210811745A CN115168953B CN 115168953 B CN115168953 B CN 115168953B CN 202210811745 A CN202210811745 A CN 202210811745A CN 115168953 B CN115168953 B CN 115168953B
Authority
CN
China
Prior art keywords
grid
section
control point
model
numerical calculation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210811745.0A
Other languages
English (en)
Other versions
CN115168953A (zh
Inventor
孟祥海
魏鹏儒
张明扬
张瑜荣
阿格交力·沙丁
郑敏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202210811745.0A priority Critical patent/CN115168953B/zh
Publication of CN115168953A publication Critical patent/CN115168953A/zh
Application granted granted Critical
Publication of CN115168953B publication Critical patent/CN115168953B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0635—Risk analysis of enterprise or organisation activities
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Human Resources & Organizations (AREA)
  • Computer Hardware Design (AREA)
  • Strategic Management (AREA)
  • Entrepreneurship & Innovation (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Economics (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Game Theory and Decision Science (AREA)
  • Educational Administration (AREA)
  • Computer Graphics (AREA)
  • Tourism & Hospitality (AREA)
  • Marketing (AREA)
  • Development Economics (AREA)
  • Quality & Reliability (AREA)
  • General Business, Economics & Management (AREA)
  • Architecture (AREA)
  • Civil Engineering (AREA)
  • Structural Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Road Paving Structures (AREA)

Abstract

一种基于边坡稳定性的山区公路线位优化方法，本发明涉及基于边坡稳定性的山区公路线位优化方法。本发明目的是为了解决目前山区公路优化方法并未充分考虑边坡路段对山区公路线位优化的影响，致使山区公路边坡路段的滑坡风险高，使用寿命低，维修频率高，易发生交通事故，且缺少一种综合化的山区公路线位优化方法的问题。过程为：1：确定边坡和山区公路横断面几何参数；2：得到网格化的模型；3：确定边坡的物理力学参数；4：对网格化的模型进行弹性计算；5：计算模型的动能与弹性剪切应变能；6：计算动剪比；7：重复1至6，直至取遍公路横断面嵌入边坡深度的所有数值；8：确定合理的嵌入深度。本发明用于山区公路设计领域。

Description

一种基于边坡稳定性的山区公路线位优化方法
技术领域
本发明属于山区公路设计领域，具体涉及基于边坡稳定性的山区公路线位优化方法。
背景技术
在山区公路设计选线阶段就应当充分考虑如滑坡等地质灾害的影响，尤其为了解决山区复杂艰险地区的道路线路安全问题，需对山区公路通过边坡路段的线位优化开展研究。
目前山区公路通过边坡路段的线位优化方法研究不多，现有研究又主要集中在公路线形优化或山区公路边坡治理等方面，并未充分考虑边坡路段对山区公路线位优化的影响，所以开展基于边坡稳定性的山区公路线位优化方法的研究有较好的理论意义与工程价值。
发明内容
本发明目的是为了解决目前山区公路优化方法主要集中在公路线形优化或山区公路边坡治理等方面，并未充分考虑边坡路段对山区公路线位优化的影响，致使山区公路边坡路段的滑坡风险高，使用寿命低，维修频率高，易发生交通事故，且缺少一种综合化的山区公路线位优化方法的问题，从而提出一种基于边坡稳定性的山区公路线位优化方法。
一种基于边坡稳定性的山区公路线位优化方法具体过程为：
步骤1：确定边坡和山区公路横断面几何参数；
步骤2：基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的模型；
步骤3：确定边坡的物理力学参数；
步骤4：基于步骤2和步骤3，将网格化的模型导入FLAC 3D软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算，得到地应力平衡后的三维数值计算模型；
步骤5：在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型参数中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1×10-5，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
步骤6：基于最后一次计算收敛后模型的动能与弹性剪切应变能数值，计算“动剪比”系数；
步骤7：改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值；
步骤8：对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度。
本发明的有益效果为：
本发明一种基于边坡稳定性的山区公路线位优化方法，首先，确定边坡和山区公路横断面几何参数；基于几何参数，确定边坡与公路横断面的数值计算模型中的控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的模型；基于“莫尔—库伦”准则，确定边坡的物理力学参数；
然后，将网格化的模型导入FLAC 3D软件，对网格化的模型赋值、设置边界条件并施加重力场以进行弹性计算；在FLAC 3D软件中采用强度折减法，对作归零平衡处理后模型参数中的粘聚力和内摩擦角进行强度折减，设置收敛精度条件1×10-5，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；基于最后一次计算收敛后模型的动能与弹性剪切应变能数值，计算“动剪比”系数；
最后，改变公路横断面嵌入边坡深度的数值，重复上述的过程，直至取遍公路横断面嵌入边坡深度的所有数值；对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度。
解决了目前山区公路优化方法主要集中在公路线形优化或山区公路边坡治理等方面，并未充分考虑边坡路段对山区公路线位优化的影响，致使山区公路边坡路段的滑坡风险高，使用寿命低，维修频率高，易发生交通事故，且缺少一种综合化的山区公路线位优化方法的问题；
本发明充分考虑边坡路段对山区公路线位优化的影响，提高了山区公路使用寿命，降低了维修频率，提高了通行效率，降低了交通事故发生率，且能根据山区地质情况和公路线形参数动态优化山区公路线位。
附图说明
图1为本发明一种基于边坡稳定性的山区线位优化方法的流程图；
图2为本发明边坡与公路横断面二维计算模型示意图；
图3为本发明边坡与公路横断面计算模型示意图；
图4为本发明XZ平面的网格划分数量示意图；
图5为本发明三个复合滑动面示意图。
具体实施方式
具体实施方式一：本实施方式一种基于边坡稳定性的山区公路线位优化方法具体过程为：
步骤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为整数(向下取整，比如3.4取3，3.9取3)；
基于边坡坡度α确定未作处理的边坡；
基于路堤填方坡度β确定公路路堤；
基于边坡相对高度b确定未作处理的边坡、公路路堤、公路路面、公路路堑；
沿未作处理的边坡段前端向前延长30m，形成第一平台段；
在第一平台段的最前端向下延长20m，形成稳定基岩段；
接着从稳定基岩段最底端向右延长L，形成数值计算模型的底部边界；
最后从底部边界最右端向上垂直延长至边坡最高的位置，形成数值计算模型的右侧边界，右侧边界最顶端与公路路堑最顶端形成了第二平台段；
与数值计算模型的底部边界对应的顶部从左到右依次为：第一平台段、未作处理的边坡、公路路堤、公路路面、公路路堑、第二平台段。
几何参数及单位如图2和表1所示。
表1边坡与公路横断面几何参数表
参数符号 参数含义 单位
α 边坡坡度 度(°)
β 路堤填方坡度 度(°)
b 边坡相对高度 米(m)
s 公路横断面宽度(取整) 米(m)
a 公路横断面嵌入边坡内侧深度(取整) 米(m)
L 计算模型长度 米(m)
其它步骤及参数与具体实施方式一相同。
具体实施方式三：本实施方式与具体实施方式一或二不同的是，所述步骤2中基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，建立三维计算模型，确定网格划分数量，得到网格化的三维网格数值计算模型；具体过程如下：
步骤21、基于几何参数，在XZ平面内确定数值计算模型的平面控制点坐标及其表达式；
步骤22、基于XZ平面内控制点形成的平面，沿Y轴方向拓展1m，形成三维数值计算模型；
步骤23：进行三维数值计算模型的网格划分。
其它步骤及参数与具体实施方式一至二之一相同。
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是，所述步骤21中基于几何参数，在XZ平面内确定数值计算模型的平面控制点坐标及其表达式，如图3和表2所示：
在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在底部边界的投影点(控制点B向底部边界做垂线，在底部边界的交点)设置为控制点X1；
将控制点C在底部边界的投影点设置为控制点X2；
将控制点D在底部边界的投影点设置为控制点X3；
将控制点E在底部边界的投影点设置为控制点X4；
将控制点F在底部边界的投影点设置为控制点X5；
将控制点G在底部边界的投影点设置为控制点X6，控制点X6和控制点H为同一控制点；
将控制点A在右侧边界的投影点(控制点A向右侧边界做垂线，在右侧边界的交点)设置为控制点Z1；
将控制点B在右侧边界的投影点设置为控制点Z1；
将控制点C在右侧边界的投影点设置为控制点Z2；
将控制点D在右侧边界的投影点设置为控制点Z3；
将控制点E在右侧边界的投影点设置为控制点Z3；
将控制点F在右侧边界的投影点设置为控制点Z4，控制点Z4和控制点G为同一控制点；
控制点O的X轴坐标为0，Z轴坐标为0；
控制点A的X轴坐标为0，Z轴坐标为Z1；
控制点B的X轴坐标为X1，Z轴坐标为Z1；
控制点C的X轴坐标为X2，Z轴坐标为Z2；
控制点D的X轴坐标为X3，Z轴坐标为Z3；
控制点E的X轴坐标为X4，Z轴坐标为Z3；
控制点F的X轴坐标为X5，Z轴坐标为Z4；
控制点G的X轴坐标为X6，Z轴坐标为Z4；
控制点H的X轴坐标为X6，Z轴坐标为0；
表2数值计算模型控制点坐标表达式
控制点 X坐标 Z坐标
O 0 0
A 0 Z1
B X1 Z1
C X2 Z2
D X3 Z3
E X4 Z3
F X5 Z4
G X6 Z4
H X6 0
由几何关系可得各坐标的表达式如下：
Figure SMS_1
Figure SMS_2
Figure SMS_3
其它步骤及参数与具体实施方式一至三之一相同。
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是，所述步骤22中在XZ平面内基于控制点形成的平面，沿Y轴方向拓展1m，形成三维数值计算模型，在此过程中拓展方向应为朝着Y轴正方向。
其它步骤及参数与具体实施方式一至四之一相同。
具体实施方式六：本实施方式与具体实施方式一至五之一不同的是，所述步骤23中进行三维数值计算模型的网格划分，具体过程为：
网格划分单元应为八面体，每边尺寸不宜超过0.5m，XZ平面的网格划分数量如图4所示：
由几何关系和坐标可计算各边划分的网格数量，int()表示取整数：
将三维数值计算模型进行网格划分，划分为八个网格体；
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段划分的网格数量为w1，网格体1的OA段划分的网格数量为w4，网格体1的OX1段划分的网格数量为w1，网格体1的X1B段划分的网格数量为w4；
网格体2的X1X4段划分的网格数量为w2，网格体2的Bδ段划分的网格数量为w2，网格体2的X1B段划分的网格数量为w4，网格体2的X4δ段划分的网格数量为w4；
网格体3的X4X6段划分的网格数量为w3，网格体3的X4δ段划分的网格数量为w4，网格体3的δZ1段划分的网格数量为w3，网格体3的Z1X6段划分的网格数量为w4，
网格体4的Bδ段划分的网格数量为w2，网格体4的BC段划分的网格数量为w5，网格体4的Cφ段划分的网格数量为w2，网格体4的φδ段划分的网格数量为w5；
网格体5的φδ段划分的网格数量为w5，网格体5的φZ2段划分的网格数量为w3，网格体5的Z2Z1段划分的网格数量为w5，网格体5的Z1δ段划分的网格数量为w3；
网格体6的CD段划分的网格数量为w6，网格体6的DE段划分的网格数量为w2，网格体6的Eφ段划分的网格数量为w6，网格体6的φC段划分的网格数量为w2；
网格体7的Eφ段划分的网格数量为w6，网格体7的EZ3段划分的网格数量为w3，网格体7的Z3Z2段划分的网格数量为w6，网格体7的Z2φ段划分的网格数量为w3；
网格体8的EF段划分的网格数量为w7，网格体8的FG段划分的网格数量为w3，网格体8的GZ3段划分的网格数量为w7，网格体8的Z3E段划分的网格数量为w3；
其中
Figure SMS_4
Figure SMS_5
Figure SMS_6
沿Y轴方向的网格划分数量wy＝2，网格尺寸为0.5m(Y轴方向拓展1m，划分为2，每个网格Y轴方向尺寸为0.5m)，int()表示取整数。
其它步骤及参数与具体实施方式一至五之一相同。
具体实施方式七：本实施方式与具体实施方式一至六之一不同的是，所述步骤3中基于“莫尔—库伦”准则，确定边坡的物理力学参数，如表3所示：
边坡的物理力学参数包括：弹性模量E1、泊松比υ、土体密度ρ、内摩擦角
Figure SMS_7
粘聚力c。
表3边坡的物理力学参数
Figure SMS_8
Figure SMS_9
其它步骤及参数与具体实施方式一至六之一相同。
具体实施方式八：本实施方式与具体实施方式一至七之一不同的是，所述步骤4中基于步骤2和步骤3，将网格化的模型导入FLAC 3D软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算，得到地应力平衡后的三维数值计算模型；具体过程为：
步骤41、基于步骤3中的弹性模量E1和泊松比υ计算体积模量K1和切变模量G1，将网格化的模型导入FLAC 3D软件中，并将体积模量K1、切变模量G1、土体密度ρ、内摩擦角
Figure SMS_10
和粘聚力c的参数值赋予FLAC 3D软件中的网格模型；
其中体积模量K1和剪切模量G1计算公式如下：
Figure SMS_11
Figure SMS_12
式中：K1为体积模量，单位为Pa；G1为切变模量，单位为Pa；E1为弹性模量，单位为Pa；v为泊松比；
步骤42、针对三维网格数值计算模型设置边界条件，对模型的左侧边界(AO段及拓展形成的面，即稳定基岩段)和右侧边界(GH段及拓展形成的面)进行左右约束，即边界上的网格单元沿X轴方向不能移动，对底部边界(OH段及拓展形成的面)进行全约束，即底部边界的网格单元完全固定；对三维网格数值计算模型前后两个面进行拓展方向约束，即前后两个面沿Y轴方向不能移动；
步骤43、对模型施加重力场，取重力加速度g＝-9.8m/s2(沿Z轴负方向)，计算收敛精度设置为1×10-5，即表示计算过程中当网格体系的最大不平衡力与典型内力的比率小于1×10-5时计算结束(所谓体系最大不平衡力，是指每一个计算时步中，外力通过网格节点传递分配到体系各节点时，所有节点的外力与内力之差中的最大值；所谓典型内力，则是指计算模型所有网格点力的平均值)，然后将已经赋值且约束完成的三维网格数值计算模型基于“莫尔—库伦”准则进行弹性求解(“莫尔—库伦”准则是FLAC 3D软件中自带的，此处弹性求解即为地应力平衡)，得到速度、位移和塑性区等，对速度、位移和塑性区作归零平衡处理，得到地应力平衡后的三维数值计算模型。
其它步骤及参数与具体实施方式一至七之一相同。
具体实施方式九：本实施方式与具体实施方式一至八之一不同的是，所述步骤5中在 FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1×10-5，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；具体过程为：
步骤51、在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型中的粘聚力和内摩擦角进行强度折减(强度折减后内摩擦角
Figure SMS_13
与粘聚力c改变)，将未改变的体积模量K1、切变模量G1、土体密度ρ的参数值和改变后内摩擦角/>
Figure SMS_14
和粘聚力c的参数值赋予FLAC 3D软件中的网格模型(三维数值计算模型)，再次基于“莫尔-库伦”准则进行弹性求解，设置收敛精度为1×10-5，即表示计算过程中当网格体系的最大不平衡力与典型内力的比率小于等于1×10-5时计算结束(所谓体系最大不平衡力，是指每一个计算时步中，外力通过网格节点传递分配到体系各节点时，所有节点的外力与内力之差中的最大值；所谓典型内力，则是指计算模型所有网格点力的平均值)，保存计算模型与结果，之后再对内摩擦角
Figure SMS_15
与粘聚力c继续折减，重复步骤51步骤中的内容；
强度折减法的原理是同步等比例折减边坡土体的内摩擦角
Figure SMS_16
与粘聚力c，方法如下：
Figure SMS_17
Figure SMS_18
式中：Fi为强度折减系数，且Fi＞1；
ci、
Figure SMS_19
分别为折减后的土体粘聚力(单位：Pa)与内摩擦角(单位：°)；
c、
Figure SMS_20
分别为折减前的土体粘聚力(单位：Pa)与内摩擦角(单位：°)；
步骤52、强度折减过程中，观察数值计算收敛精度的变化，直至数值计算不收敛为止，即FLAC 3D软件中模型网格体系最大不平衡力与典型内力的比率一直大于1×10-5，(FLAC 3D软件中模型网格体系最大不平衡力与典型内力的比率一直大于1×10-5，表示收敛)，即数值计算精度长时间无法收敛至1×10-5以内，完成数值计算；
步骤53、数值计算完成后，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
所述步骤53中计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能，计算方法如下：
基于最后一次折减计算收敛后的模型，提取第j个网格单元的质点速度
Figure SMS_21
(j为网格单元编号)，计算第j个网格单元的动能：
Figure SMS_22
式中：
Figure SMS_23
为第j个网格单元的动能，单位为J；ρj为第j个网格单元的密度，等于土体密度ρ，即ρj＝ρ，单位为kg/m3；/>
Figure SMS_24
为第j个网格单元的质点速度，单位为m/s；
逐一提取所有网格单元的质点速度并计算所有网格单元的动能，累计所有网格单元的动能，求得整个三维网格数值计算模型的动能，计算公式如下：
Figure SMS_25
式中：Uk为三维网格数值计算模型的动能，单位为J；n为三维网格数值计算模型的网格单元数量；
Figure SMS_26
为第j个网格单元的动能，单位为J；ρj为第j个网格单元的密度，等于土体密度ρ，即ρj＝ρ，单位为kg/m3；/>
Figure SMS_27
为第j个网格单元的质点速度，单位为m/s；
所述步骤53中计算并提取最后一次折减计算收敛后三维网格数值计算模型的弹性剪切应变能，其定义和计算方法如下：
针对岩土材料，若从能量角度解释其屈服行为，可在三个主应力构成的莫尔库伦圆中，将其代表性体积单元(模型中的网格单元)的内部看作存在三个复合滑动面(CMP)，即图5中
Figure SMS_28
和/>
Figure SMS_29
的作用面，同时假设岩土材料屈服时三个复合滑动面的弹性剪切应变能之和为常数。
基于能量准则假设，三条莫尔库伦圆的切线在横轴上相交于三轴等拉强度应力点A，假定
Figure SMS_30
是土体的内摩擦角/>
Figure SMS_31
和/>
Figure SMS_32
可由几何关系分别求得；然后，基于几何关系得出三个复合滑动面上法向应力和剪切应力，接着可求得岩石材料屈服时三个复合滑动面的弹性剪切应变能P13、P12和P23，最终得到弹性剪切应变能之和P。
基于网格单元的第一主应力(最大主应力)和第二主应力(中间主应力)确定第三半圆，基于网格单元的第一主应力(最大主应力)和第三主应力(最小主应力)确定第一半圆，基于网格单元的第二主应力(中间主应力)和第三主应力(最小主应力)确定第二半圆；
假定
Figure SMS_33
是土体的内摩擦角/>
Figure SMS_34
第一半圆的切线在横轴上相交于点A，点A与第一半圆的切线在横轴上的交点夹角为/>
Figure SMS_35
令第二半圆的切线在横轴上的交点位于A点，确定A与第二半圆的切线在横轴上的交点夹角为
Figure SMS_36
令第三半圆的切线在横轴上的交点位于A点，确定A与第三半圆的切线在横轴上的交点夹角为
Figure SMS_37
切线与半圆相交的点为复合滑动面；
步骤531、基于最后一次折减计算收敛后的模型，提取第j个网格单元的三个主应力σ1-j、σ2-j、σ3-j，计算第j个网格单元的Lode角正切值：
Figure SMS_38
式中：tanθj为第j个网格单元的Lode角正切值；σ1-j、σ2-j、σ3-j分别为第j个网格单元的第一主应力(最大主应力)、第二主应力(中间主应力)和第三主应力(最小主应力)，单位为Pa；第一主应力(最大主应力)、第二主应力(中间主应力)和第三主应力(最小主应力)是上面自动保存时每个网格就有的，“莫尔—库伦”准则得到的；
步骤532、基于三个复合滑动面准则，计算第j个网格单元的三个复合滑动面的几何角度
Figure SMS_39
和/>
Figure SMS_40
的正弦值：
Figure SMS_41
Figure SMS_42
Figure SMS_43
式中：
Figure SMS_44
分别为第j个网格单元三个复合滑动面的几何角度/>
Figure SMS_45
的正弦值；tanθj为第j个网格单元的Lode角正切值；/>
Figure SMS_46
为第j个网格单元的内摩擦角，等于土体的内摩擦角/>
Figure SMS_47
即/>
Figure SMS_48
步骤533、计算第j个网格单元三个复合滑动面上的法向应力和剪切应力：
Figure SMS_49
Figure SMS_50
Figure SMS_51
Figure SMS_52
Figure SMS_53
Figure SMS_54
式中：σ12-j、τ12-j分别为第j个网格单元第一个复合滑动面的法向应力和剪切应力，单位为Pa；σ23-j、τ23-j分别为第j个网格单元第二个复合滑动面的法向应力和剪切应力，单位为Pa；σ13-j、τ13-j分别为第j个网格单元第三个复合滑动面的法向应力和剪切应力，单位为Pa；
步骤534、分别计算第j个网格单元三个复合滑动面的弹性剪切应变能及其之和：
Figure SMS_55
Figure SMS_56
Figure SMS_57
Pj＝P12-j+P23-j+P13-j (20)
式中：Pj为第j个网格单元的总弹性剪切应变能，单位为J；P13-j、P12-j和P23-j分别为第j个网格单元中三个复合滑动面的弹性剪切应变能，单位为J；G1-j为第j个网格单元的切变模量，单位为Pa，等于土体的切变模量，即G1-j＝G1；
步骤535、逐一计算所有网格单元的总弹性剪切应变能并累计求和，得到三维网格数值计算模型总的弹性剪切应变能，计算公式如下：
Figure SMS_58
式中：P为模型总的弹性剪切应变能，单位为J；n为模型的网格单元数量；Pj为第 j个网格单元的总弹性剪切应变能，单位为J。
其它步骤及参数与具体实施方式一至九之一相同。
具体实施方式十：本实施方式与具体实施方式一至十之一不同的是，所述步骤6中基于最后一次计算收敛后模型的动能与弹性剪切应变能，计算“动剪比”系数up，其计算公式如下：
Figure SMS_59
式中：up为“动剪比”系数；Uk为整个模型的动能，单位为J；P为模型总的弹性剪切应变能，单位为J；
所述步骤7中改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值，0≤a≤s；
公路横断面嵌入边坡深度a取整数(向下取整数，比如5.2取5；5.9取5)，公路横断面s取整数(向下取整数，比如5.2取5；5.9取5)，a的取值为0，1，2，·····，s，合计s+1种取值方式。
特别指出，公路横断面嵌入边坡深度a和公路横断面s的数值都取整数，所以a的取值为0，1，2，·····，s(合计s+1种取值方式)。
所述步骤8中对比分析不同公路横断面嵌入边坡深度的“动剪比”系数，以“动剪比”系数小于0.2为标准，确定合理的嵌入深度；具体过程为：
根据步骤7中公路横断面嵌入边坡深度a的s+1种取值，求得s+1个“动剪比”系数，以“动剪比”系数小于0.2为标准(up＜0.2)确定合理的嵌入深度。
其它步骤及参数与具体实施方式一至九之一相同。
采用以下实施例验证本发明的有益效果：
实施例一：
下面结合附图以及具体实施例子对本发明的技术方案做出详细说明。如图1所示，为本发明的一种基于边坡稳定性的山区线位优化方法的流程图，包括以下步骤：
步骤1：确定边坡和山区公路横断面几何参数；如表4和图2所示
表4边坡与公路横断面几何参数表
参数符号 参数含义 取值
α 边坡坡度 45°
β 路堤填方坡度 56.3°
b 边坡相对高度 60m
s 公路横断面宽度(取整) 10m
a 公路横断面嵌入边坡内侧深度(取整) 5m
L 计算模型长度 140m
步骤2：基于步骤1中的几何参数，确定边坡与公路横断面的数值计算模型控制点坐标，如表5和图3所示，建立三维计算模型，确定网格划分数量，得到网格化的模型；
Figure SMS_60
Figure SMS_61
Figure SMS_62
表5数值计算模型控制点坐标
Figure SMS_63
Figure SMS_64
基于控制点形成的平面，沿Y轴方向拓展1m，形成三维数值计算模型，在此过程中拓展方向应为朝着Y轴正方向。
网格划分单元应为八面体，每边尺寸不宜超过0.5m，XZ平面的网格划分如图4所示：
由几何关系和坐标可计算各边划分的网格数量，int()表示取整数：
Figure SMS_65
Figure SMS_66
/>
Figure SMS_67
沿Y轴方向的网格划分数量wy＝2，网格尺寸为0.5m。
步骤3：基于“莫尔—库伦”准则，确定边坡的物理力学参数，如表6所示：
表6边坡的物理力学参数取值
Figure SMS_68
步骤4：基于步骤2和步骤3，将网格化的模型导入FLAC 3D软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算；
步骤41、基于步骤3中的弹性模量E1和泊松比v换算成体积模量K1和切变模量G1，将网格化的模型导入FLAC 3D软件中，并将体积模量K1、切变模量G1、土体密度ρ、内摩擦角
Figure SMS_69
和粘聚力c的参数值赋予FLAC 3D软件中的网格模型，其中体积模量K1和剪切模量G1计算如下：
Figure SMS_70
Figure SMS_71
式中：K1—体积模量；G1—切变模量；E1—弹性模量；v—泊松比。
步骤42、针对三维网格数值计算模型设置边界条件，对模型的左侧边界(AO段及拓展形成的面，即稳定基岩段)和右侧边界(GH段及拓展形成的面)进行左右约束，即边界上的网格单元沿X轴方向不能移动，对底部边界(OH段及拓展形成的面)进行全约束，即底部边界的网格单元完全固定；对三维网格数值计算模型前后两个面进行拓展方向约束，即前后两个面沿Y轴方向不能移动；
步骤43、对模型施加重力场，取重力加速度g＝-9.8m/s2(沿Z轴负方向)，计算收敛精度设置为1×10-5，然后将已经赋值且约束完成的三维网格数值计算模型基于“莫尔—库伦”准则进行弹性求解(“莫尔—库伦”准则是FLAC 3D软件中自带的)，得到速度、位移和塑性区，对速度、位移和塑性区作归零平衡处理，得到地应力平衡后的三维数值计算模型。
步骤5：在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型参数中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1×10-5，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
步骤51、在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型参数中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1×10-5，当FLAC 3D软件中体系最大不平衡力与典型力的比值小于等于1×10-5时表示收敛，对每次折减计算收敛后的三维网格数值计算模型进行自动保存。强度折减法的原理是同步等比例折减边坡土体的内摩擦角
Figure SMS_72
与粘聚力c，方法如下：
Figure SMS_73
Figure SMS_74
式中：Fi为强度折减系数，且Fi＞1；
ci、
Figure SMS_75
分别为折减后的土体粘聚力(单位：Pa)与内摩擦角(单位：°)；
c、
Figure SMS_76
分别为折减前的土体粘聚力(单位：Pa)与内摩擦角(单位：°)。
步骤52、强度折减过程中，观察数值计算收敛精度的变化，直至数值计算不收敛为止，即FLAC 3D软件中体系最大不平衡力与典型力的比值一直大于等于1×10-5，完成数值计算。
步骤53、数值计算完成后，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
根据数值计算结果，最后一次折减计算收敛时的折减系数为1.321，计算和提取得到的动能和弹性剪切应变能分别：
动能：Uk＝7.03923×10-7J；弹性剪切应变能：P＝4.59857×107J。
步骤6：基于最后一次计算收敛后模型的动能与弹性剪切应变能数值，计算“动剪比”系数：
Figure SMS_77
式中：up—“动剪比”系数；
Uk—整个模型的动能(单位：J)；
P—模型总的弹性剪切应变能(单位：J)。
步骤7：改变公路横断面嵌入边坡深度的数值，重复步骤1至步骤6的过程，直至取遍公路横断面嵌入边坡深度的所有数值(a＝1m，2m，3m，4m，5m，6m，7m，8m，9m 和10m)。
步骤8：对比分析公路横断面嵌入边坡不同深度的“动剪比”系数，如表7所示，以“动剪比”系数小于0.2为标准(up＜0.2)确定合理的嵌入深度。
表7公路横断面嵌入边坡深度与“动剪比”系数
Figure SMS_78
Figure SMS_79
根据表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)。
本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (1)

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，
Figure FDA0004228518440000031
Figure FDA0004228518440000032
X6＝L；
Z1＝20，
Figure FDA0004228518440000033
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段划分的网格数量为w1，网格体1的OA段划分的网格数量为w4，网格体1的OX1段划分的网格数量为w1，网格体1的X1B段划分的网格数量为w4；
网格体2的X1X4段划分的网格数量为w2，网格体2的Bδ段划分的网格数量为w2，网格体2的X1B段划分的网格数量为w4，网格体2的X4δ段划分的网格数量为w4；
网格体3的X4X6段划分的网格数量为w3，网格体3的X4δ段划分的网格数量为w4，网格体3的δZ1段划分的网格数量为w3，网格体3的Z1X6段划分的网格数量为w4，
网格体4的Bδ段划分的网格数量为w2，网格体4的BC段划分的网格数量为w5，网格体4的Cφ段划分的网格数量为w2，网格体4的φδ段划分的网格数量为w5；
网格体5的φδ段划分的网格数量为w5，网格体5的φZ2段划分的网格数量为w3，网格体5的Z2Z1段划分的网格数量为w5，网格体5的Z1δ段划分的网格数量为w3；
网格体6的CD段划分的网格数量为w6，网格体6的DE段划分的网格数量为w2，网格体6的Eφ段划分的网格数量为w6，网格体6的φC段划分的网格数量为w2；
网格体7的Eφ段划分的网格数量为w6，网格体7的EZ3段划分的网格数量为w3，网格体7的Z3Z2段划分的网格数量为w6，网格体7的Z2φ段划分的网格数量为w3；
网格体8的EF段划分的网格数量为w7，网格体8的FG段划分的网格数量为w3，网格体8的GZ3段划分的网格数量为w7，网格体8的Z3E段划分的网格数量为w3；
其中
w1＝60，
Figure FDA0004228518440000041
w3＝2L-60-w2，w4＝40，
Figure FDA0004228518440000042
Figure FDA0004228518440000043
沿Y轴方向的网格划分数量wy＝2，网格尺寸为0.5m，int()表示取整数；
所述步骤3中确定边坡的物理力学参数；具体过程为：
边坡的物理力学参数包括：弹性模量E1、泊松比υ、土体密度ρ、内摩擦角
Figure FDA0004228518440000044
粘聚力/>
Figure FDA0004228518440000051
；
所述步骤4中基于步骤2和步骤3，将网格化的模型导入FLAC 3D软件，对网格化的模型赋值、设置边界条件并施加重力场进行弹性计算，得到地应力平衡后的三维数值计算模型；具体过程为：
步骤41、基于步骤3中的弹性模量E1和泊松比υ计算体积模量K1和切变模量G1，将网格化的模型导入FLAC 3D软件中，并将体积模量K1、切变模量G1、土体密度ρ、内摩擦角
Figure FDA0004228518440000054
和粘聚力c的参数值赋予FLAC 3D软件中的网格模型；
其中体积模量K1和剪切模量G1计算公式如下：
Figure FDA0004228518440000052
Figure FDA0004228518440000053
式中：K1为体积模量，单位为Pa；G1为切变模量，单位为Pa；E1为弹性模量，单位为Pa；v为泊松比；
步骤42、针对三维网格数值计算模型设置边界条件，对模型的左侧边界和右侧边界进行左右约束，即边界上的网格单元沿X轴方向不能移动，对底部边界进行全约束，即底部边界的网格单元完全固定；对三维网格数值计算模型前后两个面进行拓展方向约束，即前后两个面沿Y轴方向不能移动；
步骤43、对模型施加重力场，取重力加速度g＝-9.8m/s2，计算收敛精度设置为1×10-5，即表示计算过程中当网格体系的最大不平衡力与典型内力的比率小于1×10-5时计算结束，然后将已经赋值且约束完成的三维网格数值计算模型基于“莫尔—库伦”准则进行弹性求解，得到速度、位移和塑性区，对速度、位移和塑性区作归零平衡处理，得到地应力平衡后的三维数值计算模型；
所述步骤5中在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型中的粘聚力和内摩擦角进行强度折减，同时设置收敛精度条件为1×10-5，对每次折减计算收敛后的三维网格数值计算模型进行自动保存，直至数值计算不收敛为止，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；具体过程为：
步骤51、在FLAC 3D软件中采用强度折减法，对地应力平衡后的三维数值计算模型中的粘聚力和内摩擦角进行强度折减，将未改变的体积模量K1、切变模量G1、土体密度ρ的参数值和改变后内摩擦角
Figure FDA0004228518440000061
和粘聚力c的参数值赋予FLAC 3D软件中的网格模型，再次基于“莫尔-库伦”准则进行弹性求解，设置收敛精度为1×10-5，即表示计算过程中当网格体系的最大不平衡力与典型内力的比率小于等于1×10-5时计算结束，保存计算模型与结果，之后再对内摩擦角/>
Figure FDA0004228518440000062
与粘聚力c继续折减，重复步骤51步骤中的内容；
步骤52、强度折减过程中，观察数值计算收敛精度的变化，直至数值计算不收敛为止，即FLAC 3D软件中模型网格体系最大不平衡力与典型内力的比率一直大于1×10-5，完成数值计算；
步骤53、数值计算完成后，计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能与弹性剪切应变能；
所述步骤53中计算并提取最后一次折减计算收敛后三维网格数值计算模型的动能，计算方法如下：
基于最后一次折减计算收敛后的模型，提取第j个网格单元的质点速度
Figure FDA0004228518440000063
计算第j个网格单元的动能：
Figure FDA0004228518440000064
式中：
Figure FDA0004228518440000065
为第j个网格单元的动能，单位为J；ρj为第j个网格单元的密度，等于土体密度ρ，即ρj＝ρ，单位为kg/m3；/>
Figure FDA0004228518440000066
为第j个网格单元的质点速度，单位为m/s；
逐一提取所有网格单元的质点速度并计算所有网格单元的动能，累计所有网格单元的动能，求得整个三维网格数值计算模型的动能，计算公式如下：
Figure FDA0004228518440000067
式中：Uk为三维网格数值计算模型的动能，单位为J；n为三维网格数值计算模型的网格单元数量；
Figure FDA0004228518440000068
为第j个网格单元的动能，单位为J；ρj为第j个网格单元的密度，等于土体密度ρ，即ρj＝ρ，单位为kg/m3；/>
Figure FDA0004228518440000069
为第j个网格单元的质点速度，单位为m/s；
所述步骤53中计算并提取最后一次折减计算收敛后三维网格数值计算模型的弹性剪切应变能，计算方法如下：
基于网格单元的第一主应力和第二主应力确定第三半圆，基于网格单元的第一主应力和第三主应力确定第一半圆，基于网格单元的第二主应力和第三主应力确定第二半圆；
假定
Figure FDA0004228518440000071
是土体的内摩擦角/>
Figure FDA0004228518440000072
第一半圆的切线在横轴上相交于点A，点A与第一半圆的切线在横轴上的交点夹角为/>
Figure FDA0004228518440000073
令第二半圆的切线在横轴上的交点位于A点，确定A与第二半圆的切线在横轴上的交点夹角为
Figure FDA0004228518440000074
令第三半圆的切线在横轴上的交点位于A点，确定A与第三半圆的切线在横轴上的交点夹角为
Figure FDA0004228518440000075
切线与半圆相交的点为复合滑动面；
步骤531、基于最后一次折减计算收敛后的模型，提取第j个网格单元的三个主应力σ1-j、σ2-j、σ3-j，计算第j个网格单元的Lode角正切值：
Figure FDA0004228518440000076
式中：tanθj为第j个网格单元的Lode角正切值；σ1-j、σ2-j、σ3-j分别为第j个网格单元的第一主应力、第二主应力和第三主应力，单位为Pa；
步骤532、计算第j个网格单元的三个复合滑动面的几何角度
Figure FDA0004228518440000077
和/>
Figure FDA0004228518440000078
的正弦值：
Figure FDA0004228518440000079
Figure FDA00042285184400000710
Figure FDA00042285184400000711
式中：
Figure FDA00042285184400000712
分别为第j个网格单元三个复合滑动面的几何角度
Figure FDA00042285184400000713
和/>
Figure FDA00042285184400000714
的正弦值；tanθj为第j个网格单元的Lode角正切值；/>
Figure FDA00042285184400000715
为第j个网格单元的内摩擦角，等于土体的内摩擦角/>
Figure FDA00042285184400000716
即/>
Figure FDA00042285184400000717
步骤533、计算第j个网格单元三个复合滑动面上的法向应力和剪切应力：
Figure FDA0004228518440000081
Figure FDA0004228518440000082
Figure FDA0004228518440000083
Figure FDA0004228518440000084
Figure FDA0004228518440000085
Figure FDA0004228518440000086
式中：σ12-j、τ12-j分别为第j个网格单元第一个复合滑动面的法向应力和剪切应力，单位为Pa；σ23-j、τ23-j分别为第j个网格单元第二个复合滑动面的法向应力和剪切应力，单位为Pa；σ13-j、τ13-j分别为第j个网格单元第三个复合滑动面的法向应力和剪切应力，单位为Pa；
步骤534、分别计算第j个网格单元三个复合滑动面的弹性剪切应变能及其之和：
Figure FDA0004228518440000087
Figure FDA0004228518440000088
Figure FDA0004228518440000089
Pj＝P12-j+P23-j+P13-j (20)
式中：Pj为第j个网格单元的总弹性剪切应变能，单位为J；P13-j、P12-j和P23-j分别为第j个网格单元中三个复合滑动面的弹性剪切应变能，单位为J；G1-j为第j个网格单元的切变模量，单位为Pa，等于土体的切变模量，即G1-j＝G1；
步骤535、逐一计算所有网格单元的总弹性剪切应变能并累计求和，得到三维网格数值计算模型总的弹性剪切应变能，计算公式如下：
Figure FDA0004228518440000091
式中：P为模型总的弹性剪切应变能，单位为J；n为模型的网格单元数量；Pj为第j个网格单元的总弹性剪切应变能，单位为J；
所述步骤6中基于最后一次计算收敛后模型的动能与弹性剪切应变能，计算“动剪比”系数up，其计算公式如下：
Figure FDA0004228518440000092
式中：up为“动剪比”系数；Uk为整个模型的动能，单位为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为标准确定合理的嵌入深度。
CN202210811745.0A 2022-07-11 2022-07-11 一种基于边坡稳定性的山区公路线位优化方法 Active CN115168953B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210811745.0A CN115168953B (zh) 2022-07-11 2022-07-11 一种基于边坡稳定性的山区公路线位优化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210811745.0A CN115168953B (zh) 2022-07-11 2022-07-11 一种基于边坡稳定性的山区公路线位优化方法

Publications (2)

Publication Number Publication Date
CN115168953A CN115168953A (zh) 2022-10-11
CN115168953B true CN115168953B (zh) 2023-06-27

Family

ID=83494092

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210811745.0A Active CN115168953B (zh) 2022-07-11 2022-07-11 一种基于边坡稳定性的山区公路线位优化方法

Country Status (1)

Country Link
CN (1) CN115168953B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103485353A (zh) * 2013-09-24 2014-01-01 昆明理工大学 基于全局最优化的边坡稳定性分析条分法
WO2020186507A1 (zh) * 2019-03-20 2020-09-24 东北大学 一种基于动态强度折减dda法的边坡稳定性分析系统
CN112257140A (zh) * 2020-09-16 2021-01-22 南京工业大学 一种海底边坡稳定性的安全系数计算方法
CN113536414A (zh) * 2021-06-09 2021-10-22 桂林理工大学 基于三维建模的岩质边坡稳定性分析方法、系统及介质
CN114353876A (zh) * 2022-01-07 2022-04-15 兰州大学 一种黄土公路边坡健康监测方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101090583B1 (ko) * 2009-09-30 2011-12-08 경남과학기술대학교 산학협력단 임도비탈면의 복원을 위한 식생기반재 돌망태의 안정성 분석 방법
CN108254782B (zh) * 2018-02-09 2019-11-05 中国地质大学（北京） 一种边坡地震破坏失稳概率的获取方法及系统
CN110263483B (zh) * 2019-07-02 2022-12-13 四川农业大学 基于flac3d软件的抗滑桩设计推力曲线计算方法
CN111324942B (zh) * 2019-12-27 2023-04-07 昆明理工大学 一种考虑滑面动力渐进破坏的地震边坡稳定性分析方法
CN111287225B (zh) * 2020-02-20 2021-05-07 中南大学 锚固型边坡的锚杆应力监测及重建边坡稳定性评价方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103485353A (zh) * 2013-09-24 2014-01-01 昆明理工大学 基于全局最优化的边坡稳定性分析条分法
WO2020186507A1 (zh) * 2019-03-20 2020-09-24 东北大学 一种基于动态强度折减dda法的边坡稳定性分析系统
CN112257140A (zh) * 2020-09-16 2021-01-22 南京工业大学 一种海底边坡稳定性的安全系数计算方法
CN113536414A (zh) * 2021-06-09 2021-10-22 桂林理工大学 基于三维建模的岩质边坡稳定性分析方法、系统及介质
CN114353876A (zh) * 2022-01-07 2022-04-15 兰州大学 一种黄土公路边坡健康监测方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
FLAC~(3D)在公路边坡工程中的应用;胡义生;张玉峰;;西部探矿工程(10);全文 *
基于FLAC 3D的复合土钉在软土基坑中的应用;杨保全;周磊;吴伟;;北华大学学报(自然科学版)(01);全文 *
基于FLAC~(3D)的边坡稳定安全系数的确定;万保安;殷妮芳;熊茂东;;四川建筑科学研究(06);全文 *
基于强度折减法的岩质边坡稳定性分析;严利娥;孙元习;;广西大学学报(自然科学版)(03);全文 *
高填方路基稳定性的FLAC-3D分析;秦敏;白群立;范文斌;朱璐璐;段凯;;中外建筑(06);全文 *

Also Published As

Publication number Publication date
CN115168953A (zh) 2022-10-11

Similar Documents

Publication Publication Date Title
CN106443783B (zh) 一种基于断层活动性的多期次裂缝定量预测方法
CN102494667A (zh) 一种表征地面沉降的方法
CN102129712A (zh) 基于多地层及三维的土石方数量的三角网模型构建方法
CN105404758B (zh) 一种基于有限单元法的固体连续介质变形的数值模拟方法
CN101009024A (zh) 一种滑坡灾害可视化的实现方法
CN114297864B (zh) 一种受陡缓倾角控制的碎裂松动岩体边坡稳定性分析方法
CN105354394A (zh) 一种基于三维可视化的拱坝坝肩边坡稳定判断方法
CN102521882A (zh) 基于离散高程和自适应混合加权得到海床地形数据的方法
CN113420515B (zh) 一种基于降雨数据的滑坡泥石流形成演化模拟方法
CN109446573B (zh) 一种构建多维路面仿真模型的方法
CN115906256A (zh) 一种水库滑坡涌浪数值模拟方法及系统
CN115168953B (zh) 一种基于边坡稳定性的山区公路线位优化方法
CN112184902A (zh) 一种面向越界开采识别的地下开采面反演方法
CN115238553A (zh) 一种地埋管线渗漏浸蚀的危险区域划分方法和系统
CN101702181B (zh) 基于最近距离优先的渠道土石方迁移方法
CN112948931B (zh) 新建和既有地铁隧道双线叠交工况下盾构施工合理夹角与净距的确定方法
CN112883614A (zh) 一种基于数值模拟的岩溶地层盾构隧道溶洞处理范围判断方法
CN109215124A (zh) 一种复杂地质条件下大型地下工程3d网格模型的构建方法
CN106407569A (zh) 一种厚松散层薄基岩条件下地表下沉值计算方法
CN106846481B (zh) 一种地质剖面图的生成方法
CN115392032A (zh) 一种gis-mpm无缝集成的动态三维地质模型构建方法
CN1797035A (zh) 三维复杂地质体重构的技术
Dong et al. Analysis of development process of landslide accumulation body based on discrete element method
Wang et al. Study on steep slope stability of coal mine under open-pit and underground mining
CN109492285B (zh) 一种可变形三维任意圆化凸多面体块体离散单元法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant