CN116364218B - 基于近场动力学岩石材料率效应本构模型构建与实施方法 - Google Patents
基于近场动力学岩石材料率效应本构模型构建与实施方法 Download PDFInfo
- Publication number
- CN116364218B CN116364218B CN202310324881.1A CN202310324881A CN116364218B CN 116364218 B CN116364218 B CN 116364218B CN 202310324881 A CN202310324881 A CN 202310324881A CN 116364218 B CN116364218 B CN 116364218B
- Authority
- CN
- China
- Prior art keywords
- elongation
- break
- field
- tensile
- bond
- 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
Links
- 239000000463 material Substances 0.000 title claims abstract description 94
- 239000011435 rock Substances 0.000 title claims abstract description 58
- 230000000694 effects Effects 0.000 title claims abstract description 54
- 238000000034 method Methods 0.000 title claims abstract description 38
- 238000010276 construction Methods 0.000 title claims abstract description 10
- 239000000126 substance Substances 0.000 claims abstract description 43
- 230000006378 damage Effects 0.000 claims abstract description 30
- 238000009826 distribution Methods 0.000 claims abstract description 21
- 238000013016 damping Methods 0.000 claims abstract description 20
- 230000009471 action Effects 0.000 claims abstract description 11
- 230000008859 change Effects 0.000 claims abstract description 8
- 230000006835 compression Effects 0.000 claims description 49
- 238000007906 compression Methods 0.000 claims description 49
- 230000003993 interaction Effects 0.000 claims description 11
- 238000006073 displacement reaction Methods 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 6
- 230000010354 integration Effects 0.000 claims description 6
- 230000001133 acceleration Effects 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 5
- 230000006872 improvement Effects 0.000 claims description 3
- 230000007547 defect Effects 0.000 abstract description 3
- 239000006185 dispersion Substances 0.000 description 7
- 239000002245 particle Substances 0.000 description 5
- 230000006399 behavior Effects 0.000 description 3
- 230000000977 initiatory effect Effects 0.000 description 3
- 238000009417 prefabrication Methods 0.000 description 3
- 238000005336 cracking Methods 0.000 description 2
- 239000013013 elastic material Substances 0.000 description 2
- 238000005728 strengthening Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 208000027418 Wounds and injury Diseases 0.000 description 1
- 238000009412 basement excavation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005422 blasting Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 230000001687 destabilization Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000005489 elastic deformation Effects 0.000 description 1
- 208000014674 injury Diseases 0.000 description 1
- 238000012821 model calculation Methods 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000012887 quadratic function Methods 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
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
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- 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)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computing Systems (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
本发明涉及一种基于近场动力学岩石材料率效应本构模型构建与实施方法，构建方法步骤为：引入能够反映物质点作用强度随物质点间距变化的核函数，对近场动力学本构力函数进行改进；求解改进后本构力函数的应变能密度；令近场动力学中的应变能密度和传统连续介质力学的应变能密度相等，求得修正后的微模量常数；引入键的损伤函数；引入非局部粘性阻尼，构造近场动力学率效应本构模型；将微模量常数以及损伤函数代入近场动力学率效应本构模型，完成岩石材料率效应本构模型的构建。本发明通过考虑物质点的拉伸断裂伸长率和压缩断裂伸长率服从威布尔分布来描述岩石材料的非均质性，弥补了传统键基近场动力学无法反映岩石材料异质性的不足。
Description
技术领域
本发明涉及岩土工程领域，具体涉及一种基于近场动力学岩石材料率效应本构模型构建与实施方法。
背景技术
岩石作为一种天然介质体，内部存在大量方向各异的微裂纹，岩石的破裂通常是由预先存在的裂缝起裂、扩展及相互贯穿导致的。在隧道开挖、巷道掘进、油气开采等岩体工程实践中，爆破、地震波、矿震等各种动态荷载会对地下岩体结构产生影响，进而诱发地下岩体结构的失稳破坏，影响岩体的稳定性。因此，研究岩石动态荷载作用下的力学特性、变形特征和断裂行为，对岩体工程建设中动力灾害防治、围岩支护设计均具有重要的指导意义。
近些年来，一种新兴的基于非局部思想的无网格方法—近场动力学，被广泛用于模拟裂纹扩展问题。近场动力学通过积分形式求解物质点的运动方程，避免了不连续处局部空间导数不存在的问题，同时无需外部定义的准则，就可以自发地模拟裂纹起裂、扩展与分叉，广泛的应用于材料断裂行为的研究。
近场动力学按照本构力的求解形式可分为“键”基、常规“态”基和非常规“态”基三种。其中，键基近场动力学的微观弹脆性模型是应用最广泛的一种，该模型描述的是脆性断裂，当物质点之间键的伸长率超过临界伸长率，物质点之间的相互作用力就消失，不适合模拟岩石材料的破坏。此外，传统键基近场动力学不能反映物质点内部尺寸对物质点作用强度的影响。岩石材料在冲击荷载作用下的变形特征表现为明显的率效应，传统近场动力学不能准确描述应变率对材料力学特性及变形特征的影响。目前，部分学者采用不同的损伤函数和核函数对传统近场动力学模型进行改进，但是并未考虑应变率对岩石材料力学特性的影响，也并未考虑岩石材料的异质性，导致传统近场动力学不适合模拟冲击荷载作用下岩石的动态断裂过程。
发明内容
针对现有技术的不足，本发明的目的是提供一种基于近场动力学岩石材料率效应本构模型构建与实施方法，能够有效模拟岩石材料在冲击荷载作用下裂纹起裂、扩展、贯穿的全过程。
为了实现上述目的，本发明采用如下技术方案实现：
本发明一方面提出一种近场动力学岩石材料率效应本构模型的构建方法，步骤为：
S11引入能够反映物质点作用强度随物质点间距变化的核函数，对近场动力学本构力函数进行改进；
S12求解改进后本构力函数的应变能密度WPD(X)；
S13令近场动力学中的应变能密度WPD(X)和传统连续介质力学的应变能密度相等，求得修正后的微模量常数c(0,δ)；
S14引入键的损伤函数D(s)；
S15引入非局部粘性阻尼，构造近场动力学率效应本构模型；
S16将步骤S13中求得的微模量常数c(0,δ)以及步骤S14中键的损伤函数D(s)代入步骤S15中的近场动力学率效应本构模型，完成岩石材料率效应本构模型的构建。
进一步的，步骤S11中改进后本构力函数为：
其中核函数的表达形式为：
其中，g(ξ,δ)为反映物质点作用强度随物质点间距变化的核函数，δ为领域半径，ξ为物质点之间的初始距离。
进一步的，所述步骤S12中应变能密度WPD(X)为：
对于三维问题，应变能密度表示为
其中，δ为领域半径，c(0,δ)为微模量常数，s表示键的伸长率，ξ为物质点之间的初始距离；
对于二维问题，应变能密度表示为
其中，h为厚度，δ为领域半径，c(0,δ)为微模量常数，s表示键的伸长率；对于一维问题，应变能密度表示为
其中，h为厚度，δ为领域半径，A为面积，c(0,δ)为微模量常数，s表示键的伸长率。
进一步的，步骤S13中的微模量为：
式中，E表示弹性模量，δ为领域半径，h表示厚度，A表示面积。
进一步的，步骤S14中键的损伤函数可表示为：
其中，sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率，s表示键的伸长率。
进一步的，步骤S15中所述的率效应本构模型可表示为
其中，υ表示粘性阻尼系数，表示键的变形率，c(0,δ)为微模量常数；
式中，表示当前构型键的速度，η+ξ表示当前构型的键，||η+ξ||表示当前构型键的长度，||ξ||表示初始构型键的长度。
进一步的，步骤S16中所述的岩石材料率效应本构模型可表示为：
对于三维问题，
其中，E表示弹性模量；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量；/>表示当前构型键方向；
对于平面应力问题，
其中，E表示弹性模量；h表示厚度；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量；
对于平面应变问题，
其中，E表示弹性模量；h表示厚度；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量；/>表示当前构型键方向。
本发明另一方面提出一种近场动力学岩石材料率效应本构模型的实施方法，步骤为：
A1：数值模型参数初始化，输入材料力学参数，确定数值模型尺寸；
A2：确定近场动力学中物质点间距，领域半径，离散求解域，生成物质点坐标；
A3：搜索物质点领域范围内的点并编号，初始化所有物质点在其领域范围内的键；
A4：对每个物质点拉伸断裂伸长率和压缩断裂伸长率按威布尔分布赋值；
A5：根据实际岩石试样，判断步骤A1中的数值模型是否需要预制裂纹；
A6：对数值模型施加初始边界条件；
A7：确定数值模型的时间步长；
A8：计算所有物质点的本构力；计算每个物质点领域范围内所有键的本构力并求和；
A9：判断伸长率是否大于拉伸断裂伸长率或小于压缩断裂伸长率，如果“是”，则物质点之间的键断裂，物质点之间的相互作用力消失；如果“否”，则按步骤A8计算物质点之间的本构力；
A10：通过物质点本构力计算出物质点的加速度，利用显式向前和向后差分公式进行时域积分，求出位移和速度。
A11：输出计算结果，绘制损伤云图。
进一步的，步骤A4中所述服从威布尔分布的拉伸断裂伸长率可表示为：
其中，st表示拉伸断裂伸长率，表示拉伸断裂伸长率的平均值，m为拉伸断裂伸长率服从威布尔分布的形状参数，r(x)表示(0,1]之间均匀分布的随机数。与现有技术相比，本发明的有益效果是：
(1)提出一种可以反映物质点内部长度对非局部长程力影响的核函数，对传统键基近场动力学进行改进，可以减少波的数值色散，提高计算精度。
(2)传统键基近场动力学描述的是脆性断裂，当物质点之间键的伸长超过临界伸长率，物质点之间的相互作用力就消失，不适合模拟岩石材料的破坏，通过引入键的损伤函数，来描述键的非线性力学行为，更符合岩石材料的力学特性。
(3)岩石材料在冲击荷载作用下的变形特征表现为明显的率效应，传统键基近场动力学不适合模拟冲击荷载作用下岩石材料的变形破坏。通过在物质点之间引入非局部粘性阻尼，实现率效应本构模型。物质点之间本构力和键的伸长量关系可以看成是由损伤体和粘性体并联而成，粘性体不具有损伤特性。可以模拟岩石材料在冲击荷载作用下的变形破坏特征。
(4)通过考虑物质点的拉伸断裂伸长率和压缩断裂伸长率服从威布尔分布来描述岩石材料的非均质性，弥补了传统键基近场动力学无法反映岩石材料异质性的不足。
附图说明
图1为近场动力学物质点运动原理图；
图2为PMB材料和岩石材料本构力函数图；
图3为角频率和波数的频散曲线图；
图4为率效应本构模型图；
图5为岩石材料率效应本构模型计算流程图；
图6为单边切口三点弯梁模型图；
图7为拉伸断裂率和压缩断裂率随机分布图；
图8为冲击荷载作用下单边切口三点弯梁裂纹扩展图。
具体实施方式
以下结合本发明实施例中的附图来对本发明的技术方案进行完整清晰的描述。
实施例1
如图1所示，在近场动力学中，将连续体离散成有限个带有质量、体积、密度的物质点，每个物质点和其领域范围内的物质点存在相互作用，和其领域范围外的物质点不存在相互作用。物质点之间的相互作用力用本构力函数f来表示。根据牛顿第二定律，任一时刻物质点的运动方程可表示为：
其中，x′为x近场范围内的物质点；u(x′,t)和u(x,t)分别表示物质点x′和x在t时刻的位移；ρ(x)和分别表示物质点的密度和加速度；dV′x为物质点x′的体积微元；H＝{x′∈R,||x′-x||＜δ}为物质点x近场范围内物质点x′的集合；b(x,t)为外荷载密度；f表示与材料本身属性相关的本构力函数，对于弹脆性材料，其表达式为：
其中，s0表示临界伸长率。
如图2所示，传统键基近场动力学中的本构力函数描述的是脆性材料的破坏，当物质点之间的伸长率超过临界伸长率，物质点之间的键断裂，物质点之间的相互作用消失，不适合模拟准脆性材料的破坏。此外，岩石材料在冲击荷载作用下的力学特性和变形特征表现为明显的率效应。传统弹脆性材料并未考虑应变率的变化对键变形特征的影响。据此，本发明提出了一种基于近场动力学岩石材料率效应本构模型构建方法，包括以下步骤：
S11引入能够反映物质点作用强度随物质点间距变化的核函数，对近场动力学本构力函数进行改进；
S12推导改进后本构力函数的应变能密度；
S13令近场动力学中的应变能密度和传统连续介质力学的应变能密度相等，求得修正后的微模量常数；
S14引入键的损伤函数；
S15引入非局部粘性阻尼，构造近场动力学率效应本构模型；
S16将步骤S13中求得的微模量常数以及步骤S14中键的损伤函数代入步骤S15中的近场动力学率效应本构模型，完成岩石材料率效应本构模型的构建。
步骤S11中所述的核函数的表达形式为：
其中，δ为领域半径，ξ为物质点之间的初始距离。
步骤S11中所述的改进型本构力函数可以减少波的数值色散，具体包括以下步骤：
S11.1表示出平面波振动方程的一般形式：
u＝u0ei(kx-ωt) (4)
其中，u0表示振幅，k表示波速，ω表示角频率。
S11.2表示出近场动力学物质点的运动方程：
S11.3近场动力学本构力函数可表示成
S11.4不考虑损伤和体力，步骤S11.2所述的物质点运动方程可表示成：
S11.5将步骤S11.1所述平面波振动方程一般形式带入S11.4所述的物质点运动方程，可表示为：
S11.6进一步地，物质点运动方程可以化简为：
S11.7对步骤S11.6中所述物质点运动方程两边同乘e-ikx可表示为：
进一步地，步骤S11.7所述物质点运动方程可表示为
S11.8对步骤S11.7中所述的物质点运动方程中的指数函数用欧拉公式展开，可表示为
进一步地，因为g(ξ,δ)关于积分区间Hx是偶函数，sin(k(x′-x))关于积分区间Hx是奇函数，所以
S11.9将步骤S11.8中所述的物质点运动方程表示成波速和角频率的关系：
S11.10对步骤S11.9中所述的波速和角频率的关系进行数值离散，令δ＝m|Δx|，ξ＝j|Δx|，则角频率和波数的关系可表示为
图3为应用本发明计算出的角频率和波数的频散曲线图，其特点为基于本发明所提出的核函数，考虑了物质点间距变化对物质点作用强度的影响，对传统键基近场动力学进行改进。图中I为采用本发明方法改进型近场动力学，m＝2的曲线图；图中II为采用现有近场动力学，m＝2的曲线图；图中III为采用本发明方法改进型近场动力学，m＝3的曲线图；图中IV为采用现有近场动力学，m＝3的曲线图；图中V为采用本发明方法改进型近场动力学，m＝4的曲线图；图中VI为采用现有近场动力学，m＝4的曲线图；在本实施例中，Δx＝0.001m，两个数值模型的数值频散随着领域半径的增加而加剧，数值频散随着波数的增加、波长的减少和角频率的增加而加剧。同时也可以发现对于相同的领域半径，改进型近场动力学比传统近场动力学更接近理论解，本实施例能够有效减少数值频散。
所述步骤S12包括以下步骤：
S12.1确定改进后近场动力学本构力函数标量的表达形式：
S12.2确定改进后近场动力学的微势能函数：
S12.3对步骤S12.2中本构力函数进行积分，求得应变能密度WPD(X)。
对于三维问题，应变能密度表示为
其中，δ为领域半径，c(0,δ)为微模量常数。
对于二维问题，应变能密度表示为
其中，h为厚度，δ为领域半径。
对于一维问题，应变能密度表示为
其中，h为厚度，δ为领域半径。
步骤S13中的微模量常数通过推导可得
本构力函数非线性段是一个二次函数，将本构力函数表示成
f＝cs[1-D(s)]
步骤S14中所述的键的损伤函数可表示为：
其中，sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率。
借鉴粘弹性力学中的Kelvin-Voigt模型，在PMB材料模型中引入非局部粘性阻尼，则步骤S15中所述的率效应本构模型可表示为
其中，υ表示粘性阻尼系数，表示键的变形率。
在步骤S16所述的本构力函数中，当s1c≤s≤s1t时，键处于线弹性变形阶段，点对力与伸长率的关系为线性关系，键无损伤；当s1t＜s＜st或s1c＜s＜sc时，键处于非线性变形阶段。对于非线性变形阶段又可分为非线性强化变形阶段和非线性软化变形阶段，当s1t＜s＜s2t或s1c＜s＜s2c时，键处于非线性强化变形阶段，键开始损伤，微裂纹开始产生，点对力随着伸长率增加非线性增加；当s2t＜s＜st或s2c＜s＜sc时，键处于非线性软化阶段，点对力随着伸长率增加非线性减小；当s＞st或s＜sc时，键的损伤为1，键发生断裂，物质点之间不存在相互作用。
如图4所示，步骤S16中所述的岩石材料率效应本构模型通过在物质点之间引入非局部粘性阻尼，实现率效应本构模型。物质点之间本构力和键伸长量关系可以看成是由损伤体和粘性体并联而成，粘性体不具有损伤特性。可以模拟岩石材料在冲击荷载作用下的变形破坏特征。
对于三维问题，步骤S16中所述的岩石材料率效应本构模型可表示为：
其中，E表示弹性模量；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量。
对于平面应力问题，步骤S16中所述的岩石材料率效应本构模型可表示为
其中，E表示弹性模量；h表示厚度；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量。
对于平面应变问题，步骤S16中所述的岩石材料率效应本构模型可表示为：
其中，E表示弹性模量；h表示厚度；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量。
实施例2
如图5所示，本发明还提供一种基于近场动力学岩石材料率效应本构模型实施方法，包括以下步骤：
A1：数值模型参数初始化，输入材料力学参数，确定数值模型尺寸；
以平面应力问题为例，根据实际情况，确定数值模型(试件模型)长度L、宽度W、弹性模量E，密度ρ，泊松比μ，压缩断裂伸长率sc；压缩临界伸长率s2c；压缩弹性伸长率s1c；拉伸断裂伸长率st；拉伸临界伸长率s2t；拉伸弹性伸长率s1t，粘性阻尼系数υ；
A2：确定近场动力学中物质点间距，领域半径，离散求解域，生成物质点坐标；
数值模型长度为L，宽度为W，沿长度方向分布的物质点个数为n，沿宽度方向分布物质点个数为m，物质点间距物质点体积Vx＝(Δx)3，领域半径δ＝3.015Δx。
假设沿长度方向为x轴方向，沿宽度方向为y轴方向，则物质点的坐标可表示为：
其中，i＝1，2，…，n，表示沿x方向排列晶格的个数。j＝1，2，…，m，表示沿y方向排列晶格的个数。
A3：搜索物质点领域范围内的点并编号，初始化所有物质点在其领域范围内的键；
遍历所有物质点，计算物质点x和物质点x′之间的距离，对满足H＝{x′∈R,||x′-x||＜δ}的物质点编号。构建所有物质点近场领域范围的键并对其初始化，μ(x,t,ξ)＝1，
其中，μ(x,t,ξ)是一个标量函数，用来判断键是否断裂，u(x,t,ξ)＝1，键完好，u(x,t,ξ)＝0，键断裂。
物质点的损伤程度用表示，根据其领域内断裂键个数和原总键数所占的比值来计算：
A4：对每个物质点拉伸断裂伸长率和压缩断裂伸长率按威布尔分布赋值；
假设拉伸断裂伸长率服从威布尔分布：
其中，st表示拉伸断裂伸长率，表示拉伸断裂伸长率的平均值，m为拉伸断裂伸长率服从威布尔分布的形状参数。
则服从威布尔分布的拉伸断裂伸长率可表示为：
其中，st表示拉伸断裂伸长率，表示拉伸断裂伸长率的平均值，m为拉伸断裂伸长率服从威布尔分布的形状参数。
假设压缩断裂伸长率服从威布尔分布：
其中，sc表示压缩断裂伸长率，表示压缩断裂伸长率的平均值，n为压缩断裂伸长率服从威布尔分布的形状参数。
则服从威布尔分布的压缩断裂伸长率可表示为：
其中，sc表示拉伸断裂伸长率，表示拉伸断裂伸长率的平均值，n为拉伸断裂伸长率服从威布尔分布的形状参数。
A5：根据实际岩石试样，判断步骤A1中的数值模型是否需要预制裂纹；
根据实际模型进行预制裂纹，比如模拟预制裂隙岩石试样，就需要预制裂隙。如果需要预制裂纹，确定裂纹中心线所在坐标，将连接裂纹中心线的键断裂，确定裂纹表面。
A6：对数值模型施加初始边界条件；
A7：确定数值模型的时间步长；
时间步长按下式进行计算：
其中，ρ为密度，c为微模量。
A8：计算所有物质点的本构力；计算每个物质点领域范围内所有键的本构力并求和；
以平面应力问题为例，任意时刻物质点之间键的本构力可表示为：
任意物质点总的本构力可表示为
/>
其中，k表示物质点x在其领域范围内物质点的总个数。
A9：判断伸长率是否大于拉伸断裂伸长率或小于压缩断裂伸长率，如果“是”，则物质点之间的键断裂，物质点之间的相互作用力消失；如果“否”，则按步骤A8计算物质点之间的本构力；
所述键的伸长率可通过下式计算：
所述拉伸断裂伸长率可通过取两点平均值计算：
同理，所述压缩断裂伸长率也可通过取平均值求得：
A10：通过物质点本构力计算出物质点的加速度，利用显式向前和向后差分公式进行时域积分，求出位移和速度；
所述物质点的加速度可表示为：
所述位移和速度可表示为：
A11：输出计算结果，绘制损伤云图。
实施例3
所述计算效果的具体实施步骤包括：
B1模型参数初始化，输入材料力学参数，确定模型尺寸。
在本实施例中，如图6所示，模型的几何尺寸为228mm×76mm，弹性模量E＝31.37GPa，密度ρ＝2400kg/m3，泊松比预制裂隙长度为L＝19.05mm。
在本实例中，将模型沿x轴方向划分229个物质点，沿y轴方向划分77个物质点，物质点间距Δx＝1mm，领域尺寸δ＝3.015mm。模型物质点总数为17633。拉伸弹性伸长率s1t＝0.000315；拉伸临界伸长率s2t＝0.00035；拉伸断裂伸长率st＝0.0035；压缩弹性伸长率s1c＝-0.0021；压缩临界伸长率s2c＝-0.0296；压缩断裂伸长率sc＝-0.021。
B2确定物质点间距，领域半径，离散求解域，生成物质点坐标。
本实施例中，按从左到右排序，每个物质点沿x方向的坐标可表示为：
x＝-0.114+(i-1) (44)
其中，i＝1，2，…，229，表示沿x方向排列晶格的个数。
按从下到上排序，每个物质点沿y方向的坐标可表示为；
y＝-0.038+(j-1) (45)
其中，j＝1，2，…，77，表示沿y方向排列晶格的个数。
B3搜索物质点领域范围内的点并编号，初始化所有物质点在其领域范围内的键；
搜索每个物质点在其领域范围内H＝{x′∈R,||x′-x||＜3.015mm}的其他物质点并编号。对领域内所有的键进行初始化，将标量函数μ(x,t,ξ)、分别设为1和0。
B4对每个物质点拉伸断裂伸长率和压缩断裂伸长率按威布尔分布赋值。
如图7所示，在本实施例中，拉伸断裂伸长率的平均值服从威布尔分布拉伸断裂伸长率的形状参数m为5。压缩断裂伸长率的平均值/>服从威布尔分布压缩断裂伸长率的形状参数n为5
B5预制裂隙。
把裂纹所在中心线的键断开，即裂纹中心线上的键满足u(x,t,ξ)＝0。裂纹所在中心线的表达式为：
B6施加边界条件。
在模型上端中点位置处，模型下端左右两侧据中点位置101.6mm处设置虚拟边界层。对模型上端中点虚拟层施加冲击速度，对模型下端左侧的虚拟边界层施加位移边界条件，固定其x方向、y方向的位移，对模型下端右侧的虚拟边界层位移条件，固定其y方向的位移。
所述冲击速度可表示为：
其中，V0＝0.065m/s，t1＝1.96×10-4s。
B7计算物质点的本构力，求解物质点的速度，位移。
计算每个物质点近场领域范围内所有键的本构力并求和，即为单个物质点的本构力。利用显式向前和向后差分公式进行时域积分，求解速度、位移，时间步长为2.0×10-7s。
B8计算物质点的损伤。
计算每个物质点领域内键伸长率，如果键的伸长率s大于0.0035或者s小于-0.021，则物质点间的键断裂，即u(x,t,ξ)＝0，记录断键的个数，通过断键个数与近场领域内键总个数的比值求得物质点损伤。
图8为应用本发明仿真冲击荷载作用下单边切口三点弯梁的计算效果图，其特点为基于本发明所提出的岩石材料率效应本构模型的实施方法，建立了冲击荷载作用下单边切口三点弯梁的数值模型。应用本公开构建的岩石材料率效应本构模型构建方法可以模拟单边切口三点弯梁在冲击荷载作用下裂纹起裂、扩展、贯穿的全过程。裂纹首先从切口处起裂，之后沿着加载方向扩展，最后贯穿整个梁。
Claims (5)
1.一种近场动力学岩石材料率效应本构模型的构建方法，其特征在于：步骤为：
S11引入能够反映物质点作用强度随物质点间距变化的核函数，对近场动力学本构力函数进行改进；
改进后本构力函数为：
其中核函数的表达形式为：
其中，g(ξ,δ)为反映物质点作用强度随物质点间距变化的核函数，δ为领域半径，ξ为物质点之间的初始距离；
S12求解改进后本构力函数的应变能密度WPD(X)；
所述应变能密度WPD(X)为：
对于三维问题，应变能密度表示为
其中，δ为领域半径，c(0,δ)为微模量常数，s表示键的伸长率，ξ为物质点之间的初始距离；
对于二维问题，应变能密度表示为
其中，h为厚度，δ为领域半径，c(0,δ)为微模量常数，s表示键的伸长率；对于一维问题，应变能密度表示为
其中，δ为领域半径，A为面积，c(0,δ)为微模量常数，s表示键的伸长率；S13令近场动力学中的应变能密度WPD(X)和传统连续介质力学的应变能密度相等，求得修正后的微模量常数c(0,δ)；
微模量常数c(0,δ)为：
式中，E表示弹性模量，δ为领域半径，h表示厚度，A表示面积；
S14引入键的损伤函数D(s)；
所述键的损伤函数D(s)可表示为：
其中，sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率，s表示键的伸长率；
S15引入非局部粘性阻尼，构造近场动力学率效应本构模型；
S16将步骤S13中求得的微模量常数c(0,δ)以及步骤S14中键的损伤函数D(s)代入步骤S15中的近场动力学率效应本构模型，完成岩石材料率效应本构模型的构建。
2.根据权利要求1所述的近场动力学岩石材料率效应本构模型的构建方法，其特征在于：步骤S15中所述的率效应本构模型可表示为
其中，υ表示粘性阻尼系数，表示键的变形率，c(0,δ)为微模量常数，s表示键的伸长率；
式中，表示当前构型键的速度，η+ξ表示当前构型的键，||η+ξ||表示当前构型键的长度，||ξ||表示初始构型键的长度。
3.根据权利要求1所述的近场动力学岩石材料率效应本构模型的构建方法，其特征在于：步骤S16中所述的岩石材料率效应本构模型可表示为：
对于三维问题，
其中，E表示弹性模量；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量；/>表示当前构型键方向；
对于平面应力问题，
其中，E表示弹性模量；h表示厚度；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量；/>表示当前构型键方向；
对于平面应变问题，
其中，E表示弹性模量；h表示厚度；表示键的变形率；υ表示粘性阻尼系数；δ表示领域半径；sc表示压缩断裂伸长率；s2c表示压缩临界伸长率；s1c表示压缩弹性伸长率；st表示拉伸断裂伸长率；s2t表示拉伸临界伸长率；s1t表示拉伸弹性伸长率；ξ表示物质点之间的初始距离；s表示键的伸长量；/>表示当前构型键方向。
4.一种如权利要求1所述的近场动力学岩石材料率效应本构模型的实施方法，其特征在于：步骤为：
A1：数值模型参数初始化，输入材料力学参数，确定数值模型尺寸；
A2：确定近场动力学中物质点间距，领域半径，离散求解域，生成物质点坐标；
A3：搜索物质点领域范围内的点并编号，初始化所有物质点在其领域范围内的键；
A4：对每个物质点拉伸断裂伸长率和压缩断裂伸长率按威布尔分布赋值；
A5：根据实际岩石试样，判断步骤A1中的数值模型是否需要预制裂纹；
A6：对数值模型施加初始边界条件；
A7：确定数值模型的时间步长；
A8：计算所有物质点的本构力；计算每个物质点领域范围内所有键的本构力并求和；
A9：判断伸长率是否大于拉伸断裂伸长率或小于压缩断裂伸长率，如果“是”，则物质点之间的键断裂，物质点之间的相互作用力消失；如果“否”，则按步骤A8计算物质点之间的本构力；
A10：通过物质点本构力计算出物质点的加速度，利用显式向前和向后差分公式进行时域积分，求出位移和速度；
A11：输出计算结果，绘制损伤云图。
5.根据权利要求4所述的近场动力学岩石材料率效应本构模型的实施方法，其特征在于：步骤A4中所述服从威布尔分布的拉伸断裂伸长率可表示为：
其中，st表示拉伸断裂伸长率，表示拉伸断裂伸长率的平均值，m为拉伸断裂伸长率服从威布尔分布的形状参数，r(x)表示(0,1]之间均匀分布的随机数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310324881.1A CN116364218B (zh) | 2023-03-30 | 2023-03-30 | 基于近场动力学岩石材料率效应本构模型构建与实施方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310324881.1A CN116364218B (zh) | 2023-03-30 | 2023-03-30 | 基于近场动力学岩石材料率效应本构模型构建与实施方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116364218A CN116364218A (zh) | 2023-06-30 |
| CN116364218B true CN116364218B (zh) | 2023-09-08 |
Family
ID=86907047
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310324881.1A Active CN116364218B (zh) | 2023-03-30 | 2023-03-30 | 基于近场动力学岩石材料率效应本构模型构建与实施方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116364218B (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110414167A (zh) * | 2019-08-02 | 2019-11-05 | 贵州大学 | 准脆性材料的近场动力学本构力函数建模方法 |
| CN112131709A (zh) * | 2020-08-25 | 2020-12-25 | 山东大学 | 基于近场动力学本构模型的节理岩体力学仿真方法及系统 |
| CN114186456A (zh) * | 2021-12-02 | 2022-03-15 | 大连理工大学 | 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法 |
| CN114818404A (zh) * | 2022-02-07 | 2022-07-29 | 上海大学 | 一种脆性材料裂纹扩展的预测方法 |
| CN114896757A (zh) * | 2022-04-01 | 2022-08-12 | 首都师范大学 | 一种基于近场动力学理论的地面沉降建模方法 |
| WO2023011029A1 (zh) * | 2021-08-03 | 2023-02-09 | 大连理工大学 | 结构损伤分析的近场有限元法及在商用软件中的实现方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6631647B2 (en) * | 2001-04-26 | 2003-10-14 | Joseph B. Seale | System and method for quantifying material properties |
-
2023
- 2023-03-30 CN CN202310324881.1A patent/CN116364218B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110414167A (zh) * | 2019-08-02 | 2019-11-05 | 贵州大学 | 准脆性材料的近场动力学本构力函数建模方法 |
| CN112131709A (zh) * | 2020-08-25 | 2020-12-25 | 山东大学 | 基于近场动力学本构模型的节理岩体力学仿真方法及系统 |
| WO2023011029A1 (zh) * | 2021-08-03 | 2023-02-09 | 大连理工大学 | 结构损伤分析的近场有限元法及在商用软件中的实现方法 |
| CN114186456A (zh) * | 2021-12-02 | 2022-03-15 | 大连理工大学 | 结构冲击弹塑性断裂分析的时间间断态基近场动力学方法 |
| CN114818404A (zh) * | 2022-02-07 | 2022-07-29 | 上海大学 | 一种脆性材料裂纹扩展的预测方法 |
| CN114896757A (zh) * | 2022-04-01 | 2022-08-12 | 首都师范大学 | 一种基于近场动力学理论的地面沉降建模方法 |
Non-Patent Citations (1)
| Title |
|---|
| 《基于改进型近场动力学方法的多裂纹扩展分析》;秦洪远;《工程力学》;第31-38页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116364218A (zh) | 2023-06-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Saharan et al. | Numerical procedure for dynamic simulation of discrete fractures due to blasting | |
| CN113569442B (zh) | 一种基于rkpm-pd耦合算法的岩石裂纹扩展预测方法 | |
| Menouillard et al. | Smoothed nodal forces for improved dynamic crack propagation modeling in XFEM | |
| Sun | Damage zone prediction for rock blasting | |
| Liu et al. | Peridynamic modelling of impact damage in three-point bending beam with offset notch | |
| Das et al. | A mesh-free approach for fracture modelling of gravity dams under earthquake | |
| Han et al. | Lattice discrete particle modeling of size effect in slab scratch tests | |
| Hu et al. | Numerical study on crack propagation by using softening model under blasting | |
| Yu et al. | Explicit finite element modeling of static crack propagation in reinforced concrete | |
| Jiang et al. | Dynamic responses and damage mechanism of rock with discontinuity subjected to confining stresses and blasting loads | |
| Chen et al. | A multiscale method coupling peridynamic and boundary element models for dynamic problems | |
| Wang et al. | Influence of material heterogeneity on the blast-induced crack initiation and propagation in brittle rock | |
| CN116364218B (zh) | 基于近场动力学岩石材料率效应本构模型构建与实施方法 | |
| Wang et al. | Dynamic strain field for granite specimen under SHPB impact tests based on stress wave propagation | |
| CN116110520B (zh) | 基于近场动力学岩石类材料本构模型的建立与仿真方法 | |
| Yahaghi et al. | Development of a three-dimensional grain-based combined finite-discrete element method to model the failure process of fine-grained sandstones | |
| Liong et al. | Finite element multibody simulation of a breathing crack in a rotor with a cohesive zone model | |
| Ye et al. | Attenuation characteristics of shock waves in drilling and blasting based on viscoelastic wave theory | |
| Hu et al. | Numerical simulations of arbitrary evolving cracks in geotechnical structures using the nonlinear augmented finite element method (N-AFEM) | |
| Yang et al. | Behavior and law of crack propagation in the dynamic-static superimposed stress field | |
| Uenishi | Rupture, waves and earthquakes | |
| Fan et al. | Peridynamic modeling of stochastic fractures in bolted glass plates | |
| Ma et al. | Particle Flow Code Simulation of the Characteristics of Crack Evolution in Rock‐Like Materials with Bent Cracks | |
| Cao et al. | Study on the guiding mechanism of multiple and empty holes under explosion load | |
| CN117854638A (zh) | 适用于准脆性材料力学行为的键基近场动力学本构方法 |
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 |