土石混合体混合表征与并发式求解方法

本发明涉及岩土工程，尤其是一种土石混合体混合表征与并发式求解方法。

背景技术：

1、土石混合体是岩土工程中一种广泛存在的材料，它的非均质性和非线性导致了土石混合体具有非常复杂的力学性质。准确的再现和预测土石混合体的力学性质对于泥石流、碎屑流、滑坡等自然灾害的预警具有重要意义。通常，常见的求解土石混合体问题主要有两种方法，分别为连续类方法、非连续类方法。在连续类方法中，如有限元、有限差分、光滑粒子流体动力学等，土石混合体被当作连续介质进行计算，其力学行为通过本构模型进行描述，这可以很大程度上简化问题并减小计算量，但是土体与石块之间的相互作用不能纳入考虑，一些微观机制并不能再现；在非连续类方法中，如离散元，每个土体颗粒和石块都可以被很具体的建模出来，这可以在一定程度上考虑微观机制，但是面对大型工程问题，计算量会大大增加。另一方面，一些微观参数是难以标定的，这也是制约离散元法在实际工程中应用的重要因素之一。

技术实现思路

1、本发明的主要目的是提出一种土石混合体混合表征与并发式求解方法，旨在解决现有的土石混合体求解方法不能再现微观机制、计算量较大以及通用度较低的问题。

2、为实现上述目的，本发明提出一种土石混合体混合表征与并发式求解方法，所述土石混合体混合表征与并发式求解方法包括：

3、在确定所述目标土石混合体的土体和石块后，对所述土体以物质点法进行建模，并生成与所述土体对应的物质点，对所述石块以离散元法进行建模，并生成与所述石块对应的离散单元；

4、对所述土体和所述石块之间的作用关系以预设接触模型进行描述，并构造所述目标土石混合体的mpm-dem表征方案；

5、计算各所述物质点的虚半径；

6、根据所述虚半径，确定与所述离散单元有接触的物质点，并标记为虚粒子；

7、采用离散元接触模型计算所述虚粒子与所述离散单元的耦合接触力；

8、根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度；

9、根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时速度和位置等参数。

10、可选地，计算各所述物质点的虚半径，包括：

11、采用公式计算各所述物质点的虚半径，其中，为第p个物质点的虚半径，为第p个物质点的体积，kp为第p个物质点的孔隙率。

12、可选地，采用离散元法计算所述虚粒子与所述离散单元的耦合接触力，包括：

13、采用公式计算法向合力，其中，为第p个离散单元受到的来自于对应物质点的法向合力，为第q个物质点作用于第p个离散单元上的法向力，pcontact为与第p个离散单元相接触的粒子组成的邻域列表；

14、采用公式计算切向合力，其中，为第p个离散单元受到的来自于对应物质点的切向合力，为第q个物质点作用于第p个离散单元上的切向力，pcontact为与第p个离散单元相接触的粒子组成的邻域列表；

15、将所述法向合力和所述切向合力合成所述耦合接触力。

16、可选地，根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，包括：

17、构造物质点法平衡方程其中，为物质点法中第i个节点上的合内力，为物质点法中第i个节点上的合外力，为物质点法中第i个节点的质量，为物质点法中第i个节点的加速度；

18、求解所述物质点法平衡方程，获得物质点法中各背景网格节点上的加速度。

19、可选地，求解所述物质点法平衡方程，获得物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，还包括：

20、采用公式计算所述合内力，其中，nip为形函数，为外界施加到物质点上的体力，为物质点法中第i个节点上所受到的面力，为第p个离散单元受到的来自于对应物质点的法向合力，为第p个离散单元受到的来自于对应物质点的切向合力；

21、采用公式计算所述合外力，其中，为第p个物质点上的应力，为第p个物质点的体积，bip为形函数梯度；

22、根据所述合内力和所述合外力计算各背景网格节点上的加速度。

23、可选地，根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，还包括：

24、构造离散元法平衡方程以及其中，为第p个离散单元上受到的来自其他离散单元的法向合力，为第p个离散单元上受到的来自其他离散单元的切向合力，为第p个离散单元上受到的体力，为第p个离散单元上受到的转动阻尼力，为第p个离散单元的质量，为第p个离散单元的加速度，ip为第p个离散单元的惯性矩，为第p个离散单元的角加速度，为第p个离散单元的粒径，npq为接触点的单位法向量，δn,pq为p个离散单元与第q个离散单元(或物质点)之间的法向堆叠量，为第p个离散单元上所受到的来自第q个离散单元的力；

25、求解所述离散元法平衡方程，获得各所述离散单元上的加速度和角加速度。

