一种岩土颗粒接触蠕变变形离散元数值模拟方法

本发明涉及岩土工程，具体是一种岩土颗粒接触蠕变变形离散元数值模拟方法。

背景技术：

1、现如今工程领域中大坝的建造主要还以土石坝与堆石坝为主，此类坝体由粗粒土作为填料，然而在土石坝的实际运行后期坝体会发生预料之外的沉降，且这种变形通常持续时间较久。引起坝体的长期变形的重要原因是其主要建筑填料，也就是岩土粗颗粒所产生的长期变形。

2、对于粗颗粒堆石料等无粘性土，其宏观力学性质主要受微观层次的颗粒连接和排列等因素的影响。近年来，学者们开始逐渐深入研究颗粒接触力学性质，以更全面地从微观视角解释粗颗粒土的宏观力学现象。这种类型的研究有望深入挖掘粗颗粒土的宏观物理规律，值得注意的是，目前仍然是初级阶段，研究成果有限，但这一领域的研究具有极高的潜在价值，以上对颗粒接触力学性质的研究以球颗粒和立方体颗粒作为主要研究对象，将颗粒接触的形式简化为点-点接触，通过颗粒接触试验来探究接触力学性质。但该领域尚未有关于颗粒接触蠕变的研究成果。

3、离散元软件pfc3d(particle flow code in three dimensions)是一款专门针对散体材料的性质进行研究的数值模拟软件，可以很好的契合颗粒接触试验，其中对于无粘性散粒材料的接触蠕变提供burgers模型，burgers模型作为maxwell与kelvin的组合，结合了两种模型的优点，可以很好的描述颗粒接触的每个蠕变阶段。但在计算过程中，burgers模型需要确定参数，在现有的技术上都是参数凑配的方法，并且都是输入一个固定参数值，缺乏科学严谨性。如cn202010075697.4中为了简化计算，使burgers模型中kelvin体模型的4个离散元非线性接触流变模型参数无穷大，忽略了burgers模型中的kelvin体，也就是说将原本的burgers模型的10个参数简化成了6个，这样带来的后果就是相当于直接使用maxwell模型来描述粗颗粒接触蠕变，准确性较低。

技术实现思路

1、针对现有技术的不足，本发明的目的在于，提供一种岩土颗粒接触蠕变变形离散元数值模拟方法，该方法基于法向颗粒接触蠕变试验和切向颗粒接触蠕变试验确定burgers模型的10个参数，可准确地模拟颗粒接触试验中颗粒接触变形的变形量。

2、为实现上述目的，本发明采用的技术方案是：

3、一种岩土颗粒接触蠕变变形离散元数值模拟方法，所述方法用于对粗颗粒之间的长时间接触蠕变进行模拟，包括以下步骤：

4、法向粗颗粒接触蠕变和切向粗颗粒接触蠕变的burgers模型分别为公式(1)和公式(6)，

5、

6、

7、其中，un为岩土颗粒的法向位移，fn是试样所受的法向荷载，kn2、ηn2是maxwell体的法向弹簧刚度和法向粘滞系数，kn1、ηn1是kelvin体的法向弹簧刚度和法向粘滞系数；角标n代表法向；角标s代表切向；t为时间；us为岩土颗粒的切向位移，fs是试样所受的切向荷载，ks2、ηs2是maxwell体的切向弹簧刚度和切向粘滞系数，ks1、ηs1是kelvin体的切向弹簧刚度和切向粘滞系数；

8、分别确定两种burgers模型中的四个参数；

9、第一步、进行法向颗粒接触蠕变试验，获得粗颗粒分别在不同法向荷载条件下的粗颗粒法向时间-位移曲线；一次法向颗粒接触蠕变试验能确定出一组kn2、ηn2、kn1、ηn1，这四个值均与该次试验中的法向荷载fn条件所对应；

10、在法向荷载fn已知的情况下，记录法向颗粒接触蠕变试验时，时间t＝0所对应的瞬时位移un|t＝0，根据un|t＝0和fn的值带入un|t＝0＝fn/kn2，获得当前法向荷载下的参数kn2的值；

11、在确定参数kn2的值后，获得法向荷载下的时间与法向位移数据，在公式(1)的基础上利用oringin软件进行函数拟合，拟合后带入法向荷载及对应的kn2，然后根据时间t的系数项反解出ηn2、kn1、ηn1；

12、获得所有不同法向荷载下对应的四个参数，分别以法向载荷fn为自变量，四个参数为因变量进行函数拟合，分别获得法向本构参数kn2、ηn2、kn1、ηn1与法向荷载fn关系的函数表达式；

13、第二步、进行切向颗粒接触蠕变试验，在切向颗粒接触蠕变试验中，先进行法向加载，待稳定后施加切向力，单次切颗粒接触蠕变试验能确定出一组ks2、ηs2、ks1、ηns1，这四个值均与该次试验中的法向荷载fn和切向荷载fs所对应；

14、在切向荷载fs已知的情况下，记录切向颗粒接触蠕变试验时，时间t＝0所对应的瞬时位移us|t＝0，根据us|t＝0和fs的值带入us|t＝0＝fs/ks2，获得当前切向荷载下的参数kn2的值；

