一种基于持续同调方法的颗粒材料液化预测方法

技术领域

本发明属于颗粒物质力学研究领域，涉及颗粒材料液化现象的预测方法，具体指一种基于持续同调方法的颗粒材料液化预测方法。

背景技术

颗粒材料是大量离散的固体颗粒相互作用而组成的复杂体系，其普遍存在于工程建设、工业生产和自然界的方方面面，自然界中广泛存在的泥石流及堰塞体、土木工程中所用的块石状材料、工业生产中的煤炭及砂石等均涉及颗粒材料的堆积和运动等过程。颗粒材料存在从离散到连续、从微观到宏观、从无序到有序和从流动到堵塞等特性，是物理、材料、力学等基础学科中的重点研究对象之一。与传统意义上的固体和液体不同，颗粒体系由于其复杂的微观结构在不同剪切状态会呈现出类似气体、液体或固体的复杂运动模式。著名学术期刊《Science》曾在2005年公布了125个最具挑战性的世纪科学难题，其中之一便是“能否发展关于湍流动力学和颗粒材料运动学的综合理论”。

现有研究发现，受到剪切作用的颗粒体系存在两种破坏模式，分别称为局部剪切带破坏和整体的分散性破坏。其中，局部剪切带破坏是指体系的变形主要集中在一个带状区域，而分散性破坏发生时无局部集中变形的出现，在试验和工程中难以观察到，因此这方面的研究相对单薄。近年来针对等体积剪切颗粒体系的研究发现，当体系孔隙比达到一定值，也即体系达到一定疏松度时，体系会发生剧烈的失稳现象，被称为颗粒材料的液化现象，这一现象属于分散性破坏的一个特例。液化现象的发生极快，颗粒体系会在极快的时间内瞬间失去承载能力，变成“类液态”。此外，发生液化时体系的应力状态未达到塑性极限破坏强度，这无疑对基于传统建模方法预测液化现象的发生带来了极大的挑战，急需发展新型的液化预测方法。通过对大量的自然界地质灾害和水利、岩土工程破坏失效案例分析，人们逐渐认识到颗粒材料的液化行为是山体滑坡、泥石流、坝体溃决、土石填方失稳等地质和工程灾害的孕育机制之一。因此，研究颗粒材料液化现象的预测及建模方法对地质灾害防灾减灾和岩土工程安全运行具有重要的理论和实际应用价值。

另一方面，颗粒材料本质是大量固体颗粒的无序堆积，由此引起的结构无序性导致内部颗粒接触呈现典型的空间异质性，也即接触力的取向与量级在空间上存在显著的非均匀性。本世纪初，一些学者通过光弹等光学技术首次在实验中观测到空间异质分布的颗粒接触力网络，随后的大量研究认为颗粒体系的接触力网络是颗粒体系最本质的特征，其结构与演化在很大程度上决定了体系的宏观力学响应。为进一步研究接触力网络的结构与演化，大量学者通过人为定义力链、力环等介观结构来反映体系接触力网络的演化，相关研究显著推动了颗粒物质力学的发展。然而，上述介观结构主要是从局部的角度定义，而由颗粒集合体堆积形成的接触力网络是一个整体，特别是不同空间尺度下呈现的异构拓扑组织和复杂行为模式，所以应该从系统的角度来整体研究接触力网络的拓补结构与演化。

接触力网络本质上属于复杂网络的范畴，近年来，图论及其相关方法被广泛用于研究颗粒材料接触力网络的拓补结构与演化，主要包括通过社区检测方法提取体系的力链结构，预测颗粒间接触力，探测体系的失效路径等。作为图论及拓补数据分析中的一个基本工具，持久同调方法通过量化复杂网络中拓扑不变量的持久性，为定量研究复杂网络的动力学行为提供了一个基于代数拓扑的理论框架。其中，拓扑不变量的拓扑性质不随网络结构的连续变形而改变，被视为持久同调的关键特征，主要包括连通组件、环和空洞等结构。近年来，持续同调及其相关方法在材料科学、生物学和自然语言处理等不同领域展现了广阔的应用前景。

因此，有必要设计一种基于持续同调方法的颗粒材料液化预测方法，以克服上述问题。

发明内容

本发明提供一种基于持续同调方法的颗粒材料液化预测方法，以解决相关技术中仅从局部的角度研究介观结构的演化，无法从整体角度定量表征接触力网络的动力学行为和预测颗粒材料液化现象的问题。

为了解决上述技术问题，本发明采用如下技术方案：

一种基于持续同调方法的颗粒材料液化预测方法，包括如下步骤：

步骤S1.根据颗粒体系中颗粒间的相互接触关系，将原始颗粒体系离散为接触力网络；

步骤S2.基于持续同调方法，通过多尺度拓补分析，识别接触力网络中每个拓补不变量产生与消亡的接触力阈值，构建对应的持续图，实现对关键拓补结构信息的提取与接触力网络的降维；

步骤S3.计算相邻两个时刻持续图间的距离，量化接触力网络随剪切过程的演化速率；

