本发明属于道路工程，涉及一种基于路基土黏弹塑性力学模型的路基结构永久变形计算方法。
背景技术：
1、交通基础设施建设发展迅猛，道路运输的需求量增长、客货运车辆的吨位加重以及道路设计和建设标准不断提高，均对道路系统的服役水平提出了严格要求。路基作为路面的基础，直接承受着上部传递而来的大周数循环剪切作用而产生永久变形。由于路基永久变形产生，导致路基顶部与路面结构底部出现了一定深度的脱空区域，使得行车动荷载经由路面结构传递至路基时受阻，路面各结构层将产生弯拉应力与弯拉应变，进而导致路面结构出现疲劳开裂等多种病害，将大幅制约道路系统发挥应有的安全、舒适、耐久等功能。然而，我国路基工程建设在设计、施工及验收阶段均缺乏对路基永久变形的全面考虑，即在设计阶段以干湿/冻融折减系数的形式考虑了气候条件对路基土回弹模量的影响并以之作为路基设计指标，施工阶段以击实试验得到的路基土最佳含水率与压实度作为路基压实质量的控制指标，填筑完成后以路基顶面弯沉值作为验收指标，这一整套流程都没有涉及长期循环荷载下路基的永久变形。因此，准确计算长期循环荷载下的路基永久变形对于补充我国路基设计规范的技术要求、建造耐久性路基以及构建长寿命道路结构技术体系具有重要的科学价值和工程意义。
2、路基土永久变形计算模型的不同类型，决定了进行路基结构永久变形计算时选取的计算方法。目前，路基结构永久变形计算方法分为数值计算法和分层总和法两类。
3、数值计算法的基本思路是通过设置一定数量的增量步，将每次加载分成多个增量步予以分析。此类方法的核心是基于经典土力学的土体本构模型，可以准确模拟循环荷载下路基结构的变形行为，并能够改变参数赋值进而揭示不同因素对路基结构永久变形的影响规律。尽管有研究尝试了这一方式，但一致反映的事实是需消耗大量的计算资源。具体表现在以土样几何尺寸建立的计算模型，虽然网格单元数量较少，但也仅获得了几十次加载的变形结果。可想而知，若计算模型扩大到实际路基的大型尺寸，则基于经典土力学本构模型的数值计算法远不能满足千万级别加载次数下路基结构永久变形的分析需求。究其原因在于，数值计算法在追踪每一次循环荷载的应力路径时经典土力学本构模型必须记忆全部套叠面的形状和它们与初始屈服面间的相对位置，一旦本构方程的维数和计算模型的网格单元过多，需要消耗大量的计算资源和迭代时间，计算的效率和收敛性会大幅降低，且本构模型的部分参数很难通过试验标定。例如，文献“饱和黄土动力本构模型及地铁隧道周围土层变形分析[d].西安建筑科技大学博士学位论文,2015”曾尝试根据饱和黄土的动力本构模型计算西安地铁路基结构的永久变形，计算过程中因计算机死机和结果存储磁盘空间不足导致了计算被动停止。
4、尽管数值计算法已取得一定进展，但分层总和法仍然是确定路基结构永久变形最常用的方法。分层总和法的基本思路是将路基视为线弹性半无限空间体，根据弹性体系理论一次性计算得到交通荷载所致的路基内部动应力分布水平，然后将动力影响区内的路基分为若干薄层，最后利用事先建立的力学-经验预估模型计算各薄层的永久变形并累加得到整个路基结构的永久变形。然而，路基永久变形是路基在长期循环荷载下的塑性变形累积，可见分层总和法的基本理论存在基于路基线弹性假定得到的路基内应力场用于路基塑性变形计算这一本质矛盾。虽然有学者指出即使各个薄层的划分厚度较小，但也无法体现动应力沿路基深度方向逐渐减小的现象，即每个薄层内各点的应力水平在计算时被认为是相同的，于是采用“微元”的思想构建了路基内任一极小厚度下单元体的永久变形表达式，以此表达式为被积函数沿路基动力影响区的深度积分得到了整个路基的永久变形。尽管积分法对路基内动应力非线性分布的考虑更合理，但本质仍是分层总和法，仍存在弹性假定用于塑性变形计算的矛盾。此外，现有研究根据分层总和法的计算结果分析了路基永久变形随不同因素的演变规律，但对影响因素的设置并不全面，远没有全面考虑状态变量、荷载特征、路面性质对路基永久变形的影响，且建立的路基永久变形预估模型靠经验设置，缺乏合理的逻辑支撑。
技术实现思路
1、本发明实施例的目的在于提供一种基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，以解决现有用于路基永久变形计算的数值计算法的经典土力学本构模型参数难获取、计算效率低、消耗计算资源大、计算耗时长、计算难收敛的问题，分层总和法将路基弹性假定下的应力分布用于塑性变形计算、基本理论存在本质矛盾造成计算结果不准确的问题，以及路基永久变形预估模型考虑因素不全面、缺乏理论支撑、工程应用意义有限的问题。
2、本发明实施例所采用的技术方案是：基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，包括以下步骤：
3、确定循环加载次数n、状态变量、荷载特征、路面性质；
4、基于循环加载次数n、状态变量、荷载特征、路面性质，采用路基结构永久变形最大值计算公式，计算进行路基结构永久变形最大值。
5、进一步的，所述状态变量为路基湿度w，荷载特征为行车速度v与单轮胎压m，路面性质为面层厚度hpave、基层厚度hbase、底基层厚度hsubb；
6、采用式(14)所示的路基结构永久变形最大值计算公式，计算路基结构永久变形最大值：
7、
8、式中：εpmax为路基结构永久变形最大值，n为循环加载次数，omc为最佳含水率，w/omc为路基的湿化程度。
9、进一步的，所述状态变量为路基湿度w，荷载特征为行车速度v与单轮胎压m，路面性质为面层厚度hpave、基层厚度hbase、底基层厚度hsubb；
10、采用式(13)所示的路基结构永久变形最大值计算公式，计算路基结构永久变形最大值：
11、εpmax＝f(n)·f(w)·f(v)·f(m)·f(hpave)·f(hbase)·f(hsubb) (13)
12、其中，εpmax为路基结构永久变形最大值，f(n)为循环加载次数的影响函数，f(w)为路基湿度w的影响函数，f(v)为行车速度v的影响函数，f(m)为单轮胎压m的影响函数，f(hpave)为面层厚度hpave的影响函数，f(hbase)为基层厚度hbase的影响函数，f(hsubb)为底基层厚度hsubb的影响函数，按照下表计算：
13、
14、进一步的，所述状态变量为路基湿度w，荷载特征为行车速度v与单轮胎压m，路面性质为面层材料模量epave、基层材料模量ebase、底基层材料模量esubb、面层厚度hpave、基层厚度hbase、底基层厚度hsubb；
15、采用式(12)所示的路基结构永久变形最大值计算公式，计算路基结构永久变形最大值：
16、εpmax＝f(n)·f(w)·f(v)·f(m)·f(epave)·f(ebase)·f(esubb)·f(hpave)·f(hbase)·f(hsubb) (12)
17、其中，εpmax为路基结构永久变形最大值，f(n)为循环加载次数的影响函数，f(w)为路基湿度w的影响函数，f(v)为行车速度v的影响函数，f(m)为单轮胎压m的影响函数，f(epave)为面层材料模量epave的影响函数，f(ebase)为基层材料模量ebase的影响函数，f(esubb)为底基层材料模量esubb的影响函数，f(hpave)为面层厚度hpave的影响函数，f(hbase)为基层厚度hbase的影响函数，f(hsubb)为底基层厚度hsubb的影响函数，按照下表计算：
18、
19、进一步的，所述路基结构永久变形最大值计算公式，按照下述过程拟合：
20、步骤s1、确定计算工况与道路有限元仿真模型的尺寸；
21、步骤s2、基于路基土黏弹塑性力学模型，进行路基结构永久变形的有限元计算；
22、步骤s3、基于不同工况的路基结构永久变形的有限元计算结果，确定路基结构永久变形影响因素及对应的路基结构永久变形演化规律；
23、步骤s4、基于路基结构永久变形影响因素及对应的路基结构永久变形演化规律，拟合路基结构永久变形最大值计算公式。
24、进一步的，所述步骤s1中，设计计算工况时：
25、路面部分：将底基层厚度范围调整至15～40cm，基层材料模量的变化范围确定在7000～14000mpa、底基层材料模量的变化范围确定在5000～10000mpa、沥青混合料面层的材料模量变化范围确定在5000～13500mpa，将面层、基层、底基层视为线弹性体，面层、基层、底基层材料的泊松比分别取0.3、0.25、0.25，材料密度分别取2.35g/cm3、2.10g/cm3、2.10g/cm3；
26、路基部分：将计算所需的路基湿度变化范围确定在1.0omc～1.6omc，其中omc为最佳含水率；
27、荷载部分：将行车速度的变化范围确定在40～140km/h，将单轮接地压力的变化范围确定在0.7～1.4mpa，将循环加载次数设置为5×107次；
28、设计了下表所示的计算工况，各组工况考察的因素记为x，其他因素按“w-v-m-epave-ebase-esubb-hpave-hbase-hsubb”的顺序控制变量，其中，w为路基湿度、v为行车速度、m为单轮胎压、epave为面层材料模量、ebase为基层材料模量、esubb为底基层材料模量、hpave为面层厚度、hbase为基层厚度、hsubb为底基层厚度：
29、
30、
31、进一步的，所述步骤s2的路基土黏弹塑性力学模型为：
32、
33、
34、
35、
36、
37、式中：εp是路基土的永久变形，σmax是循环应力幅值，g1是黏弹性元件在三维情况下的黏弹性剪切模量，η1为黏弹性元件的黏滞系数，t是循环荷载周期，t1是每个循环周期内的荷载作用时间，t2是每个循环周期内的荷载间歇时间，t＝t1+t2，n是循环荷载作用次数，β是黏塑性元件的非整数阶阶数，c是土体的总黏聚力，是土体的总内摩擦角，π是圆周率，η2为黏塑性元件的黏滞系数，γ为gamma函数，w为实际含水率，omc为最佳含水率，k为压实度，σ3是围压，abs表示计算绝对值。
38、进一步的，所述步骤s2采用comsol的蠕变模块公式(1)～(5)开展路基结构永久变形有限元计算，以式(6)所示蠕变率方程作为计算控制方程：
39、
40、其中，表示式(1)对时间求导。
41、进一步的，采用comsol的蠕变模块公式(1)～(5)开展路基结构永久变形有限元计算时：
42、对于循环加载次数n，通过表达式floor(t/t)+1实现，其中t为蠕变总时间，t为循环荷载的周期，floor为comsol中的取整函数；
43、对于含水率w和压实度k，通过预定义变量的方式直接赋值调用；
44、对于土体的总黏聚力c和土体的总内摩擦角通过式(7)～(8)所示的土体的总黏聚力c和土体的总内摩擦角与含水率、压实度的函数关系写入comsol：
45、c＝32.18(12.38k-10.38)(-18.52w+4.57) (7)
46、
47、对于围压σ3，通过上覆应力换算；
48、对于循环荷载，通过表达式loadintensity*(pulse(x-a)+pulse(x-b))施加，其中loadintensity为荷载幅值，a、b对应为左右两侧双轮组中心点的位置；
49、对于荷载作用时间t1与荷载间歇时间t2，按照式(9)～(10)所示写入：
50、t1＝107.22[ln(e+v)]-3.52(h+1)0.36 (9)
51、t2＝0.50(h+1)-0.11v0.40 (10)
52、式中：v为车速，h为路基顶面往下的深度；
53、对于循环应力幅值σmax，将循环应力幅值转换为八面体剪应力，在弱形式偏微分方程有限元法中创建一个附加因变量来定义八面体剪应力并赋予初始值，具体是：
54、在弱形式偏微分方程模块中创建一个新的物理场v，并定义独立变量v1，v1即为八面体剪应力的值，定义v1的弱表达形式为-(v1-sqrt(abs(solid.ii2s)*2/3))*test(v1)，其中：test()表示试算函数，solid.ii2s为应力偏量第二不变量，test()和solid.ii2s均为comsol中的预定义变量，abs表示计算绝对值。
55、进一步的，所述步骤s3中确定得到：
56、路基顶面永久变形在横向具有中间大、两侧小的盆式分布特点，且路基顶面永久变形最大值所在位置与车辆荷载的中心位置一致；
57、所述步骤s4基于确定的路基顶面永久变形的盆式分布特点以及路基顶面永久变形最大值所在位置，拟合路基结构永久变形最大值计算公式。
58、本发明实施例的有益效果是：
59、1.将推导的路基土永久变形黏弹塑性力学模型视为总时间t＝nt的等效蠕变方程，将其对时间求导得到与路基土永久变形等效的路基土蠕变的蠕变率方程，将此蠕变率方程通过用户自定义接口写入仿真平台comsol的蠕变模块，通过弱形式偏微分方程实现循环应力幅值与路基永久变形的相互作用，实现了基于路基土黏弹塑性力学模型的路基结构永久变形有限元计算；计算能力方面，大大提高了计算效率，缩减了计算时间，约6个小时即可完成5000万次循环加载下的路基结构永久变形计算，解决了数值计算法计算耗时长、消耗资源大，且在有限的计算时间里获得的几十次加载结果远远不能满足千万级别加载次数下路基结构永久变形分析需求的计算耗时长、计算难收敛的问题；计算根基方面，因comsol的蠕变模块中内嵌了von-mises塑性屈服函数以实时更新每个增量步施加后的应力分布，使得循环荷载下塑性的路基结构永久变形时刻影响着路基内的应力分布，同时路基内的应力分布又时刻反作用于塑性的路基结构永久变形，即实现了路基永久变形与路基应力分布时刻耦合影响，解决了分层总和法计算结果不准确的问题；
60、2.根据基于路基土黏弹塑性力学模型的路基结构永久变形的有限元计算结果，确定了路基永久变形最大点的所在位置，明确了路基结构永久变形最大值计算公式建立时的理论逻辑，全面分析了加载次数、状态变量(路基湿度)、荷载特征(超载程度、行车速度)、路面性质(各结构层厚度、各结构层材料模量)对路基结构永久变形最大值的影响规律与影响程度，选取了对路基结构永久变形最大值影响显著的因素建立了多元计算公式，建立的路基结构永久变形最大值计算公式逻辑清晰、理论明确、考虑因素全面、适用范围更广，解决了现有路基结构永久变形预估模型更多靠经验设置，影响因素选择及描述方法缺乏理论支撑，考虑因素不全面，难以在工程设计时推广应用的限制。
1.基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，所述状态变量为路基湿度w，荷载特征为行车速度v与单轮胎压m，路面性质为面层厚度hpave、基层厚度hbase、底基层厚度hsubb；
3.根据权利要求1所述的基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，所述状态变量为路基湿度w，荷载特征为行车速度v与单轮胎压m，路面性质为面层厚度hpave、基层厚度hbase、底基层厚度hsubb；
4.根据权利要求1所述的基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，所述状态变量为路基湿度w，荷载特征为行车速度v与单轮胎压m，路面性质为面层材料模量epave、基层材料模量ebase、底基层材料模量esubb、面层厚度hpave、基层厚度hbase、底基层厚度hsubb；
5.基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，所述路基结构永久变形最大值计算公式，按照下述过程拟合：
6.根据权利要求5所述的基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，所述步骤s1中，设计计算工况时：
7.根据权利要求5所述的基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，所述步骤s2的路基土黏弹塑性力学模型为：
8.根据权利要求7所述的基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，所述步骤s2采用comsol的蠕变模块公式(1)～(5)开展路基结构永久变形有限元计算，以式(6)所示蠕变率方程作为计算控制方程：
9.根据权利要求7或8所述的基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，采用comsol的蠕变模块公式(1)～(5)开展路基结构永久变形有限元计算时：
10.根据权利要求5～8任一项所述的基于路基土黏弹塑性力学模型的路基结构永久变形计算方法，其特征在于，所述步骤s3中确定得到：