15、在确定参数ks2的值后，获得法向荷载及切向荷载下时间与切向位移的数据，在公式(6)的基础上利用oringin软件进行函数拟合，拟合后带入切向荷载及对应的ks2，然后根据时间t的系数项反解出ηs2、ks1、ηs1；

16、获得所有不同法向荷载、不同切向荷载下对应的四个参数，以法向载荷fn为自变量，以ks2为因变量进行函数拟合，获得ks2与法向荷载fn的函数表达式；以法向载荷fn和切向荷载fs为自变量，以ks1为因变量进行函数拟合，获得ks1与法向荷载fn、切向荷载fs的函数表达式；以切向荷载fs为自变量，以ηs2、ηs1分别为因变量进行函数拟合，获得ηs2、ηs1分别与切向荷载fs的函数表达式；

17、第三步、将八个参数与荷载的函数表达式嵌入到pfc3d软件中的burgers模型中，将pfc3d软件中的burgers模型的各个参数改成上述的函数表达式；

18、当输入的法向荷载与切向荷载确定后，能够动态调整burgers模型中的参数，进而能获得任意荷载条件与时间下的颗粒接触蠕变位移。

19、对于试验试样分别进行不同湿化程度下的法向颗粒接触蠕变试验和切向颗粒接触蠕变试验；

20、获得不同湿化程度下，法向本构参数kn2、ηn2、kn1、ηn1与法向荷载fn关系的函数表达式，以及切向本构参数ks2与法向荷载fn关系的函数表达式，ks1与法向荷载fn、切向荷载fs的函数表达式；ηs2、ηs1分别与切向荷载fs的函数表达式；

21、以自然含水率工况下的函数表达式为基准，确定与不同湿化程度与自然状态下的倍数关系，此倍数关系记为湿化系数；

22、在pfc3d软件中增加湿化系数的输入口，使burgers模型的参数根据输入的荷载与工况进行动态调整，获得任意荷载条件及工况与时间下的颗粒接触蠕变位移。

23、在pfc3d软件中burgers模型包括十个参数：

24、kelvin体的法向弹簧刚度kn1，对应软件中的bur_knk

25、maxwell体的法向弹簧刚度kn2，对应软件中的bur_cnk

26、kelvin体的法向粘滞系数ηn1，对应软件中的bur_knm

27、maxwell体的法向粘滞系数ηn2，对应软件中的bur_knm

28、kelvin体的切向弹簧刚度ks1，对应软件中的bur_ksk

29、maxwell体的切向弹簧刚度ks2，对应软件中的bur_csk

30、kelvin体的切向粘滞系数ηs1，对应软件中的bur_ksm

31、maxwell体的切向粘滞系数ηs2，对应软件中的bur_csm

32、利用fish语言将上述八个参数的函数表达式记录并输入pfc软件的burgers模型中对颗粒接触试验的变形进行计算；

33、摩擦系数fric取极值；

34、拉力开关bur_notension，设置为1。

35、试验对象为粒径5cm的立方体和锥体岩土颗粒，以点-面的形式进行接触。

36、与现有技术相比，本发明有益效果在于：

37、本发明针对无粘性土、散粒颗粒进行接触蠕变研究，颗粒之间不会承受拉力，分别进行了法向加载和切向加载条件下粗颗粒接触蠕变试验，选用burgers模型来描述拟合蠕变曲线，通过蠕变试验获得参数(获得法向参数与法向应力的关系，切向参数与法向应力和切向应力的关系)，利用fish语言对pfc3d软件中的burgers模型进行二次开发，将各个参数改成相应的函数关系式。更进一步地，获得湿化程度与自然含水率工况下参数的关系，使burgers模型的参数能根据输入的荷载与工况(可以以1表示自然含水率状态，以不同的数字表示不同工况)在计算过程中动态调整，而不是一个固定的值，以便计算任意荷载条件与时间下的颗粒接触蠕变位移。

38、本发明与发明人的在先申请cn202010075697.4burger模型的构建方法及其应用相比，cn202010075697.4专利中为了简化计算，使burgers模型中kelvin体模型的4个离散元非线性接触流变模型参数无穷大，忽略了burgers模型中的kelvin体，也就是说将原本的burgers模型的10个参数简化成了6个，这样带来的后果就是相当于直接使用maxwell模型来描述粗颗粒接触蠕变，严格来说该专利并不是burgers模型参数的确定方法和应用计算，而是maxwell模型参数的确定方法和应用计算。而本专利经过大量和全面的试验，严格按照pfc数值计算理论中的burgers给定的10个参数来对粗颗粒接触蠕变试验进行模拟计算，以离散元软件pfc3d建立颗粒接触的数值模型来进行数值模拟，研究在不同应力状态下的力学特性。因此，本技术突破了室内试验的局限，实现对实际工程中粗粒料的流变变形的量化预测，采用pfc3d离散元软件对粗粒料颗粒接触蠕变的变形进行计算。