步骤S4.探测上述接触力网络演化速率曲线是否存在突变点，如果检测存在突变点，则说明接触力网络已开始加速崩溃，预示着液化现象的发生。

优选地，所述步骤S1具体为：

对于一个确定的颗粒体系，首先确定体系内每个颗粒与其邻接颗粒的接触关系，包括是否接触与接触力量级信息；根据上述接触关系，将体系中每个颗粒作为接触力网络中的一个节点，两个颗粒间如果存在接触，则在对应的两个节点中添加一条边，该边包含着接触力量级的信息。

优选地，以边的宽度来代表不同接触力的量级，边的宽度越大代表对应两个颗粒间的接触力越大。

优选地，所述步骤S2中包括如下子步骤：

步骤S2.1.首先选定体系所有颗粒接触力中的最大值作为接触力阈值，选取所有边中接触力大于该接触力阈值的部分节点与边构造接触力子网络，再从该接触力子网络中提取各个拓补不变量，包括连通组件与接触环；

步骤S2.2.不断降低上述的接触力阈值，并重复上述操作；在这一过程中，可观测到各个拓补不变量随接触力阈值降低而出现或消失；

步骤S2.3.同时记录下每个拓补不变量出现或消失时的接触力阈值,基于此，构建对应的持续图，实现对关键拓补结构信息的提取与接触力网络的降维。

优选地，所述接触力阈值均采用体系所有接触力的平均值作为接触力阈值的基本单位

优选地，所述步骤S3具体包括：

在颗粒体系剪切过程中，于每个时刻分别构建所述的持续图，采用WassersteinDistance计算每两个相邻时刻的持续图之间的距离，此距离量化了接触力网络在这两个相邻时刻之间的演化速率，得到接触力网络随剪切过程的演化速率曲线。

优选地，计算相邻两个时刻持续图之间的距离，具体计算公式为：

其中，PD和PD′代表前后两个时刻的持续图，WD(PD,PD′)为两个相邻时刻持续图之间的距离，φ(PD,PD′)代表由持续图PD到持续图PD′的点对点映射，对于持续图PD中的一点p，φ(p)代表其经过映射φ在持续图PD′中的对应点，d(p,φ(p))代表两点之间距离。

优选地，两点之间距离d(p,φ(p))的计算公式为：

