CN112329290A - 可用于施工过程模拟的有限元离散元多尺度耦合计算方法 - Google Patents

可用于施工过程模拟的有限元离散元多尺度耦合计算方法 Download PDF

Info

Publication number
CN112329290A
CN112329290A CN202011169620.XA CN202011169620A CN112329290A CN 112329290 A CN112329290 A CN 112329290A CN 202011169620 A CN202011169620 A CN 202011169620A CN 112329290 A CN112329290 A CN 112329290A
Authority
CN
China
Prior art keywords
discrete element
finite element
dem
construction process
discrete
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
Application number
CN202011169620.XA
Other languages
English (en)
Inventor
孟庆祥
王环玲
徐卫亚
赵兰浩
毛佳
曹亚军
闫龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202011169620.XA priority Critical patent/CN112329290A/zh
Publication of CN112329290A publication Critical patent/CN112329290A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • 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
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • 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
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了可用于施工过程模拟的有限元离散元多尺度耦合计算方法，属于岩土体工程技术领域，本发明革新了现有的有限元离散元多尺度耦合计算框架，构建了可变的离散元子结构序列Sdem，使之与施工过程中宏观有限元网格的减少或者增加对应起来，从而实现施工过程的模拟。与现有技术相比，本发明扩展了有限元离散元多尺度耦合计算方法的使用范围，可以更好地预测复杂岩土体材料在开挖、堆载等过程中的力学行为，为相关岩土石力学研究服务和工程安全性评价提供技术支撑。

Description

