运用多尺度离散实体元解决连续体结构力学问题的仿真方法与流程

1.本发明涉及一种多尺度离散实体单元法力学仿真方法，用于解决结构或构件的大变形、损伤断裂和倒塌破坏等强非线性力学问题，属于离散元分析技术领域。


背景技术：

2.离散单元法是美国学者cundall于1971年提出来的一种非连续性数值计算方法
[1
‑
2]
，最初被用于分析岩石边坡运动的分析计算。([1]cundall p a,strack o d l.a discrete numerical model for granular assemblies[j].geotechnique,1979,29(1):47
‑
65.[2]cundall p a,hart r d.numerical modelling of discontinua[j].engineering computations,1992,9(2):101
‑
113.)离散单元法的基本原理为：把研究对象划分成一个个离散的单元，在结构受力变形、运动过程中，单元可以分离，即一个单元与其邻近单元可以接触，也可以分开
[3
‑
4]
。([3]邢继波,王泳嘉.离散元法的改进及其在颗粒介质研究中的应用[j].岩土工程学报,1990,12(5):51
‑
57.[4]王泳嘉,邢纪波.离散单元法及其在岩土力学中的应用.东北工学院出版社[m].1991.)传统离散单元法主要研究对象为散粒体，对连续体结构不适用。近年来开始有学者尝试将离散单元法用于钢结构倒塌、稳定、冲击破坏问题的数值计算。叶继红课题组
[5
‑
7]
提出了杆系结构的离散元模型，推导了模型中弹簧的刚度系数，建立了适用于杆系结构的离散单元法塑性铰模型、考虑截面塑性开展的纤维模型、断裂准则和单元碰撞模型。([5]齐念,叶继红.基于离散元法的杆系结构几何非线性大变形分析[j].东南大学学报(自然科学版),2013(5):917
‑
922.[6]齐念.dem/fem耦合计算方法研究及其在网壳倒塌破坏模拟中的应用[d].南京:东南大学,2016.[7]覃亚男.基于离散单元法的简单结构静动力响应数值模拟研究[d].南京:东南大学,2015.)。
[0003]
发明人所在的冯若强课题组
[8
‑
10]
提出了离散实体单元法，改进了球元间接触弹簧的设置规则，在球元立方体排列模型的棱边和面对角线上采用弹簧进行连接，球元间弹簧包括法向弹簧和两个方向切线弹簧，使其能够准确地体现连续体结构的力学性能，同时可以考虑材料的泊松比效应。([8]zhu b,feng r,wang x.3d discrete solid
‑
element method for elastoplastic problems of continuity[j].journal of engineering mechanics,2018,144(7):04018051.[9]zhu b,feng r.discrete solid element model applied to plasticity and dynamic crack propagation in continuous medium[j].computational particle mechanics,2019:1
‑
17.[10]zhu b,feng r.investigation of a boundary simulation of continuity using the discrete solid element method[j].advances in mechanical engineering,2019,11(1):1687814018822397.)
[0004]
综上所述，现有的离散实体单元法已可以对简单的连续体结构进行力学分析，但在解决受力复杂的复杂结构还存在球元数目大、计算效率低等问题。因此，迫切需要多尺度离散元计算方法以高效模拟受力复杂的复杂结构力学问题。


技术实现要素：

