CN118094872A - 一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法 - Google Patents
一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法 Download PDFInfo
- Publication number
- CN118094872A CN118094872A CN202410052913.1A CN202410052913A CN118094872A CN 118094872 A CN118094872 A CN 118094872A CN 202410052913 A CN202410052913 A CN 202410052913A CN 118094872 A CN118094872 A CN 118094872A
- Authority
- CN
- China
- Prior art keywords
- particle
- calculation
- particles
- stress
- scale
- 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.)
- Pending
Links
- 239000002245 particle Substances 0.000 title claims abstract description 153
- 238000004364 calculation method Methods 0.000 title claims abstract description 52
- 238000000034 method Methods 0.000 claims abstract description 56
- 238000005516 engineering process Methods 0.000 claims abstract description 6
- 239000007787 solid Substances 0.000 claims description 9
- 230000001133 acceleration Effects 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 3
- 238000013016 damping Methods 0.000 claims description 2
- 238000006073 displacement reaction Methods 0.000 claims description 2
- 238000011002 quantification Methods 0.000 claims description 2
- 238000010187 selection method Methods 0.000 claims description 2
- 239000002689 soil Substances 0.000 claims description 2
- 238000011160 research Methods 0.000 abstract description 6
- 230000008901 benefit Effects 0.000 abstract description 5
- 230000007547 defect Effects 0.000 abstract description 5
- 230000009286 beneficial effect Effects 0.000 abstract description 2
- 238000004088 simulation Methods 0.000 abstract description 2
- 238000012935 Averaging Methods 0.000 description 26
- 230000005484 gravity Effects 0.000 description 6
- 230000003993 interaction Effects 0.000 description 5
- 230000002265 prevention Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 230000000737 periodic effect Effects 0.000 description 3
- 230000009471 action Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000005206 flow analysis Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000012876 topography Methods 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000004141 dimensional analysis Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000011236 particulate material Substances 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000000087 stabilizing effect Effects 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种跨尺度应力计算方法，其通过融合泰森多边形计算和光滑粒子近似技术，实现了高效准确实现跨尺度应力计算。本计算通过引入泰森多边形计算，将颗粒间作用力直接转化为微观颗粒应力张量，进一步基于光滑粒子近似技术实现了宏观应力的准确计算，克服了传统分层平均法计算精度低的难题。本发明的有益效果是：首先，克服了分层平均法计算精度过度依赖代表性单元尺寸和时间平均的弊端；其次，本方法在保证计算效率的前提下，显著提高了计算精度，并具有易于编程实现的优势，极大方便了跨尺度计算。本发明的公开将极大推动工程领域的跨尺度计算与模拟研究，为提高我国高性能计算水平奠定基础。
Description
技术领域
本发明涉及一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法，属于岩土力学、地质灾害防治等工程技术领域。
背景技术
颗粒材料在自然界中广泛存在，全球频发的地质灾害(滑坡、崩塌流等)本质是颗粒流动行为。因此，准确描述颗粒的流变学行为，是提高地质灾害预测精确度及防灾减灾水平的关键。基于颗粒介质力学的离散单元法，从颗粒间微观相互作用力出发，实现了颗粒流动行为的准确描述。然而，如何高效准确的将微观颗粒特性转化为应力、应变等宏观力学特征是制约跨尺度计算的难题。
目前，最常用的跨尺度计算方法是分层(块)平均法。该方法在拟求解位置构建一定体积的代表性单元，通过将代表性单元的颗粒信息进行统计平均，获取宏观力学特征。然而该方法计算精度依赖于代表性单元尺寸，同时需进行长时间的时间平均处理，才能获取稳定准确的宏观变量。因此，分层(块)平均法存在计算准确性低、效率低等缺陷，同时其依赖时间平均的特性，导致其无法应用于非稳态颗粒流动分析中。
发明内容
本发明所描述的基于泰森多边形和光滑粒子近似技术的新型跨尺度应力计算方法，可以弥补目前跨尺度应力计算方法的不足，实现了准确、高效的宏观应力计算，为研究地质灾害致灾机理和防治措施奠定基础。
本发明的目的是为了弥补传统跨尺度应力计算中准确性低及计算效率低的缺陷，提供了一种基于泰森多边形和光滑粒子近似技术的跨尺度应力计算方法，可促进地质灾害的跨尺度研究实施。
为了达到上述目的，本发明采用如下技术方案：
一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法，其特征在于，具体步骤如下：
步骤1，颗粒体系中，根据颗粒空间位置构建泰森多边形，并计算每个颗粒的泰森多边形体积Vi。
进一步，所述泰森多边形构建方法：将颗粒与其临近颗粒球心相连，并基于球心连线绘制垂直平分面，绘制所有临近颗粒构成的垂直平分面，上述垂直平分面组成的封闭空间即为该颗粒的泰森多边形。
步骤2，基于离散单元法，当颗粒接触时，根据如下公式(1)和(2)计算颗粒微观接触力Fij：
其中，j代表所有与颗粒i接触的颗粒，和/>分别为微观法向和切向接触力，δ为重合量，d为颗粒直径；f(x)选择方法如下：当颗粒间重合量比较小时，选择f(x)＝1，为线弹性接触力模型，在地质灾害领域，岩土体颗粒相互作用时重合量很有限，故选用此函数；当颗粒塑性变形比较大时，选择/>为Hertzian型非线性接触力模型。kn和kt分别为法向和切向刚度，δij为两颗粒重叠量，nij为两颗粒中心点连线单位向量，/>为弹性剪切位移，γn和γt分别为法向和切向阻尼系数，/>和/>分别为法向和切向相对速度，meff为有效质量。
步骤3，根据步骤1中计算的泰森多边形体积Vi和步骤2计算的计算颗粒微观接触力Fij，采用如下公式计算微观颗粒应力张量σi：
步骤4，计算颗粒集合体宏观应力σ、速度u及固体体积分数具体步骤如下：
基于统计力学，颗粒体系在空间点坐标为r、时间点t时的固体体积分数可表示为：
进一步，根据运动学中的动量矩概念和统计力学理论，颗粒体系在空间点坐标为r、时间点t时的动量矩密度M(r，t)可表示为：
M(r，t)＝∑mkupkW(r-rk(t))
动量矩密度M(r，t)除以颗粒体系在空间点坐标为r、时间点t时的密度可以计算空间点坐标为r、时间点t时的速度u(r，t)：
对公式(5)速度u(r，t)对时间求导可以得到计算空间点坐标为r、时间点t时的加速度a(r，t)：
根据牛顿第二定量，在公式(6)两侧乘以密度ρp，即可以获取计算空间点坐标为r、时间点t时的宏观应力σ(r，t)，其中σk＝apkρp，由此可以获得宏观应力计算公式如下：
其中r为空间点坐标，t为时间点，rk(t)、Vpk、ρp、mk、upk、σk分别为粒子k在t时刻位置、体积、密度、质量、速度及接触力合力，W为权函数(可选用高斯函数、三次钟形函数等)。
步骤5，输出并绘制颗粒体系宏观应力、速度、固体体积分数等曲线。
与现有技术相比，本发明有益效果如下：
(a)克服了分层(块)平均法计算精度依赖代表性单元尺寸和时间平均的不足。本方法通过融合泰森多边形和和光滑粒子近似技术，仅通过空间平均，而无需时间平均即可获取准确结果；同时降低了结果对代表性单元尺寸的依赖性。
(b)本方法具有计算效率快及精度高的优点。与仅进行空间平均而不进行时间平均的方法相比，本方法在不降低计算效率的前提下，具有精度高的优点；虽然，长时间的分层平均法可以提高精度，但计算效率将显著低于本方法。
(c)本方法为地质灾害中的非稳态流动分析提供了跨尺度分析工具。由于受到地形地貌的影响，地质灾害通常属于非稳态流动，即在同一坐标点的物理量随时间不断变化，采用长时间的分层平均法将无法刻画这种随时间剧烈变化的非稳态流动。本方法无需时间平均，仅采用空间平均即可获取准确结果，从而提供了一种高效准确的跨尺度应力计算方法，对于推进地质灾害的多尺度研究具有重要作用，为准确模拟地质灾害致灾过程和提高防灾能力奠定基础。
附图说明
图1为本发明处理流程图
图2为本发明实施例斜面颗粒流计算模型
图3为本发明实施例泰森多边形构建及颗粒微观受力计算示意图
图4为本发明实施例跨尺度应力计算示意图
图5为本发明实施例固体体积分数沿模型高度分布图
图6为本发明实施例流动速度沿模型高度分布图
图7为本发明实施例应变率沿模型高度分布图
图8为本发明实施例正应力沿模型高度分布图
图9为本发明实施例剪应力沿模型高度分布图
具体实施方式
下面结合实施例对本发明方法进行说明，以更清楚地表达本发明的目的及内容。本实施例为本发明优选实施方式，但不作为本发明的限定，其他凡是其原理和基本结构或实现方法与本实施例相同或近似的，均在本发明保护范围之内。
实施例1：稳定斜面颗粒流
斜面颗粒流(如图2所示)是颗粒集合体在重力作用下沿底部粗糙斜面的颗粒流动模型。它通过在流动方向(x轴)和垂直流动方向(y轴)设置周期边界，实现了采用有限颗粒模拟半无限空间的颗粒流动，有效解决了传统计算模型存在的计算效率低及边界效应明显等问题，被广泛应用于研究斜坡地质灾害。当斜面颗粒流到达稳态时，颗粒系统处于受力平衡状态，可精确获得其流动速度、正应力及剪应力的理论解，十分适合检验本发明方法的准确性及有效性。因此，本项目采用稳定斜面颗粒流作为实施例。
1)采用离散单元法建立斜面颗粒流模型，具体包括颗粒生成、周期边界设置、底部粗糙边界形成及重力施加等。如图2所示，定义尺寸为40d×20d×80d的盒子，d为颗粒直径，在顶部区域生成粒子，粒子在重力作用下向下运动，到达盒子底部开始堆积，持续生成粒子，直至达到计划的粒子数目。在盒子四周设置周期边界，所谓周期边界是当粒子到达一侧边界时，自动变换到对面边界，从对边边界以相同的速度进入盒子。采用周期边界可以用有限的区域模拟无限长的区域。将底部一定厚度的颗粒(如图2黑色粒子，厚度为4d)设为固定边界，计算时该部分粒子保持不动，形成粗糙的底部边界。倾斜盒子24°，颗粒体系在重力作用下开始流动，保持计算至总动能不随时间变化，确认斜面颗粒流进入稳定流动状态；
2)根据颗粒位置，将颗粒与其临近颗粒球心相连，并基于球心连线绘制垂直平分面，绘制所有临近颗粒构成的垂直平分面，上述垂直平分面组成的封闭空间即为该颗粒的泰森多边形，计算该颗粒的泰森多边形体积。进一步，根据颗粒位置计算颗粒间重叠量，基于颗粒微观接触力(公式1-2)计算颗粒间相互作用力，将所有与该颗粒接触的粒子的颗粒剂相互作用力按照公式(3)计算该粒子的颗粒间接触力合力张量；进一步，在此基础上，用颗粒间接触力合力张量除以泰森多边形体积，获取并输出每个颗粒的微观颗粒应力张量。对每个颗粒重复上述步骤，获得所有颗粒的微观颗粒应力张量。
以二维圆盘颗粒为例说明上述过程，三维计算步骤类似。如图3，以计算颗粒0的微观受力和泰森多边形举例，颗粒1、颗粒2和颗粒3与颗粒0接触，产生颗粒间相互作用力，颗粒4与颗粒0无接触，无颗粒间相互作用力。
将临近颗粒与其临近颗粒球心相连(图3中虚线)，并基于球心连线绘制垂直平分面(图3中点划线)，绘制所有临近颗粒构成的垂直平分面，上述垂直平分面组成的封闭空间即为该颗粒的泰森多边形(点划线围成的多边形)，计算该颗粒的泰森多边形体积V0。
根据颗粒位置计算颗粒间重叠量，基于颗粒间接触模型(公式1-2)计算颗粒间相互作用力(即F01,F02和F03)。进一步，根据公式(3)计算颗粒0的颗粒间接触力合力张量从而获取颗粒的微观颗粒应力张量。对每个颗粒重复上述步骤，获得所有颗粒的微观颗粒应力张量。
3)根据颗粒位置、速度及2)中计算的微观颗粒应力张量确定宏观固体体积分数、宏观速度及宏观应力沿高度的分布曲线。以图4虚线位置处的宏观应力计算为例说明光滑粒子近似计算过程。权函数在本实施例中选用高斯函数为例，截断距离为10d。即与计算位置的距离在截断距离之内的粒子，会对计算位置的宏观变量计算有贡献。根据粒子空间位置确定计算位置(虚线位置)的权重，后根据公式(7)考虑截断距离内的所有粒子的贡献，计算出计算位置(虚线位置)的宏观应力。对不同坐标重复上述过程，即可获得宏观应力分布曲线。
对于体积分数和速度分布曲线，同样按照上述流程，但须注意需采用公式(4)方可确定体积分数分布曲线，采用公式(5)方可确认速度分布曲线。
本方法仅通过空间平均，无需进行时间平均，获取斜面颗粒流的宏观运动特征及应力分布曲线，如图5-9中的红色实线。
稳定斜面流速度vx(z)、应变率正应力σzz(z)及剪应力σxz(z)的理论解(绿色虚线)为：
其中ABag为Bagnold系数、h为稳态时斜面颗粒流高度、g为重力加速度、ρ为颗粒密度、固体体积分数，θ为斜面倾角。图中结果均进行无量纲化处理，即将颗粒直径d，颗粒密度ρ和重力加速度g选为控制变量，通过量纲分析确定不同物理物理量的单位，如高度单位为颗粒直径d，速度单位为/>应力单位为gρd。无量纲处理的优势在于易于发现物理规律，并建立相关数学关系，是科学研究中的常用手段。理论解与本发明计算结果吻合较好，仅在靠近边界部分出现误差，证明了本发明的准确性及有效性。
为进一步说明本发明的性能，将基于本方法(红色实线)与仅进行空间平均的分层平均法(蓝色圆圈)进行了对比研究。结果如图5-9所示，其中蓝色圆圈为采用仅进行空间平均的分层平均法的计算结果，红色实线为采用本发明的计算结果，绿色虚线为理论解。相比于本发明，基于仅进行空间平均的分层平均法计算的速度与本发明方法均可与理论解较好吻合，但是分层平均法计算的应变率、正应力、剪应力与理论解差别较大，精确度较低。
同时，本案例数据处理在同一仿真环境(硬件为主频2.7GHz i9-12900CPU,4GBRAM)执行，处理时间对比：仅进行空间平均的分层平均法耗时29s，本发明提出的方法耗时30s，两种计算效率相近，但是本方法的计算精度较高；同时，鉴于该实施例为稳态流动，可进行时间平均，为达到与本方法计算结果类似的精度，需对分层平均法进行200时间步的时间平均，需耗时3200s，分层平均法计算效率显著降低。此外，对于地质灾害常见的非稳态流动问题，由于流动行为随时间急剧变化，无法进行时间平均，因此分层平均法将无法获取准确结果；但由于本方法无需采用时间平均，因此适用于此类非稳态流动问题的跨尺度分析。通过上述对比研究可知，本发明无需进行时间平均即能将微观接触力转化为宏观应力，证明了本发明的高效性、准确性及对于处理非稳态流动的适用性。
上述描述仅是对本申请较佳实施例的描述，并非是对本申请范围的任何限定。任何熟悉该领域的普通技术人员根据上述揭示的技术内容做出的任何变更或修饰均应当视为等同的有效实施例，均属于本申请技术方案保护的范围。
Claims (2)
1.一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法，其特征在于，具体步骤如下：
步骤1，颗粒体系中，根据颗粒空间位置构建泰森多边形，并计算每个颗粒的泰森多边形体积Vi；
步骤2，基于离散单元法，当颗粒接触时，根据公式(1)和(2)计算颗粒微观接触力Fij：
其中，j代表所有与颗粒i接触的颗粒，和/>分别为微观法向和切向接触力，δ为重合量，d为颗粒直径；f(x)选择方法如下：当颗粒间重合量比较小时，选择f(x)＝1，为线弹性接触力模型，在地质灾害领域，岩土体颗粒相互作用时重合量很有限，故选用此函数；当颗粒塑性变形比较大时，选择/>为Hertzian型非线性接触力模型；kn和kt分别为法向和切向刚度，δij为两颗粒重叠量，nij为两颗粒中心点连线单位向量，/>为弹性剪切位移，γn和γt分别为法向和切向阻尼系数，/>和/>分别为法向和切向相对速度，meff为有效质量；
步骤3，根据步骤1中计算的泰森多边形体积Vi和步骤2计算的计算颗粒微观接触力Fij，采用公式(3)计算微观颗粒应力张量σi：
步骤4，基于光滑粒子近似技术计算颗粒集合体宏观应力σ、速度u及固体体积分数具体步骤如下：
基于统计力学，颗粒体系在空间点坐标为r、时间点t时的固体体积分数表示为：
进一步，根据运动学中的动量矩概念和统计力学理论，颗粒体系在空间点坐标为r、时间点t时的动量矩密度M(r,t)表示为：
M(r,t)＝∑mkupkW(r-rk(t))
动量矩密度M(r,t)除以颗粒体系在空间点坐标为r、时间点t时的密度计算空间点坐标为r、时间点t时的速度u(r,t)：
对公式(5)速度u(r,t)对时间求导得到计算空间点坐标为r、时间点t时的加速度a(r,t)：
根据牛顿第二定量，在公式(6)两侧乘以密度ρp，即可以获取计算空间点坐标为r、时间点t时的宏观应力σ(r,t)，其中σk＝apkρp，由此可以获得宏观应力计算公式如下：
其中r为空间点坐标，t为时间点，rk(t)、Vpk、ρp、mk、upk、σk分别为粒子k在t时刻位置、体积、密度、质量、速度及接触力合力，W为权函数；
步骤5，输出并绘制颗粒体系宏观应力、速度、固体体积分数等曲线。
2.如权利要求1所述一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法，其特征在于，步骤1所述泰森多边形构建方法：将颗粒与其临近颗粒球心相连，并基于球心连线绘制垂直平分面，绘制所有临近颗粒构成的垂直平分面，上述垂直平分面组成的封闭空间即为该颗粒的泰森多边形。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410052913.1A CN118094872A (zh) | 2024-01-14 | 2024-01-14 | 一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410052913.1A CN118094872A (zh) | 2024-01-14 | 2024-01-14 | 一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN118094872A true CN118094872A (zh) | 2024-05-28 |
Family
ID=91141125
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410052913.1A Pending CN118094872A (zh) | 2024-01-14 | 2024-01-14 | 一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN118094872A (zh) |
-
2024
- 2024-01-14 CN CN202410052913.1A patent/CN118094872A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Bian et al. | Micromechanical particle interactions in railway ballast through DEM simulations of direct shear tests | |
| Zhao et al. | A composite particle model for non-spherical particles in DEM simulations | |
| Izadi et al. | Simulating direct shear tests with the Bullet physics library: A validation study | |
| Zhu et al. | Animating Sand as a Surface Flow. | |
| CN113221200B (zh) | 一种适用于堆芯颗粒分布不确定性分析的三维高效随机排布方法 | |
| CN112258643A (zh) | 岩土边坡任意形状落石运动轨迹三维分析方法 | |
| CN109992830B (zh) | 基于物质点方法的山体滑坡灾害场景模拟方法 | |
| Li et al. | Fluid-solid interaction simulation for particles and walls of arbitrary polygonal shapes with a coupled LBM-IMB-DEM method | |
| Benabbou et al. | Geometrical modeling of granular structures in two and three dimensions. application to nanostructures | |
| CN112241603B (zh) | 高位滑坡冲击铲刮及下垫层汇入过程的数值仿真方法 | |
| CN118094872A (zh) | 一种基于泰森多边形和光滑粒子近似的跨尺度应力计算法 | |
| CN113033060A (zh) | 一种用于预测复杂煤层开采结构的优化方法 | |
| CN115758766A (zh) | 连续-非连续数值模型连接方法 | |
| CN114154384A (zh) | 一种三维立方体空间的球形颗粒随机填充算法 | |
| Frery et al. | Stochastic particle packing with specified granulometry and porosity | |
| Alonso-Marroquín et al. | The incremental response of soils. An investigation using a discrete-element model | |
| CN114357717A (zh) | 一种应用于双色币压印成形仿真中的基于改进接触算法的物质点法 | |
| He et al. | Comparing realistic particle simulation using discrete element method and physics engine | |
| Sun et al. | A resolved SPH-DEM coupling method for analysing the interaction of polyhedral granular materials with fluid | |
| CN109492285B (zh) | 一种可变形三维任意圆化凸多面体块体离散单元法 | |
| Alwan et al. | Investigation of Wave Forces on Fixed Monopile Foundation of Offshore Wind Turbine | |
| Unger | Characterization of static and dynamic structures in granular materials | |
| Ahmadian et al. | Simulating the fluid–solid interaction of irregularly shaped particles using the LBM-DEM coupling method | |
| CN109408973B (zh) | 一种基于距离势函数二维可变形凸多边形块体离散单元法 | |
| Ahmadian | Simulation of irregularly shaped particles using coupling method of lattice Boltzmann and discrete element modelling |
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 |