留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于计算接触力学的粗颗粒土体材料细观性质模拟

潘洪武 王伟 张丙印

潘洪武, 王伟, 张丙印. 基于计算接触力学的粗颗粒土体材料细观性质模拟[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.08.0490
引用本文: 潘洪武, 王伟, 张丙印. 基于计算接触力学的粗颗粒土体材料细观性质模拟[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.08.0490
Hong-wu PAN, Wei WANG, Bing-yin ZHANG. SIMULATION ON MESO-MECHANICAL PROPERTY OF COARSE-GRAINED SOIL MATERIALS BASED ON COMPUTATIONAL CONTACT METHOD[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.08.0490
Citation: Hong-wu PAN, Wei WANG, Bing-yin ZHANG. SIMULATION ON MESO-MECHANICAL PROPERTY OF COARSE-GRAINED SOIL MATERIALS BASED ON COMPUTATIONAL CONTACT METHOD[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.08.0490

基于计算接触力学的粗颗粒土体材料细观性质模拟

doi: 10.6052/j.issn.1000-4750.2019.08.0490
基金项目: 国家重点研发计划课题项目(2017YFC0404802);国家自然科学基金项目(51979143)
详细信息
    作者简介:

    潘洪武(1993−),男,湖南人,博士生,主要从事岩土工程方面的研究(E-mail: panhw11@foxmail.com)

    王 伟(1989−),男,安徽人,博士生,主要从事岩土工程方面的研究(E-mail: wang1900w@163.com)

    通讯作者: 张丙印(1963−),男,河北人,教授,博士,博导,主要从事岩土工程方面的研究(E-mail: byzhang@tsinghua.edu.cn)
  • 中图分类号: TU43

SIMULATION ON MESO-MECHANICAL PROPERTY OF COARSE-GRAINED SOIL MATERIALS BASED ON COMPUTATIONAL CONTACT METHOD

  • 摘要: 采用离散元等数值模拟方法研究粗颗粒土体材料的细观力学行为,是近年来岩土力学的热点研究课题。提出了一种基于计算接触力学的粗颗粒土体材料细观力学行为模拟分析新方法。该法将土体颗粒剖分成一定数量的单元,通过计算接触力学方法模拟颗粒间的多体接触行为,是一种基于有限元变分原理的隐式数值求解方法,具有严格的理论基础。和离散元法方法相比,该方法在描述颗粒本身的力学特性方面具有优势,可以计算得到颗粒内部的应力分布情况,从而为颗粒破碎等细观行为的模拟计算提供依据。利用所发展的粗颗粒土体多体接触有限元计算程序,进行了不同围压的常规三轴数值试验,模拟计算结果符合粗颗粒土体材料三轴试验的一般规律,初步验证了所提出的计算方法的适用性。
  • 图  1  两个颗粒的接触问题示意图

    Figure  1.  Schematic diagram of two particle contact problem

    图  2  四球堆积算例

    Figure  2.  Four sphere accumulation example

    图  3  计算模型示意图

    Figure  3.  Schematic diagram of simulation model

    图  4  试验过程中的试样侧视图(σ3=600 kPa)

    Figure  4.  Side view of specimen during experiment (σ3=600 kPa)

    图  5  不同围压数值试验的计算结果

    Figure  5.  Computation results of numerical experiment at different confining pressures

    图  6  试样力链分布

    Figure  6.  Distribution of the force chain

    图  7  颗粒间法向接触力大小统计分布

    Figure  7.  Distribution of normal contact forces between particles

    图  8  颗粒单元大主应力统计分布

    Figure  8.  Distribution of large principle stress in a particle

    图  9  大主应力剖视图  /MPa

    Figure  9.  Section view of large principal stress

    表  1  接触力计算结果的对比

    Table  1.   Comparison between calculation and theoretical values

    接触力计算值/N理论值/N相对误差/(%)
    Fn113.43113.4340.022
    Ft14.5354.5320.065
    Fn253.72353.7230.000
    Ft24.5354.5410.129
    下载: 导出CSV

    表  2  两种模拟方法对比

    Table  2.   Comparison between two simulation methods

    基本原理求解方法颗粒模拟接触条件
    离散元方法牛顿第二
    定律
    显示求解一般为刚体
    块体
    通过罚函数法引入法向和切向弹簧
    计算接触力学方法有限元变分原理隐式求解有限元剖分的变形体通过约束变分引入KKT条件
    下载: 导出CSV
  • [1] Cundall P A. A computer model for simulating progressive large scale movements in blocky systems [C]// Proceedings of the Symposium of the International Society of Rock Mechanics, Nancy, France, International Society for Rock Mechanics, 1971, 1: 8 − 12.
    [2] Cundall P A, Hart R D. Numerical modelling of discontinua [J]. Engineering Computations, 1992, 9(2): 101 − 113. doi:  10.1108/eb023851
    [3] 王蕴嘉, 周梦佳, 宋二祥. 考虑颗粒破碎的堆石料湿化变形特性离散元模拟研究[J]. 工程力学, 2018, 35(增刊1): 217 − 222. doi:  10.6052/j.issn.1000-4750.2017.05.S043

    Wang Yunjia, Zhou Mengjia, Song Erxiang. DEM simulation of wetting deformation characteristics of rockfill considering particle breakage [J]. Engineering Mechanics, 2018, 35(增刊1): 217 − 222. (in Chinese) doi:  10.6052/j.issn.1000-4750.2017.05.S043
    [4] 杨智勇, 曹子君, 李典庆, 等. 颗粒接触摩擦系数空间变异性对颗粒流双轴数值试验的影响[J]. 工程力学, 2017, 34(5): 235 − 246.

    Yang Zhiyong, Cao Zijun, Li Dianqing, et al. Effect of spatially variable friction coefficient of granular materials on its macro-mechanical behaviors using biaxial compression numerical simulation [J]. Engineering Mechanics, 2017, 34(5): 235 − 246. (in Chinese)
    [5] 蒋明镜. 现代土力学研究的新视野−破折号之间断线宏微观土力学[J]. 岩土工程学报, 2019, 41(2): 195 − 254.

    Jiang Mingjing. New paradigm for modern soil mechanics: Geomechanics from micro to macro [J]. Chinese Journal of Geotechnical Engineering, 2019, 41(2): 195 − 254. (in Chinese)
    [6] 狄少丞, 王庆, 薛彦卓, 等. 破冰船冰区操纵性能离散元分析[J]. 工程力学, 2018, 35(11): 249 − 256. doi:  10.6052/j.issn.1000-4750.2017.09.0698

    Di Shaochen, Wang Qing, Xue Yanzhuo, et al. Manoeuvrability analysis of an icebreaker based on discrete element method [J]. Engineering Mechanics, 2018, 35(11): 249 − 256. (in Chinese) doi:  10.6052/j.issn.1000-4750.2017.09.0698
    [7] 严颖, 赵金凤, 季顺迎, 等. 块石含量和空间分布对土石混合体抗剪强度影响的离散元分析[J]. 工程力学, 2017, 34(6): 146 − 156.

    Yan Ying, Zhao Jinfeng, Ji Shunying, et al. Discrete element analysis of the influence of rock content and rock spatial distribution on shear strength of rock-soil mixtures [J]. Engineering Mechanics, 2017, 34(6): 146 − 156. (in Chinese)
    [8] Zhou W, Wu W, Ma G, et al. Undrained behavior of binary granular mixtures with different fines contents [J]. Powder Technology, 2018, 34: 139 − 153.
    [9] Munjiza A. The combined finite-discrete element method [M]. London: Wiley, 2004.
    [10] 马刚, 周创兵, 常晓林, 等. 岩石破坏全过程的连续-离散耦合分析方法[J]. 岩石力学与工程学报, 2011, 30(12): 2444 − 2455.

    Ma Gang, Zhou Chuangbin, Chang Xiaolin, et al. Continuous-discontinuous coupling analysis for whole failure process of rock [J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(12): 2444 − 2455. (in Chinese)
    [11] 王永亮, 鞠杨, 陈佳亮, 等. 自适应有限元-离散元算法、ELFEN软件及页岩体积压裂应用[J]. 工程力学, 2018, 35(9): 17 − 25. doi:  10.6052/j.issn.1000-4750.2017.06.0421

    Wang Yongliang, Ju Yang, Chen Jialiang, et al. Adaptive finite element-discrete element algorithm, software ELFEN and application in stimulated reservoir volume of shale [J]. Engineering Mechanics, 2018, 35(9): 17 − 25. (in Chinese) doi:  10.6052/j.issn.1000-4750.2017.06.0421
    [12] Moslem N, Hossein G. An efficient design tool based on FEM for evaluating effects of components properties and operating conditions on interaction of tire with rigid road [J]. Petroleum, 2015, 1(04): 388 − 396. doi:  10.1016/j.petlm.2015.10.012
    [13] Han C J, Yu C, Li Y, et al. Mechanical performance analysis of hollow cylindrical roller bearing of cone bit by FEM [J]. Journal of Yangtze River Scientific Research Institue, 2017, 34(3): 80 − 84.
    [14] Gregor H, Mitsunori T, Bojan D. Justification for a 2D versus 3D fingertip finite element model during static contact simulations [J]. Computer Methods in Biomechanics and Biomedical Engineering, 2016, 19(13): 1409 − 1417. doi:  10.1080/10255842.2016.1146712
    [15] Alexander S, Philipp F, Johannes K, et al. Isogeometric dual mortar methods for computational contact mechanics [J]. Computer Methods in Applied Mechanics and Engineering, 2016, 301: 259 − 280. doi:  10.1016/j.cma.2015.12.018
    [16] Zhou M Z, Zhang B Y, Peng C. Numerical evaluation of soft-slab joint in concrete-faced rockfill dam with dual mortar element method [J]. International Journal for Numerical and Analytical methods in Geomechanics, 2018, 42(5): 781 − 805. doi:  10.1002/nag.2768
    [17] Wang W, Zhou M Z, Zhang B Y, et al. A dual mortar contact method for porous media and its application to clay-core rockfill dams [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2019, 43(9): 1744 − 1769. doi:  10.1002/nag.2930
    [18] 王伟, 殷殷, 潘洪武, 等. 基于接触力学的堰塞坝防渗墙应力变形分析[J]. 水力发电学报, 2018, 37(12): 94 − 101. doi:  10.11660/slfdxb.20181210

    Wang Wei, Yin Yin, Pan Hongwu, et al. Stress-deformation analysis of diaphragm wall in landslide dam based on contact mechanics [J]. Journal of Hydroelectric Engineering, 2018, 37(12): 94 − 101. (in Chinese) doi:  10.11660/slfdxb.20181210
    [19] Tur M, Fuenmayor F J, Wriggers P. A mortar-based frictional contact formulation for large deformations using Lagrange multipliers [J]. Computer Methods in Applied Mechanics and Engineering, 2009, 198(37-40): 2860 − 2873. doi:  10.1016/j.cma.2009.04.007
    [20] Wohlmuth B I. A mortar finite element method using dual spaces for the Lagrange multiplier [J]. SIAM Journal on Numerical Analysis, 2000, 38(3): 989 − 1012.
    [21] Popp A, Gee M W, Wall W A. A finite deformation mortar contact formulation using a primal-dual active set strategy [J]. International Journal for Numerical Methods in Engineering, 2009, 79(11): 1354 − 1391. doi:  10.1002/nme.2614
    [22] 刘凯欣, 高凌天. 离散元法研究的评述[J]. 力学进展, 2003, 33(4): 483 − 490. doi:  10.3321/j.issn:1000-0992.2003.04.005

    Liu Kaixin, Gao Lingtian. A review on the discrete element method [J]. Advances in Mechanics, 2003, 33(4): 483 − 490. (in Chinese) doi:  10.3321/j.issn:1000-0992.2003.04.005
  • [1] 陶慕轩, 赵继之.  采用通用有限元程序的弥散裂缝模型和分层壳单元模拟钢筋混凝土构件裂缝宽度 . 工程力学, doi: 10.6052/j.issn.1000-4750.2019.07.0342
    [2] 栗志杰, 由小川, 柳占立, 庄茁, 杨策.  基于三维头部数值模型的颅脑碰撞损伤机理研究 . 工程力学, doi: 10.6052/j.issn.1000-4750.2018.04.0254
    [3] 王文杰, 张磊, 林峰.  预应力钢丝缠绕“正交预紧机架”刚度比计算方法研究 . 工程力学, doi: 10.6052/j.issn.1000-4750.2014.07.0647
    [4] 井国庆, 王子杰, 施晓毅.  多围压下三轴压缩试验与不可破裂颗粒离散元法分析 . 工程力学, doi: 10.6052/j.issn.1000-4750.2014.03.0205
    [5] 王明强, 杨军.  HISS模型在ABAQUS中的三维实现、验证和应用 . 工程力学, doi: 10.6052/j.issn.1000-4750.2012.09.0666
    [6] 王大奎, 张 军, 王春艳, 郭 璇.  车轮与高速道岔长短心轨的接触分析 . 工程力学, doi: 10.6052/j.issn.1000-4750.2011.07.0471
    [7] 赵秋, 吴冲.  U肋加劲板焊接残余应力的一种简化计算方法 . 工程力学, doi: 10.6052/j.issn.1000-4750.2010.12.0936
    [8] 周素霞, 肖 楠, 谢基龙.  多轴载荷下车轮辐板裂纹扩展特性研究 . 工程力学,
    [9] 张建国, 王 芳, 薛 强.  后碰撞中人体颈部动力学响应的有限元分析 . 工程力学,
    [10] 杨允表, 吕忠达.  大跨度斜拉桥索塔锚固区钢-混凝土结构竖向受力机理的有限元法 . 工程力学,
    [11] 魏红卫, 喻泽红, 尹华伟.  土工合成材料加筋粘性土的三轴试验研究 . 工程力学,
    [12] 童育强, 向天宇, 赵人达.  基于退化理论的空间梁单元有限元分析 . 工程力学,
    [13] 王琥, 李光耀, 钟志华.  有限元并行计算中网格自动分区的优化 . 工程力学,
    [14] 邱文亮, 姜萌, 张哲.  钢-混凝土组合梁收缩徐变分析的有限元方法 . 工程力学,
    [15] 田玉基, 蓝宗建, 杨庆山.  巨型框架多功能减振结构的有限元计算方法 . 工程力学,
    [16] 郭然, 施惠基, 郁向东.  多层圆柱结构接触预紧力失效分析 . 工程力学,
    [17] 徐刚, 任文敏, 郑兆昌.  水下运动体的三维动力特性分析 . 工程力学,
    [18] 王爱俊, 乔新, 厉蕾, 陈宇宏.  层合透明材料抗冲击有限元数值模拟 . 工程力学,
    [19] 王成国, 梁国平, 刘金朝.  多柔性体系统动力学的有限元方法(MUFEM) . 工程力学,
    [20] 宋二祥.  软化材料有限元分析的一种非局部方法 . 工程力学,
  • 加载中
计量
  • 文章访问数:  49
  • HTML全文浏览量:  34
  • PDF下载量:  1
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-08-28
  • 修回日期:  2020-01-17
  • 网络出版日期:  2020-06-02

基于计算接触力学的粗颗粒土体材料细观性质模拟

doi: 10.6052/j.issn.1000-4750.2019.08.0490
    基金项目:  国家重点研发计划课题项目(2017YFC0404802);国家自然科学基金项目(51979143)
    作者简介:

    潘洪武(1993−),男,湖南人,博士生,主要从事岩土工程方面的研究(E-mail: panhw11@foxmail.com)

    王 伟(1989−),男,安徽人,博士生,主要从事岩土工程方面的研究(E-mail: wang1900w@163.com)

    通讯作者: 张丙印(1963−),男,河北人,教授,博士,博导,主要从事岩土工程方面的研究(E-mail: byzhang@tsinghua.edu.cn)
  • 中图分类号: TU43

摘要: 采用离散元等数值模拟方法研究粗颗粒土体材料的细观力学行为,是近年来岩土力学的热点研究课题。提出了一种基于计算接触力学的粗颗粒土体材料细观力学行为模拟分析新方法。该法将土体颗粒剖分成一定数量的单元,通过计算接触力学方法模拟颗粒间的多体接触行为,是一种基于有限元变分原理的隐式数值求解方法,具有严格的理论基础。和离散元法方法相比,该方法在描述颗粒本身的力学特性方面具有优势,可以计算得到颗粒内部的应力分布情况,从而为颗粒破碎等细观行为的模拟计算提供依据。利用所发展的粗颗粒土体多体接触有限元计算程序,进行了不同围压的常规三轴数值试验,模拟计算结果符合粗颗粒土体材料三轴试验的一般规律,初步验证了所提出的计算方法的适用性。

English Abstract

潘洪武, 王伟, 张丙印. 基于计算接触力学的粗颗粒土体材料细观性质模拟[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.08.0490
引用本文: 潘洪武, 王伟, 张丙印. 基于计算接触力学的粗颗粒土体材料细观性质模拟[J]. 工程力学. doi: 10.6052/j.issn.1000-4750.2019.08.0490
Hong-wu PAN, Wei WANG, Bing-yin ZHANG. SIMULATION ON MESO-MECHANICAL PROPERTY OF COARSE-GRAINED SOIL MATERIALS BASED ON COMPUTATIONAL CONTACT METHOD[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.08.0490
Citation: Hong-wu PAN, Wei WANG, Bing-yin ZHANG. SIMULATION ON MESO-MECHANICAL PROPERTY OF COARSE-GRAINED SOIL MATERIALS BASED ON COMPUTATIONAL CONTACT METHOD[J]. Engineering Mechanics. doi: 10.6052/j.issn.1000-4750.2019.08.0490
  • 颗粒物质是离散和相互作用的大量固体颗粒所组成的复杂体系。颗粒物质广泛存在于自然界,与人类的日常生活和生产密切相关。但是,颗粒体系的很多静力学和动力学行为尚不能用一般的连续介质力学理论很好地解释。Science杂志于2005年把颗粒物质与湍流并列为100个科学难题之一。

    土体由离散的颗粒组成。颗粒间以强耗散的接触摩擦为主,其宏观力学特性主要取决于细观层次的颗粒联结和排列等。近年来,人们逐渐开始重视土颗粒体系细观力学行为的研究。由于在物理试验中存在对颗粒细观行为测量技术的难题,因而各种数值模拟方法逐步成为主要的研究手段。对此,学者们已经发展出了离散元法、非连续变形方法和连续-离散耦合分析方法等多种数值计算方法。

    离散元法是Cundall[1]于1971年提出的一种基于牛顿第二运动定律的非连续变形数值方法。该法将土体看成是离散单元的集合,通过给定颗粒间的接触关系[2],来反映颗粒材料的宏细观力学响应。离散元法具有很强的描述颗粒间不连续变形的能力,受到人们广泛的青睐和重视,目前是进行颗粒材料细观力学特性分析的最主要的方法,取得了丰硕的研究成果[3-8]

    由于计算方法的限制,离散元法对颗粒本身细观力学特性的描述(如颗粒内部应力分布等)相对较为粗糙,其细观数值模拟难以定量地再现颗粒材料的宏细观力学特性。为此,Munjiza、马刚、王永亮等[9-11]发展了连续-离散耦合分析方法(FDEM)。

    有限元方法是目前应用最广泛的数值模拟计算方法。但一般认为,有限元法主要适合于连续介质的模拟计算分析,其在不连续变形的描述方面存在缺陷。近年来,基于有限元方法的计算接触力学算法[12-15]得到了迅速发展,使得有限元方法在模拟不连续变形和多体接触问题方面的能力大幅提高。

    土体是由离散的土颗粒组成的,其细观力学性质的研究为一个典型的多体接触问题。从原理上讲,计算接触力学方法是适用于进行这种多体接触分析的。基于上述的思路,本文开展了基于计算接触力学的土颗粒材料细观力学特性模拟的探索性研究工作。通过对粗颗粒土体材料进行三轴试验的数值计算分析,探讨了基于有限元方法的计算接触力学算法对粗颗粒土体材料多体接触问题的适用性。

    • 计算接触力学方法来自固体力学领域,是非线性数值分析中最具挑战性的课题之一,提出以来得到了来自计算数学和计算力学等领域诸多学者的研究。计算接触力学研究内容繁多,涉及范围很广。主要研究的问题包括接触界面离散形式(点对点、点对面和面对面等)、接触非线性求解、接触本构方程和接触约束施加方法等。计算接触力学方法被广泛应用在机械、航空和生物力学等领域。在岩土工程领域近年来才开始得到应用[16-18]

    • 直观地讲,接触问题反映的是不同物体之间的相互作用。对于本文考虑的粗颗粒土体材料,不失一般性,考虑如图1所示的可变形颗粒${\Omega ^1}$${\Omega ^2}$的接触问题。将颗粒${\Omega ^i}$的边界记作$\partial {\Omega ^i}(i{\rm{ = }}1,2)$,每个边界可分为Dirichlet边界$\varGamma _{\rm{D}}^i$、Neumann边界$\varGamma _{\rm{N}}^i$和接触边界$\varGamma _{\rm{C}}^i$,三者组成整体边界且相互无重叠:

      $$\varGamma _{\rm{D}}^i \cup \varGamma _{\rm{N}}^i \cup \varGamma _{\rm{C}}^i = \partial {\Omega ^i}\quad\quad\quad\quad\;\;$$ (1)
      $$\varGamma _{\rm{D}}^i \cap \varGamma _{\rm{N}}^i{\rm{ = }}\varGamma _{\rm{D}}^i \cap \varGamma _{\rm{C}}^i{\rm{ = }}\varGamma _{\rm{N}}^i \cap \varGamma _{\rm{C}}^i = \emptyset$$ (2)

      图  1  两个颗粒的接触问题示意图

      Figure 1.  Schematic diagram of two particle contact problem

      接触边界上,两物体被记作从面和主面。将物体表面的法向和切向单位向量分别记为${n}$${\\text{τ}}_1$${{\\text{τ}}_2}$,接触区域的法向间隙${g_{\rm{n}}}$可表示为:

      $${g_{\rm{n}}}({{x}^{\rm{A}}}) = - {n} \cdot ({{x}^{\rm{A}}} - {{\hat x}^{\rm{B}}})$$ (3)

      式中:${{x}^{\rm{A}}}$为从面点的坐标;${{\hat x}^{\rm{B}}}$为点${{x}^{\rm{A}}}$$\varGamma _{\rm{C}}^{\rm{B}}$上的投影点。接触区域的物体接触力${{t}_{\rm{C}}}$可记作:

      $${{t}_{\rm{C}}} = {p_{\rm{n}}}{n} + {t_{\tau 1}}{{\\text{τ}}_1} + {t_{\tau 2}}{{\\text{τ}}_2}$$ (4)

      式中:${p_{\rm{n}}}$为接触力法向分量;${t_{\tau 1}}$${t_{\tau 2}}$为接触力的切向分量。

      在计算接触力学中,一般将上述接触条件表示为如下的互补KKT条件。

      法向:${g_{\rm{n}}} \geqslant 0,\;{p_{\rm{n}}} \geqslant 0,\;{g_{\rm{n}}}{p_{\rm{n}}} = 0$

      切向:$\phi : = \left\| {{{{t}}_{\tau}}} \right\| - \mu \left| {{p_{\rm{n}}}} \right| {\leqslant} 0$

      $$ {{v}_{{\tau} ,\;{\rm{rel}}}} + \beta {{t}_{\tau} } = {0},\;\beta \geqslant 0,\;\phi \beta = 0 $$ (5)

      式中:$\mu $为摩擦系数;${{v}_{\tau ,{\rm{rel}}}}$为切向相对速度;${p_{\rm{n}}}$为法向压力;${{t}_{\tau} } = {t_{{\tau} 1}}{{\\text{τ}}_1} + {t_{{\tau} 2}}{{\\text{τ}}_2}$为切向接触力;$\phi $$\beta $为中间参数。需要说明的是,式(5)中切向采用的是摩擦定律。实际上,也完全可以采用其他更加复杂的接触本构方程(如包括粘结等)。

    • 考虑两物体间的接触问题,将接触面边界视为面力边界,则该问题的虚功原理可表示为:

      $$\delta W{\rm{ = }}\delta {W_{{\rm{int,ext}}}} + \delta {W_{\rm{C}}} = 0$$ (6)

      式中:$\delta {W_{{\rm{int,ext}}}}$为不考虑接触时物理的内、外力虚功之和;$\delta {W_{\rm{C}}}$为接触力虚功。

      在接触边界上,由牛顿第三定律可知,物体间接触力满足${t}_{\rm{C}}^{\rm{s}} = - {t}_{\rm{C}}^{\rm{m}}$,则接触力虚功可记作:

      $$\delta {W_{\rm{C}}} = - \int_{\varGamma _{\rm{C}}^i} {{t}_{\rm{C}}^{\rm{s}}} (\delta {{u}^{\rm{s}}} - \delta {{u}^{\rm{m}}}){\rm{d}}\Gamma $$ (7)

      在接触从面上引入拉格朗日乘子:

      $$ {\\text{λ}} = {{\lambda }_n}{{n}} + {{\lambda }_{\tau 1}}{{{\\text{τ}}}_1} + {{\lambda }_{\tau 2}}{{{\\text{τ}}}_2} $$ (8)

      式中,${\lambda _n}$为乘子法向分量,${\lambda _{\tau 1}}$${\lambda _{\tau 2}}$为乘子的切向分量。考虑到拉格朗日乘子和接触力等值法向。接触虚功可以记作:

      $$ \delta {W_{\rm{C}}} = \int_{\varGamma _{\rm{C}}^i} {\\text{λ}} \left(\delta {{u}^{\rm{s}}} - \delta {{u}^{\rm{m}}}\right){\rm{d}}\Gamma $$ (9)

      式(5)中的法向不可贯入条件和切向摩擦条件的弱形式可表示为:

      $$\int_{\Gamma _{\rm{C}}^{\rm{S}}} {\delta {\lambda _n}g({{x}^{\rm{s}}}){\rm{d}}} \Gamma \geqslant 0\quad\quad\quad$$ (10)
      $$\int_{\Gamma _{\rm{C}}^{\rm{S}}} {\delta {{\\text{λ}} _{\tau} } \cdot [{{v}_{{\tau} ,{\rm{rel}}}} + \beta {{t}_{\tau} }]{\rm{d}}} \Gamma = 0$$ (11)

      式中,${{\\text{λ}}_{\tau} }={\lambda _{\tau 1} }{\tau}_1+{\lambda _{\tau 2} }{\tau}_2$为乘子${\\text{λ}}$的切向分量。式(9)~式(11)表示了两体接触问题的混合变分格式。

    • 为了在接触界面上计算接触虚功,需要对接触界面进行离散,常用的界面离散方法有点点方法、点面方法和面面方法等,其中,面面方法具有计算精度高、能通过接触分片试验等优点,是国内外学者们目前研究的热点,也是本文数值研究采用的方法。

      在接触界面上,引入拉格朗日乘子,通过节点量插值记为:

      $${{\\text{λ}}^{\rm{s}}} = \sum\limits_{i = 1}^{{n_{\rm{s}}}} {{\Psi _i}} {{\\text{λ}}_i}$$ (12)

      式中:${n_{\rm{s}}}$为主面和从面节点数目;${{\\text{λ}}^{\rm{s}}}$为分布于从面上的拉格朗日乘子;${\Psi _i}$为Lagrange乘子的形函数。基于上述空间离散格式,接触力对应的虚功项可采用矩阵形式表示为:

      $$\delta {W_{\rm{c}}} = \delta {u}_{\rm{s}}^{\rm{T}}{D^{\rm{T}}}{\\text{λ}} - \delta {u}_{\rm{m}}^{\rm{T}}{M^{\rm{T}}}{\\text{λ}}$$ (13)

      式中:${\\text{λ}}{\rm{ = }}\left\{ {{\lambda _1},{\lambda _1}, \cdots ,{\lambda _n}} \right\}$为从面节点拉格朗日乘子向量;矩阵DM为接触矩阵,表示从面和主点间相关物理量的关联系数,其3×3子块分别为:

      $${D_{jk}} = \int_{\Gamma _{\rm{C}}^{\rm{S}}} {{\Psi _j}{N_k}{\rm{d}}} \Gamma {{{I}}_3}$$ (14)
      $${M_{jl}} = \int_{\Gamma _{\rm{C}}^{\rm{S}}} {{\Psi _j}{N_l}{\rm{d}}} \Gamma {{{I}}_3}$$ (15)

      式中,${{{I}}_3}$为3×3单位矩阵。

      在计算接触力学方法中,拉格朗日乘子可取不同分布形式。在Mortar元方法[19]中,乘子形函数和有限元单元形函数相同,即${\Psi _i} = {N_i}$。在Dual mortar元方法[20]中,乘子分布服从对偶形函数,即乘子形函数满足所谓双正交条件:

      $$\int_{\Gamma _{\rm{C}}^{\rm{S}}} {{\Psi _j}{N_l}{\rm{d}}} \Gamma {\rm{ = }}{\delta _{jl}}\int_{\Gamma _{\rm{C}}^{\rm{S}}} {{N_l}{\rm{d}}} \Gamma ,j,l = 1,2, \cdots ,{n_{\rm{s}}}$$ (16)

      式中,${\delta _{jl}}$为克罗内克函数。此时,接触矩阵D为对角阵,可实现计算格式中乘子的高效凝聚。

    • 由于接触区域和接触状态的双重不确定性,接触问题具有极强的非线性。计算接触力学中常采用PDASS (primal dual active-set strategy)方法[21]对接触问题进行非线性迭代求解。其基本思路为:首先将接触界面上的从点划分为脱开、滑移和粘结等不同状态的点集,并分别引入相应的接触边界条件;然后再根据相应互补的接触条件进行校核,判别从点的接触状态是否正确;对接触状态不正确的从点重新划分其接触状态。如此迭代直至接触界面从点的状态不再发生改变为止。

    • 本文在界面离散上采用Dual mortar元方法,其得到的D矩阵为对角矩阵,可以高效地进行求逆操作。在激活节点上引入相对位移:

      $$\bar d_j^{\rm{s}} = R_{jj}^{\rm{T}}\left( {d_j^{\rm{s}} - \sum\limits_{l \in {\rm{M}}} {D_{jj}^{ - 1}} {M_{jl}}{d_l}} \right),\forall j \in {\rm{A}}$$ (17)

      式中:$\bar d_j^{\rm{s}} = \{ \bar d_j^{{\tau _1}},\bar d_j^{{\tau _2}},\bar d_j^n\} $为局部坐标系中从面和主面的相对位移;R为旋转坐标。

      $${R_j}_k = {\delta _j}_k{R_j}_j,{R_j}_j = \{ {({{\\text{τ}}_1})_j},{({{\\text{τ}}_2})_j},{{n}_j}\} $$ (18)

      有限元计算格式可表示为:

      $$\tilde {{K}}\Delta \tilde {{d}} = \tilde {{r}}$$ (19)

      式中:

      $$\tilde K = {{{P}}^{\rm{T}}}K{{P}}\quad\quad\quad\quad\quad\quad\quad$$ (20)
      $${{P}} = \left[ {\begin{array}{*{20}{c}} {{I_{\rm{N}}}}&0&0&0 \\ 0&{{I_{\rm{M}}}}&0&0 \\ 0&0&{{I_{\rm{I}}}}&0 \\ 0&0&{{D^{ - 1}}M}&R \end{array}} \right]\;\;$$ (21)
      $$\tilde {{r}} = \left[ {\begin{array}{*{20}{l}} {F_{\rm{N}}^{{\rm{ext}},{\rm{int}}}}\\ {F_{\rm{M}}^{{\rm{ext}},{\rm{int}}}{\rm{ + }}({D^{{\rm{ - }}1}}M)F_{\rm{A}}^{{\rm{ext}},{\rm{int}}}}\\ {F_{\rm{I}}^{{\rm{ext}},{\rm{int}}}}\\ {{R^{\rm{T}}}F_{\rm{A}}^{{\rm{ext}},{\rm{int}}} - D{\\text{λ}}} \end{array}} \right]$$ (22)

      将有限元网格节点划为不同点集[21],其约束条件分别为:

      $${\lambda _j} = 0,\forall j \in I\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$$ (23)
      $${({\lambda _{\tau i}})_j} = - \mu \left| {{\lambda _n}} \right|{({{\\text{τ}}_i})_j} \cdot \frac{{{\bar v}}}{{\left\| {{\bar v}} \right\|}},i = 1,2,\forall j \in {\rm{L}}$$ (24)
      $$\Delta \bar d_j^{\rm{s}} = (0,0, - {(\bar d_{\rm{n}}^{})_j}),\forall j \in T\;\;\;\quad\quad\quad\quad\quad$$ (25)
      $$\Delta \bar d_{\rm{n}}^{\rm{s}} = - {(\bar d_{\rm{n}}^{})_j},\forall j \in L\quad\quad\quad\quad\quad\quad\quad\quad$$ (26)
    • 基于课题组发展的多体接触有限元计算程序Contana,针对颗粒集合体的特点,发展了相应的土体颗粒多体接触有限元计算程序。

      采用四球堆积算例对所发展的有限元程序进行验证。如图3所示,四个球体的半径R相等。球体B2B3B4位于平面W1上,其球心连线为边长等于l=0.105 m的等边三角形。球体B1对称放置在三个球体的顶部,四个球体的球心连线形成一个四面体。球体半径R=0.05 m,密度ρ=7800 kg/m3,弹性模型E=100 GPa,泊松比ν=0.3。本算例中,采用六面体单元对球体进行建模,每个球体含有1734个节点、1152个单元。

      当球体间以及球体和底面间的摩擦系数μ1μ2足够大时,四个球体处于稳定状态。记球体之间以及球体和底面之间接触力的法向和切向分量分别为Fn1Ft1Fn2Ft2。根据静力平衡条件,可以求得不考虑球体变形条件下各接触力的大小。由于在算例情况下球体变形非常小,下文称该解答为理论解。

      图  2  四球堆积算例

      Figure 2.  Four sphere accumulation example

      表1给出了有限元接触力计算结果与理论值的对比情况。其中,理论解扣除了有限元网格划分所引起的质量误差。从表1中可以看出,接触力的计算值和相应理论解的最大相对误差为0.129%,数值解和理论解吻合较好。

      表 1  接触力计算结果的对比

      Table 1.  Comparison between calculation and theoretical values

      接触力计算值/N理论值/N相对误差/(%)
      Fn113.43113.4340.022
      Ft14.5354.5320.065
      Fn253.72353.7230.000
      Ft24.5354.5410.129
    • 利用所发展的土体颗粒多体接触有限元计算程序,进行了不同围压的常规三轴数值试验,探讨了计算接触力学方法对土体颗粒系统细观力学特性分析的适用性。

    • 三轴试样由粒径分别为21 mm、19 mm、17 mm、15 mm和13 mm的460个球颗粒组成。上述不同粒径的颗粒数目分别为70、95、130、95和70。

      图3所示,计算采用10 cm×10 cm×20 cm的柱体试样。试样的四个侧面和两个顶面共设置6块加载板。加载板使用六面体单元离散。在加载板之间不设任何的接触关系,它们在计算中不会发生相互作用,可独立施加所需的应力或位移。每个球颗粒均被剖分为510个四面体单元、133个节点。整个计算模型共包含有78748个节点、245564个单元。

      图  3  计算模型示意图

      Figure 3.  Schematic diagram of simulation model

      计算中,颗粒和加载板均采用线弹性模型模拟,弹性模量分别取E=1 GPa和25 GPa,泊松比均取ν=0.3。不考虑重力影响。颗粒-颗粒之间以及颗粒-加载板之间均采用库仑摩擦定律模拟,摩擦系数分别取0.4和0.2。

    • 为了能与土体三轴试验的结果进行定性对比,参照土工三轴物理试验的流程,将本文数值试验的计算流程也划分为制样、固结和加载三个过程。数值试验的试样参照目前离散元方法所采用的三轴试样生成方法进行生成[5]

      采用生成的同一个数值试验试样,分别进行了围压σ3=200 kPa、400 kPa和600 kPa的常规三轴数值压缩试验。计算时,首先模拟施加围压的固结过程。采用Δσ3=50 kPa分级施加围压力σ3

      对每一个围压的试验,固结完成后,保持侧向围压力σ3大小不变,采用应变控制的加载方式进行剪切加载。具体方式为,每个加载步使加载顶板向下位移0.01 cm,对应试样发生0.05%的轴向应变。

      计算在个人计算机在上进行。计算机CPU的型号为Intel® CoreTM i7-2600K(3.40 GHz),每个围压试验的计算耗时约为40 h。

    • 在本文的计算结果中,应力和应变均以压为正。如前所述,采用计算接触力学的方法,可以对颗粒本身划分有限元网格。因此,所得数值试验的结果不仅能反映试样的整体宏观力学特性,并且也能反映颗粒本身的细部力学性质。

    • 图4给出了σ3=600 kPa数值试验过程中试样的侧视图。可以看出,在试验过程中,颗粒本身不发生大的变形,试样的变形主要由颗粒在空间相对位置的变化所致,即主要由颗粒间的滑移和颗粒转动等所导致的颗粒集合体排列方式的改变所产生。

      图  4  试验过程中的试样侧视图(σ3=600 kPa)

      Figure 4.  Side view of specimen during experiment (σ3=600 kPa)

      图5给出了不同围压的应力-应变曲线和体变计算曲线。由图5(a)所示的(σ1σ3)−ε1曲线可以看到,各围压的应力-应变曲线均表现出一定的应变软化特征。偏差应力随轴向应变的增加逐步增大,在开始阶段增加相对较快,之后逐步变得较为平缓,并达到峰值。峰值后随轴向应变增加,偏差应力会发生一定的下降。对不同围压的试验,对应相同的轴向应变,围压大对应的偏差应力也越高,表明试样的刚度或强度随围压逐步增大,表现出了颗粒材料压硬性的特点。这些特点均符合土体常规三轴压缩试验应力-应变曲线的一般特征。

      图  5  不同围压数值试验的计算结果

      Figure 5.  Computation results of numerical experiment at different confining pressures

      图5(b)所示的εv~ε1曲线可以看出,所得体变曲线总体符合粗粒土三轴试验中的体变特性。对一个围压的数值试验,试样在初始先发生体缩。当轴向应变较大时,体积变形逐步转变为剪胀。对不同围压的数值试验,围压较小时,试验初始段的体缩较小,后段的剪胀较为显著。随着围压的提高,试样的剪缩性总体增大,剪胀性逐步减小。根据试验结果绘制莫尔圆可得试样内摩擦角约为44°。

    • 图6给出围压σ3=600 kPa试验试样内部的力链分布情况。由图6可见,在试验过程中,水平方向接触力小于竖直方向接触力,荷载主要沿竖直方向传递。对比不同轴向应变下试样力链分布可知,颗粒间接触力随着轴向应变增加而增大。

      图  6  试样力链分布

      Figure 6.  Distribution of the force chain

      图7进一步给出了围压σ3=600 kPa试验中颗粒间法向接触力大小的分布情况。由图7可以看出,当试样固结完成尚未进行剪切时,颗粒间法向接触力总体相对较小,绝大多数分布在小于4 kN以下的范围。开始加载剪切后,颗粒间法向接触力逐步增大。随轴向应变的增大,小法向力接触对数目逐渐降低,大法向力接触对数目逐渐增加。当轴向应变达到14%时,一些颗粒的接触应力可达16 kN以上。

      如前所述,采用计算接触力学方法可以对颗粒本身划分所需的有限元的单元,从而可以计算得到其详细的细观应力和变形状态。这是计算接触力学方法相比离散元方法的一个优势。在本次计算中,每个球颗粒均包含133个节点、510个单元,因而可以得到任意时刻所有颗粒详细的应力状态。

      本文的计算模型中,颗粒采用四面体常应变单元进行网格剖分,且单元尺寸相近,因此可根据颗粒内部单元应力的统计分布来评估颗粒的总体应力状态。图8给出了σ3=600 kPa围压数值试验所得某个颗粒中单元大主应力大小的分布情况。当施加围压力之后但尚未进行剪切之前,试样处于等向压缩应力状态。此时,假如试样是连续的均质材料,则试样中所有单元的三个主应力都应为600 kPa。但由于试样由土体颗粒组成,因而计算所得单元应力的大小并不相等。其中,少数单元处于受拉状态,计算所得大主应力的最大值约为1.5 MPa。

      图  7  颗粒间法向接触力大小统计分布

      Figure 7.  Distribution of normal contact forces between particles

      图  8  颗粒单元大主应力统计分布

      Figure 8.  Distribution of large principle stress in a particle

      当进行剪切后,由图8可见,计算所得颗粒单元的大主应力总体增大,承受高应力的单元数目显著增加。少数单元仍会处于受拉状态。当轴向应变达到14%时,一些颗粒单元计算所得大主应力的最大值增大约为3.0 MPa。

      图9σ3=600 kPa围压的三轴数值试验,当轴向应变ε1=10%时,某典型颗粒中的大主应力分布的计算结果。由图9可以看出,计算所得单元的大主应力在颗粒中的分布非常的不均匀。其中,在颗粒中的大部分区域,计算所得到的应力值相对均较低。但在颗粒接触点附近的局部小区域,则会计算得到相对很高的应力分布。对于颗粒材料,其所受到的作用力是通过颗粒接触点进行传递的。对于本文所采用的球形颗粒,在接触点部位会作用集中的粒间力的作用。计算所得上述的应力分布特征符合赫兹接触问题的一般规律。

      图  9  大主应力剖视图  /MPa

      Figure 9.  Section view of large principal stress

    • 由本文前面的分析可知,计算接触力学方法是基于有限元方法的一种进行多体接触问题的数值求解方法。与离散元方法相比,在基本原理、求解方法、颗粒本身的模拟以及接触条件的模拟等方面都存在显著的不同。表2总结列出了两种方法在模拟颗粒材料特性方面的主要差别。

      表 2  两种模拟方法对比

      Table 2.  Comparison between two simulation methods

      基本原理求解方法颗粒模拟接触条件
      离散元方法牛顿第二
      定律
      显示求解一般为刚体
      块体
      通过罚函数法引入法向和切向弹簧
      计算接触力学方法有限元变分原理隐式求解有限元剖分的变形体通过约束变分引入KKT条件

      总体而言,将计算接触力学方法引入土体颗粒集合体细观应力变形特性的模拟分析,可能会具有如下的优势:

      1)在基本原理和求解方面,计算接触力学方法是一种基于变分原理的隐式数值求解方法,具有严格的理论基础。离散元方法的基本原理是牛顿第二定律,根据颗粒受到的合力计算颗粒在各时间步的位移,是一种显示求解方法。自提出至今,离散元及其相关方法的研究已有近50年历史,众多学者发表了大量学术论文。但是,离散元方法缺乏理论的严密性。刘凯欣和高凌天[22]曾评述“研究离散元的理论和算法的文章却很少”“自诞生的那天起就带有缺乏理论严密性的先天不足”。

      2)在颗粒本身的模拟方面,采用计算接触力学方法可以对颗粒本身划分所需的有限单元,从而可以得到其详细的细观应力和变形状态,可为颗粒破碎等的模拟计算提供基础。传统离散元方法的模型一般由刚性块体组成,此时无法考虑颗粒体本身变形,也无法计算得到颗粒体内部应力的分布。

      3)在颗粒间的接触特性模拟方面,计算接触力学方法通过拉格朗日乘子法引入KKT接触条件,来模拟颗粒接触界面的接触特性,是一种高精度的接触问题求解方法。拉格朗日乘子的物理意义为接触应力。该方法无须引入接触弹簧,可避免罚模量选取对人工经验的依赖等。离散元方法采用罚函数法,通过在接触点引入法向和切向弹簧模拟颗粒间的接触行为,存在弹簧刚度以及罚模量选取对人工经验的依赖等。

      但是,大规模的土体颗粒集合体是一个大型强非连续颗粒系统,将海量的颗粒间接触条件引入有限元方程后,会造成系统刚度矩阵的高度病态。为此,研发针对大型复杂颗粒集合体的高效接触搜索算法、高度病态矩阵预处理技术和具有较强鲁棒性的迭代算法等,是该种方法能否成功应用于土体颗粒集合体细观力学行为模拟分析的关键。这也是该法亟待在未来研究中解决的关键技术难题。

    • 本文开展了基于计算接触力学的粗颗粒土体材料细观力学行为模拟分析的探索性研究工作,主要取得如下的结论:

      (1)提出了一种基于计算接触力学的粗颗粒土体材料细观力学行为模拟分析新方法。该法将颗粒剖分成一定数量的有限元单元,通过计算接触力学方法模拟颗粒间的接触行为。

      (2)和离散元法方法相比,本文提出的模拟计算方法是一种基于变分原理的隐式数值求解方法,具有严格的理论基础。此外,该方法在描述颗粒本身的力学特性方面具有优势,既可以计算分析粗颗粒土体材料的宏细观力学特性,也可以计算得到颗粒内部的应力分布情况,可为颗粒破碎等的细观行为的模拟计算提供依据。

      (3)利用所发展的土体颗粒多体接触有限元计算程序,进行了不同围压的常规三轴数值试验,模拟计算结果符合土体三轴试验的一般规律,初步验证了所提出的计算方法的适用性。

目录

    /

    返回文章
    返回