[0005]
发明目的：针对传统离散单元法不能很好模拟连续介质，现有的离散实体单元法已可以对简单的连续体结构进行力学分析，但在解决受力复杂的复杂结构还存在球元数目大、计算效率低等问题。本发明提供了一种多尺度离散实体单元法，基于能量等效原则，构建了多尺度离散实体元模型，推导了细观模型参数(弹簧刚度)与宏观弹性常数(弹性模量和泊松比)之间的解析关系，从而解决受力复杂的复杂结构采用传统离散实体元法计算中球元数目大、计算效率低等问题。
[0006]
技术方案：为了实现上述发明目的，本发明提出的多尺度离散实体单元法力学仿真方法，步骤如下：
[0007]
步骤1，根据仿真需求，将结构或构件离散为多尺度离散元模型，建立模型的球元为规则排列，材料的变形完全由球元间的接触弹簧来存储和表示；
[0008]
多尺度离散实体元模型包括大球区域、小球区域和过渡区域，球元之间通过弹簧系统连接，弹簧系统包括轴向弹簧和两个方向的切向弹簧。
[0009]
步骤2，t＝0，t≥0时刻，对各颗粒单元的外力、内力、位移、速度和加速度赋初值。
[0010]
步骤3，对所有颗粒单元：计算t时刻的外力及阻尼力，应用牛顿第二定律求解运动控制方程，得到t δt时刻的颗粒单元的位置与速度。
[0011]
步骤4，计算位移增量δu，并计算试探接触力增量从而得到t δt的接触力f，根据屈服准则判断当前的颗粒单元的弹塑性状态。
[0012]
步骤5，根据弹塑性接触本构方程，计算δt时步内颗粒单元的接触力增量δf，并更新当前时刻的颗粒单元间的接触力。
[0013]
步骤6，根据平衡方程，计算各颗粒单元的内力、阻尼力以及外力，为t 2δt时刻的计算准备，如此循环往复，直至模型达到稳定状态。
[0014]
具体来说，包括以下具体步骤：
[0015]
(1)初步估计仿真对象受力复杂的区域
[0016]
根据力学概念、经验估计仿真对象受力复杂的区域，确定大球、小球球元的半径，将受力复杂区域划分为小球区域，将受力较简单的区域划分为大球区域，小球区域和大球区域之间为过渡区域。
[0017]
(2)将仿真对象离散为多尺度离散元模型
[0018]
根据步骤1确定的大球区域、小球区域和过渡区域，建立仿真对象的多尺度离散元模型，大球区域、小球区域和过渡区域球元按各自区域的规律排列，材料的变形完全由单元的接触弹簧来存储和表示，输入模型参数以及外荷载信息，确定边界条件。
[0019]
大球区域弹簧刚度可表示为：
[0020]
法向弹簧刚度：
[0021]
切向弹簧刚度：
[0022]
其中，为大球区域棱弹簧法向弹簧刚度，为大球区域面对角弹簧法向弹簧刚度，度，分别为大球区域棱弹簧两个方向切向弹簧刚度，分别为大球区域面对角
弹簧两个方向切向弹簧刚度，r为大球半径，e为弹性模量，μ为材料的泊松比。
[0023]
小球区域弹簧刚度可表示为：
[0024]
法向弹簧刚度：
[0025]
切向弹簧刚度：
[0026]
其中，为小球区域棱弹簧法向弹簧刚度，为小球区域面对角弹簧法向弹簧刚度，度，分别为小球区域棱弹簧两个方向切向弹簧刚度，分别为小球区域面对角弹簧两个方向切向弹簧刚度，r为小球半径，e为弹性模量，μ为材料的泊松比。
[0027]
过渡区域弹簧刚度可表示为：
[0028]
法向弹簧刚度：
[0029]
切向弹簧刚度：切向弹簧刚度：
[0030]
其中，分别为过渡区域与大球区域、小球区域相连部分棱弹簧轴向刚度；分别为过渡区域与大球区域、小球区域相连部分面对角弹簧轴向刚度；分别为过渡区域棱弹簧和面对角弹簧轴向刚度；分别为过渡区域与大球区域、小球区域相连部分棱弹簧两个方向切向刚度；分别为过渡区域与大球区域、小球区域相连部分面对角弹簧两个方向切向刚度；分别为过渡区域棱弹簧、面对角弹簧两个方向切向刚度；r为小球半径，r为小球半径，e为弹性模量，μ为材料的泊松比。
[0031]
(3)t＝0时刻，对球元的速度、加速度等赋初值，对弹簧系统的内力赋初值；
[0032]
(4)对所有球元，由t时刻的外力及阻尼力，应用牛顿第二定律，求得t δt时刻球元的位置与速度；
[0033]
任取一个球元i，设球元i上作用了n个外力，有m个球元与其相邻，作用在球元i上的外力为相邻球元对其作用力之和为球元i阻尼力为f
ic
。根据牛顿第二定律，其运动控制方程如公式(1)所示。
[0034][0035]
其中，m
i
为球元i的质量；为球元i的加速度；为球元i上第j个外力；为球元i上第j个接触内力；f
ic
为球元i的阻尼力。
[0036]
在公式(1)的基础上，运用中心有限差分算法确定球元球心的加速度，进而确定球
元球心的位移。计算球元i在t
△
t时刻的位移如公式(2)所示。
[0037][0038]
其中，u
i
(t δt)为球元i在t δt时刻的位移；u
i
(t)为球元i在t时刻的位移；为球元i在时刻的加速度。
[0039]
(5)计算时步δt内，弹簧法向和切向接触力增量与法向和切向位移增量的关系如公式(3)所示。采用坐标转换矩阵将弹簧的法向和切向接触力增量转换为整体坐标下接触力增量，如公式(4)所示。采用增量的形式计算球元间的接触力，t δt时刻与t时刻整体坐标下弹簧接触力的关系如公式(5)所示。
[0040][0041][0042][0043]
(6)实体单元法计算模型中，一个球元与若干个球元通过弹簧连接，弹簧接触力和外力之和即为该球元在t δt时刻的不平衡力，这个不平衡力将用于求解下一计算时步t 2δt球元的运动计算，进而求得球元的加速度、速度和位移。
[0044]
有益效果：
[0045]
(1)本发明基于连续介质力学理论以能量等效为原则，得到细观模型参数(弹簧刚度)和宏观弹性常数(弹性模量和泊松比)之间的解析关系。有效减少球元和弹簧的数量，提高计算效率，使离散实体元法提供具有更多可能性，适用于结构或构件的大变形、损伤断裂和倒塌破坏等复杂力学仿真。
[0046]
(2)在多尺度离散实体单元法中，无需组集刚度矩阵，进行迭代求解，计算流程简单清晰。结构越复杂，单元数量就越多，相应的计算量也越大，但只是循环重复计算，不会给计算带来实质性困难。
[0047]
(3)传统方法相比，在大型结构和复杂力学行为的模拟中离散实体单元法将表现出更大的优势。构建的三维离散实体元模型能够体现泊松比效应，实现不同尺度球元的连接。
附图说明
[0048]
图1为本发明的流程图；
[0049]
图2为多尺度离散实体单元法物理模型示意；
[0050]1‑
大球区域面弹簧、2
‑
大球区域棱弹簧、3
‑
过渡区域棱弹簧、4
‑
小球区域面弹簧、5
‑
小球区域棱弹簧、6
‑
过渡区域棱弹簧、7
‑
过渡区域面弹簧；
[0051]
a
‑
大球区域、b
‑
过渡区域、c
‑
小球区域；
[0052]
a
‑
端部固定约束、b
‑
端部大球固定约束；
[0053]
图3为多尺度离散实体单元法物理模型分区示意；
[0054]
图4为离散实体元模型；
[0055]
图5为多尺度离散实体元模型；
[0056]
图6为离散实体元模型变形图；
[0057]
图7为多尺度离散实体元模型变形图；
[0058]
图8为离散实体元(dsem)、多尺度离散实体元(mdsem)与有限元模型(fem)变形对比图。
具体实施方式
[0059]
以下将结合实施例具体说明本发明的技术方案：
[0060]
参照图1所示，本发明运用多尺度离散实体元解决连续体结构力学问题的仿真方法，步骤如下：
[0061]
(1)初步估计仿真对象受力复杂的区域，根据力学概念、经验估计仿真对象受力复杂的区域，确定大球、小球球元的半径，将受力复杂区域划分为小球区域，将受力较简单的区域划分为大球区域，小球区域和大球区域之间为过渡区域。
[0062]
(2)将仿真对象离散为多尺度离散元模型，建立仿真对象的多尺度离散元模型，大球区域、小球区域和过渡区域球元按各自区域的规律排列，材料的变形完全由单元的接触弹簧来存储和表示，输入模型参数以及外荷载信息，确定边界条件，参照图2、3所示。
[0063]
(3)t＝0时刻，对球元的速度、加速度等赋初值，对弹簧系统的内力赋初值；
[0064]
(4)对所有球元，由t时刻的外力及阻尼力，应用牛顿第二定律，求得t δt时刻球元的位置与速度；
[0065]
(5)计算时步δt内法向和切向位移增量，采用坐标转换矩阵将弹簧的法向和切向接触力增量转换为整体坐标下接触力增量，计算球元间的接触力。
[0066]
(6)实体单元法计算模型中，一个球元与若干个球元通过弹簧连接，弹簧接触力和外力之和即为该球元在t δt时刻的不平衡力，这个不平衡力将用于求解下一计算时步t 2δt球元的运动计算，进而求得球元的加速度、速度和位移。
[0067]
实施例1
[0068]
下面用算例验证多尺度离散实体元在模拟构件大变形的优越性。
[0069]
悬臂梁长160mm，截面尺寸为30mm
×
30mm，弹性模量为2.06
×
105mpa，泊松比为0.2，密度7850kg/m3，梁端部施加剪切力1000mpa。
[0070]
离散实体元模型球元半径为2.5mm，球元数量1617，连接弹簧数量12092
×
3，如图4所示。多尺度离散实体元模型大球区域长80mm，大球直径10mm，小球区域长70mm，小球直径5mm，过渡区域长10mm，球元数量879，连接弹簧数量6332
×
3，如图5所示。本实施例考察了悬臂梁端部在弯矩m作用下y向位移的位移变化情况。
[0071]
图6和图7为离散实体元模型和多尺度离散实体元模型剪力作用下变形。离散实体元、多尺度离散实体元与有限元计算结果对比如表1所示。
[0072]
图8为离散实体元、多尺度离散实体元与有限元模型变形对比图。比较离散实体元
模型，多尺度离散实体元模型的球元数量和连接弹簧数量下降约50％，计算效率提高50％，但计算精度下降约2％，表明多尺度离散实体元在满足精度的前提下，计算效率有显著提高。
[0073]
表1悬臂梁受剪力(y＝0，z＝0)直线上各点y向位移计算结果
[0074]