d(p,φ(p))＝max(|f

其中，f

优选地，所述步骤S4具体包括：采用Bernaola Galvan分割算法检测曲线中是否存在突变点，若所述曲线存在突变点，说明接触力网络发生加速崩溃，预示着液化现象的发生；反之说明颗粒材料不会发生液化现象。

与现有技术相比，本申请具有如下有益效果：

(1)本发明提供了一种基于持续同调方法的颗粒材料液化预测方法；

(2)本发明所提供的方法可以实现接触力网络中各种拓补不变量的全面有效提取，可有效描述接触力网络的拓补结构，并为其他非晶体系的相关研究提供复杂体系的拓补特征提取方式；

(3)本发明所提供的方法可以从系统和整体的角度定量研究接触力网络的动力学演化，为颗粒材料液化现象的预测提供了一条崭新且有效的途径；

(4)由于颗粒材料接触力网络的高维复杂性，以往构建颗粒材料本构模型的研究难以考虑接触力网络的演化，因此，本发明所提供的方法可为建立考虑接触力网络演化的颗粒材料本构模型提供参考。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种基于持续同调方法的颗粒材料液化预测方法；

图2是本发明实施例提供的原始颗粒体系的示意图；

图3是本发明实施例提供的颗粒剪切实验得到的宏观应力应变曲线图；

图4是本发明实施例提供的颗粒体系接触力网络示意图的示意图；

图5是本发明实施例提供的接触力网络随接触力阈值的变化图；

图6是本发明实施例提供的基于持续同调方法获取的接触力网络持续图；

图7是本发明实施例提供的接触力网络演化速率随剪切过程的曲线图。

具体实施方式

以下结合附图和实施例详细说明本发明的技术方案。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供一种基于持续同调方法的颗粒材料液化预测方法，其能解决相关技术中仅从局部的角度研究介观结构的演化，无法从整体角度定量表征接触力网络的动力学行为和预测颗粒材料液化现象的问题。

参见图1所示，为本发明实施例提供的一种基于持续同调方法的颗粒材料液化预测方法，其可以包括以下步骤：

步骤S1.根据颗粒体系中颗粒间的相互接触关系，将原始颗粒体系离散为接触力网络；

步骤S2.基于持续同调方法，通过多尺度拓补分析，识别接触力网络中每个拓补不变量产生与消亡的接触力阈值，构建对应的持续图，实现对关键拓补结构信息的提取与接触力网络的降维；

步骤S3.计算相邻两个时刻持续图间的距离，量化接触力网络随剪切过程的演化速率；

步骤S4.探测上述接触力网络演化速率曲线是否存在突变点，如果检测存在突变点，则说明接触力网络已开始加速崩溃，预示着液化现象的发生。其具体过程如下：

1.进行颗粒材料的离散元数值模拟或物理实验。模拟或实验参数视具体情况确定即可。如附图2所示，给定一原始颗粒体系，对其施加一定的剪切作用，如双轴剪切、单剪、等体积剪切等不同的剪切形式。在剪切过程中，首先记录体系宏观偏应力随宏观剪切应变的变化趋势。如附图3所示，本发明实施例采用的颗粒体系在剪切过程中，宏观偏应力逐渐减小至0，表明体系在该处彻底失去承载能力，也即发生了液化行为。此外，还需记录每个时刻体系中每个颗粒的位置信息及其与其他颗粒的接触关系。其中，位置信息是指每个颗粒在空间中的坐标位置，接触关系主要包括每个颗粒与其周围颗粒是否接触和接触力量级等信息；

2.构造颗粒体系接触力网络。如附图4所示，将颗粒体系中的每个颗粒都作为接触力网络中的一个节点。根据上述记录的颗粒位置信息与接触关系，在存在接触的两个节点之间添加一条边，同时这条边包含着接触力量级的信息，附图4中以边的宽度来代表不同接触力的量级，边的宽度越大代表对应两个颗粒间的接触力越大。基于此，可将原始颗粒体系离散为接触力网络的形式；

3.基于持续同调方法，构造接触力网络的持续图。如附图4所示，颗粒材料的接触力网络在空间中呈现高度异质性的分布特点。为了量化这一复杂网络的拓补结构及其演化，基于持续同调方法的思想，重点关注接触力网络中拓补不变量随不同接触力阈值的诞生与消亡。

具体步骤为：

a.首先选定体系所有颗粒接触力中的最大值作为接触力阈值，再选取接触力网络所有边中接触力大于该接触力阈值的部分节点与边构造接触力子网络，再从该接触力子网络中提取各个拓补不变量，包括连通组件与接触环。其中连通组件是指由若干节点通过边连接而成的连通结构，接触环是指由若干节点通过边连接而形成的闭环结构；

b.不断降低上述的接触力阈值，并重复上述操作。如附图5所示，采用体系所有接触力的平均值F

c.根据所记录的每个拓补不变量诞生或消亡时的接触力阈值，以诞生时的接触力阈值为X轴数据，消亡时的接触力阈值为Y轴数据的接触力网络持续图。如附图6所示，持续图中的每个点均代表一个拓补不变量诞生和消亡时对应的接触力阈值。

基于此，可实现基于持续同调方法，构造接触力网络的持续图；

4.计算相邻两个时刻持续图间的距离，量化接触力网络的演化速率。在颗粒体系剪切过程中，在每个时刻均构造上述接触力网络的持续图。持续图编码了每个拓补不变量诞生和消亡的信息，代表着接触力网络的最本质特征。因此，量化接触力网络持续图在相邻两个时刻之间的距离，等同于量化接触力网络在对应时刻的演化速率。采用Wassersteindistance计算相邻两个时刻持续图之间的距离，具体计算公式为：

其中，PD和PD′代表前后两个时刻的持续图，WD(PD,PD′)为上述两个持续图之间的距离，φ(PD,PD′)代表由持续图PD到持续图PD′的点对点映射，由于可能存在着多种可能的点对点映射组合，最终选定的距离是所有可能映射中计算结果的最小值。对于持续图PD中的一点p，φ(p)代表其经过映射φ在持续图PD′中的对应点，该两点之间距离d(p,φ(p))的计算公式为：

d(p,φ(p))＝max(|f

其中，f

基于上述计算过程，可以获得如附图7所示的接触力网络随剪切应变增加的演化速率曲线图；

5.探测接触力网络的演化曲线是否存在突变点。根据上述的接触力网络随剪切应变增加的演化速率曲线图，采用Bernaola-Galvan算法探测其中是否存在突变点。Bernaola-Galvan算法基于启发式分割思想，可从非平稳时间序列中快速获取突变点等信息。突变点的存在说明接触力网络经历着加速崩溃的过程，预示着颗粒材料液化现象的发生；反之，不存在突变点说明接触力网络在经历平稳演化，体系没有发生液化现象的风险。如附图7中的箭头所示，采用Bernaola-Galvan算法对前述接触力网络演化速率曲线图进行探测后，确实发现演化曲线中存在两个突变点。说明本发明例采用的颗粒体系在随后的剪切过程中会发生接触力网络的加速崩溃与液化行为，这一点可以从附图3中得到佐证，证明了本发明具有较好的实施效果，并确能实现定量描述颗粒材料接触力网络的动力学演化与宏观液化行为的预测。

在本发明的描述中，需要说明的是，术语“上”、“下”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

需要说明的是，在本发明中，诸如“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

以上所述仅是本发明的具体实施方式，使本领域技术人员能够理解或实现本发明。对这些实施例的多种修改对本领域的技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所申请的原理和新颖特点相一致的最宽的范围。