CN114613458B - 一种基于常规态型近场动力学的颗粒破碎模拟方法及系统 - Google Patents

一种基于常规态型近场动力学的颗粒破碎模拟方法及系统 Download PDF

Info

Publication number
CN114613458B
CN114613458B CN202210411866.6A CN202210411866A CN114613458B CN 114613458 B CN114613458 B CN 114613458B CN 202210411866 A CN202210411866 A CN 202210411866A CN 114613458 B CN114613458 B CN 114613458B
Authority
CN
China
Prior art keywords
particle
point
bond
value
key
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210411866.6A
Other languages
English (en)
Other versions
CN114613458A (zh
Inventor
肖源杰
华文俊
李志勇
陈宇亮
谢峰
王小明
李婷玉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hunan Communications Research Institute Co ltd
Central South University
Original Assignee
Hunan Communications Research Institute Co ltd
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hunan Communications Research Institute Co ltd, Central South University filed Critical Hunan Communications Research Institute Co ltd
Priority to CN202210411866.6A priority Critical patent/CN114613458B/zh
Publication of CN114613458A publication Critical patent/CN114613458A/zh
Application granted granted Critical
Publication of CN114613458B publication Critical patent/CN114613458B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C60/00—Computational materials science, i.e. ICT specially adapted for investigating the physical or chemical properties of materials or phenomena associated with their design, synthesis, processing, characterisation or utilisation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Theoretical Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于常规态型近场动力学的颗粒破碎模拟方法及系统，包括以下步骤：进行室内单颗粒压碎试验，测定单颗粒破碎强度试验值；建立单颗粒物质点模型；率定排斥力常数等模型参数；进行单颗粒破碎数值试验近场动力学模拟；获取单颗粒破碎强度模拟值与碎片分布。本发明方法解决了单颗粒破碎室内试验离散性大、成本高、周期长等问题。与传统的离散元方法相比，该方法在计算上可以更方便地考虑不规则颗粒的形状，且更易于校准细观力学参数，并用于获取准确的脆性颗粒破碎强度值，为进一步认识单颗粒破碎的细观力学机制和可破碎颗粒材料的宏观力学行为提供了有效技术手段。

Description

一种基于常规态型近场动力学的颗粒破碎模拟方法及系统
技术领域
本发明属于岩土工程技术领域，特别是涉及一种基于常规态型近场动力学的颗粒破碎模拟方法及系统。
背景技术
颗粒材料中的颗粒破碎行为会改变颗粒级配和颗粒间的接触关系，从而影响颗粒材料的宏细观力学行为。由于室内试验往往存在离散性大、成本高、周期长等问题，通过数值模拟研究单个颗粒在单轴压缩荷载作用下的破碎条件和破碎模式对准确理解可破碎颗粒材料的宏观力学行为（如变形、屈服、流动和堆积特性）具有重要意义。目前，研究者们多利用离散元方法来模拟颗粒材料中多颗粒破碎行为，包括颗粒聚团方法和碎片替换方法，前一种方法计算规模特别大，而后一种方法无法真实反映单个颗粒的破碎行为，且离散元方法涉及的内部参数较难确定。相比之下，近场动力学方法采用单一的准则来捕获裂纹的萌生和扩展，并且可以自然地处理具有任意形状的颗粒。因此，有必要借助这一新兴的数值模拟方法开展单颗粒破碎机制研究。
发明内容
本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于常规态型近场动力学的颗粒破碎模拟方法及系统，借助近场动力学数值模拟，该方法可以得到不规则形状颗粒的破碎强度值，并能准确地表征颗粒在破碎过程中的损伤演化规律以及颗粒破碎后所产生的碎片分布情况，具有计算精度和计算效率高、适用性广等特点。
本发明的目的可以通过以下技术方案来实现：本发明提供一种基于常规态型近场动力学的颗粒破碎模拟方法及系统，包括步骤：
S1，进行室内单颗粒压碎试验，测定单颗粒破碎强度试验值；
S2，建立单颗粒物质点模型：根据室内试验所用单颗粒构建三维表面形貌模型，并根据所述三维表面形貌模型构建单颗粒物质点模型，确定所述单颗粒物质点模型的名义粒径；
S3，建立近场动力学模型：在所述单颗粒物质点模型的基础上，布置加载板物质点，定义边界条件，构建基于键的链表并计算键的有关参数；
S4，进行单颗粒破碎数值试验：采用显式时间积分方案模拟所述近场动力学模型中所述单颗粒物质点模型在单轴压缩荷载作用下的破碎过程，监测所述加载板物质点中接触力的最大值随所述加载板物质点位移的变化情况；
S5，确定所述近场动力学模型的参数：以步骤S1所述的单颗粒破碎室内试验结果为基础对排斥力常数
Figure 915304DEST_PATH_IMAGE001
进行率定；
S6，获取单颗粒破碎强度模拟值：确定上、下两块加载板与颗粒之间的接触力，并将所述接触力的平均值除以所述名义粒径的平方值作为颗粒破碎强度模拟值；
S7，获取颗粒破碎后的碎片分布：采用碎片提取算法识别颗粒破碎后所产生的碎片，并对碎片的数量、尺寸和形状进行量化分析。
本发明还公开一种基于常规态型近场动力学的颗粒破碎模拟系统，该系统包括上述基于常规态型近场动力学的颗粒破碎模拟方法。
与现有技术相比，本发明的有益效果体现在：
（1）本发明采用体素化离散方法对单颗粒三维表面形貌模型进行空间离散，可以方便地考虑不规则颗粒形状对颗粒破碎强度和破碎模式的影响。通过设定不同大小的物质点尺寸,该方法还可以在计算效率与几何模拟精度之间取得适当的平衡。
（2）本发明采用Kd树算法构建近邻物质点集合，通过建立有效的索引结构可以大大加快检索的速度，在物质点数量庞大时具有明显的优势。相较于时间复杂度为
Figure 424783DEST_PATH_IMAGE002
的全配对搜索法，该算法的时间复杂度仅为/>
Figure 411194DEST_PATH_IMAGE003
。
（3）本发明采用常规态型近场动力学理论计算物质点之间的长程作用力，可准以确地模拟不同材料类型的颗粒在单轴压缩荷载作用下的破碎特性，克服了键基近场动力学仅能模拟固定泊松比材料的缺陷。
（4）本发明所建立的数值模型无需定义内聚区单元或着引入额外的节点自由度和局部强化函数，可以方便地模拟材料内部损伤演化与宏观裂纹的萌生和扩展，除采用近场动力学理论进行力学计算外，还可推广应用于分子动力学、物质点法、光滑粒子流体动力学等无网格方法。
附图说明
图1为本发明一实施例中的球形颗粒单轴压缩模型示意图；
图2为本发明一实施例中的常规态型近场动力学模型示意图；
图3为本发明一实施例中的近邻物质点树形搜索结构示意图；
图4为本发明一实施例中的键的体积修正示意图；
图5为本发明一实施例中的键的刚度修正示意图；
图6为本发明一实施例中的短程接触力模型示意图；
图7为本发明一实施例中的键的断裂准则示意图；
图8为本发明一实施例中的球形颗粒单轴压缩破碎模拟中的荷载-位移曲线；
图9为本发明一实施例中的常规态型近场动力学模拟球形颗粒单轴压缩破碎过程损伤云图；
图10为本发明一实施例中的球形颗粒破碎形态及碎片分布图。
其中，1-单颗粒物质点，2-加载板物质点，3-近邻物质点集合，4-单颗粒物质点之间的键，5-单颗粒物质点与加载板物质点之间的键，6-加载板物质点之间的键。
具体实施方式
下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。
请一并结合图1-10，以下以一个具体的实例和参数来说明本发明的一种基于常规态型近场动力学的颗粒破碎模拟方法及系统，包括步骤：
S1，进行室内单颗粒压碎试验，测定单颗粒破碎强度试验值。在本发明中，采用直径为10 mm，形状近似为球形的石英砂颗粒进行单颗粒破碎试验；
具体的，将试验得到的荷载-位移曲线的最大峰值点或第一个典型峰值点定义为破碎点，将所述破碎点对应的荷载除以所述颗粒的名义粒径的平方值作为颗粒破碎强度试验值，为后续率定数值模型的排斥力常数等细观力学参数提供试验数据支撑。
S2，建立单颗粒物质点模型：根据室内试验所用单颗粒构建球形颗粒的三维表面形貌模型，并根据所述三维表面形貌模型构建单颗粒物质点模型，确定所述单颗粒物质点模型的名义粒径为10 mm；
所述单颗粒物质点模型为由所述单颗粒三维表面形貌模型通过体素化离散而成的物质点集合，物质点的尺寸设定为0.2 mm，满足小于单颗粒直径值的十分之一的条件。
具体的，体素化离散方法是把空间均匀地分割为小立方体的三维栅格，如图1所示。这些小立方体通常称作物质点（material point）或体素（voxel），与坐标轴对齐。考虑计算量与颗粒特征尺寸的平衡，选取物质点的尺寸为0.2 mm，则所述单颗粒物质点模型被离散为69694个方形物质点。
S3，根据表1所示的参数建立近场动力学模型，其中近场范围m值为近场范围尺寸
Figure 842175DEST_PATH_IMAGE004
与物质点尺寸/>
Figure 939444DEST_PATH_IMAGE005
的比值,参考已有研究将该值设定为3.015：在所述单颗粒物质点模型的基础上，布置加载板物质点，定义边界条件，构建基于键的链表并计算键的有关参数：
表1 近场动力学模型参数
参数 数值
颗粒的密度（kg/m3） 2650
颗粒的弹性模量（GPa） 100
颗粒的泊松比 0.15
颗粒的临界能量释放率（J·m-2） 10
加载板的密度（kg/m3） 3850
加载板的弹性模量（GPa） 380
加载板的泊松比 0.22
加载板的临界能量释放率（J·m-2） 1×106（足够大）
物质点的近场范围（m值） 3.015
物质点的尺寸（mm） 0.2
接触临界距离（mm） 0.19
排斥力常数 -5×1011
颗粒与加载板之间的摩擦系数 0.5
具体的，近场动力学是一种非局部的、基于连续体的无网格方法。如图2所示，该方法将一个连续的物质体离散为多个具有特定体积和质量的物质点，每个物质点通过“键”与它“近场范围”内的其他物质点相互作用，比如物质点
Figure 506692DEST_PATH_IMAGE006
通过键/>
Figure 664003DEST_PATH_IMAGE007
与物质点/>
Figure 815237DEST_PATH_IMAGE008
相互作用，物质点/>
Figure 450618DEST_PATH_IMAGE006
和物质点/>
Figure 137951DEST_PATH_IMAGE008
上分别作用有键力/>
Figure 262902DEST_PATH_IMAGE009
和键力/>
Figure 606158DEST_PATH_IMAGE010
。近场范围/>
Figure 45230DEST_PATH_IMAGE011
的大小由参数“近场范围尺寸/>
Figure 587070DEST_PATH_IMAGE004
”定义，该范围表示为图2所示的半径。物质点/>
Figure 86184DEST_PATH_IMAGE012
的近场范围内的所有物质点称为物质点/>
Figure 11677DEST_PATH_IMAGE012
的近邻物质点。对于未受损伤破坏的完整材料，每个物质点都与它的近邻物质点存在键的作用。
建立所述近场动力学模型包括步骤：
S31，在所述单颗粒物质点模型的顶部和底部对称布置10
Figure 254440DEST_PATH_IMAGE013
10/>
Figure 119627DEST_PATH_IMAGE013
0.8 mm（长度/>
Figure 586381DEST_PATH_IMAGE013
宽度/>
Figure 169809DEST_PATH_IMAGE013
高度）的长方体块作为加载板，则单块加载板由50/>
Figure 216262DEST_PATH_IMAGE014
50/>
Figure 732694DEST_PATH_IMAGE014
4 mm（长度/>
Figure 72147DEST_PATH_IMAGE014
宽度/>
Figure 408450DEST_PATH_IMAGE014
高度）个物质点组成，加载板的面积大于所述单颗粒物质点模型在水平面上的投影面积，加载板的厚度大于近场范围尺寸/>
Figure 789753DEST_PATH_IMAGE004
，可确保所施加的荷载和约束条件能在所述单颗粒物质点模型中得到准确地反映；
S32，定义一个大小为
Figure 363954DEST_PATH_IMAGE015
的二维逻辑数组作为位移约束标记/>
Figure 375772DEST_PATH_IMAGE016
，其长度设定为等于所述近场动力学模型中物质点的总数/>
Figure 261688DEST_PATH_IMAGE017
，其第一列、第二列和第三列数值分别标识物质点在X方向、Y方向和Z方向上的位移约束情况，其中逻辑值“0”代表所述物质点在所述方向上的位移不受约束，逻辑值“1”代表所述物质点在所述方向上的位移受到限制。将所述加载板物质点的所述位移约束标记/>
Figure 384365DEST_PATH_IMAGE016
的第一列和第二列数值设定为“0”，以确保模拟过程中加载板不发生水平面内的相对滑移。并将下加载板物质点的所述位移约束标记/>
Figure 376854DEST_PATH_IMAGE016
的第三列数值设定为“0”，以模拟试验过程中下加载板始终保持静止的情况；
S33，采用树形搜索算法确定所述近场动力学模型中所有物质点的近邻物质点集合
Figure 559574DEST_PATH_IMAGE018
，并创建键列表/>
Figure 870469DEST_PATH_IMAGE019
；所述近场动力学模型中所有物质点的近邻物质点集合/>
Figure 99DEST_PATH_IMAGE018
的确定方法为Kd树算法，该方法遵循一个有效的物质点的近邻物质点集合的存储程序，可以克服存储空间有限的问题。该方法采用两个不同的数组：数组#1将所有物质点的近邻物质点按顺序存储为一列；数组#2用作所述数组#1的指针，以便从所述数组#1中提取特定物质点的近邻物质点。所述数组#2中的每个元素都对应于数组#1中一个特定物质点的第一个近邻物质点的索引；所述键列表L是大小为/>
Figure 142368DEST_PATH_IMAGE020
的二维索引数组，其长度设定为等于所述近场动力学模型中键的总数/>
Figure 699251DEST_PATH_IMAGE021
，所述键列表L的第一列和第二列数值分别对应存在长程力作用的两个物质点的索引；
具体的，树形搜索算法建立在树形数据结构基础上，根据计算域
Figure 497443DEST_PATH_IMAGE022
中物质点的空间位置，将整个计算域/>
Figure 227501DEST_PATH_IMAGE023
递归分割成多个挂限，直到每个挂限有且仅有一个物质点，其中的每一个树节点都对应x、y、z维度上的一个区间。在构建好树形数据结构后，可以高效地求解物质点/>
Figure 734530DEST_PATH_IMAGE025
的近邻物质点。
考虑图3，该图中有11个物质点（编号为1-11）待添加进Kd树。以物质点
Figure 259052DEST_PATH_IMAGE026
为中心，建立一个边长为/>
Figure 544540DEST_PATH_IMAGE027
的方形近场范围，然后检测物质点/>
Figure 78290DEST_PATH_IMAGE026
所在层次的其他物质点的挂限是否与上述方向近场范围有重合。若重合，则判定该物质点为物质点/>
Figure 867254DEST_PATH_IMAGE026
的近邻物质点；若不重合，则进入下一个树节点所在的层次，针对其层次内的物质点重复上述检测和判断过程，直至所有物质点都被放置到这颗树中。尽管这是个二维的例子，但其可以直接扩展至三维。
S34，计算键的有关参数，所述键的有关参数包括未变形前键的长度
Figure 562678DEST_PATH_IMAGE028
、体积修正系数/>
Figure 335462DEST_PATH_IMAGE029
、键的类型以及刚度修正系数/>
Figure 407323DEST_PATH_IMAGE030
；根据参考构型中键两端的物质点
Figure 817838DEST_PATH_IMAGE031
和/>
Figure 949742DEST_PATH_IMAGE032
的空间坐标，按照下式计算所述未变形前键的长度/>
Figure 944243DEST_PATH_IMAGE028
：
Figure 819795DEST_PATH_IMAGE034
/>
式中：
Figure 380089DEST_PATH_IMAGE035
和/>
Figure 682895DEST_PATH_IMAGE036
分别为参考构型中物质点/>
Figure 164692DEST_PATH_IMAGE037
和/>
Figure 76890DEST_PATH_IMAGE038
的横坐标值、纵坐标值和竖坐标值。
积分域在理论上是圆形或球形区域，而在实际计算中，计算区域被离散为正方形或立方体子域，导致在圆形或球形边界上的积分不准确。
因此，为了减少邻域截断误差的影响，按照下式计算所述体积修正系数
Figure 960533DEST_PATH_IMAGE040
：
Figure 434240DEST_PATH_IMAGE041
式中：
Figure 137753DEST_PATH_IMAGE042
为未变形前键的长度；/>
Figure 355108DEST_PATH_IMAGE043
为近场范围尺寸；/>
Figure 93257DEST_PATH_IMAGE044
为物质点的直径。
考虑图4，图中每一个正方形子域的边长等于物质点的尺寸
Figure 737865DEST_PATH_IMAGE044
，灰色（包括浅灰色和深灰色）子域中的物质点为需要积分的物质点，白色子域中的物质点则不与中心物质点
Figure 194254DEST_PATH_IMAGE045
存在相互作用。/>
Figure 451185DEST_PATH_IMAGE046
为与中心物质点/>
Figure 309420DEST_PATH_IMAGE045
通过键/>
Figure 859350DEST_PATH_IMAGE047
作用的近邻物质点，键/>
Figure 803035DEST_PATH_IMAGE047
的体积修正系数值为1；/>
Figure 362192DEST_PATH_IMAGE048
为与中心物质点/>
Figure 74933DEST_PATH_IMAGE045
通过键/>
Figure 61344DEST_PATH_IMAGE049
作用的近邻物质点，键/>
Figure 226746DEST_PATH_IMAGE049
的体积修正系数值为
Figure 619288DEST_PATH_IMAGE050
。尽管这是个二维的粒子，但其可以直接扩展至三维。
根据所述键的两端物质点的材料类型确定所述键的类型，所述键的类型包括颗粒-颗粒键、加载板-加载板键以及颗粒-加载板键。
具体的，当所述键两端的物质点均为单颗粒物质点时，则所述键的类型为颗粒-颗粒键；当所述键两端的物质点均为加载板物质点时，则所述键的类型为加载板-加载板键；当所述键一端的物质点为单颗粒物质点，而另一端的物质点为加载板物质点时，则所述键的类型为颗粒-加载板键。
在近场动力学模型中，键通常用于描述连续体内部不同物质点间的长程力作用。对于分属不同物体的物质点（比如单颗粒物质点和加载板物质点），则应该采用短程力来描述它们由于相互接触所产生的局部作用。因此在模型求解过程中，需要计算颗粒与加载板之间的接触力，而无需计算颗粒-加载板键的键力。
同样的，理论上假设近场域为完整的圆形或球形区域，而在实际计算过程中，模型表面附近不可避免地存在非圆形（或非球形）的近场域。这些近场域内部的物质点具有更少的相互作用（键），导致模型表面附近材料刚度的降低（即所谓的表面效应），如图5所示。针对靠近模型表面的近场域内的键，应按照下式计算所述刚度修正系数
Figure 452115DEST_PATH_IMAGE051
：
Figure 609427DEST_PATH_IMAGE053
Figure 465387DEST_PATH_IMAGE055
式中：
Figure 366347DEST_PATH_IMAGE056
为近场范围的体积；/>
Figure 53680DEST_PATH_IMAGE057
和/>
Figure 178631DEST_PATH_IMAGE058
分别为所述物质点/>
Figure 521888DEST_PATH_IMAGE031
和/>
Figure 462424DEST_PATH_IMAGE032
的近邻物质点的体积；/>
Figure 4264DEST_PATH_IMAGE027
为近场范围尺寸；/>
Figure 503379DEST_PATH_IMAGE059
为圆周率。
S4，进行单颗粒破碎数值试验：采用显式时间积分方案模拟所述近场动力学模型中所述单颗粒物质点模型在单轴压缩荷载作用下的破碎过程，监测所述加载板物质点中接触力的最大值随所述加载板物质点位移的变化情况：
步骤S4中，采用所述显式时间积分方案模拟所述近场动力学模型中所述单颗粒物质点模型在单轴压缩荷载作用下的破碎过程包括步骤：
S41，初始化所述近场动力学模型中所有物质点的信息并开始时程积分。所述信息包括：键的断裂状态
Figure 927407DEST_PATH_IMAGE060
、变形后键的长度/>
Figure 904590DEST_PATH_IMAGE061
、键的伸长率/>
Figure 97674DEST_PATH_IMAGE062
、键的软化系数/>
Figure 266225DEST_PATH_IMAGE063
、软化标记/>
Figure 380811DEST_PATH_IMAGE064
、节点力/>
Figure 427265DEST_PATH_IMAGE065
、物质点的速度/>
Figure 678117DEST_PATH_IMAGE066
、物质点的的位移/>
Figure 519035DEST_PATH_IMAGE067
；
S42，施加边界条件。在恒位移加载条件下，通过控制所述加载板物质点在Z方向上的运动速度以及在X方向和Y方向上的位移约束来模拟所述单颗粒破碎室内试验中加载板的压缩作用；
具体的，根据所述步骤S32中所述位移约束标记
Figure 120917DEST_PATH_IMAGE068
将所述加载板物质点相应方向上的位移值设定为零，并对上加载板物质点施加沿Z方向向下的大小为/>
Figure 439903DEST_PATH_IMAGE069
的位移。
S43，计算变形后键的长度
Figure 810842DEST_PATH_IMAGE070
和键的伸长率/>
Figure 324125DEST_PATH_IMAGE071
。根据现时构型中键两端的物质点的空间坐标计算所述变形后键的长度/>
Figure 413303DEST_PATH_IMAGE070
，并按照下式计算所述键的伸长率/>
Figure 535980DEST_PATH_IMAGE071
：
Figure 558163DEST_PATH_IMAGE072
式中：
Figure 6462DEST_PATH_IMAGE073
为键未变形前的长度；/>
Figure 317357DEST_PATH_IMAGE074
为键变形后的长度。
S44，更新软化标记
Figure 446987DEST_PATH_IMAGE075
并计算键的软化系数/>
Figure 792518DEST_PATH_IMAGE076
。所述软化标记/>
Figure 453131DEST_PATH_IMAGE077
是大小为/>
Figure 251322DEST_PATH_IMAGE078
的一维逻辑数组，其长度设定为等于所述近场动力学模型中键的总数
Figure 981381DEST_PATH_IMAGE079
，其数值标识键的软化状态，其中逻辑值“0”代表所述键的伸长率/>
Figure 915839DEST_PATH_IMAGE080
大于弹性伸长率/>
Figure 705940DEST_PATH_IMAGE081
，键的变形进入软化阶段，逻辑值“1”代表所述键的伸长率/>
Figure 991428DEST_PATH_IMAGE080
小于弹性极限/>
Figure 259598DEST_PATH_IMAGE081
，键的变形仍处于弹性阶段。为模拟准脆性材料的渐进损伤，按照下式计算所述键的软化系数
Figure 815607DEST_PATH_IMAGE082
：
Figure 776610DEST_PATH_IMAGE083
Figure 283814DEST_PATH_IMAGE084
式中：
Figure 90096DEST_PATH_IMAGE080
为键的伸长率；/>
Figure 264726DEST_PATH_IMAGE085
为键的临界伸长率；/>
Figure 927788DEST_PATH_IMAGE081
为键的弹性伸长率；/>
Figure 187868DEST_PATH_IMAGE086
为材料的单轴抗拉强度；/>
Figure 296376DEST_PATH_IMAGE087
为材料的弹性模量。
具体的，对于所述步骤S1中的所述石英砂颗粒，其在侧限压力较小的压缩试验中多呈现脆性破坏特性。因此，为准确模拟所述单颗粒在单轴压缩荷载作用下的破碎行为，模型求解过程中键的软化系数
Figure 325512DEST_PATH_IMAGE088
始终设定为零，即不考虑键的软化现象。
S45，计算键力
Figure 362738DEST_PATH_IMAGE089
、接触力/>
Figure 844535DEST_PATH_IMAGE090
以及节点力/>
Figure 258199DEST_PATH_IMAGE091
。采用典型的的线性近场动力学固体（linear peridynamic solid）材料本构模拟物质点之间的非局部长程作用力，按照下式计算所述键力/>
Figure 141842DEST_PATH_IMAGE092
：
Figure 615548DEST_PATH_IMAGE094
/>
Figure 882844DEST_PATH_IMAGE095
Figure 834619DEST_PATH_IMAGE096
Figure 572768DEST_PATH_IMAGE097
Figure 217376DEST_PATH_IMAGE098
Figure 673765DEST_PATH_IMAGE099
式中：
Figure 225969DEST_PATH_IMAGE100
为键的某一端物质点/>
Figure 818625DEST_PATH_IMAGE101
的加权体积；/>
Figure 867090DEST_PATH_IMAGE102
为键的某一端物质点/>
Figure 607513DEST_PATH_IMAGE101
的体积膨胀率；/>
Figure 166670DEST_PATH_IMAGE103
为键的某一端物质点/>
Figure 879411DEST_PATH_IMAGE101
的近场范围；/>
Figure 600242DEST_PATH_IMAGE104
为键的某一端物质点/>
Figure 765645DEST_PATH_IMAGE101
的体积；/>
Figure 862914DEST_PATH_IMAGE105
为键的伸长量，分为各向同性部分/>
Figure 993943DEST_PATH_IMAGE106
和偏态部分/>
Figure 885676DEST_PATH_IMAGE107
；/>
Figure 538374DEST_PATH_IMAGE108
和/>
Figure 236071DEST_PATH_IMAGE109
分别为未变形前键的长度和变形后键的长度；/>
Figure 392246DEST_PATH_IMAGE110
和/>
Figure 15732DEST_PATH_IMAGE111
分别为刚度修正系数和体积修正系数；/>
Figure 155726DEST_PATH_IMAGE112
和/>
Figure 860377DEST_PATH_IMAGE113
分别代表体积模量和剪切模量。
采用短程排斥力模型和库伦定律描述所述单颗粒物质点与所述加载板物质点之间的相互作用。当分属不同物体（比如单颗粒和加载板）的两个物质点之间的距离小于临界距离时，认为这两个物质点相互接触，沿着两个物质点的连线方向作用有大小相等，方向相反的法向接触力，如图6所示。当两个相接触的物质点在垂直于法线方向上的相对运动有非零分量时，则作用有沿切线方向的摩擦力。按照下式计算法向接触力
Figure 402217DEST_PATH_IMAGE114
和切向接触力
Figure 901332DEST_PATH_IMAGE115
：
Figure 263043DEST_PATH_IMAGE117
Figure 505805DEST_PATH_IMAGE119
式中：
Figure 902152DEST_PATH_IMAGE120
和/>
Figure 73632DEST_PATH_IMAGE121
分别为现时构型中物质点/>
Figure 188219DEST_PATH_IMAGE122
和/>
Figure 969093DEST_PATH_IMAGE123
的位置矢量；/>
Figure 282263DEST_PATH_IMAGE124
为物质点/>
Figure 123180DEST_PATH_IMAGE122
和/>
Figure 725062DEST_PATH_IMAGE123
之间的短程排斥力作用的临界距离，采用均匀离散形式时，可取/>
Figure 309627DEST_PATH_IMAGE125
；/>
Figure 190820DEST_PATH_IMAGE126
为排斥力常数，仅与颗粒材料类型与特征尺寸有关，可通过与试验对比方法率定；/>
Figure 937059DEST_PATH_IMAGE127
为摩擦系数，设定为为0.5。
具体的，本发明中接触临界距离
Figure 760658DEST_PATH_IMAGE124
设定为0.19 mm，略小于所述物质点的尺寸/>
Figure 680073DEST_PATH_IMAGE128
，以允许互相接触的物质点间有微小的重叠，此举可在一定程度上解决浮点数运算的精度问题
遍历所有键分别作用在两端的物质点处的键力
Figure 905518DEST_PATH_IMAGE129
和作用在相接触的物质点处的接触力/>
Figure 353817DEST_PATH_IMAGE130
，按照下式计算作用在所有物质点处的节点力/>
Figure 930291DEST_PATH_IMAGE131
：
Figure 92545DEST_PATH_IMAGE132
式中：
Figure 172496DEST_PATH_IMAGE133
为键/>
Figure 791696DEST_PATH_IMAGE134
作用在某端物质点/>
Figure 386626DEST_PATH_IMAGE135
处的键力；/>
Figure 319947DEST_PATH_IMAGE136
为与物质点/>
Figure 51142DEST_PATH_IMAGE135
相连的键/>
Figure 841244DEST_PATH_IMAGE134
的断裂状态;/>
Figure 359688DEST_PATH_IMAGE137
为作用在物质点/>
Figure 424596DEST_PATH_IMAGE135
处的接触力。
S46，更新键的断裂状态
Figure 479139DEST_PATH_IMAGE138
。若变形后所述键/>
Figure 440142DEST_PATH_IMAGE134
的伸长率/>
Figure 947347DEST_PATH_IMAGE139
大于临界伸长率/>
Figure 19208DEST_PATH_IMAGE140
，则两个物质点之间的键断开，且断开后不再承载，所述键/>
Figure 928258DEST_PATH_IMAGE134
的断裂状态/>
Figure 296048DEST_PATH_IMAGE141
；若变形后所述键/>
Figure 556128DEST_PATH_IMAGE134
的伸长率/>
Figure 166101DEST_PATH_IMAGE139
小于临界伸长率/>
Figure 991974DEST_PATH_IMAGE140
，则两个物质点之间的键是完整的，所述键/>
Figure 29200DEST_PATH_IMAGE134
的断裂状态/>
Figure 510997DEST_PATH_IMAGE142
，如图7所示。其中，键的临界伸长率与创建新的材料表面所需要的能量有关。按照下式计算所述键的临界伸长率/>
Figure 924661DEST_PATH_IMAGE140
：/>
Figure 103576DEST_PATH_IMAGE143
式中：
Figure 780545DEST_PATH_IMAGE144
为临界能量释放率；/>
Figure 546376DEST_PATH_IMAGE145
和/>
Figure 498152DEST_PATH_IMAGE146
分别代表体积模量和剪切模量；/>
Figure 501880DEST_PATH_IMAGE147
为近场范围尺寸。
S47，计算物质点的损伤值
Figure 615329DEST_PATH_IMAGE148
。物质点的损伤定义为物质点近场范围内断键数目与键总数之比，按照下式计算所述物质点的损伤值/>
Figure 337297DEST_PATH_IMAGE148
：
Figure 328650DEST_PATH_IMAGE149
式中：
Figure 186884DEST_PATH_IMAGE150
，/>
Figure 2393DEST_PATH_IMAGE151
表示材料未损伤，/>
Figure 211658DEST_PATH_IMAGE152
表示所述物质点/>
Figure 36394DEST_PATH_IMAGE153
完全损伤；
Figure 749135DEST_PATH_IMAGE136
为与物质点/>
Figure 469967DEST_PATH_IMAGE153
相连的键/>
Figure 399483DEST_PATH_IMAGE134
的断裂状态；/>
Figure 496752DEST_PATH_IMAGE103
为物质点/>
Figure 64000DEST_PATH_IMAGE153
的近场范围；/>
Figure 18049DEST_PATH_IMAGE104
为物质点/>
Figure 139589DEST_PATH_IMAGE153
的体积。
S48，更新物质点的速度
Figure 571707DEST_PATH_IMAGE154
和位移/>
Figure 462303DEST_PATH_IMAGE155
。时间步进格式采用维尔莱速度（Verlet-Velocity）差分格式，按照下式计算所述物质点的速度/>
Figure 88719DEST_PATH_IMAGE154
和位移/>
Figure 431975DEST_PATH_IMAGE155
：
Figure 136626DEST_PATH_IMAGE156
Figure 475204DEST_PATH_IMAGE157
式中：
Figure 974318DEST_PATH_IMAGE158
、/>
Figure 601609DEST_PATH_IMAGE159
和/>
Figure 578792DEST_PATH_IMAGE160
分别为当前时间步/>
Figure 240717DEST_PATH_IMAGE161
下物质点的速度、位移和所受到的节点力；/>
Figure 155408DEST_PATH_IMAGE162
、/>
Figure 269994DEST_PATH_IMAGE163
和/>
Figure 50869DEST_PATH_IMAGE164
分别为下一个时间步/>
Figure 364038DEST_PATH_IMAGE165
内物质点的速度、位移和所受到的节点力；/>
Figure 204955DEST_PATH_IMAGE166
为计算时间步长；/>
Figure 10100DEST_PATH_IMAGE167
为物质点的密度。
具体的，计算时间步长
Figure 391403DEST_PATH_IMAGE166
应小于稳定时间步长/>
Figure 263806DEST_PATH_IMAGE168
，稳定时间步长/>
Figure 275625DEST_PATH_IMAGE168
为保证计算收敛的前提下所允许的最大时间步长。稳定时间步长/>
Figure 833645DEST_PATH_IMAGE168
可由柯郎-弗里德里希斯-列维（Courant-Friedrichs-Lewy）条件（简称CFL条件）计算得到。
按照下式计算CFL条件对应的稳定时间步长
Figure 487480DEST_PATH_IMAGE168
：
Figure 447346DEST_PATH_IMAGE169
式中：
Figure 630066DEST_PATH_IMAGE170
为物质点的尺寸；/>
Figure 206540DEST_PATH_IMAGE167
为材料的密度；/>
Figure 428181DEST_PATH_IMAGE171
为材料的弹性模量。特别地，当所述物质点模型采用非均匀离散方法时，/>
Figure 508132DEST_PATH_IMAGE170
应取为物质点的最小尺寸；当所述物质点模型包含多种材料的物质点时，/>
Figure 330595DEST_PATH_IMAGE172
应取为不同种材料属性中的最小值。本发明中，计算时间步长设定为/>
Figure 925524DEST_PATH_IMAGE173
，共计算50000个时间步。
S49，重复所述步骤S42到步骤S47，直到所有计算时间步结束。
S5，确定所述近场动力学模型的参数：以步骤S1所述的单颗粒破碎室内试验结果为基础对排斥力常数
Figure 390004DEST_PATH_IMAGE174
进行率定：
所述排斥力常数
Figure 590041DEST_PATH_IMAGE174
的率定方法包括步骤：
S51，按照下式计算所述排斥力常数
Figure 380142DEST_PATH_IMAGE174
的可能取值范围：/>
Figure 400051DEST_PATH_IMAGE175
式中：
Figure 169686DEST_PATH_IMAGE174
为排斥力常数，/>
Figure 224230DEST_PATH_IMAGE176
；/>
Figure 185233DEST_PATH_IMAGE167
为颗粒材料的密度；/>
Figure 692437DEST_PATH_IMAGE170
为所述物质点的尺寸；
Figure 561036DEST_PATH_IMAGE177
为所述计算时间步长。
S52，根据模拟过程中所述单颗粒物质点模型的破碎形态进一步缩小所述排斥力常数
Figure 470086DEST_PATH_IMAGE178
的取值范围。若所述单颗粒物质点模型呈现溃散性破坏，应减小所述排斥力常数/>
Figure 601990DEST_PATH_IMAGE174
的取值；若所述单颗粒物质点模型与所述加载板物质点发生穿透现象，则应增大所述排斥力常数/>
Figure 95026DEST_PATH_IMAGE174
的取值，使得所述单颗粒物质点模型表现出合理的破碎形态；
S53，比较数值计算的所述加载板物质点的接触力的最大值-所述加载板物质点位移关系曲线与物理试验测定的荷载-位移曲线，采用试错法（trial and error）准确率定排斥力常数
Figure 704999DEST_PATH_IMAGE174
。
如图8所示。可以看出，数值模拟与室内试验得到的荷载-位移曲线较为吻合，说明了本发明提出的一种基于常规态型近场动力学的颗粒破碎模拟方法及系统的准确性。荷载-位移曲线在加载初期（a点~b点）基本呈线性增长，当加载板的位移增大到0.0385 mm（d点）时，颗粒与加载板之间的接触力达到了峰值，随后迅速减小到零，表现出明显的弹脆性。
此外，图9还给出了球形颗粒在整个破碎过程中的损伤云图。当荷载-位移曲线到达b点时，加载端出开始出现方形损伤区域；在荷载-位移曲线的b点~c点这段过程中，方形损伤区域的损伤程度不断加深；当荷载-位移曲线到达d点时，方形损伤区域的四个角点处分别开始出现向下发展的四条裂纹，裂纹不断延伸直至与下加载点发展出来的四条裂纹交汇；当荷载-位移曲线到达e点时，球形颗粒破碎成均等的四部分，且裂纹区域的损伤值仍在不断加深。
S6，获取单颗粒破碎强度模拟值：确定上、下两块加载板与颗粒之间的接触力，并将所述接触力的平均值除以所述名义粒径的平方值作为颗粒破碎强度模拟值：
步骤S6中，所述加载板与颗粒之间的接触力为加载板物质点的接触力的最大值，所述名义粒径为等效直径
Figure 734135DEST_PATH_IMAGE179
或筛析粒径/>
Figure 771361DEST_PATH_IMAGE180
，按照下式计算所述名义粒径：
Figure 784317DEST_PATH_IMAGE181
Figure 197980DEST_PATH_IMAGE182
式中：
Figure 816044DEST_PATH_IMAGE183
为所述单颗粒的体积；/>
Figure 24171DEST_PATH_IMAGE184
和/>
Figure 494729DEST_PATH_IMAGE185
分别为所述单颗粒的中轴长度和长轴长度；
Figure 712084DEST_PATH_IMAGE186
为圆周率。
具体的，所述球形单颗粒物质点模型的等效直径
Figure 450232DEST_PATH_IMAGE179
等于筛析粒径/>
Figure 94840DEST_PATH_IMAGE187
，所述名义粒径为10 mm。
S7，获取颗粒破碎后的碎片分布：采用碎片提取算法识别颗粒破碎后所产生的碎片，并对碎片的数量、尺寸和形状进行量化分析：
所述碎片提取算法通过为每个所述单颗粒物质点指定标签值来识别碎片，属于同一块碎片的物质点将共享同一个标签值。所述碎片提取算法采用了和广度优先搜索（breadth first search）类似的思想，包括以下步骤：
S71，将所有单颗粒物质点的所述标签值初始化为“0”；
S72，遍历所有单颗粒物质点，查找第一个损伤值
Figure 347967DEST_PATH_IMAGE188
小于指定损伤阈值/>
Figure 103434DEST_PATH_IMAGE189
的物质点/>
Figure 961668DEST_PATH_IMAGE190
，将所述物质点/>
Figure 10133DEST_PATH_IMAGE190
的标签值设置为当前标签值“1”，并将所述物质点/>
Figure 688239DEST_PATH_IMAGE190
压入堆栈
Figure 247397DEST_PATH_IMAGE191
；
S73，遍历与步骤S72中所述物质点的所有相连的物质点
Figure DEST_PATH_IMAGE192
，若所述相连的物质点
Figure 288034DEST_PATH_IMAGE192
的损伤值/>
Figure 274444DEST_PATH_IMAGE188
小于指定损伤阈值/>
Figure 705426DEST_PATH_IMAGE193
，则将所述相连的物质点/>
Figure 38580DEST_PATH_IMAGE192
的标签值设置为当前标签值“1”，并将所述相连的物质点/>
Figure 871407DEST_PATH_IMAGE192
压入堆栈/>
Figure 28719DEST_PATH_IMAGE191
，同时从堆栈/>
Figure 415838DEST_PATH_IMAGE191
中弹出步骤S72中所述物质点/>
Figure 316798DEST_PATH_IMAGE190
；
S74，针对步骤S73中所有所述相连的物质点
Figure 4131DEST_PATH_IMAGE192
，重复所述步骤S73，直到堆栈/>
Figure 332344DEST_PATH_IMAGE191
变空，所有所述相连的物质点/>
Figure 472339DEST_PATH_IMAGE192
的所有近邻物质点/>
Figure DEST_PATH_IMAGE194
均已被访问，并将当前标签值设置为“2”；
S75，针对尚未被访问且损伤值
Figure 483981DEST_PATH_IMAGE188
小于指定损伤阈值/>
Figure 25821DEST_PATH_IMAGE189
的物质点，重复所述步骤S72到步骤S74，直到所有单颗粒物质点都已被访问并被指定标签值。
具体的，所述相连的物质点是指物质点之间通过完整的键相连（即所述键的损伤状态
Figure 524935DEST_PATH_IMAGE195
），且所述键的未变形前的长度/>
Figure DEST_PATH_IMAGE196
小于指定键的长度阈值/>
Figure 683384DEST_PATH_IMAGE197
。本发明中，键的长度阈值/>
Figure DEST_PATH_IMAGE198
设定为物质点的尺寸的1.5倍（即/>
Figure 755508DEST_PATH_IMAGE199
），损伤阈值/>
Figure DEST_PATH_IMAGE200
设定为0.2。
图10所示为模拟得到的球形石英砂颗粒在单轴压缩荷载作用下的最终破坏形态和根据所述碎片提取算法得到的四块颗粒碎片。可以看出，四块颗粒碎片呈“蒜瓣型”分布，且四块颗粒碎片体积相近，与试验观测到的碎裂模式十分吻合，进一步说明了本发明提出的一种基于常规态型近场动力学的颗粒破碎模拟方法及系统的准确性和可靠性。
通过上述步骤S1-S7，本发明方法解决了单颗粒破碎室内试验离散性大、成本高、周期长等问题。该方法在计算上可以方便地考虑不规则颗粒的形状，所涉及的细观力学参数易于校准。此外，该方法还可用于获取准确的脆性颗粒破碎强度值，并识别颗粒破碎后所产生的碎片。进一步的，借助本发明方法开展大量不同材料、粒径和形状的颗粒的单颗粒破碎数值模拟，可以揭示潜在的单颗粒破碎细观机理。
以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (9)

1.一种基于常规态型近场动力学的颗粒破碎模拟方法，其特征在于：包括步骤：
S1，进行室内单颗粒压碎试验，测定单颗粒破碎强度试验值；
S2，建立单颗粒物质点模型：根据室内试验所用单颗粒构建三维表面形貌模型，并根据所述三维表面形貌模型构建单颗粒物质点模型，确定所述单颗粒物质点模型的名义粒径；
S3，建立近场动力学模型：在所述单颗粒物质点模型的基础上，布置加载板物质点，定义边界条件，构建基于键的链表并计算键的有关参数；
S4，进行单颗粒破碎数值试验：采用显式时间积分方案模拟所述近场动力学模型中所述单颗粒物质点模型在单轴压缩荷载作用下的破碎过程，监测所述加载板物质点中接触力的最大值随所述加载板物质点位移的变化情况；所述破碎过程包括步骤：
S41，初始化所述近场动力学模型中所有物质点的信息并开始时程积分；所述信息包括：键的断裂状态
Figure QLYQS_3
、变形后键的长度/>
Figure QLYQS_5
、键的伸长率/>
Figure QLYQS_7
、键的软化系数/>
Figure QLYQS_2
、软化标记
Figure QLYQS_4
、节点力/>
Figure QLYQS_6
、物质点的速度/>
Figure QLYQS_8
、物质点的位移/>
Figure QLYQS_1
；
S42，施加边界条件：在恒位移加载条件下，通过控制所述加载板物质点在Z方向上的运动速度以及在X方向和Y方向上的位移约束来模拟所述单颗粒破碎室内试验中加载板的压缩作用；
S43，计算变形后键的长度
Figure QLYQS_9
和键的伸长率/>
Figure QLYQS_10
：根据现时构型中键两端的物质点/>
Figure QLYQS_11
和
Figure QLYQS_12
的空间坐标计算所述变形后键的长度/>
Figure QLYQS_13
，并按照下式计算所述键的伸长率/>
Figure QLYQS_14
：
Figure QLYQS_15
；
式中：
Figure QLYQS_16
为键变形前的长度；/>
Figure QLYQS_17
为键变形后的长度；
S44，更新软化标记
Figure QLYQS_19
并计算键的软化系数/>
Figure QLYQS_23
；所述软化标记/>
Figure QLYQS_26
是大小为
Figure QLYQS_20
的一维逻辑数组，其长度设定为等于所述近场动力学模型中键的总数/>
Figure QLYQS_22
，其数值标识键的软化状态，其中逻辑值“0”代表所述键的伸长率/>
Figure QLYQS_25
大于弹性伸长率/>
Figure QLYQS_27
，键的变形进入软化阶段，逻辑值“1”代表所述键的伸长率/>
Figure QLYQS_18
小于弹性极限/>
Figure QLYQS_21
，键的变形仍处于弹性阶段；为模拟准脆性材料的渐进损伤，按照下式计算所述键的软化系数/>
Figure QLYQS_24
：
Figure QLYQS_28
；
Figure QLYQS_29
；
式中：
Figure QLYQS_30
为键的伸长率；/>
Figure QLYQS_31
为键的临界伸长率；/>
Figure QLYQS_32
为键的弹性伸长率；/>
Figure QLYQS_33
为材料的单轴抗拉强度；/>
Figure QLYQS_34
为材料的弹性模量；
S45，计算键力
Figure QLYQS_35
、接触力/>
Figure QLYQS_36
以及节点力/>
Figure QLYQS_37
：采用典型的线性近场动力学固体材料本构模拟物质点之间的非局部长程作用力，按照下式计算所述键力/>
Figure QLYQS_38
：
Figure QLYQS_39
；
Figure QLYQS_40
；
Figure QLYQS_41
；/>
Figure QLYQS_42
；
Figure QLYQS_43
；
Figure QLYQS_44
；
式中：
Figure QLYQS_50
为键的某一端物质点/>
Figure QLYQS_45
处的加权体积；/>
Figure QLYQS_58
为键的某一端物质点/>
Figure QLYQS_49
的体积膨胀率；/>
Figure QLYQS_56
为键的某一端物质点/>
Figure QLYQS_53
的近场范围；/>
Figure QLYQS_61
为键的某一端物质点/>
Figure QLYQS_51
的体积；/>
Figure QLYQS_54
为键的伸长量，分为各向同性部分/>
Figure QLYQS_46
和偏态部分/>
Figure QLYQS_59
；/>
Figure QLYQS_48
和/>
Figure QLYQS_57
分别为未变形前键的长度和变形后键的长度；/>
Figure QLYQS_52
和/>
Figure QLYQS_60
分别为刚度修正系数和体积修正系数；/>
Figure QLYQS_47
和/>
Figure QLYQS_55
分别代表体积模量和剪切模量；
采用短程排斥力模型和库伦定律描述所述单颗粒物质点与所述加载板物质点之间的相互作用；当分属不同物体的两个物质点之间的距离小于临界距离时，认为这两个物质点相互接触，沿着两个物质点的连线方向作用有大小相等，方向相反的法向接触力；当两个相接触的物质点在垂直于法线方向上的相对运动有非零分量时，则作用有沿切线方向的摩擦力；按照下式计算法向接触力
Figure QLYQS_62
和切向接触力/>
Figure QLYQS_63
：
Figure QLYQS_64
；
Figure QLYQS_65
；
式中：
Figure QLYQS_66
和/>
Figure QLYQS_69
分别为现时构型中物质点/>
Figure QLYQS_72
和/>
Figure QLYQS_67
的位置矢量；/>
Figure QLYQS_70
为物质点/>
Figure QLYQS_73
和/>
Figure QLYQS_75
之间的短程排斥力作用的临界距离，采用均匀离散形式时，可取/>
Figure QLYQS_68
；/>
Figure QLYQS_71
为排斥力常数，仅与颗粒材料类型与特征尺寸有关，可通过与试验对比方法率定；/>
Figure QLYQS_74
为摩擦系数；
遍历所有键分别作用在两端的物质点处的键力
Figure QLYQS_76
和作用在相接触的物质点处的接触力/>
Figure QLYQS_77
，按照下式计算作用在所有物质点处的节点力/>
Figure QLYQS_78
：
Figure QLYQS_79
；
式中：
Figure QLYQS_82
为键/>
Figure QLYQS_83
作用在某端物质点/>
Figure QLYQS_85
处的键力；/>
Figure QLYQS_81
为与物质点/>
Figure QLYQS_84
相连的键/>
Figure QLYQS_86
的断裂状态;/>
Figure QLYQS_87
为作用在物质点/>
Figure QLYQS_80
处的接触力；
S46，更新键的断裂状态
Figure QLYQS_88
：若变形后所述键/>
Figure QLYQS_95
的伸长率/>
Figure QLYQS_98
大于临界伸长率/>
Figure QLYQS_90
，则两个物质点之间的键断开，且断开后不再承载，所述键/>
Figure QLYQS_92
的断裂状态/>
Figure QLYQS_96
；若变形后所述键/>
Figure QLYQS_99
的伸长率/>
Figure QLYQS_91
小于临界伸长率/>
Figure QLYQS_93
，则两个物质点之间的键是完整的，所述键/>
Figure QLYQS_97
的断裂状态/>
Figure QLYQS_100
；其中，键的临界伸长率与创建新的材料表面所需要的能量有关；按照下式计算所述键/>
Figure QLYQS_89
的临界伸长率/>
Figure QLYQS_94
：
Figure QLYQS_101
；
式中：
Figure QLYQS_102
为临界能量释放率；/>
Figure QLYQS_103
为近场范围尺寸；
S47，计算物质点的损伤值
Figure QLYQS_104
；物质点的损伤定义为物质点近场范围内断键数目与键总数之比，按照下式计算所述物质点的损伤值/>
Figure QLYQS_105
：
Figure QLYQS_106
；
式中：
Figure QLYQS_107
，/>
Figure QLYQS_111
表示材料未损伤，/>
Figure QLYQS_114
表示所述物质点/>
Figure QLYQS_108
完全损伤；/>
Figure QLYQS_110
为与物质点/>
Figure QLYQS_113
相连的键/>
Figure QLYQS_116
的断裂状态；/>
Figure QLYQS_109
为物质点/>
Figure QLYQS_112
的近场范围；/>
Figure QLYQS_115
为物质点/>
Figure QLYQS_117
的体积；
S48，更新物质点的速度
Figure QLYQS_118
和位移/>
Figure QLYQS_119
；时间步进格式采用维尔莱速度差分格式，按照下式计算所述物质点的速度/>
Figure QLYQS_120
和位移/>
Figure QLYQS_121
：
Figure QLYQS_122
；
Figure QLYQS_123
；
式中：
Figure QLYQS_125
、/>
Figure QLYQS_127
和/>
Figure QLYQS_130
分别为当前时间步/>
Figure QLYQS_126
下物质点的位移、速度和所受到的节点力；/>
Figure QLYQS_129
、/>
Figure QLYQS_132
和/>
Figure QLYQS_133
分别为下一个时间步/>
Figure QLYQS_124
内物质点的速度、位移和所受到的节点力；/>
Figure QLYQS_128
为计算时间步长；/>
Figure QLYQS_131
为物质点的密度；
S49，重复所述步骤S42到步骤S47，直到所有计算时间步结束；
S5，确定所述近场动力学模型的参数：以步骤S1所述的单颗粒破碎室内试验结果为基础对排斥力常数
Figure QLYQS_134
进行率定；
S6，获取单颗粒破碎强度模拟值：确定上、下两块加载板与颗粒之间的接触力，并将所述接触力的平均值除以所述名义粒径的平方值作为颗粒破碎强度模拟值；
S7，获取颗粒破碎后的碎片分布：采用碎片提取算法识别颗粒破碎后所产生的碎片，并对碎片的数量、尺寸和形状进行量化分析。
2.根据权利要求1所述的一种基于常规态型近场动力学的颗粒破碎模拟方法，其特征在于：步骤S2中，单颗粒三维表面形貌通过三维扫描方式获取。
3.根据权利要求1所述的一种基于常规态型近场动力学的颗粒破碎模拟方法，其特征在于：步骤S2中，所述单颗粒物质点模型为由所述单颗粒三维表面形貌模型通过体素化离散而成的物质点集合，物质点的尺寸小于单颗粒直径值的十分之一。
4.根据权利要求1所述的一种基于常规态型近场动力学的颗粒破碎模拟方法，其特征在于：步骤S3中，建立所述近场动力学模型包括步骤：
S31，在所述单颗粒物质点模型的顶部和底部对称布置一定数量的加载板物质点，加载板的面积应大于所述单颗粒物质点模型在水平面上的投影面积，加载板的厚度应大于近场范围尺寸
Figure QLYQS_135
，以确保所施加的荷载和约束条件能在所述单颗粒物质点模型中得到准确地反映；
S32，定义一个大小为
Figure QLYQS_136
的二维逻辑数组作为位移约束标记/>
Figure QLYQS_137
，其长度设定为等于所述近场动力学模型中物质点的总数/>
Figure QLYQS_138
，其第一列、第二列和第三列数值分别标识物质点在X方向、Y方向和Z方向上的位移约束情况，其中逻辑值“0”代表所述物质点在所述方向上的位移不受约束，逻辑值“1”代表所述物质点在所述方向上的位移受到限制；将所述加载板物质点的所述位移约束标记/>
Figure QLYQS_139
的第一列和第二列数值设定为“0”，以确保模拟过程中加载板不发生水平面内的相对滑移；
S33，采用树形搜索算法确定所述近场动力学模型中所有物质点的近邻物质点集合
Figure QLYQS_140
，并创建列表/>
Figure QLYQS_141
；
S34，计算键的有关参数，所述键的有关参数包括未变形前键的长度
Figure QLYQS_142
、体积修正系数/>
Figure QLYQS_143
、键的类型以及刚度修正系数/>
Figure QLYQS_144
。
5.根据权利要求1所述的一种基于常规态型近场动力学的颗粒破碎模拟方法，其特征在于：步骤S5中，所述排斥力常数
Figure QLYQS_145
的率定方法包括步骤：
S51，按照下式计算所述排斥力常数
Figure QLYQS_146
的可能取值范围：
Figure QLYQS_147
；
式中：
Figure QLYQS_148
为排斥力常数，/>
Figure QLYQS_149
；/>
Figure QLYQS_150
为颗粒材料的密度；/>
Figure QLYQS_151
为所述物质点的尺寸；/>
Figure QLYQS_152
为所述计算时间步长；
S52，根据模拟过程中所述单颗粒物质点模型的破碎形态进一步缩小所述排斥力常数
Figure QLYQS_153
的取值范围；若所述单颗粒物质点模型呈现溃散性破坏，应减小所述排斥力常数/>
Figure QLYQS_154
的取值；若所述单颗粒物质点模型与所述加载板物质点发生穿透现象，则应增大所述排斥力常数/>
Figure QLYQS_155
的取值，使得所述单颗粒物质点模型表现出合理的破碎形态；
S53，比较数值计算的所述加载板物质点的接触力的最大值-所述加载板物质点位移关系曲线与物理试验测定的荷载-位移曲线，采用试错法准确率定排斥力常数
Figure QLYQS_156
。
6.根据权利要求1所述的一种基于常规态型近场动力学的颗粒破碎模拟方法，其特征在于：步骤S6中，所述加载板与颗粒之间的接触力为加载板物质点的接触力的最大值，所述名义粒径为等效直径
Figure QLYQS_157
或筛析粒径/>
Figure QLYQS_158
，按照下式计算所述名义粒径：
Figure QLYQS_159
；
Figure QLYQS_160
；
式中：
Figure QLYQS_161
为所述单颗粒的体积；/>
Figure QLYQS_162
和/>
Figure QLYQS_163
分别为所述单颗粒的中轴长度和长轴长度；/>
Figure QLYQS_164
为圆周率。
7.根据权利要求1所述的一种基于常规态型近场动力学的颗粒破碎模拟方法，其特征在于：步骤S7中，所述碎片提取算法通过为每个所述单颗粒物质点指定标签值来识别碎片，属于同一块碎片的物质点将共享同一个标签值；所述碎片提取算法包括以下步骤：
S71，将所有单颗粒物质点的所述标签值初始化为“0”；
S72，遍历所有单颗粒物质点，查找第一个损伤值
Figure QLYQS_165
小于指定损伤阈值/>
Figure QLYQS_166
的物质点/>
Figure QLYQS_167
，将所述物质点/>
Figure QLYQS_168
的标签值设置为当前标签值“1”，并将所述物质点/>
Figure QLYQS_169
压入堆栈/>
Figure QLYQS_170
；
S73，遍历与步骤S72中所述物质点的所有相连的物质点
Figure QLYQS_172
，若所述相连的物质点/>
Figure QLYQS_176
的损伤值/>
Figure QLYQS_178
小于指定损伤阈值/>
Figure QLYQS_173
，则将所述相连的物质点/>
Figure QLYQS_175
的标签值设置为当前标签值“1”，并将所述相连的物质点/>
Figure QLYQS_177
压入堆栈/>
Figure QLYQS_179
，同时从堆栈/>
Figure QLYQS_171
中弹出步骤S72中所述物质点/>
Figure QLYQS_174
；
S74，针对步骤S73中所有所述相连的物质点
Figure QLYQS_180
，重复所述步骤S73，直到堆栈/>
Figure QLYQS_181
变空，所有所述相连的物质点/>
Figure QLYQS_182
的所有近邻物质点/>
Figure QLYQS_183
均已被访问，并将当前标签值设置为“2”；
S75，针对尚未被访问且损伤值
Figure QLYQS_184
小于指定损伤阈值/>
Figure QLYQS_185
的物质点，重复所述步骤S72到步骤S74，直到所有单颗粒物质点都已被访问并被指定标签值。
8.根据权利要求1所述的一种基于常规态型近场动力学的颗粒破碎模拟方法，其特征在于：步骤S34中，
根据参考构型中键两端的物质点的空间坐标，按照下式计算所述未变形前键的长度
Figure QLYQS_186
；/>
Figure QLYQS_187
；
式中：
Figure QLYQS_189
、/>
Figure QLYQS_191
、/>
Figure QLYQS_193
和/>
Figure QLYQS_190
、/>
Figure QLYQS_192
、/>
Figure QLYQS_194
、分别为参考构型中物质点/>
Figure QLYQS_195
和/>
Figure QLYQS_188
的横坐标值、纵坐标值和竖坐标值；
为了减少邻域截断误差的影响，按照下式计算所述体积修正系数
Figure QLYQS_196
：
Figure QLYQS_197
；
式中：
Figure QLYQS_198
为未变形前键的长度；/>
Figure QLYQS_199
为近场范围尺寸；/>
Figure QLYQS_200
为物质点的直径；
按照下式计算所述刚度修正系数
Figure QLYQS_201
：
Figure QLYQS_202
；
Figure QLYQS_203
；
式中：
Figure QLYQS_204
为近场范围的体积；/>
Figure QLYQS_205
和/>
Figure QLYQS_206
分别为所述物质点/>
Figure QLYQS_207
和/>
Figure QLYQS_208
的近邻物质点的体积；/>
Figure QLYQS_209
为圆周率。
9.一种基于常规态型近场动力学的颗粒破碎模拟系统，其特征为：所述系统包括权利要求1所述的一种基于常规态型近场动力学的颗粒破碎模拟方法。
CN202210411866.6A 2022-04-19 2022-04-19 一种基于常规态型近场动力学的颗粒破碎模拟方法及系统 Active CN114613458B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210411866.6A CN114613458B (zh) 2022-04-19 2022-04-19 一种基于常规态型近场动力学的颗粒破碎模拟方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210411866.6A CN114613458B (zh) 2022-04-19 2022-04-19 一种基于常规态型近场动力学的颗粒破碎模拟方法及系统

Publications (2)

Publication Number Publication Date
CN114613458A CN114613458A (zh) 2022-06-10
CN114613458B true CN114613458B (zh) 2023-06-02

Family

ID=81870168

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210411866.6A Active CN114613458B (zh) 2022-04-19 2022-04-19 一种基于常规态型近场动力学的颗粒破碎模拟方法及系统

Country Status (1)

Country Link
CN (1) CN114613458B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116705200B (zh) * 2023-04-27 2024-05-03 安徽理工大学 一种基于分形理论的聚变堆球床颗粒材料破碎仿真计算方法
CN117291060B (zh) * 2023-11-23 2024-02-27 成都理工大学 一种考虑动力碎裂效应的岩崩运动过程三维模拟预测方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109766603A (zh) * 2018-12-26 2019-05-17 哈尔滨工程大学 一种冰船接触过程数值计算方法
CN111153689A (zh) * 2019-12-23 2020-05-15 西南交通大学 一种用于轨道交通轮轨界面增黏的增黏颗粒及制备方法
CN112417715A (zh) * 2020-10-26 2021-02-26 山东大学 一种真三轴伺服加载状态下的岩体破裂模拟方法及系统

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109655342B (zh) * 2019-01-29 2021-04-06 西南交通大学 一种微小颗粒试样破碎过程研究实验装置及其实验方法
CN111125920A (zh) * 2019-12-26 2020-05-08 湘潭大学 一种基于蒙特卡罗模拟的颗粒破碎强度确定方法
CN111339594B (zh) * 2020-02-26 2022-12-16 河北工业大学 基于dic技术的近场动力学参数实验反演系统及使用方法
CN112380736B (zh) * 2020-12-14 2022-03-25 中南大学 一种实现单颗粒破碎强度weibull分布精准控制的离散元模型构建方法
CN113049413A (zh) * 2021-03-25 2021-06-29 中国矿业大学 一种基于分形理论与落锤实验表征岩石坚固性的方法
CN113821926B (zh) * 2021-09-18 2024-01-26 塔里木大学 一种坚果果壳破碎的离散元数值模拟方法
CN114186456B (zh) * 2021-12-02 2022-06-14 大连理工大学 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109766603A (zh) * 2018-12-26 2019-05-17 哈尔滨工程大学 一种冰船接触过程数值计算方法
CN111153689A (zh) * 2019-12-23 2020-05-15 西南交通大学 一种用于轨道交通轮轨界面增黏的增黏颗粒及制备方法
CN112417715A (zh) * 2020-10-26 2021-02-26 山东大学 一种真三轴伺服加载状态下的岩体破裂模拟方法及系统

Also Published As

Publication number Publication date
CN114613458A (zh) 2022-06-10

Similar Documents

Publication Publication Date Title
CN114613458B (zh) 一种基于常规态型近场动力学的颗粒破碎模拟方法及系统
Khatir et al. Structural health monitoring using modal strain energy damage indicator coupled with teaching-learning-based optimization algorithm and isogoemetric analysis
Sun et al. Topological configuration analysis and design for foam filled multi-cell tubes
TWI344611B (en) Method of simulating deformation of rubber material
Liu et al. DEM investigation of the effect of intermediate principle stress on particle breakage of granular materials
Jin et al. Stability prediction for asphalt mixture based on evolutional characterization of aggregate skeleton
Lum et al. Structural optimization for flexure-based parallel mechanisms–Towards achieving optimal dynamic and stiffness properties
JP4329848B2 (ja) 構造解析装置および構造解析方法
Jha et al. Peridynamics-based discrete element method (PeriDEM) model of granular systems involving breakage of arbitrarily shaped particles
Wang et al. A 3D Voronoi clump based model for simulating failure behavior of brittle rock
Müller et al. A hyperelastic model for simulating cells in flow
Peng et al. Performance analysis and multi-objective optimization of bionic dendritic furcal energy-absorbing structures for trains
CN112258643A (zh) 岩土边坡任意形状落石运动轨迹三维分析方法
CN111259593B (zh) 基于edem实现岩石破碎过程可视化的方法
Nader et al. Grain breakage under uniaxial compression using a three-dimensional discrete element method
Shin et al. On simulating impact fracture in high-strength concrete using GraFEA
Suwa et al. Cluster morphology as a test of different cosmological models
Yang et al. Fragmentation analyses of rocks under high-velocity impacts using the combined finite-discrete element simulation
Shi et al. Multi-scale analysis of shear behaviour of crushable granular sand under general stress conditions
Zhang et al. Cavitation, crazing and bond scission in chemically cross-linked polymer nanocomposites
Ujihara et al. A mechanical cell model and its application to cellular biomechanics
Bhattacharjee et al. Predicting high rate granular transition and fragment statistics at the onset of granular flow for brittle ceramics
Feng et al. A GPU based Hybrid Material point and Discrete element method (MPDEM) algorithm and validation
Jahanfekr et al. Insight to damage identification in truss-type structures using a second-order gradient-based algorithm
Silling Data‐driven nonlocal model for fragmentation in the crushing of solids

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant