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

摘要

本发明公开了一种基于离散元的混凝土结构破坏过程数值模拟方法，包括以下步骤：基于目标混凝土结构件，采集目标混凝土结构件的边界墙尺寸，构建三维刚性边界墙，生成随机颗粒，压缩颗粒至高紧凑状态，根据颗粒中心间距与两颗粒的半径之和，获取处于接触状态的第二颗粒，并在第二颗粒的内部生成接触键；基于接触键的法向单位向量、切向单位向量，获取第二颗粒的合力向量及弯矩，以及第二颗粒在第一单位时间的加速度和角加速度，根据显示中心差分法，迭代更新第二颗粒在第二单位时间的位置、速度和角速度；根据迭代更新的结果，获取第二颗粒的接触作用力；本发明提供的方法可准确、快速预测混凝土从细观到宏观结构构件或结构的破坏过程。

说明书

技术领域

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

背景技术

随着现代化建设脚步的加快，混凝土作为主导材建筑料凭借取材便捷、造价低廉、抗压强度高等优势被广泛应用到各大土木工程中。研究混凝土结构力学性质及破坏过程，对工程设计施工、安全评估、质量保障等方面具有重要意义。混凝土是水泥、细骨料和粗骨料胶结而成的多相复合材料，具有明显的非均质属性。而混凝土结构的损伤往往萌生于材料局部截面，后逐渐扩展至整体。因此，有必要从细观角度出发，对混凝土结构破坏过程中的损伤机理及受力特性进行深入研究。

数值模拟研究材料的行为来减少试验开支，加速试验试错过程，是节约研究成本的有效方法。对于混凝土结构的数值模拟大多采用均质方法，以有限元法见常。但受到有限元方法本身限制，如相邻单元边界需要满足变形协调条件、裂缝出现时计算复杂(不连续)、增量法分析非线性问题时解的稳定差。同时，材料内部相互作用总会有新的受力表面产生，导致初始裂缝分布具有一定随机性。因此不宜采用均质方法对混凝土结构断裂行为进行系统分析.因此，急需一种基于离散元的混凝土结构损伤数值模拟方法,来解决现有技术中存在的技术问题。

发明内容

针对以上问题，本发明的目的是提出一种基于离散元方法的混凝土破坏过程数值模拟方法，以克服现有数值模拟方法局限性问题，且该方法可准确、快速预测混凝土从细观到宏观结构构件或结构的破坏过程。

为了实现上述目的，本发明提供了一种基于离散元的混凝土结构破坏过程数值模拟方法，包括以下步骤：

基于目标混凝土结构件，采集目标混凝土结构件的边界墙尺寸，构建三维刚性边界墙，其中，三维刚性边界墙的墙体单元采用平面单元；

在三维刚性边界墙的墙体内部空间生成随机颗粒，压缩颗粒至高紧凑状态，根据颗粒中心间距与两颗粒的半径之和，获取处于接触状态的第二颗粒，并在第二颗粒的内部生成接触键；

基于接触键的法向单位向量、切向单位向量，获取第二颗粒的合力向量及弯矩，以及第二颗粒在第一单位时间的加速度和角加速度，根据显示中心差分法，迭代更新第二颗粒在第二单位时间的位置、速度和角速度；

根据迭代更新的结果，获取第二颗粒的接触作用力，其中，通过将接触作用力的法向力和剪切力，与目标混凝土结构件的抗拉强度、抗压强度、剪切强度进行对比，获取目标混凝土结构件的损伤情况。

优选地，在压缩颗粒至高紧凑状态的过程前，采用控制孔隙度方法生成颗粒的颗粒数量，包括：

其中，n表示生成目标粒子数，V表示给定空间总体积，C表示孔隙度，R表示平均半径；体积分数在0.615-0.735之间。

优选地，在压缩颗粒至高紧凑状态的过程中，通过弹性半无线空间压板压缩颗粒，压缩速度控制在0.0001m/s-0.0012m/s之间；

压缩结束后，颗粒与压缩板之间的接触力小于0.1N；

高紧凑状态的判断依据为：检测颗粒间的配位数是否在4.01-8.15之间，判断颗粒是否处于高紧凑状态，以及检查颗粒的重叠性。

优选地，在第二颗粒的内部生成接触键的过程中，令颗粒间的内聚力处于零应力状态，两个颗粒之间存在接触键需满足以下条件：

其中，d

优选地，在获取第二颗粒的合力向量及弯矩的过程中，法向单位向量用

切向单位向量用

其中，Δp

优选地，在获取第二颗粒的合力向量及弯矩的过程中，利用法向力和切向力的单位向量分别点乘相应力的标量，然后相加得到合力向量；

弯矩的计算公式为

优选地，法向力为：

切向力为：

颗粒间的法向刚度为：

颗粒间的切向刚度为K

优选地，在获取加速度和角加速度的过程中，根据合力向量及弯矩，通过牛顿第二定律或达朗贝尔原理，获取第二颗粒在第一单位时间的加速度和角加速度。

优选地，在根据显示中心差分法进行迭代更新的过程中，包括以下步骤：

设置half_t＝dt*(0:N)+dt/2，t＝dt*(0:N)，其中，dt为时间增量，N为大于0的正整数；

第二颗粒在第一单位时间t的加速度为：

第二颗粒在第一单位时间t的角加速度为：

第二颗粒在第二单位时间t+1的加速度：

第二颗粒在第二单位时间t+1的角速度：

优选地，目标混凝土结构件的损伤情况包括，法向力判别公式和切向力判别公式；

法向力判别公式为：

F

-F

切向力判别公式为：

F

其中，α

与现有技术相比，本发明具有以下有益效果：

本发明提供的方法可准确、快速预测混凝土从细观到宏观结构构件或结构的破坏过程。

附图说明

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

图1为本发明所述的模型截面示意图；

图2为本发明所述的模型的应力应变曲线；

图3为本发明所述的模拟方法流程图。

具体实施方式

下为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本申请的实施例的详细描述并非旨在限制要求保护的本申请的范围，而是仅仅表示本申请的选定实施例。基于本申请的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本申请保护的范围。

如图1-3所示，本发明公开了一种基于离散元的混凝土结构破坏过程数值模拟方法，包括以下步骤：

基于目标混凝土结构件，采集目标混凝土结构件的边界墙尺寸，构建三维刚性边界墙，其中，三维刚性边界墙的墙体单元采用平面单元；

在三维刚性边界墙的墙体内部空间生成随机颗粒，压缩颗粒至高紧凑状态，根据颗粒中心间距与两颗粒的半径之和，获取处于接触状态的第二颗粒，并在第二颗粒的内部生成接触键；

基于接触键的法向单位向量、切向单位向量，获取第二颗粒的合力向量及弯矩，以及第二颗粒在第一单位时间的加速度和角加速度，根据显示中心差分法，迭代更新第二颗粒在第二单位时间的位置、速度和角速度；

根据迭代更新的结果，获取第二颗粒的接触作用力，其中，通过将接触作用力的法向力和剪切力，与目标混凝土结构件的抗拉强度、抗压强度、剪切强度进行对比，获取目标混凝土结构件的损伤情况。

进一步优选地，在压缩颗粒至高紧凑状态的过程前，采用控制孔隙度方法生成颗粒的颗粒数量，包括：

其中，n表示生成目标粒子数，V表示给定空间总体积，C表示孔隙度，R表示平均半径；体积分数在0.615-0.735之间。

进一步优选地，在压缩颗粒至高紧凑状态的过程中，通过弹性半无线空间压板压缩颗粒，压缩速度控制在0.0001m/s-0.0012m/s之间；

压缩结束后，颗粒与压缩板之间的接触力小于0.1N；

高紧凑状态的判断依据为：检测颗粒间的配位数是否在4.01-8.15之间，判断颗粒是否处于高紧凑状态，以及检查颗粒的重叠性。

进一步优选地，在第二颗粒的内部生成接触键的过程中，令颗粒间的内聚力处于零应力状态，两个颗粒之间存在接触键需满足以下条件：

其中，d

进一步优选地，在获取第二颗粒的合力向量及弯矩的过程中，法向单位向量用

切向单位向量用

其中，Δp

进一步优选地，在获取第二颗粒的合力向量及弯矩的过程中，利用法向力和切向力的单位向量分别点乘相应力的标量，然后相加得到合力向量；

弯矩的计算公式为

进一步优选地，法向力为：

切向力为：

颗粒间的法向刚度为：

颗粒间的切向刚度为K

进一步优选地，在获取加速度和角加速度的过程中，根据合力向量及弯矩，通过牛顿第二定律或达朗贝尔原理，获取第二颗粒在第一单位时间的加速度和角加速度。

进一步优选地，在根据显示中心差分法进行迭代更新的过程中，包括以下步骤：

设置half_t＝dt*(0：N)+dt/2，t＝dt*(0：N)，其中，dt为时间增量，N为大于0的正整数；

第二颗粒在第一单位时间t的加速度为：

第二颗粒在第一单位时间t的角加速度为：

第二颗粒在第二单位时间t+1的加速度：

第二颗粒在第二单位时间t+1的角速度：

进一步优选地，目标混凝土结构件的损伤情况包括，法向力判别公式和切向力判别公式；

法向力判别公式为：

F

-F

切向力判别公式为：

F

其中，α

本发明还包括一种基于离散元的混凝土结构破坏过程数值模拟系统，用于实现上述模拟方法，包括：

数据采集与处理模块，用于基于目标混凝土结构件，采集目标混凝土结构件的边界墙尺寸，构建三维刚性边界墙，其中，三维刚性边界墙的墙体单元采用平面单元；

接触键生成模块，用于在三维刚性边界墙的墙体内部空间生成随机颗粒，压缩颗粒至高紧凑状态，根据颗粒中心间距与两颗粒的半径之和，获取处于接触状态的第二颗粒，并在第二颗粒的内部生成接触键；

数据运算模块，用于基于接触键的法向单位向量、切向单位向量，获取第二颗粒的合力向量及弯矩，以及第二颗粒在第一单位时间的加速度和角加速度，根据显示中心差分法，迭代更新第二颗粒在第二单位时间的位置、速度和角速度；

损伤评估模块，用于根据迭代更新的结果，获取第二颗粒的接触作用力，其中，通过将接触作用力的法向力和剪切力，与目标混凝土结构件的抗拉强度、抗压强度、剪切强度进行对比，获取目标混凝土结构件的损伤情况。

实施例1：本发明的本发明的具体实现方法，包括以下步骤：

设置几何和物理两组参数。几何参数包括迭代计算次数N、时间增量dt、颗粒质量、密度、半径、转动惯量等。物理参数包括弹性模量、泊松比、摩擦角、粘聚力、局部(全局)阻尼系数等。其中，微观参数法相刚度K

确定混凝土结构构件整体尺寸，确定边界墙尺寸，建立三维刚性边界墙，墙体单元采用平面单元。

在墙体内部空间生成随机颗粒，采用控制孔隙度方法生成颗粒数量，计算公式为

内部生成颗粒后，通过弹性半无线空间压板压缩颗粒，压缩速度控制在0.0001m/s-0.0012m/s之间。压缩结束时颗粒与压缩板之间的接触力控制在0.1N以下。待模型达到平衡状态，检查颗粒的两个指标。第一，检测颗粒间的配位数是否在4.01-8.15之间，判断颗粒是否处于高紧凑状态；第二，检查颗粒的重叠性。若以上满足条件，记录所有颗粒的位置坐标及初始速度和角速度。若不满足条件，返回上一步调整。

确定接触检测方法，根据颗粒中心间距与两颗粒的半径之和做比较，若小于半径之和，则认定为两个颗粒接触；反之不接触。利用布尔矩阵记录判断结果，被判定为接触的颗粒，可以进行下一步计算接触力，被判定为不接触的颗粒不需要计算接触力。

确定接触模型。确定好颗粒的数量和排列后，根据以下方程生成颗粒内部的接触键。在创建初始时，令颗粒间的内聚力处于零应力状态。两个颗粒之间存在接触键需满足以下条件：

法向力和切向力的具体求解方法。利用材料的拉伸、压缩和剪切刚度等相关损伤折减系数定义法向力F

颗粒利用法向和切的单位向量分别点乘相应力的标量，然后相加得到合力向量。弯矩计算公式为

获得力和弯矩后，利用牛顿第二定律或达朗贝尔原理，获得在经过时间dt后颗粒的加速度

利用显示中心差分法，通过以下公式迭代更新颗粒的位置、速度和角速度。设half_t＝dt*(0：N)+dt/2，t＝dt*(0：N)其中dt为时间增量，N为大于0的正整数。则有，在t时刻颗粒的速度和角速度为：

t+1时刻的速度和角速度为：

对t+1时刻的速度和角速度进行积分，得到该时刻的颗粒位置坐标。此过程反复迭代N次。

导入损伤模型。通过以上方程计算任意颗粒间的接触作用力，根据迭代计算的法向力和剪切力与材料的抗拉强度、抗压强度、剪切强度对比，判断材料的损伤情况。法向力和切向力的判别公式如下：

F

-F

F

其中，α

记录模型的相关模拟数据结果，用于后续具体分析。

实施例2：

采用本发明方法对无筋混凝土柱压缩进行数值模拟。

根据设计尺寸建立刚性单元边界墙。梁宽度统一设置为200mm，长高比例关系用a表示。利用孔隙度法在墙体内部生成平均半径约为4.5mm的颗粒堆积体模型。利用重力及弹性半无线压力板使颗粒堆积到一起，待模型平衡后，发生接触颗粒创建接触键。在梁中心位置施加线荷载。颗粒的物理参数如表1所示，模型颗粒数和接触键数如表2所示。

表1

表2

模型截面示意图如图1所示。模型的应力应变曲线如图2所示。

本发明提出混凝土破坏过程数值模拟方法，克服了现有数值模拟方法局限性问题，且该方法可准确、快速预测混凝土从细观到宏观结构构件或结构的破坏过程。