一种基于离散元的混凝土结构破坏过程数值模拟方法

1.本技术涉及土木工程仿真领域，具体而言，涉及一种基于离散元的混凝土结构破坏过程数值模拟方法。


背景技术：

2.随着现代化建设脚步的加快，混凝土作为主导材建筑料凭借取材便捷、造价低廉、抗压强度高等优势被广泛应用到各大土木工程中。研究混凝土结构力学性质及破坏过程，对工程设计施工、安全评估、质量保障等方面具有重要意义。混凝土是水泥、细骨料和粗骨料胶结而成的多相复合材料，具有明显的非均质属性。而混凝土结构的损伤往往萌生于材料局部截面，后逐渐扩展至整体。因此，有必要从细观角度出发，对混凝土结构破坏过程中的损伤机理及受力特性进行深入研究。
3.数值模拟研究材料的行为来减少试验开支，加速试验试错过程，是节约研究成本的有效方法。对于混凝土结构的数值模拟大多采用均质方法，以有限元法见常。但受到有限元方法本身限制，如相邻单元边界需要满足变形协调条件、裂缝出现时计算复杂(不连续)、增量法分析非线性问题时解的稳定差。同时，材料内部相互作用总会有新的受力表面产生，导致初始裂缝分布具有一定随机性。因此不宜采用均质方法对混凝土结构断裂行为进行系统分析.因此，急需一种基于离散元的混凝土结构损伤数值模拟方法,来解决现有技术中存在的技术问题。


技术实现要素：

4.针对以上问题，本发明的目的是提出一种基于离散元方法的混凝土破坏过程数值模拟方法，以克服现有数值模拟方法局限性问题，且该方法可准确、快速预测混凝土从细观到宏观结构构件或结构的破坏过程。
5.为了实现上述目的，本发明提供了一种基于离散元的混凝土结构破坏过程数值模拟方法，包括以下步骤：
6.基于目标混凝土结构件，采集目标混凝土结构件的边界墙尺寸，构建三维刚性边界墙，其中，三维刚性边界墙的墙体单元采用平面单元；
7.在三维刚性边界墙的墙体内部空间生成随机颗粒，压缩颗粒至高紧凑状态，根据颗粒中心间距与两颗粒的半径之和，获取处于接触状态的第二颗粒，并在第二颗粒的内部生成接触键；
8.基于接触键的法向单位向量、切向单位向量，获取第二颗粒的合力向量及弯矩，以及第二颗粒在第一单位时间的加速度和角加速度，根据显示中心差分法，迭代更新第二颗粒在第二单位时间的位置、速度和角速度；
9.根据迭代更新的结果，获取第二颗粒的接触作用力，其中，通过将接触作用力的法向力和剪切力，与目标混凝土结构件的抗拉强度、抗压强度、剪切强度进行对比，获取目标混凝土结构件的损伤情况。
10.优选地，在压缩颗粒至高紧凑状态的过程前，采用控制孔隙度方法生成颗粒的颗粒数量，包括：
[0011][0012]
其中，n表示生成目标粒子数，v表示给定空间总体积，c表示孔隙度，r表示平均半径；体积分数在0.615-0.735之间。
[0013]
优选地，在压缩颗粒至高紧凑状态的过程中，通过弹性半无线空间压板压缩颗粒，压缩速度控制在0.0001m/s-0.0012m/s之间；
[0014]
压缩结束后，颗粒与压缩板之间的接触力小于0.1n；
[0015]
高紧凑状态的判断依据为：检测颗粒间的配位数是否在4.01-8.15之间，判断颗粒是否处于高紧凑状态，以及检查颗粒的重叠性。
[0016]
优选地，在第二颗粒的内部生成接触键的过程中，令颗粒间的内聚力处于零应力状态，两个颗粒之间存在接触键需满足以下条件：
[0017][0018]
其中，d
ij
表示两颗粒中心间距，ri表示第i个颗粒的半径，rj表示第j个颗粒的半径，α表示允许不完全接触的颗粒存在接触键。
[0019]
优选地，在获取第二颗粒的合力向量及弯矩的过程中，法向单位向量用表示，
[0020]
切向单位向量用表示：
[0021][0022]
其中，δp
c，rel
＝δp
cj-δp
ci
，δp
cj
表示第i个颗粒的接触点的位移增量，δp
ci
表示第j个颗粒的接触点的位移增量，位移增量等于颗粒的位移增量与角度增量之和。
[0023]
优选地，在获取第二颗粒的合力向量及弯矩的过程中，利用法向力和切向力的单位向量分别点乘相应力的标量，然后相加得到合力向量；
[0024]
弯矩的计算公式为其中kr表示旋转刚度矩阵，表示旋转角度增量，ri表示第i个颗粒的半径，表示第i个颗粒的和第j个颗粒之间的作用力。
[0025]
优选地，法向力为：
[0026]
切向力为：其中，其中当颗粒间相对位移为正时，βn表示拉伸折减系数；当颗粒间相对位移为负时，βn表示压缩折减系数；βx表示剪切强度折减系数；
[0027]
颗粒间的法向刚度为：其中，ki表示第i个颗粒的材料刚度；kj表示第j个颗粒的材料刚度；
[0028]
颗粒间的切向刚度为ks＝μkn，其中，μ表示摩擦系数。
[0029]
优选地，在获取加速度和角加速度的过程中，根据合力向量及弯矩，通过牛顿第二定律或达朗贝尔原理，获取第二颗粒在第一单位时间的加速度和角加速度。
[0030]
优选地，在根据显示中心差分法进行迭代更新的过程中，包括以下步骤：
[0031]
设置half_t＝dt*(0:n) dt/2，t＝dt*(0:n)，其中，dt为时间增量，n为大于0的正整数；
[0032]
第二颗粒在第一单位时间t的加速度为：
[0033][0034]
第二颗粒在第一单位时间t的角加速度为：
[0035][0036]
第二颗粒在第二单位时间t 1的加速度：
[0037][0038]
第二颗粒在第二单位时间t 1的角速度：
[0039][0040]
优选地，目标混凝土结构件的损伤情况包括，法向力判别公式和切向力判别公式；
[0041]
法向力判别公式为：
[0042]fn
≤αffr[0043]-fn≤αcfc[0044]
切向力判别公式为：
[0045]fs
≤α
sfs-μfn[0046]
其中，α
t
、αc、αs分别为拉伸、压缩和剪切强度作用参数，fn表示颗粒间的法向力，fs表示颗粒间的切向力，f
t
表示颗粒材料的抗拉强度，fc表示材料的抗压强度。
[0047]
与现有技术相比，本发明具有以下有益效果：
[0048]
本发明提供的方法可准确、快速预测混凝土从细观到宏观结构构件或结构的破坏过程。
附图说明
[0049]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0050]
图1为本发明所述的模型截面示意图；
[0051]
图2为本发明所述的模型的应力应变曲线；
[0052]
图3为本发明所述的模拟方法流程图。
具体实施方式
[0053]
下为使本技术实施例的目的、技术方案和优点更加清楚，下面将结合本技术实施
例中附图，对本技术实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本技术一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本技术实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本技术的实施例的详细描述并非旨在限制要求保护的本技术的范围，而是仅仅表示本技术的选定实施例。基于本技术的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本技术保护的范围。
[0054]
如图1-3所示，本发明公开了一种基于离散元的混凝土结构破坏过程数值模拟方法，包括以下步骤：
[0055]
基于目标混凝土结构件，采集目标混凝土结构件的边界墙尺寸，构建三维刚性边界墙，其中，三维刚性边界墙的墙体单元采用平面单元；
[0056]
在三维刚性边界墙的墙体内部空间生成随机颗粒，压缩颗粒至高紧凑状态，根据颗粒中心间距与两颗粒的半径之和，获取处于接触状态的第二颗粒，并在第二颗粒的内部生成接触键；
[0057]
基于接触键的法向单位向量、切向单位向量，获取第二颗粒的合力向量及弯矩，以及第二颗粒在第一单位时间的加速度和角加速度，根据显示中心差分法，迭代更新第二颗粒在第二单位时间的位置、速度和角速度；
[0058]
根据迭代更新的结果，获取第二颗粒的接触作用力，其中，通过将接触作用力的法向力和剪切力，与目标混凝土结构件的抗拉强度、抗压强度、剪切强度进行对比，获取目标混凝土结构件的损伤情况。
[0059]
进一步优选地，在压缩颗粒至高紧凑状态的过程前，采用控制孔隙度方法生成颗粒的颗粒数量，包括：
[0060][0061]
其中，n表示生成目标粒子数，v表示给定空间总体积，c表示孔隙度，r表示平均半径；体积分数在0.615-0.735之间。
[0062]
进一步优选地，在压缩颗粒至高紧凑状态的过程中，通过弹性半无线空间压板压缩颗粒，压缩速度控制在0.0001m/s-0.0012m/s之间；
[0063]
压缩结束后，颗粒与压缩板之间的接触力小于0.1n；
[0064]
高紧凑状态的判断依据为：检测颗粒间的配位数是否在4.01-8.15之间，判断颗粒是否处于高紧凑状态，以及检查颗粒的重叠性。
[0065]
进一步优选地，在第二颗粒的内部生成接触键的过程中，令颗粒间的内聚力处于零应力状态，两个颗粒之间存在接触键需满足以下条件：
[0066][0067]
其中，d
ij
表示两颗粒中心间距，ri表示第i个颗粒的半径，rj表示第j个颗粒的半径，α表示允许不完全接触的颗粒存在接触键。
[0068]
进一步优选地，在获取第二颗粒的合力向量及弯矩的过程中，法向单位向量用表示，
[0069]
切向单位向量用表示：
[0070][0071]
其中，δp
c，rel
＝δp
cj-δp
ci
，δp
cj
表示第i个颗粒的接触点的位移增量，δp
ci
表示第j个颗粒的接触点的位移增量，位移增量等于颗粒的位移增量与角度增量之和。
[0072]
进一步优选地，在获取第二颗粒的合力向量及弯矩的过程中，利用法向力和切向力的单位向量分别点乘相应力的标量，然后相加得到合力向量；
[0073]
弯矩的计算公式为其中kr表示旋转刚度矩阵，表示旋转角度增量，ri表示第i个颗粒的半径，表示第i个颗粒的和第j个颗粒之间的作用力。
[0074]
进一步优选地，法向力为：
[0075]
切向力为：其中，其中当颗粒间相对位移为正时，βn表示拉伸折减系数；当颗粒间相对位移为负时，βn表示压缩折减系数；βs表示剪切强度折减系数；
[0076]
颗粒间的法向刚度为：其中，ki表示第i个颗粒的材料刚度；kj表示第j个颗粒的材料刚度；
[0077]
颗粒间的切向刚度为ks＝μkn，其中，μ表示摩擦系数。
[0078]
进一步优选地，在获取加速度和角加速度的过程中，根据合力向量及弯矩，通过牛顿第二定律或达朗贝尔原理，获取第二颗粒在第一单位时间的加速度和角加速度。
[0079]
进一步优选地，在根据显示中心差分法进行迭代更新的过程中，包括以下步骤：
[0080]
设置half_t＝dt*(0：n) dt/2，t＝dt*(0：n)，其中，dt为时间增量，n为大于0的正整数；
[0081]
第二颗粒在第一单位时间t的加速度为：
[0082][0083]
第二颗粒在第一单位时间t的角加速度为：
[0084][0085]
第二颗粒在第二单位时间t 1的加速度：
[0086][0087]
第二颗粒在第二单位时间t 1的角速度：
[0088][0089]
进一步优选地，目标混凝土结构件的损伤情况包括，法向力判别公式和切向力判别公式；
[0090]
法向力判别公式为：
[0091]fn
≤α
tft
[0092]-fn≤αcfc[0093]
切向力判别公式为：
[0094]fs
≤α
sfs-μfn[0095]
其中，α
t
、αc、αs分别为拉伸、压缩和剪切强度作用参数，fn表示颗粒间的法向力，fs表示颗粒间的切向力，f
t
表示颗粒材料的抗拉强度，fc表示材料的抗压强度。
[0096]
本发明还包括一种基于离散元的混凝土结构破坏过程数值模拟系统，用于实现上述模拟方法，包括：
[0097]
数据采集与处理模块，用于基于目标混凝土结构件，采集目标混凝土结构件的边界墙尺寸，构建三维刚性边界墙，其中，三维刚性边界墙的墙体单元采用平面单元；
[0098]
接触键生成模块，用于在三维刚性边界墙的墙体内部空间生成随机颗粒，压缩颗粒至高紧凑状态，根据颗粒中心间距与两颗粒的半径之和，获取处于接触状态的第二颗粒，并在第二颗粒的内部生成接触键；
[0099]
数据运算模块，用于基于接触键的法向单位向量、切向单位向量，获取第二颗粒的合力向量及弯矩，以及第二颗粒在第一单位时间的加速度和角加速度，根据显示中心差分法，迭代更新第二颗粒在第二单位时间的位置、速度和角速度；
[0100]
损伤评估模块，用于根据迭代更新的结果，获取第二颗粒的接触作用力，其中，通过将接触作用力的法向力和剪切力，与目标混凝土结构件的抗拉强度、抗压强度、剪切强度进行对比，获取目标混凝土结构件的损伤情况。
[0101]
实施例1：本发明的本发明的具体实现方法，包括以下步骤：
[0102]
设置几何和物理两组参数。几何参数包括迭代计算次数n、时间增量dt、颗粒质量、密度、半径、转动惯量等。物理参数包括弹性模量、泊松比、摩擦角、粘聚力、局部(全局)阻尼系数等。其中，微观参数法相刚度kn、切向刚度ks，需要依据颗粒的宏观弹性模量和泊松比确定。同时，利用材料局部塑性准则和断裂准则模拟混凝土的非线性行为，进一步标定材料的微观影响系数。
[0103]
确定混凝土结构构件整体尺寸，确定边界墙尺寸，建立三维刚性边界墙，墙体单元采用平面单元。
[0104]
在墙体内部空间生成随机颗粒，采用控制孔隙度方法生成颗粒数量，计算公式为其中n表示生成目标粒子数；v表示给定空间总体积；c表示孔隙度；r表示平均半径。为避免细化程度过高影响计算效率，控制体积分数在0.615-0.735之间。
[0105]
内部生成颗粒后，通过弹性半无线空间压板压缩颗粒，压缩速度控制在0.0001m/s-0.0012m/s之间。压缩结束时颗粒与压缩板之间的接触力控制在0.1n以下。待模型达到平衡状态，检查颗粒的两个指标。第一，检测颗粒间的配位数是否在4.01-8.15之间，判断颗粒是否处于高紧凑状态；第二，检查颗粒的重叠性。若以上满足条件，记录所有颗粒的位置坐标及初始速度和角速度。若不满足条件，返回上一步调整。
[0106]
确定接触检测方法，根据颗粒中心间距与两颗粒的半径之和做比较，若小于半径之和，则认定为两个颗粒接触；反之不接触。利用布尔矩阵记录判断结果，被判定为接触的颗粒，可以进行下一步计算接触力，被判定为不接触的颗粒不需要计算接触力。
[0107]
确定接触模型。确定好颗粒的数量和排列后，根据以下方程生成颗粒内部的接触键。在创建初始时，令颗粒间的内聚力处于零应力状态。两个颗粒之间存在接触键需满足以
下条件：其中d
ij
表示两颗粒中心间距；ri表示第i个颗粒的半径；rj表示第j个颗粒的半径；表示允许不完全接触的颗粒存在接触键。第i个颗粒的质心与接触点c之间的距离用d
ci
表示，接触点的位移增量等于颗粒的位移增量与角度增量之和。两个接触点的相对位移增量为：δp
c，rel
＝δp
cj-δp
ci
，其中δp
ci
表示第i个颗粒的接触点的位移增量；δp
cj
表示第j个颗粒的接触点的位移增量。接触键的法向单位向量用表示，接触键的切向单位向量用表示根据上一步得到的接触检测判断结果，分别计算在第i个颗粒周围的颗粒对其产生的接触作用力向量，将所有接触力向量相加视为第i个颗粒的合力。
[0108]
法向力和切向力的具体求解方法。利用材料的拉伸、压缩和剪切刚度等相关损伤折减系数定义法向力fn、切向力及弯矩m
ij
。法向力及切向力其中当颗粒间相对位移为正时，βn表示拉伸折减系数；当颗粒间相对位移为负时，βn表示压缩折减系数；βs表示剪切强度折减系数。颗粒间的法向刚度计算公式为其中ki表示第i个颗粒的材料刚度；kj表示第j个颗粒的材料刚度。ki及kj可通过宏观物理参数计算获得。颗粒间的切向刚度计算公式为ks＝μkn，其中μ表示摩擦系数。
[0109]
颗粒利用法向和切的单位向量分别点乘相应力的标量，然后相加得到合力向量。弯矩计算公式为其中kr表示旋转刚度矩阵，表示旋转角度增量。
[0110]
获得力和弯矩后，利用牛顿第二定律或达朗贝尔原理，获得在经过时间dt后颗粒的加速度和角加速度
[0111]
利用显示中心差分法，通过以下公式迭代更新颗粒的位置、速度和角速度。设half_t＝dt*(0：n) dt/2，t＝dt*(0：n)其中dt为时间增量，n为大于0的正整数。则有，在t时刻颗粒的速度和角速度为：
[0112][0113][0114]
t 1时刻的速度和角速度为：
[0115][0116][0117]
对t 1时刻的速度和角速度进行积分，得到该时刻的颗粒位置坐标。此过程反复迭代n次。
[0118]
导入损伤模型。通过以上方程计算任意颗粒间的接触作用力，根据迭代计算的法
向力和剪切力与材料的抗拉强度、抗压强度、剪切强度对比，判断材料的损伤情况。法向力和切向力的判别公式如下：
[0119]fn
≤α
tft
(法向拉力作用)
[0120]-fn≤αcfc(法向压力作用)
[0121]fs
≤α
sf5-μfn(剪切力作用)
[0122]
其中，α
t
、αc、αs分别为拉伸、压缩和剪切强度作用参数，与材料的线性变形极值和变形极值有关。两个颗粒间作用力超出以上允许范围内时，则可以认定接触键断裂。
[0123]
记录模型的相关模拟数据结果，用于后续具体分析。
[0124]
实施例2：
[0125]
采用本发明方法对无筋混凝土柱压缩进行数值模拟。
[0126]
根据设计尺寸建立刚性单元边界墙。梁宽度统一设置为200mm，长高比例关系用a表示。利用孔隙度法在墙体内部生成平均半径约为4.5mm的颗粒堆积体模型。利用重力及弹性半无线压力板使颗粒堆积到一起，待模型平衡后，发生接触颗粒创建接触键。在梁中心位置施加线荷载。颗粒的物理参数如表1所示，模型颗粒数和接触键数如表2所示。
[0127]
表1
[0128]
名称骨料颗粒砂浆颗粒杨氏模量(gpa)5050泊松比0.20.2摩擦角(
°
)2020密度(kg/m3)24002100内聚强度(mpa)2020抗拉强度(mpa)5.257.2阻尼系数0.10.1
[0129]
表2
[0130][0131]
模型截面示意图如图1所示。模型的应力应变曲线如图2所示。
[0132]
本发明提出混凝土破坏过程数值模拟方法，克服了现有数值模拟方法局限性问题，且该方法可准确、快速预测混凝土从细观到宏观结构构件或结构的破坏过程。