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

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


背景技术：

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


技术实现要素：

3.本发明的主要目的是提出一种土石混合体混合表征与并发式求解方法，旨在解决现有的土石混合体求解方法不能再现微观机制、计算量较大以及通用度较低的问题。
4.为实现上述目的，本发明提出一种土石混合体混合表征与并发式求解方法，所述土石混合体混合表征与并发式求解方法包括：
5.在确定所述目标土石混合体的土体和石块后，对所述土体以物质点法进行建模，并生成与所述土体对应的物质点，对所述石块以离散元法进行建模，并生成与所述石块对应的离散单元；
6.对所述土体和所述石块之间的作用关系以预设接触模型进行描述，并构造所述目标土石混合体的mpm-dem表征方案；
7.计算各所述物质点的虚半径；
8.根据所述虚半径，确定与所述离散单元有接触的物质点，并标记为虚粒子；
9.采用离散元接触模型计算所述虚粒子与所述离散单元的耦合接触力；
10.根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度；
11.根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时速度和位置等参数。
12.可选地，计算各所述物质点的虚半径，包括：
13.采用公式计算各所述物质点的虚半径，其中，
为第p个物质点的虚半径，为第p个物质点的体积，k
p
为第p个物质点的孔隙率。
14.可选地，采用离散元法计算所述虚粒子与所述离散单元的耦合接触力，包括：
15.采用公式计算法向合力，其中，为第p个离散单元受到的来自于对应物质点的法向合力，为第q个物质点作用于第p个离散单元上的法向力，pcontact为与第p个离散单元相接触的粒子组成的邻域列表；
16.采用公式计算切向合力，其中，为第p个离散单元受到的来自于对应物质点的切向合力，为第q个物质点作用于第p个离散单元上的切向力，pcontact为与第p个离散单元相接触的粒子组成的邻域列表；
17.将所述法向合力和所述切向合力合成所述耦合接触力。
18.可选地，根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，包括：
19.构造物质点法平衡方程其中，为物质点法中第i个节点上的合内力，为物质点法中第i个节点上的合外力，为物质点法中第i个节点的质量，为物质点法中第i个节点的加速度；
20.求解所述物质点法平衡方程，获得物质点法中各背景网格节点上的加速度。
21.可选地，求解所述物质点法平衡方程，获得物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，还包括：
22.采用公式计算所述合内力，其中，n
ip
为形函数，为外界施加到物质点上的体力，为物质点法中第i个节点上所受到的面力，为第p个离散单元受到的来自于对应物质点的法向合力，为第p个离散单元受到的来自于对应物质点的切向合力；
23.采用公式计算所述合外力，其中，为第p个物质点上的应力，为第p个物质点的体积，b
ip
为形函数梯度；
24.根据所述合内力和所述合外力计算各背景网格节点上的加速度。
25.可选地，根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，还包括：
26.构造离散元法平衡方程以及其中，为第p个离散单元上受到的来自其他离散单元的法向合力，为第p个离散单元上受到的来自其他离散单元的切向合力，为第p个离散单元上受到的体力，为第
p个离散单元上受到的转动阻尼力，为第p个离散单元的质量，为第p个离散单元的加速度，i
p
为第p个离散单元的惯性矩，为第p个离散单元的角加速度，为第p个离散单元的粒径，n
pq
为接触点的单位法向量，δ
n,pq
为p个离散单元与第q个离散单元(或物质点)之间的法向堆叠量，为第p个离散单元上所受到的来自第q个离散单元的力；
27.求解所述离散元法平衡方程，获得各所述离散单元上的加速度和角加速度。
28.可选地，所述物质点的实时参数包括物质点速度，所述离散单元的实时参数包括离散单元速度；
29.根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时参数，包括：
30.采用flip速度更新格式对所述物质点速度进行实时更新，其中，所述flip速度更新格式为以及以及为第p个物质点的相邻节点所组成的邻域列表，δt为物质点的时间步长，为第p个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度，为第p个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度；
31.采用预设速度更新格式对所述离散单元速度进行实时更新，其中，所述预设速度更新格式为以及以及为第p个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度，为第p个物质点在t
n+1/2
时刻的加速度，为第p个物质点在t
n-1/2
时刻的加速度。
32.可选地，所述物质点的实时参数包括物质点位置，所述离散单元的实时参数包括离散单元位置；
33.根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时参数，包括：
34.采用更新格式对所述物质点位置进行实时更新，其中，为第p个物质点在t
n+1
时刻的物质点坐标，为第p个物质点在t
n-1
时刻的物质点坐标；
35.采用预设位置更新格式对所述离散单元位置进行实时更新，其中，所述预设位置更新格式为以及以及为第p个离散单元在t
n+1
时刻的离散单元坐标，为第p个离散单元在t
n-1
时刻的离散单元坐标，为第p个离散单元在t
n+1
时刻的离散单元坐标，为第p个离散单元在t
n-1
时刻的离散单元坐标。
36.可选地，根据求解结果，以预设更新格式分别更新所述物质点的应力、速度和位置、以及所述离散单元的速度和位置，还包括：
37.将物质点的动量参数映射至背景网格上；
38.根据所述动量参数，以预设方程求解网格节点上的速度，其中，所述预设方程为以及
39.以应变更新格式更新物质点上的应变，其中，所述应变更新格式设置为
40.以旋量更新格式更新物质点上的旋量，其中，所述旋量更新格式设置为
41.根据所述应变和所述旋量，采用预设本构模型更新物质点上的应力。
42.可选地，在确定所述目标土石混合体的土体和石块后，对所述土体以物质点法进行建模，获取与所述土体对应的物质点，对所述石块以离散元法进行建模，获取与所述石块对应的离散单元的步骤之前，还包括：
43.获取目标土石混合体的分界粒径；
44.根据所述分界粒径，确定所述目标土石混合体的土体和石块。
45.本发明的技术方案中，小于分界粒径的颗粒视为土体，采用物质点法(mpm)建模，并采用宏观本构描述其力学行为，可以适当提高物质点的粒径，以减少总粒子数，节省计算资源。对于大于分界粒径的颗粒视为石块，采用离散元法(dem)建模，以便追踪单个石块的运动和接触行为；将每个物质点赋予虚半径，用以进行接触检测，将与dem单元相互接触的物质点当作虚粒子处理，从而计算耦合接触力；将耦合接触力当作外力边界条件，分别施加到物质点和离散元各自的部分，并构建二者各自的平衡方程；求解物质点与离散元的平衡方程，并采用musl应力更新格式和flip速度更新格式分别更新物质点的应力和速度以及物质点的位置；而对于离散元，其速度和位置可以直接更新；需要说明的是，物质点法是一种结合了拉格朗日描述和欧拉描述的数值方法。相比于基于网格的数值方法，物质点法可以克服网格畸变等问题；相比于其他无网格法，物质点法也有着优秀的计算效率。离散元法在处理多体接触等方面有着明显的优势，非常适合处理土石混合体中的石块以及土体和石块之间的接触；而本发明中将物质点法和离散元法耦合，采用物质点法模拟土体，采用离散元法模拟石块，充分发挥二者的优势。鉴于物质点法的连续特性，可以采用宏观本构来描述土体的力学响应，而不必在乎单个土颗粒之间的接触行为。在保证整体力学性质不变的情况下，可以适当加大物质点的尺寸，以减少模型中的总粒子数，提高计算效率；而将每个石块当作一个单独的dem单元，不再进行精细剖分，也可减少模型中的总粒子数，提高计算效率。同时，石块与石块之间，石块与土体之间的接触得以保留，也可以在一定程度上反映微观机制。最后，本发明中所给出的物质点-离散元耦合框架中，物质点法与离散元法的计算可以并发进行，物质点法与离散元法之间除了计算耦合接触力之外，其余步骤互不影响，简化了代码实现难度，并且提高了计算效率。
附图说明
46.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
47.图1为本发明提供的土石混合体混合表征与并发式求解方法第一实施例的流程图；
48.图2为本发明提供的土石混合体的mpm-dem表征方案；
49.图3为本发明提供的基于虚粒子接触检测技术的接触力耦合方案；
50.图4为本发明提供的二元混合体的级配曲线；
51.图5为本发明提供的粒径比为1：10的二元混合体试样的建模图；
52.图6为本发明提供的粒径比为1：10的二元混合体试样计算结果与离散元法对比的曲线图；
53.图7为本发明提供的粒径比为1：5的二元混合体试样的建模图；
54.图8为本发明提供的粒径比为1：5的试样计算结果与粒径比1：10的结果对比的曲线图；
55.图9为本发明提供的室内中型三轴试验的试样的建模图；
56.图10为本发明提供的200kpa围压下室内中型三轴试验模拟与实际结果对比的曲线图；
57.图11为本发明提供的400kpa围压下室内中型三轴试验模拟与实际结果对比的曲线图；
58.图12为本发明提供的800kpa围压下室内中型三轴试验模拟与实际结果对比的曲线图。
59.本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
60.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
61.需要说明，若本发明实施例中有涉及方向性指示(诸如上、下、左、右、前、后
……
)，则该方向性指示仅用于解释在某一特定姿态(如附图所示)下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。
62.另外，若本发明实施例中有涉及“第一”、“第二”等的描述，则该“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。另外，全文中出现的“和/或”的含义，包括三个并列的方案，以“a和/或b”为例，包括a方案、或b方案、或a和b同时满足的方案。另外，各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法
实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
63.土石混合体是岩土工程中一种广泛存在的材料，它的非均质性和非线性导致了土石混合体具有非常复杂的力学性质。准确的再现和预测土石混合体的力学性质对于泥石流、碎屑流、滑坡等自然灾害的预警具有重要意义。通常，常见的求解土石混合体问题主要有两种方法，分别为连续类方法、非连续类方法。在连续类方法中，如有限元、有限差分、光滑粒子流体动力学等，土石混合体被当作连续介质进行计算，其力学行为通过本构模型进行描述，这可以很大程度上简化问题并减小计算量，但是土体与石块之间的相互作用不能纳入考虑，一些微观机制并不能再现；在非连续类方法中，如离散元，每个土体颗粒和石块都可以被很具体的建模出来，这可以在一定程度上考虑微观机制，但是面对大型工程问题，计算量会大大增加。另一方面，一些微观参数是难以标定的，这也是制约离散元法在实际工程中应用的重要因素之一。
64.鉴于此，本发明提出一种土石混合体混合表征与并发式求解方法。图1至图12为土石混合体混合表征与并发式求解方法的具体实施例。
65.请参阅图1，所述土石混合体混合表征与并发式求解方法包括：
66.s30：在确定所述目标土石混合体的土体和石块后，对所述土体以物质点法进行建模，并生成与所述土体对应的物质点，对所述石块以离散元法进行建模，并生成与所述石块对应的离散单元；
67.s40：对所述土体和所述石块之间的作用关系以预设接触模型进行描述，并构造所述目标土石混合体的mpm-dem表征方案；
68.s50：计算各所述物质点的虚半径；
69.s60：根据所述虚半径，确定与所述离散单元有接触的物质点，并标记为虚粒子；
70.s70：采用离散元接触模型计算所述虚粒子与所述离散单元的耦合接触力；
71.s80：根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度；
72.s90：根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时速度和位置等参数。
73.本发明的技术方案中，小于分界粒径的颗粒视为土体，采用物质点法(mpm)建模，并采用宏观本构描述其力学行为，可以适当提高物质点的粒径，以减少总粒子数，节省计算资源。对于大于分界粒径的颗粒视为石块，采用离散元法(dem)建模，以便追踪单个石块的运动和接触行为；将每个物质点赋予虚半径，用以进行接触检测，将与dem单元相互接触的物质点当作虚粒子处理，从而计算耦合接触力；将耦合接触力当作外力边界条件，分别施加到物质点和离散元各自的部分，并构建二者各自的平衡方程；求解物质点与离散元的平衡方程，并采用musl应力更新格式和flip速度更新格式分别更新物质点的应力和速度以及物质点的位置；而对于离散元，其速度和位置可以直接更新；需要说明的是，物质点法是一种结合了拉格朗日描述和欧拉描述的数值方法。相比于基于网格的数值方法，物质点法可以克服网格畸变等问题；相比于其他无网格法，物质点法也有着优秀的计算效率。离散元法在处理多体接触等方面有着明显的优势，非常适合处理土石混合体中的石块以及土体和石块之间的接触；而本发明中将物质点法和离散元法耦合，采用物质点法模拟土体，采用离散元法模拟石块，充分发挥二者的优势。鉴于物质点法的连续特性，可以采用宏观本构来描述土
体的力学响应，而不必在乎单个土颗粒之间的接触行为。在保证整体力学性质不变的情况下，可以适当加大物质点的尺寸，以减少模型中的总粒子数，提高计算效率；而将每个石块当作一个单独的dem单元，不再进行精细剖分，也可减少模型中的总粒子数，提高计算效率。同时，石块与石块之间，石块与土体之间的接触得以保留，也可以在一定程度上反映微观机制。最后，本发明中所给出的物质点-离散元耦合框架中，物质点法与离散元法的计算可以并发进行，物质点法与离散元法之间除了计算耦合接触力之外，其余步骤互不影响，简化了代码实现难度，并且提高了计算效率。
74.具体地，计算各所述物质点的虚半径s50，包括：
75.s501：采用公式计算各所述物质点的虚半径，其中，为第p个物质点的虚半径，为第p个物质点的体积，k
p
为第p个物质点的孔隙率。
76.在一实施例中，采用离散元法计算所述虚粒子与所述离散单元的耦合接触力s70，包括：
77.s701：采用公式计算法向合力，其中，为第p个离散单元受到的来自于对应物质点的法向合力，为第q个物质点作用于第p个离散单元上的法向力，pcontact为与第p个离散单元相接触的粒子组成的邻域列表；
78.s702：采用公式计算切向合力，其中，为第p个离散单元受到的来自于对应物质点的切向合力，为第q个物质点作用于第p个离散单元上的切向力，pcontact为与第p个离散单元相接触的粒子组成的邻域列表；
79.s703：将所述法向合力和所述切向合力合成所述耦合接触力。
80.在本实施例中，物质点法与离散元法仅在计算耦合接触力时进行交互，从而进行接触检测和接触力的求解。而二者各自的计算流程是相对独立的，因此可以用并发式的求解方法，即让物质点和离散元各自的计算各自分开进行，只在耦合接触力求解时强制同步。
81.同时，根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度s80，包括：
82.s801：构造物质点法平衡方程其中，为物质点法中第i个节点上的合内力，为物质点法中第i个节点上的合外力，为物质点法中第i个节点的质量，为物质点法中第i个节点的加速度；
83.s802：求解所述物质点法平衡方程，获得物质点法中各背景网格节点上的加速度。
84.在本实施例中，物质点用于积分的体积和用于计算接触力的体积是不同的，用于积分的体积是物质点的固体部分和孔隙部分的总体积，而用于计算接触力的体积是物质点中固体部分的体积。
85.具体地，求解所述物质点法平衡方程，获得物质点法中各背景网格节点上的加速度s802，还包括：
86.s8021：采用公式计算所述合内力，其中，n
ip
为形函数，为外界施加到物质点上的体力，为物质点法中第i个节点上所受到的面力，为第p个离散单元受到的来自于对应物质点的法向合力，为第p个离散单元受到的来自于对应物质点的切向合力；
87.s8022：采用公式计算所述合外力，其中，为第p个物质点上的应力，为第p个物质点的体积，b
ip
为形函数梯度；
88.s8023：根据所述合内力和所述合外力计算各背景网格节点上的加速度。
89.同时，根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度s80，还包括：
90.s801
′
：构造离散元法平衡方程以及其中，为第p个离散单元上受到的来自其他离散单元的法向合力，为第p个离散单元上受到的来自其他离散单元的切向合力，为第p个离散单元上受到的体力，为第p个离散单元上受到的转动阻尼力，为第p个离散单元的质量，为第p个离散单元的加速度，i
p
为第p个离散单元的惯性矩，为第p个离散单元的角加速度，为第p个离散单元的粒径，n
pq
为接触点的单位法向量，δ
n,pq
为p个离散单元与第q个离散单元(或物质点)之间的法向堆叠量，为第p个离散单元上所受到的来自第q个离散单元的力；
91.s802
′
：求解所述物质点法平衡方程，获得各所述离散单元上的加速度和角加速度。
92.在本发明中，所述物质点的实时参数包括物质点速度，所述离散单元的实时参数包括离散单元速度；
93.根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时参数s90，包括：
94.s901：采用flip速度更新格式对所述物质点速度进行实时更新，其中，所述flip速度更新格式为以及以及为第p个物质点的相邻节点所组成的邻域列表，δt为物质点的时间步长，为第p个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度，为第p个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度；
95.s902：采用预设速度更新格式对所述离散单元速度进行实时更新，其中，所述预设速度更新格式为以及以及为第p
个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度，为第p个物质点在t
n+1/2
时刻的加速度，为第p个物质点在t
n-1/2
时刻的加速度。
96.同时，所述物质点的实时参数包括物质点位置，所述离散单元的实时参数包括离散单元位置；
97.根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时参数s90，包括：
98.s901
′
：采用更新格式对所述物质点位置进行实时更新，其中，为第p个物质点在t
n+1
时刻的物质点坐标，为第p个物质点在t
n-1
时刻的物质点坐标；
99.s902
′
：采用预设位置更新格式对所述离散单元位置进行实时更新，其中，所述预设位置更新格式为以及以及为第p个离散单元在t
n+1
时刻的离散单元坐标，为第p个离散单元在t
n-1
时刻的离散单元坐标，为第p个离散单元在t
n+1
时刻的离散单元坐标，为第p个离散单元在t
n-1
时刻的离散单元坐标。
100.并且，根据求解结果，以预设更新格式分别更新所述物质点的应力、速度和位置、以及所述离散单元的速度和位置s90，还包括：
101.s901
″
：将物质点的动量参数映射至背景网格上；
102.s902
″
：根据所述动量参数，以预设方程求解网格节点上的速度，其中，所述预设方程为以及
103.s903
″
：以应变更新格式更新物质点上的应变，其中，所述应变更新格式设置为
104.s904
″
：以旋量更新格式更新物质点上的旋量，其中，所述旋量更新格式设置为
105.s905
″
：根据所述应变和所述旋量，采用预设本构模型更新物质点上的应力。
106.在本发明中，在确定所述目标土石混合体的土体和石块后，对所述土体以物质点法进行建模，获取与所述土体对应的物质点，对所述石块以离散元法进行建模，获取与所述石块对应的离散单元的步骤s30之前，还包括：
107.s10：获取目标土石混合体的分界粒径；
108.s20：根据所述分界粒径，确定所述目标土石混合体的土体和石块。
109.需要说明的是，本发明所给出的验证算例中，主要使用了基于drucker-prager屈服准则的理想弹塑性本构模型。
110.以下以具体实验进行验证：
111.实施例1：
112.将实验样品简化，即同一体系中只有两种粒径——粗颗粒与细颗粒，这种情况也被称为二元混合体。其级配曲线如图4所示，粒径比为1：10。其试样如图5所示，试样为立方体试样，边长为0.075m。
113.首先，采用纯细颗粒试样分别在50kpa、100kpa、150kpa三种围压下进行不排水三轴试验，以标定纯细颗粒的宏观参数，用于物质点-离散元耦合方法的模拟。最终模拟采用的参数见表1。
114.表1纯细颗粒的宏观参数
[0115][0116]
为了模拟出三轴试验的应变硬化和软化现象，引入了三次硬化函数和指数型软化函数，如下公式：
[0117][0118]
公式中：ε
p
是当前时刻的累计等效塑性应变；到达峰值时的等效塑性应变，这里取0.022；和是残余和峰值内摩擦角，是弹性段的内摩擦角，此处假设它和相等；a，b，c，d是硬化函数的系数，这里分别取1075.64、-407.93、16.39、0.32；h是与软化速率相关的系数，这里取50.0。
[0119]
最终，将上述参数用于物质点法代表的土体中，在50kpa围压下进行含石量为0％、20％和40％的三轴试验，获得各个试样的应力-应变曲线如图6所示。可以看出，在各个含石量的条件下，物质点-离散元耦合方法的结果与离散元的结果都十分接近，硬化和软化现象都得到了准确的反映。而且，随着含石量的增加，试样越早达到峰值，且峰值强度更高这一趋势也得到了准确再现。
[0120]
采用上述宏观参数，用一个物质点代替多个土颗粒，分别建立了粒径比为1：5的试样，共有0％，20％和40％三种含石量，其模型如图7所示。将上述试样，在50kpa的围压下进行三轴试验，最终的结果如图8所示。
[0121]
可以看出，即使减少了物质点的数目，1：10粒径比条件下的力学特性仍能得到准确反映。而通过比较试样的总粒子数，发现1：5的试样粒子数仅为1：10的1/8，这说明至少可以节约87.5％的计算量。
[0122]
实施例2：
[0123]
采用物质点-离散元耦合方法模拟了实际的中型三轴试验，试样为0.1m*0.1m*0.2m的长方体试样，含石量均为50％，但是有三种不同的粒径比，其模型如图9所示。其级配如表2所示：
[0124]
表2中型三轴试样的级配
[0125]
粒径(m)级配1级配2级配30.016
–
0.02050％0％0％
0.010
–
0.0160％500％0.005
–
0.0100％0％50％0.5e-3
–
1.0e-30.6％0.6％0.6％0.25e-3
–
0.5e-31.3％1.3％1.3％0.075e-3
–
0.25e-314.1％14.1％14.1％《0.075e-334.0％34.0％34.0％
[0126]
对于每种级配的试样，分别在200kpa、400kpa、800kpa围压下，采用物质点-离散元耦合方法进行了模拟。石块的密度、杨氏模量、泊松比分别设置为2700kg/m3、50gpa和0.2，土颗粒的微观参数难以直接获得，所以假设与石块的参数一致。土体宏观参数参数如表3，为了模拟出硬化现象，采用了线性的硬化函数。土体与石块之间、石块与石块之间的摩擦系数均设为0.3，忽略土体与石块之间的粘结。
[0127][0128]
获得最终结果如图10-12所示。可以看出，模拟的结果和实际三轴试验的应力-应变曲线非常接近，试样的持续硬化现象得到了很好的反映。同时，可以看出，石块的粒径会明显的影响试样的力学特性。特别是在高围压条件下，石块的粒径越小，试样的硬化现象就越明显，试样的峰值强度就越高。
[0129]
由上述实施可知，物质点法和离散元法相互耦合，充分发挥了二者的优势。采用物质点法来离散土体，利用其连续介质特性，不再需要关注土颗粒之间的接触，而是采用宏观本构来描述土体的力学响应。同时，在不影响到混合体的整体的力学性质的前提下，可以适当增大物质点的尺寸，以减少模型中的粒子数，节省计算资源；采用离散元法来建模石块，可以准确捕捉每个石块的运动和接触行为。同时，石块与石块之间、石块与土体之间的相互作用可以被考虑，一定程度上，保留了反映微观机制的能力。采用物质点法作为土体建模的方法，可以避免网格畸变的影响，便于模拟大变形问题。同时，物质点法在众多无网格法中也有着较高的计算效率，这使得耦合方法在效率上可以得到保证。在本发明中，物质点法与离散元法仅在计算耦合接触力时进行交互，二者各自的计算流程是相互独立的，因此可以采用并发式求解方法，实现方便，且能大幅提高计算效率。并且，采用本方法对土石混合体问题进行求解，相比于离散元法，单个土颗粒的参数不需要被精确测量，仅需给出一个大概的数值即可，其建模难度大大降低。也就是说，本发明提供的方法方便、快捷、精度高，既简化了工作流程也方便处理力学特性复杂的土石混合体问题，能够更好地指导工程实践，具有广阔的推广空间。
[0130]
以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。