可用于施工过程模拟的有限元离散元多尺度耦合计算方法
技术领域
本发明属于岩土体工程技术领域，具体涉及可用于施工过程模拟的有限元离散元多尺度耦合计算方法。
背景技术
均匀化多尺度模拟方法是通过将具有复杂结构的不均匀材料假设为均质材料，基于不均匀介质等效宏观特性来开展分析重要方法。均匀化分析方法可以分为三大类：解析、数值和计算均匀化方法。
解析均匀化通过严格的公式推导，建立反映不均匀介质材料宏观特性的数学表达。Hill最早提出RVE的概念和基于能量等效原理(Hill-Mandel Macro-homogeneityPrinciple)的复合材料弹性特征。解析多尺度需要有较强的理论基础和复杂的公式推导，仅适用于具有规则细观结构的不均匀材料，因此仅适用于内部结构较为简单的不均匀介质。
数值均匀化弥补了解析多尺度的不足，通过对不均匀介质标准体积单元(RVE)进行分析，针对表征体积单元开展均匀化分析，拟合确定现有的宏观本构模型的基本参数。目前，常见的数值均匀化实现方式较多，主要有快速傅里叶变换FFT、有限单元法FEM、离散单元法DEM和边界元BEM。通过连接细观结构和宏观本构模型，数值均匀化方法在解决大尺度工程问题方面具有重大应用潜力。
相对于解析和数值均匀化方法，计算均匀化方法不需要预先定义宏观本构模型，应力和应变关系直接通过数值计算得到。因此，该方法具有极大的灵活性，既适用于弹塑性问题，也适用于多场耦合及大变形等复杂问题。笔者基于开源离散元程序Yade和自编有限元程序，实现了FEM×DEM，开展了散体材料多尺度研究，并应用于边坡安全系数分析中。但是目前所有的研究方法网格是固定的，无法模拟开挖、堆载等施工过程。本专利克服了这一问题，提出了一种可用于施工过程模拟的有限元离散元多尺度耦合计算方法。
发明内容
发明目的：本发明的目的在于提供可用于施工过程模拟的有限元离散元多尺度耦合计算方法，
技术方案：为实现上述目的，本发明提供如下技术方案：
可用于施工过程模拟的有限元离散元多尺度耦合计算方法，包括如下步骤：
(1)建立宏观有限元网格模型与对应的离散元子模型序列Sdem；
(2)根据工程边界条件和材料重度，进行自重平衡计算初始状态；
(3)施工过程模拟。
进一步地，所述的步骤(1)中，具体步骤为：根据工程实际建立宏观有限元网格模型Mfem和表征材料力学行为的离散元子模型Mdem，给定子模型离散元的材料参数。
进一步地，针对宏观有限元模型的单元类型，确定每个单元高斯积分点的数量，将每个高斯点对应一个所述的离散元子模型Mdem，将多个离散元子模型Mdem组成1个序列Sdem，该序列根据宏观网格的变化进行更改。
进一步地，所述的步骤(2)包括如下步骤：
1)如果初始应力状态已知，将每一个高斯点处的应力作为边界施加给离散元子模型Mdem，通过伺服加载实现；
2)如果初始应力状态未知，通过如下公式(I)计算自重，转换为宏观模型上的节点力F；
F＝γ∫∫[N]Tdxdydz (I)；
其中γ为材料重度，N为形函数，T为转置运算符；
3)根据离散元子模型Mdem的初始状态，采用均匀化公式(II)得到材料的刚度，采用传统有限元按照弹性计算得到初始应力状态，然后采用步骤1)的方式把初始应力施加到离散元子模型Mdem上；
Figure BDA0002746876060000021
其中kn为法向刚度，ks为切向刚度，nc为法向接触力的法线方向，tc切向接触力的切线方向，D为切线刚度矩阵。
进一步地，所述的步骤(3)包括如下步骤：
1)随着施工过程的开展，宏观有限元网格会减少或者增加，确定对应减少或者增加网格对应的离散元子结构模拟，在原有的离散元子结构序列中增加或者删除离散元子结构模拟；
2)根据高斯点上的应变增量转换为离散元子结构的周期性边界条件，采用并行计算的方式对序列Sdem开展计算，根据公式(III)把离散元子模型计算结果得到积分点的应力；
Figure BDA0002746876060000031
其中V是离散元子模型的体积DEM，Nc是接触数量，dc是接触颗粒中心向量，fc是接触力的向量；
3)通过公式(IV)计算得到的节点反力与外力差值判断是否收敛，如果收敛输出宏观有限元的应力、应变、位移等信息和细观离散元的接触力、粘结断裂数等结果，否则进行下一个循环，从而实现施工过程的模拟；
R＝∫ΩBTσdΩ-F (IV)；
其中B为有限元单元的B矩阵，σ为均匀画应力，R为残余节点荷载。
有益效果：与现有技术相比，本发明根据宏观有限元网格会减少或者增加，对应减少或者增加网格对应的离散元子结构，扩展了有限元离散元多尺度耦合计算方法的使用范围，可以更好地预测复杂岩土体材料在开挖、堆载等过程中的力学行为，为相关岩土石力学研究服务和为工程安全性评价提供技术支撑。
附图说明
图1是本发明实施例的总体流程图；
图2是实施例中堆载三个过程的宏观有限元模型；
图3是实施例中每个高斯点对应的离散元子子结构模型；
图4是实施例中初始状态下的宏观位移场和细观颗粒断裂分布；
图5是实施例中第一步堆载后的宏观位移场和细观颗粒断裂分布；
图6是实施例中第二步堆载后的宏观位移场和细观颗粒断裂分布；
图7是实施例中第三步堆载后的宏观位移场和细观颗粒断裂分布。
具体实施方式
下面结合附图和具体实施例对本发明作更进一步的说明。
可用于施工过程模拟的有限元离散元多尺度耦合计算方法，包括如下步骤：
(1)宏观模型与对应的离散元子模型序列
1)根据工程实际建立宏观有限元网格模型Mfem和可以表征材料力学行为的离散元子模型Mdem，给定子模型离散元的材料参数；
2)针对宏观有限元模型的单元类型，确定每个单元高斯积分点的数量，将每个高斯点对应一个离散元子模型Mdem，将多个离散元子模型Mdem组成1个序列Sdem，该序列可以根据宏观网格的变化进行更改；
(2)初始状态计算
1)根据工程边界条件和材料重度，进行自重平衡计算，如果初始应力状态已知，可以直接将每一个高斯点处的应力作为边界施加给离散元子模型Mdem，通过伺服加载实现；
2)如果初始应力状态未知，通过如下公式(I)计算自重，转换为宏观模型上的节点力F；
F＝γ∫∫[N]Tdxdydz (I)；
其中γ为材料重度，N为形函数，T为转置运算符；
3)根据离散元子模型Mdem的初始状态，采用均匀化公式(II)得到材料的刚度，采用传统有限元按照弹性计算得到初始应力状态，然后采用步骤1)的方式把初始应力施加到离散元子模型Mdem上；
Figure BDA0002746876060000041
其中kn为法向刚度，ks为切向刚度，nc为法向接触力的法线方向，tc切向接触力的切线方向，D为切线刚度矩阵。
(3)施工过程模拟
1)随着施工过程的开展，宏观有限元网格会减少或者增加，确定对应减少或者增加网格对应的离散元子结构模拟，在原有的离散元子结构序列中增加或者删除离散元子结构模拟；
2)根据高斯点上的应变增量转换为离散元子结构的周期性边界条件，采用并行计算的方式对序列Sdem开展计算，根据公式(III)把离散元子模型计算结果得到积分点的应力；
Figure BDA0002746876060000042
其中V是离散元子模型的体积DEM，Nc是接触数量，dc是接触颗粒中心向量，fc是接触力的向量；
3)通过公式(IV)计算得到的节点反力与外力差值判断是否收敛，如果收敛输出宏观有限元的应力、应变、位移等信息和细观离散元的接触力、粘结断裂数等结果，否则进行下一个循环，从而实现施工过程的模拟。
R＝∫ΩBTσdΩ-F (IV)；
其中B为有限元单元的B矩阵，σ为均匀画应力，R为残余节点荷载。
实施例
如图1所示，本实施例应用Python语言开发了有限元分析程序，结合开源离散元软件Yade开发了一种可用于施工过程模拟的有限元离散元多尺度耦合计算程序，针对常见的堆载施工过程开展分析，包括如下步骤：
(1)宏观模型与对应的离散元子模型序列
1)根据工程实际，在初始长度为60m高度为5m的路基上建设堤坝，假定堤坝高度为10m，分三次填筑，第一次填筑高度为2m，第二次为5m，第三次为10m，该模型总共有350个四边形单元(如图2所示)；
2)假定每个单元有四个高斯积分点，每个高斯积分点对应一个如图3所示的离散元子结构模型；
(2)初始状态计算
1)根据工程边界条件和材料重度，进行自重平衡计算，本实施例假定初始应力状态未知，通过公式(I)计算自重应力情况下宏观模型上的节点力F，初始模型有150个单元，对应600个离散元子结构；
F＝γ∫∫[N]Tdxdydz (I)；
2)根据离散元子模型Mdem的初始状态，采用均匀化公式(II)得到材料的刚度，采用传统有限元按照弹性计算得到初始应力状态，然后把初始应力施加到离散元子模型Mdem上；
Figure BDA0002746876060000051
其中kn为法向刚度，ks为切向刚度，nc为法向接触力的法线方向，tc切向接触力的切线方向，D为切线刚度矩阵。
3)计算得到的位移场分布如图4(a)所示，细观离散元颗粒粘结断裂数如图4(b)所示。
(3)施工过程模拟
1)当第一步施工时，宏观有限元网格增加到190个，更新网格对应的离散元子结构数量为760，形成离散元子结构序列；
2)将新添加的单元产生的荷载采用公式(I)计算得到施加了宏观模型上，计算得到位移的增量和高斯点应变增量，根据高斯点上的应变增量转换为离散元子结构的周期性边界条件，采用并行计算的方式对离散元子结构序列Sdem开展计算，采用公式(III)将离散元子模型计算结果转换为积分点的应力；
Figure BDA0002746876060000061
其中V是离散元子模型的体积DEM，Nc是接触数量，dc是接触颗粒中心向量，fc是接触力的向量；
3)通过公式(IV)计算得到的节点反力与外力差值R判断是否收敛，如果收敛输出宏观有限元的应力、应变、位移等信息和细观离散元的接触力、粘结断裂数等结果，否则进行下一个循环，从而实现施工过程的模拟。
R＝∫ΩBTσdΩ-F (IV)；
其中B为有限元单元的B矩阵，σ为均匀画应力，R为残余节点荷载。
4)经过19此迭代计算后，残余节点力的最大值与上一步的比值小于1％，认为收敛，输出宏观模型的位移和细观结构的颗粒粘结断裂数如图5所示；堆载高度2m时，计算得到的位移场分布如图5(a)所示，颗粒粘结断裂数如图5(b)所示。
5)采用相同的流程，分别计算得到堆载5m的结果如图6所示，计算得到的位移场分布如图6(a)所示，颗粒粘结断裂数如图6(b)所示。堆载高度为10m的结果如图7所示，计算得到的位移场分布如图7(a)所示，颗粒粘结断裂数如图7(b)所示。可见高度为10m时，颗粒粘结的断裂数量较多，形成一条明显的局部化带，此时如果继续填筑会有一定的安全隐患。
需要说明的是，以上说明仅是本发明的优选实施方式，应当理解，对于本领域技术人员来说，在不脱离本发明技术构思的前提下还可以做出若干改变和改进，这些都包括在本发明的保护范围内。