26、可选地，所述物质点的实时参数包括物质点速度，所述离散单元的实时参数包括离散单元速度；

27、根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时参数，包括：

28、采用flip速度更新格式对所述物质点速度进行实时更新，其中，所述flip速度更新格式为以及为第p个物质点的相邻节点所组成的邻域列表，δt为物质点的时间步长，为第p个物质点在tn+1/2时刻的速度，为第p个物质点在tn-1/2时刻的速度，为第p个物质点在tn+1/2时刻的速度，为第p个物质点在tn-1/2时刻的速度；

29、采用预设速度更新格式对所述离散单元速度进行实时更新，其中，所述预设速度更新格式为以及为第p个物质点在tn+1/2时刻的速度，为第p个物质点在tn-1/2时刻的速度，为第p个物质点在tn+1/2时刻的加速度，为第p个物质点在tn-1/2时刻的加速度。

30、可选地，所述物质点的实时参数包括物质点位置，所述离散单元的实时参数包括离散单元位置；

31、根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时参数，包括：

32、采用更新格式对所述物质点位置进行实时更新，其中，为第p个物质点在tn+1时刻的物质点坐标，为第p个物质点在tn-1时刻的物质点坐标；

33、采用预设位置更新格式对所述离散单元位置进行实时更新，其中，所述预设位置更新格式为以及为第p个离散单元在tn+1时刻的离散单元坐标，为第p个离散单元在tn-1时刻的离散单元坐标，为第p个离散单元在tn+1时刻的离散单元坐标，为第p个离散单元在tn-1时刻的离散单元坐标。

34、可选地，根据求解结果，以预设更新格式分别更新所述物质点的应力、速度和位置、以及所述离散单元的速度和位置，还包括：

35、将物质点的动量参数映射至背景网格上；

36、根据所述动量参数，以预设方程求解网格节点上的速度，其中，所述预设方程为以及

37、以应变更新格式更新物质点上的应变，其中，所述应变更新格式设置为

38、以旋量更新格式更新物质点上的旋量，其中，所述旋量更新格式设置为

39、根据所述应变和所述旋量，采用预设本构模型更新物质点上的应力。

40、可选地，在确定所述目标土石混合体的土体和石块后，对所述土体以物质点法进行建模，获取与所述土体对应的物质点，对所述石块以离散元法进行建模，获取与所述石块对应的离散单元的步骤之前，还包括：

41、获取目标土石混合体的分界粒径；

42、根据所述分界粒径，确定所述目标土石混合体的土体和石块。

43、本发明的技术方案中，小于分界粒径的颗粒视为土体，采用物质点法(mpm)建模，并采用宏观本构描述其力学行为，可以适当提高物质点的粒径，以减少总粒子数，节省计算资源。对于大于分界粒径的颗粒视为石块，采用离散元法(dem)建模，以便追踪单个石块的运动和接触行为；将每个物质点赋予虚半径，用以进行接触检测，将与dem单元相互接触的物质点当作虚粒子处理，从而计算耦合接触力；将耦合接触力当作外力边界条件，分别施加到物质点和离散元各自的部分，并构建二者各自的平衡方程；求解物质点与离散元的平衡方程，并采用musl应力更新格式和flip速度更新格式分别更新物质点的应力和速度以及物质点的位置；而对于离散元，其速度和位置可以直接更新；需要说明的是，物质点法是一种结合了拉格朗日描述和欧拉描述的数值方法。相比于基于网格的数值方法，物质点法可以克服网格畸变等问题；相比于其他无网格法，物质点法也有着优秀的计算效率。离散元法在处理多体接触等方面有着明显的优势，非常适合处理土石混合体中的石块以及土体和石块之间的接触；而本发明中将物质点法和离散元法耦合，采用物质点法模拟土体，采用离散元法模拟石块，充分发挥二者的优势。鉴于物质点法的连续特性，可以采用宏观本构来描述土体的力学响应，而不必在乎单个土颗粒之间的接触行为。在保证整体力学性质不变的情况下，可以适当加大物质点的尺寸，以减少模型中的总粒子数，提高计算效率；而将每个石块当作一个单独的dem单元，不再进行精细剖分，也可减少模型中的总粒子数，提高计算效率。同时，石块与石块之间，石块与土体之间的接触得以保留，也可以在一定程度上反映微观机制。最后，本发明中所给出的物质点-离散元耦合框架中，物质点法与离散元法的计算可以并发进行，物质点法与离散元法之间除了计算耦合接触力之外，其余步骤互不影响，简化了代码实现难度，并且提高了计算效率。