技术特征：
1.一种土石混合体混合表征与并发式求解方法，其特征在于，所述土石混合体混合表征与并发式求解方法包括：在确定所述目标土石混合体的土体和石块后，对所述土体以物质点法进行建模，并生成与所述土体对应的物质点，对所述石块以离散元法进行建模，并生成与所述石块对应的离散单元；对所述土体和所述石块之间的作用关系以预设接触模型进行描述，并构造所述目标土石混合体的mpm-dem表征方案；计算各所述物质点的虚半径；根据所述虚半径，确定与所述离散单元有接触的物质点，并标记为虚粒子；采用离散元接触模型计算所述虚粒子与所述离散单元的耦合接触力；根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度；根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时速度和位置等参数。2.根据权利要求1所述的土石混合体混合表征与并发式求解方法，其特征在于，计算各所述物质点的虚半径，包括：采用公式计算各所述物质点的虚半径，其中，为第p个物质点的虚半径，为第p个物质点的体积，k
p
为第p个物质点的孔隙率。3.根据权利要求1所述的土石混合体混合表征与并发式求解方法，其特征在于，采用离散元法计算所述虚粒子与所述离散单元的耦合接触力，包括：采用公式计算法向合力，其中，为第p个离散单元受到的来自于对应物质点的法向合力，为第q个物质点作用于第p个离散单元上的法向力，pcontact为与第p个离散单元相接触的粒子组成的邻域列表；采用公式计算切向合力，其中，为第p个离散单元受到的来自于对应物质点的切向合力，为第q个物质点作用于第p个离散单元上的切向力，pcontact为与第p个离散单元相接触的粒子组成的邻域列表；将所述法向合力和所述切向合力合成所述耦合接触力。4.根据权利要求1所述的土石混合体混合表征与并发式求解方法，其特征在于，根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，包括：构造物质点法平衡方程其中，为物质点法中第i个节点上的合内力，为物质点法中第i个节点上的合外力，为物质点法中第i个节点的质量，为物质点法中第i个节点的加速度；求解所述物质点法平衡方程，获得物质点法中各背景网格节点上的加速度。
5.根据权利要求4所述的土石混合体混合表征与并发式求解方法，其特征在于，求解所述物质点法平衡方程，获得物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，还包括：采用公式计算所述合内力，其中，n
ip
为形函数，为外界施加到物质点上的体力，为物质点法中第i个节点上所受到的面力，为第p个离散单元受到的来自于对应物质点的法向合力，为第p个离散单元受到的来自于对应物质点的切向合力；采用公式计算所述合外力，其中，为第p个物质点上的应力，为第p个物质点的体积，b
ip
为形函数梯度；根据所述合内力和所述合外力计算各背景网格节点上的加速度。6.根据权利要求4所述的土石混合体混合表征与并发式求解方法，其特征在于，根据所述耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各所述离散单元上的加速度和角加速度，还包括：构造离散元法平衡方程以及其中，为第p个离散单元上受到的来自其他离散单元的法向合力，为第p个离散单元上受到的来自其他离散单元的切向合力，为第p个离散单元上受到的体力，为第p个离散单元上受到的转动阻尼力，为第p个离散单元的质量，为第p个离散单元的加速度，i
p
为第p个离散单元的惯性矩，为第p个离散单元的角加速度，为第p个离散单元的粒径，n
pq
为接触点的单位法向量，δ
n,pq
为p个离散单元与第q个离散单元(或物质点)之间的法向堆叠量，为第p个离散单元上所受到的来自第q个离散单元的力；求解所述离散元法平衡方程，获得各所述离散单元上的加速度和角加速度。7.根据权利要求1所述的土石混合体混合表征与并发式求解方法，其特征在于，所述物质点的实时参数包括物质点速度，所述离散单元的实时参数包括离散单元速度；根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时参数，包括：采用flip速度更新格式对所述物质点速度进行实时更新，其中，所述flip速度更新格式为以及以及为第p个物质点的相邻节点所组成的邻域列表，δt为物质点的时间步长，为第p个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度，为第p个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度；采用预设速度更新格式对所述离散单元速度进行实时更新，其中，所述预设速度更新格式为以及以及为第p个物质点在t
n+1/2
时刻的速度，为第p个物质点在t
n-1/2
时刻的速度，为第p个物质点在t
n+1/2
时刻的加速度，为第p个物质点在t
n-1/2
时刻的加速度。8.根据权利要求1所述的土石混合体混合表征与并发式求解方法，其特征在于，所述物质点的实时参数包括物质点位置，所述离散单元的实时参数包括离散单元位置；根据所述背景网格节点上的加速度、所述离散单元的加速度以及所述离散单元的角加速度，计算所述物质点和所述离散单元的实时参数，包括：采用更新格式对所述物质点位置进行实时更新，其中，为第p个物质点在t
n+1
时刻的物质点坐标，为第p个物质点在t
n-1
时刻的物质点坐标；采用预设位置更新格式对所述离散单元位置进行实时更新，其中，所述预设位置更新格式为以及以及为第p个离散单元在t
n+1
时刻的离散单元坐标，为第p个离散单元在t
n-1
时刻的离散单元坐标，为第p个离散单元在t
n+1
时刻的离散单元坐标，为第p个离散单元在t
n-1
时刻的离散单元坐标。9.根据权利要求1所述的土石混合体混合表征与并发式求解方法，其特征在于，根据求解结果，以预设更新格式分别更新所述物质点的应力、速度和位置、以及所述离散单元的速度和位置，还包括：将物质点的动量参数映射至背景网格上；根据所述动量参数，以预设方程求解网格节点上的速度，其中，所述预设方程为以及以应变更新格式更新物质点上的应变，其中，所述应变更新格式设置为以旋量更新格式更新物质点上的旋量，其中，所述旋量更新格式设置为根据所述应变和所述旋量，采用预设本构模型更新物质点上的应力。10.根据权利要求1所述的土石混合体混合表征与并发式求解方法，其特征在于，在确定所述目标土石混合体的土体和石块后，对所述土体以物质点法进行建模，获取与所述土体对应的物质点，对所述石块以离散元法进行建模，获取与所述石块对应的离散单元的步
骤之前，还包括：获取目标土石混合体的分界粒径；根据所述分界粒径，确定所述目标土石混合体的土体和石块。

技术总结
本发明公开一种土石混合体混合表征与并发式求解方法，包括：在确定目标土石混合体的土体和石块后，对土体以物质点法进行建模，并生成与土体对应的物质点，对石块以离散元法进行建模，并生成与石块对应的离散单元；对土体和石块之间的作用关系以预设接触模型进行描述；计算各物质点的虚半径；并确定与离散单元有接触的物质点，标记为虚粒子；采用离散元接触模型计算虚粒子与离散单元的耦合接触力；根据耦合接触力，构造平衡方程，并计算物质点法中各背景网格节点上的加速度、以及各离散单元上的加速度和角加速度；根据背景网格节点上的加速度、离散单元的加速度以及离散单元的角加速度，计算物质点和离散单元的实时速度和位移等参数。等参数。等参数。


技术研发人员：王斌 李建国 王頔
受保护的技术使用者：中国科学院武汉岩土力学研究所
技术研发日：2023.03.28
技术公布日：2023/7/25