Claims (5)

1.可用于施工过程模拟的有限元离散元多尺度耦合计算方法，其特征在于：包括如下步骤：
(1)建立宏观有限元网格模型与对应的离散元子模型序列Sdem；
(2)根据工程边界条件和材料重度，进行自重平衡计算初始状态；
(3)施工过程模拟。
2.根据权利要求1所述的可用于施工过程模拟的有限元离散元多尺度耦合计算方法，其特征在于：所述的步骤(1)中，具体步骤为：根据工程实际建立宏观有限元网格模型Mfem和表征材料力学行为的离散元子模型Mdem，给定子模型离散元的材料参数。
3.根据权利要求2所述的可用于施工过程模拟的有限元离散元多尺度耦合计算方法，其特征在于：针对宏观有限元模型的单元类型，确定每个单元高斯积分点的数量，将每个高斯点对应一个所述的离散元子模型Mdem，将多个离散元子模型Mdem组成1个序列Sdem，该序列根据宏观网格的变化进行更改。
4.根据权利要求1所述的可用于施工过程模拟的有限元离散元多尺度耦合计算方法，其特征在于：所述的步骤(2)包括如下步骤：
1)如果初始应力状态已知，将每一个高斯点处的应力作为边界施加给离散元子模型Mdem，通过伺服加载实现；
2)如果初始应力状态未知，通过如下公式(I)计算自重，转换为宏观模型上的节点力F；
F=γ∫∫[N]Tdxdydz (I)；
其中γ为材料重度，N为形函数，T为转置运算符；
3)根据离散元子模型Mdem的初始状态，采用均匀化公式(II)得到材料的刚度，采用传统有限元按照弹性计算得到初始应力状态，然后采用步骤1)的方式把初始应力施加到离散元子模型Mdem上；
Figure FDA0002746876050000011
其中kn为法向刚度，ks为切向刚度，nc为法向接触力的法线方向，tc切向接触力的切线方向，D为切线刚度矩阵。
5.根据权利要求1所述的可用于施工过程模拟的有限元离散元多尺度耦合计算方法，其特征在于：所述的步骤(3)包括如下步骤：
1)随着施工过程的开展，宏观有限元网格会减少或者增加，确定对应减少或者增加网格对应的离散元子结构模拟，在原有的离散元子结构序列中增加或者删除离散元子结构模拟；
2)根据高斯点上的应变增量转换为离散元子结构的周期性边界条件，采用并行计算的方式对序列Sdem开展计算，根据公式(III)把离散元子模型计算结果得到积分点的应力；
Figure FDA0002746876050000021
其中V是离散元子模型的体积DEM，Nc是接触数量，dc是接触颗粒中心向量，fc是接触力的向量；
3)通过公式(IV)计算得到的节点反力与外力差值判断是否收敛，如果收敛输出宏观有限元的应力、应变、位移等信息和细观离散元的接触力、粘结断裂数，否则进行下一个循环，实现施工过程的模拟；
R＝∫ΩBTσdΩ-F (IV)；
其中B为有限元单元的B矩阵，σ为均匀画应力，R为残余节点荷载。
CN202011169620.XA 2020-10-28 2020-10-28 可用于施工过程模拟的有限元离散元多尺度耦合计算方法 Pending CN112329290A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011169620.XA CN112329290A (zh) 2020-10-28 2020-10-28 可用于施工过程模拟的有限元离散元多尺度耦合计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011169620.XA CN112329290A (zh) 2020-10-28 2020-10-28 可用于施工过程模拟的有限元离散元多尺度耦合计算方法

Publications (1)

Publication Number Publication Date
CN112329290A true CN112329290A (zh) 2021-02-05

Family

ID=74296127

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011169620.XA Pending CN112329290A (zh) 2020-10-28 2020-10-28 可用于施工过程模拟的有限元离散元多尺度耦合计算方法

Country Status (1)

Country Link
CN (1) CN112329290A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113177248A (zh) * 2021-04-21 2021-07-27 武汉大学 隧道围岩破裂碎胀大变形失稳灾变过程2d-fdem数值模拟方法
CN113204902A (zh) * 2021-04-21 2021-08-03 武汉大学 一种有限元-离散元耦合(fdem)恒阻大变形锚杆隧道围岩加固数值模拟方法
CN116070790A (zh) * 2023-03-21 2023-05-05 中国建筑一局(集团)有限公司 一种减少施工现场重大安全风险隐患的预测方法和系统
CN116226982A (zh) * 2023-01-31 2023-06-06 武汉大学 一种粘性土-岩石隧道开挖耦合数值方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107169236A (zh) * 2017-06-15 2017-09-15 中国水利水电科学研究院 一种基于有限元与离散元耦合的虚拟三轴试验仿真方法
CN108984829A (zh) * 2018-06-06 2018-12-11 中国农业大学 堆石混凝土堆石体堆积过程的计算方法和系统

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107169236A (zh) * 2017-06-15 2017-09-15 中国水利水电科学研究院 一种基于有限元与离散元耦合的虚拟三轴试验仿真方法
CN108984829A (zh) * 2018-06-06 2018-12-11 中国农业大学 堆石混凝土堆石体堆积过程的计算方法和系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Q.X. MENG 等: "Multiscale Strength Reduction Method for Heterogeneous Slope Using Hierarchical FEM/DEM Modeling", 《COMPUTERS AND GEOTECHNICS》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113177248A (zh) * 2021-04-21 2021-07-27 武汉大学 隧道围岩破裂碎胀大变形失稳灾变过程2d-fdem数值模拟方法
CN113204902A (zh) * 2021-04-21 2021-08-03 武汉大学 一种有限元-离散元耦合(fdem)恒阻大变形锚杆隧道围岩加固数值模拟方法
CN113177248B (zh) * 2021-04-21 2022-07-19 武汉大学 隧道围岩破裂碎胀大变形失稳灾变过程数值模拟方法
CN113204902B (zh) * 2021-04-21 2022-08-30 武汉大学 一种恒阻大变形锚杆隧道围岩加固数值模拟方法
CN116226982A (zh) * 2023-01-31 2023-06-06 武汉大学 一种粘性土-岩石隧道开挖耦合数值方法
CN116226982B (zh) * 2023-01-31 2024-05-28 武汉大学 一种粘性土-岩石隧道开挖耦合数值方法
CN116070790A (zh) * 2023-03-21 2023-05-05 中国建筑一局(集团)有限公司 一种减少施工现场重大安全风险隐患的预测方法和系统
CN116070790B (zh) * 2023-03-21 2023-09-01 中国建筑一局(集团)有限公司 一种减少施工现场重大安全风险隐患的预测方法和系统

Similar Documents

Publication Publication Date Title
CN112329290A (zh) 可用于施工过程模拟的有限元离散元多尺度耦合计算方法
CN109684693B (zh) 一种基于有限元分析预计加筋壁板后屈曲的方法
Samir et al. Damage detection in CFRP composite beams based on vibration analysis using proper orthogonal decomposition method with radial basis functions and cuckoo search algorithm
CN110555229B (zh) 一种无网格固体力学仿真方法、电子设备及存储介质
CN109033537B (zh) 堆石混凝土浇筑过程数值模拟的计算方法和系统
Tanaka et al. Study on crack propagation simulation of surface crack in welded joint structure
US10685154B2 (en) Analytical consistent sensitivities for nonlinear equilibriums, where the only source of nonlinearities is small sliding contact constraints
CN108984829B (zh) 堆石混凝土堆石体堆积过程的计算方法和系统
CN111125963B (zh) 基于拉格朗日积分点有限元的数值仿真系统及方法
Hautefeuille et al. A multi-scale approach to model localized failure with softening
CN114970260B (zh) 一种用于模拟复合材料破坏的格构相场方法
CN103942091A (zh) Matlab自定义模型和psasp联合仿真的励磁系统仿真方法及系统
CN109446731A (zh) 一种基于abaqus的岩土工程数值模拟方法
Yu et al. Incremental sequentially linear analysis to control failure for quasi-brittle materials and structures including non-proportional loading
CN111159943A (zh) 一种动翼面封严结构的屈曲处理方法
CN108491612A (zh) 为复合管液压胀形工艺提供选材方案的有限元模拟方法
Dang et al. Variable fractional modeling and vibration analysis of variable-thickness viscoelastic circular plate
CN117275633B (zh) 一种航空复合材料结构破坏过程分析方法及计算机设备
CN113011073B (zh) 基于深度学习的一维复杂滞回关系构建与结构模拟方法
Xie et al. A life prediction method of mechanical structures based on the phase field method and neural network
Strniša et al. A meshless solution of a small-strain plasticity problem
CN104730938A (zh) 一种电池特征模拟方法及系统
CN111159946B (zh) 一种基于最小势能原理的非连续性问题分区求解方法
Li et al. Structural dynamic reanalysis method for transonic aeroelastic analysis with global structural modifications
Mitra et al. Large displacement of crossbeam structure through energy method

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