CN111611695A - 一种模拟岩土材料时离散元线性刚度参数的自动标定方法 - Google Patents

一种模拟岩土材料时离散元线性刚度参数的自动标定方法 Download PDF

Info

Publication number
CN111611695A
CN111611695A CN202010393302.5A CN202010393302A CN111611695A CN 111611695 A CN111611695 A CN 111611695A CN 202010393302 A CN202010393302 A CN 202010393302A CN 111611695 A CN111611695 A CN 111611695A
Authority
CN
China
Prior art keywords
sample
parameter
parameters
simulation
particle
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
CN202010393302.5A
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.)
Taiyuan University of Technology
Original Assignee
Taiyuan University of Technology
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 Taiyuan University of Technology filed Critical Taiyuan University of Technology
Priority to CN202010393302.5A priority Critical patent/CN111611695A/zh
Publication of CN111611695A publication Critical patent/CN111611695A/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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • 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)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种模拟岩土材料时离散元线性刚度参数的自动标定方法。该方法利用解析公式作为线性弹性接触参数的初始估计值，通过将颗粒刚度参数当作自变量，数值模拟实际结果当作因变量，实现自动训练离散元模型中的接触参数。通过本发明，能够自动训练得到目标参数，极大减少了传统人工校核方法所需的时间；将基于各向同性散体试样的宏观变形解析解作为参数训练的初始值，有效避免了参数训练过程中的局部最优陷阱和训练时间过长等问题；通过在小应变条件下试样的响应来标定弹性参数，参数标定效率高。本发明实施过程简单，适用于快速标定砂土、堆石料等各类散体材料以及沉积岩石和陶瓷等粘结散体在采用离散元进行模拟研究时的细观变形参数。

Description

一种模拟岩土材料时离散元线性刚度参数的自动标定方法
技术领域
本发明属于岩土工程领域，具体的说，本发明涉及一种模拟岩土材料时离散元线性刚度参数的自动标定方法。
背景技术
由于砂土和堆石料等土工材料具有天然的离散特性，以颗粒为基本单元的离散元方法广泛应用于岩土材料和岩土工程问题的分析和研究中。然而，由于离散元算法采用的细观尺度参数不易直接通过物理试验测得，当前大多数针对岩土工程的离散元模拟，首先都要通过模拟某类常规岩土试验(如三轴试验)，通过人工不断调试离散元细观参数，直到模拟对象在目标宏观指标上和物理试验基本一致时，所用的模拟参数才被当作是一组可靠的输入参数。传统的人工标定方法不仅费时费力，而且难以达到一个较高的精度。
线性模型和赫兹模型是广泛应用于离散元模拟的两个基本接触模型。各个模型所采用的输入参数不一。其中，线性接触模型采用法向接触刚度和切向接触刚度作为控制变形的参数；赫兹模型采用颗粒剪切模量和泊松比来作为控制变形的参数。由于细观变形机理的不同，两者在标定策略上存在差异。本发明主要解决线性模型中的细观变形参数(即：法向接触刚度和切向接触刚度)的参数标定问题。
近年来，离散元参数的标定问题得到了广泛关注。
中国发明专利(申请号：201810632432.2，专利名称：一种快速确定岩石类脆性材料离散元模型参数的方法)提供了一种确定岩石类脆性材料离散元模型参数的方法。该方法需要模拟足够数量数值试样来做参数敏感性分析，以估计参数的取值范围，并为后一步采用统计方法优化参数取值提供依据。该方法为标定离散元参数提供了一套可供参考的流程，然而本质上依然是基于足够样本经验数据的统计拟合方法。该方法的不足之处在于：(1)需要模拟足够数量的数值样本来为敏感性分析和统计拟合提供依据，计算量较大；(2)在数值输出之后，从敏感性分析确定单个宏观指标背后的主要细观影响参数，再到拉丁超立方取样和蚁群算法来优化参数取值，过于依赖统计方法而忽视散体材料在宏细观参数背后的物理联系，且实施过程相对复杂。
中国发明专利(申请号：201811477615.8，专利名称：一种离散元模拟中的细观参数标定方法)提供了一种离散元模拟中的细观参数标定方法。该方法在标定颗粒刚度参数问题上，本质还是传统的人工试错的校核方法，在标定效率上并不能做到实质提高。
中国发明专利(申请号：201910080559.2，专利名称：一种离散元材料自动训练方法)提供了一种离散元材料参数自动训练方法。该方法将基于规则排列的颗粒体宏细观参数解析关系作为材料训练初始值，采用实际宏观参数和目标宏观参数的比值乘以当前细观参数来得到下一步更新的细观参数，并不断训练直到达到目标值。该方法上实现了自动训练，相对传统方法取得了明显进步。然而也存在以下不足：(1)所有细观参数同时训练，过程中各个参数相互干扰，影响标定效率；(2)更新策略来源于经验性的假设，缺乏严格数学基础，无法保证足够的收敛效率和收敛精度；(3)输入细观初始值来源于单一粒径且规则排列颗粒的变形关系，这与实际岩土材料中颗粒组构和粒径大小都随机分布的现实情况差异较大。过于偏离实际散体的宏观响应的初始估计，不仅使得训练时间过长，而且材料在训练过程中面对局部最优陷阱的风险明显增加；(4)参数更新过程忽略了散体的宏细观参数关系，例如：法向刚度和切向刚度同时对宏观参数中的弹性模量和泊松比存在影响，根据该方法使用的参数更新策略，单独只用弹性模量或者泊松比来更新法向刚度或者切向刚度，实际上都存在不妥；(5)不区分二维模型和三维模型在参数训练策略和过程上的差异。例如：对二维和三维模型采用同样的初值估计，模型不能达到最佳收敛效率。
发明内容
本发明的目的在于，为离散元模拟砂土和堆石料等散粒体时，提供一套适用于二维和三维模拟的、兼具高效和稳定性能的线性接触刚度自动标定方法。
本发明的目的通过以下步骤实现：本发明提供了一种模拟岩土材料时离散元线性刚度参数的自动标定方法，包括步骤：
一、通过室内真三轴或者双轴试验方式，确定待模拟岩土材料的杨氏模量和泊松比；
二、根据小应变条件下推导的均匀散体宏细观变形参数解析公式，计算颗粒接触法向刚度和切向刚度的初始估计值；
三、使用计算出来的刚度参数初始估计值，建立离散元双轴或者三轴数值试验；
四、在加载之前，将试样颗粒间的摩擦系数设置为1.0或者更大的值，实施小应变三轴试验数值仿真，直到轴向应变达到一个预设阈值ε时，停止加载，得到当前模拟试样的弹性模量与泊松比；
五、构建误差函数L；
六、判断误差函数L的大小是否满足误差允许值；如满足，当前模型所用刚度参数即为标定参数，标定结束；如不满足，继续执行标定；
七、基于改进的梯度下降方法更新刚度参数；
八、采用更新后的刚度参数重新计算数值模型，更新误差函数大小并评估误差函数大小是否满足误差允许值；如不满足，继续执行步骤七和八所述过程，并直到误差函数低于误差允许值。
其中，颗粒接触法向初始估计值
Figure BDA0002486424920000021
和切向刚度的初始估计值
Figure BDA0002486424920000022
的计算公式为：
二维离散元模型：
Figure BDA0002486424920000023
三维离散元模型：
Figure BDA0002486424920000024
其中，E为试验得到的小应变杨氏模量，ν为试验得到的泊松比，S为二维试样的面积，V为三维试样的体积，Nc为试样中的接触数目；当试样为单粒径分布时，r为颗粒半径；当试样为多粒径分布时，r为试样的中值粒径。
其中，对于待模拟岩土材料散体试样，面积/体积与接触数目的比值计算方法为：
Figure BDA0002486424920000031
其中，2D和3D分别表示2维和3维模型；φ为孔隙率；
Figure BDA0002486424920000032
为颗粒体的平均配位数。
其中，采用生成颗粒时的初始摩擦系数fc来估计最终试样的孔隙率和平均配位数经验值，如下：
Figure BDA0002486424920000033
。
其中，对于多粒径散体试样，采用的中值粒径来计算模型的初始刚度。
其中，真三轴或者双轴模拟结束后，应检查试样加载过程中的应力应变曲线，确保试样在此加载期间仅发生可恢复的弹性变形。
其中，误差函数L为：
Figure BDA0002486424920000034
其中，β1和β2是两个值域为[0,1]的加权系数，且满足β1+β2＝1，E*和ν*分别是数值试样计算得到的弹性模量和泊松比。
其中，若颗粒接触刚度参数为第一次更新，按照如下方式更新：
Figure BDA0002486424920000035
其中，“:＝”为赋值符号，αn和αs为系数，取值在-0.5-0.5之间。
其中，若当前颗粒接触刚度参数为第二次或者更多次更新时，则按照差分方法计算误差函数L对kn和ks的偏导数
Figure BDA0002486424920000036
和
Figure BDA0002486424920000037
若令以kn和ks为自变量的误差函数为f，即令：
f(kn,ks)＝L(E*,ν*) (7)
则两个偏导数按照如下差分方式估计得到：
Figure BDA0002486424920000038
。
其中，若当前颗粒接触刚度参数为第二次或者更多次更新时，则按照改进的梯度下降方法更新，更新公式如下：
Figure BDA0002486424920000041
。
与现有技术相比，本发明的技术效果在于：(1)能自动快速标定离散元模拟过程中颗粒单元中的细观线性刚度参数，极大促进了离散元对岩土工程问题的分析和研究；(2)基于改进的梯度下降方法来迭代标定过程，标定收敛性好，效率高；(3)基于小应变三轴或者双轴加载试验来标定变形参数，模拟时间短且标定结果可靠。
综上所述，本发明是一种高效、稳定、实施流程清晰且编程方便的方法，可为标定离散元的线性刚度参数提供技术支撑，促进离散元技术对岩土工程问题的分析和研究；不仅适用于单一粒径颗粒散体和多粒径散体，而且对二维和三维试样均有较好的标定能力。在数个迭代步内，参数的预测精度可达到1％以内。
附图说明
附图1是所述的离散元线性接触刚度参数自动标定技术路线图；
附图2是所述的二维单粒径随机颗粒试样模型图；
附图3是所述的二维多粒径随机颗粒试样模型图；
附图4是所述的三维单粒径随机颗粒试样模型图。
具体实施方式
为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。
如图1所示，本发明提供了一种模拟岩土材料时离散元线性刚度参数的自动标定方法，包括步骤：
步骤(1)：通过室内真三轴或者双轴试验等方式，确定待模拟岩土试样的宏观参数。所用宏观参数为小应变杨氏模量E与泊松比υ。
步骤(2)：将试验得到的目标宏观参数输入小应变条件下推导的均匀散体宏细观变形参数解析公式(公式如下，推导过程见后文)，分别确定二维和三维模拟情况下，颗粒接触法向刚度和切向刚度的初始估计值
Figure BDA0002486424920000044
和
Figure BDA0002486424920000045
二维离散元模型：
Figure BDA0002486424920000042
三维离散元模型：
Figure BDA0002486424920000043
其中，S为二维试样的面积，V为三维试样的体积，Nc为试样中的参与传递荷载的接触总数目；当试样单粒径分布时，r为颗粒半径；当试样为多粒径分布时，r为试样的中值粒径。试样的面积/体积与接触数目的比值为：
Figure BDA0002486424920000051
其中，2D和3D分别表示2维和3维模型；φ为孔隙率；
Figure BDA0002486424920000052
为颗粒体的平均配位数。对于离散元散体试样，生成颗粒时的初始摩擦系数往往对颗粒体的最终密实程度有很大影响。当初始摩擦系数为fc时，最终试样的孔隙率和平均配位数的值可按下列经验公式估计：
Figure BDA0002486424920000053
步骤(3)：建立离散元双轴或者真三轴数值试验，将步骤(2)计算出的颗粒接触刚度参数作为初始参数投入模型进行数值计算，数值试验中采用和室内试验相同的围压。
步骤(4)：在加载之前，将试样颗粒间的摩擦系数设置为1.0，实施小应变三轴试验数值仿真，直到轴向应变达到一个较小值ε时，停止加载。(例如：轴向应变ε＝0.05％)。模拟结束后，应检查试样加载过程中的应力应变曲线，确保试样在此加载期间仅发生可恢复的弹性变形。将颗粒间摩擦系数设置为1.0的目的是使试样在整体发生较小轴向变形时，颗粒间不发生滑动变形(塑性变形)。计算数值试样模拟得到的弹性模量E*和泊松比ν*。计算方法为：
Figure BDA0002486424920000054
Figure BDA0002486424920000055
步骤(5)：构建误差函数L:
Figure BDA0002486424920000056
其中，β1和β2是两个值域为[0,1]的加权系数，且满足β1+β2＝1。当弹性模量和泊松比需要同时标定时，β1＝β2＝0.5；当仅以弹性模量为目标宏观参数进行标定时，β1＝1，β1＝0.
步骤(6)：判断误差函数L的大小是否满足误差允许值。如满足，当前模型所用刚度参数即为标定参数，标定结束；如不满足，继续执行标定步骤。具体过程为：将步骤(4)模拟得到的弹性模量E*和泊松比ν*代入步骤(5)构建的误差函数，判断误差函数是否小于误差允许值(例如：0.01％)。如满足，则使用的细观参数为模拟该试验的合理估计参数；如不满足，则按照步骤(7)的描述更新参数。
步骤(7)：更新刚度参数。若颗粒接触刚度参数为第一次更新，按照如下方式更新：
Figure BDA0002486424920000061
其中，“:＝”为赋值符号，αn和αs为系数，一般取值为-0.5-0.5之间。
若当前颗粒接触刚度参数为第二次或者更多次更新时，则按照改进的梯度下降方法更新，更新公式如下：
Figure BDA0002486424920000062
其中，η为学习速率，
Figure BDA0002486424920000063
和
Figure BDA0002486424920000064
分别为误差函数L对kn和ks的偏导数。若令以kn和ks为自变量的误差函数为f，即令：
f(kn,ks)＝L(E*,ν*) (7)
则两个偏导数可按照差分方式进行估计得到：
Figure BDA0002486424920000065
步骤(8)：采用更新后的刚度参数重新计算数值模型，更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足，继续执行步骤(7)和(8)所述过程，并直到误差函数低于误差允许值。具体过程为：将步骤(7)更新得到的刚度参数进行离散元双轴或者真三轴压缩试验计算，根据新计算得到的弹性模量E*和泊松比ν*计算误差函数的大小，并判断误差函数是否小于误差容许值。如小于，结束计算，当前模型的刚度参数即为标定参数；如未满足误差要求，则按照步骤(7)更新参数值并重新计算，直至误差满足要求。
附：小应变条件下的均匀散体宏细观变形参数解析公式推导过程
(1)基本假设和推导思路
基本假设：均匀散体试样在整个试样空间内变形均匀(均匀应变)，散体试样的变形为小变形(可恢复的弹性变形)。
推导思路：使用连续介质力学理论，将散体试样均质化并等效为对应连续体。颗粒***与对应的等效连续体在应变能上保持一致。
注：以下推导过程基于张量运算完成。张量计算中，笛卡尔坐标系下的x轴，y轴和z轴，通常采用下标1,2,3表示；张量x的分量记法为xi(一阶)或者xij(二阶)，其中i和j可能是1,2,3中任意一项。
(2)应变能计算：
储存在颗粒***中的应变能量U为：
Figure BDA0002486424920000066
其中：Nc为整体颗粒***中的参与传递荷载的接触总数目；
Figure BDA0002486424920000067
和
Figure BDA0002486424920000068
分别代表颗粒***中第k个接触的法向和切向相对位移，k＝1,…Nc；
Figure BDA0002486424920000069
和
Figure BDA00024864249200000610
分别为接触k传递的法向和切向接触力；
Figure BDA0002486424920000071
和
Figure BDA0002486424920000072
分别指在接触k在法向和切向上的微小变形。
对于线性接触模型而言，在接触未发生滑动前，接触力与颗粒间相对位移的关系为：
Figure BDA0002486424920000073
其中：kn和ks是颗粒的法向和切向刚度。储存在第k个线性接触中的应变能Uk为：
Figure BDA0002486424920000074
(3)颗粒位移与应变的关系
假设颗粒A和B是试样内任意两个相互接触的颗粒，令颗粒的位移与等效连续体在对应点上的位移相等，则颗粒间的相对位移
Figure BDA0002486424920000075
与对应等效连续体上的应变
Figure BDA0002486424920000076
的关系为：
Figure BDA0002486424920000077
其中：
Figure BDA0002486424920000078
和
Figure BDA0002486424920000079
分别为颗粒A和B沿xj方向的坐标，Lk为颗粒A和B之间的距离，
Figure BDA00024864249200000710
为颗粒A与B之间的等效应变。
Figure BDA00024864249200000711
为第k个接触的第j个分量。
颗粒间的法向相对位移
Figure BDA00024864249200000712
可表示为：
Figure BDA00024864249200000713
颗粒间的切向相对位移
Figure BDA00024864249200000714
可表示为：
Figure BDA00024864249200000715
将公式(14)，(16)和(17)带入公式(12)，则颗粒***的应变能可以表示为：
Figure BDA00024864249200000716
(4)颗粒***的等效弹性刚度矩阵
整个颗粒***的应变能密度可使用应变能除以颗粒***的面积(二维)或者体积(三维)来获得：
Figure BDA00024864249200000717
其中，S为二维试样的面积，V为三维试样的体积。为简化推导过程，以下推导将基于三维颗粒***完成，二维模型的结果可通过将体积V更换为面积S，并忽略各个物理张量的第3个指标得到。
基于在散体内变形均匀的假定(应变各处相等)，所以颗粒***中任意位置的局部应变等于***的总体应变：
Figure BDA00024864249200000718
根据弹性力学理论，一个连续体的应力张量可通过对应变能密度对应变张量求偏微分得到：
Figure BDA0002486424920000081
颗粒***的弹性刚度矩阵Cijmn可进一步通过使用应力张量σij对应变张量εmn取偏微分得到：
Figure BDA0002486424920000082
其中，δin为克罗内克函数。
(5)均匀颗粒***中宏细观变形参数的解析关系
在颗粒尺寸一致，内部组构各向同性的均匀颗粒***中，颗粒***的弹性刚度矩阵Cijmn可进一步简化为：
Figure BDA0002486424920000083
在弹性力学中，一般各向同性弹性连续体刚度矩阵Cijmn通常表述为：
Figure BDA0002486424920000084
其中，E和ν为一般各向同性弹性连续体的弹性模量和泊松比。
通过将颗粒***的弹性刚度矩阵(公式(23))与一般各向同性弹性连续体的刚度矩阵(公式(24))对比，则可以求出小应变情况下颗粒***宏细观变形参数间的关系。
Figure BDA0002486424920000085
Figure BDA0002486424920000086
联立公式(25)和公式(26)，可推导得出颗粒刚度参数的估计公式为：
Figure BDA0002486424920000087
Figure BDA0002486424920000088
此公式基于半径相同的各向同性均匀颗粒***推导得出，对于具有任意颗粒粒径的颗粒***，其内部组构不可避免存在一定程度的各向异性特征。因此，本解析公式在本发明中仅用于对真实颗粒刚度参数的估计值。精确的结果将由本专利提出的迭代技术方案自动训练得出。
具体实施方式一
将本发明的方法应用于二维单粒径随机颗粒试样(如图2)的标定，试验参数如下：
Figure BDA0002486424920000091
表1 二维单粒径随机排列颗粒模型参数
模拟步骤如下：
步骤(1)：通过室内试验方式，确定待模拟岩土材料的宏观参数。本实施例中岩土材料的宏观参数为：杨氏模量E为10Gpa，泊松比υ为0.2。
步骤(2)：将试验得到的目标宏观参数输入公式，分别确定二维模型中，颗粒接触法向刚度和切向刚度的初始估计值
Figure BDA0002486424920000092
和
Figure BDA0002486424920000093
Figure BDA0002486424920000094
其中，S为二维试样的面积，Nc为试样中的接触数目，r为颗粒体的中值粒径，φ为孔隙率；
Figure BDA0002486424920000096
为颗粒体的平均配位数。
步骤(3)：建立离散元双轴数值试验，将步骤(2)计算出的颗粒接触刚度性参数作为初始参数放入模型进行数值计算，数值试验中采用和室内试验相同的围压，为1MPa。
步骤(4)：在加载之前，将试样颗粒间的摩擦系数设置为1.0，实施小应变三轴试验数值仿真，直到轴向应变达到一个较小值0.05％时，停止加载。检查试样加载过程中的应力应变曲线，确认试样在此加载期间仅发生可恢复的弹性变形。计算数值试样模拟得到的弹性模量E*和泊松比ν*。
步骤(5)：构建误差函数L:
Figure BDA0002486424920000095
其中，β1＝β2＝0.5。
步骤(6)：判断误差函数L的大小是否满足误差允许值。得到第一组误差函数的大小为0.0671，不满足容许误差要求。
步骤(7)：更新刚度参数。若颗粒接触刚度参数为第一次更新，按照如下方式更新：
Figure BDA0002486424920000101
其中，“:＝”为赋值符号，αn和αs为系数，均取值为0.22。
若当前颗粒接触刚度参数为第二次或者更多次更新时，则按照改进的梯度下降方法更新，更新公式如下：
Figure BDA0002486424920000102
其中，η为学习速率，
Figure BDA0002486424920000103
和
Figure BDA0002486424920000104
分别为误差函数L对kn和ks的偏导数。若令以kn和ks为自变量的误差函数为f，即令：
f(kn,ks)＝L(E*,ν*) (7)
则两个偏导数可按照差分方式进行估计得到：
Figure BDA0002486424920000105
步骤(8)：采用更新后的刚度参数重新计算数值模型，更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足，继续执行步骤(7)和(8)所述过程，并直到误差函数低于误差允许值。具体过程为：将步骤(7)更新得到的刚度参数进行离散元双轴或者真三轴压缩试验计算，根据新计算得到的弹性模量E*和泊松比ν*计算误差函数的大小，并判断误差函数是否小于误差容许值。如小于，结束计算，当前模型的刚度参数即为标定参数；如未满足误差要求，则按照步骤(7)更新参数值并重新计算。本模型迭代至第9次时，误差函数大小为2.80E-05，误差满足要求，标定过程结束，得到最优刚度参数估计值为：kn＝9605966N/m,ks＝3687070N/m。记录标定中的参数更新过程如表2所示：
Figure BDA0002486424920000106
表2 二维单粒径随机颗粒试样标定过程
具体实施方式二
将本发明的方法应用于二维多粒径随机颗粒试样的标定(如图3)，试验参数如表3：
Figure BDA0002486424920000111
表3 二维多粒径随机排列颗粒模型参数
模拟步骤如下：
步骤(1)：通过室内试验方式，确定待模拟岩土材料的宏观参数。本实施例中岩土材料的宏观参数为：杨氏模量E为10Gpa，泊松比υ为0.2。
步骤(2)：将试验得到的目标宏观参数输入公式，分别确定二维模型中，颗粒接触法向刚度和切向刚度的初始估计值
Figure BDA0002486424920000112
和
Figure BDA0002486424920000113
Figure BDA0002486424920000114
其中，S为二维试样的面积，Nc为试样中的接触数目，r为颗粒体的中值粒径，φ为孔隙率；
Figure BDA0002486424920000115
为颗粒体的平均配位数。
步骤(3)：建立离散元双轴数值试验(如图3)，将步骤(2)计算出的颗粒接触刚度性参数作为初始参数放入模型进行数值计算，数值试验中采用和室内试验相同的围压，为1MPa。
步骤(4)：在加载之前，将试样颗粒间的摩擦系数设置为1.0，实施小应变三轴试验数值仿真，直到轴向应变达到一个较小值0.05％时，停止加载。检查试样加载过程中的应力应变曲线，确认试样在此加载期间仅发生可恢复的弹性变形。计算数值试样模拟得到的弹性模量E*和泊松比ν*。
步骤(5)：构建误差函数L:
Figure BDA0002486424920000116
其中，β1＝β2＝0.5。
步骤(6)：判断误差函数L的大小是否满足误差允许值。得到第一组误差函数的大小为0.0253，不满足容许误差要求。
步骤(7)：更新刚度参数。若颗粒接触刚度参数为第一次更新，按照如下方式更新：
Figure BDA0002486424920000117
其中，“:＝”为赋值符号，αn和αs为系数，均取值为0.2。
若当前颗粒接触刚度参数为第二次或者更多次更新时，则按照改进的梯度下降方法更新，更新公式如下：
Figure BDA0002486424920000121
其中，η为学习速率，
Figure BDA0002486424920000122
和
Figure BDA0002486424920000123
分别为误差函数L对kn和ks的偏导数。若令以kn和ks为自变量的误差函数为f，即令：
f(kn,ks)＝L(E*,ν*) (7)
则两个偏导数可按照差分方式进行估计得到：
Figure BDA0002486424920000124
步骤(8)：采用更新后的刚度参数重新计算数值模型，更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足，继续执行步骤七和八所述过程，并直到误差函数低于误差允许值。具体过程为：将步骤(7)更新得到的刚度参数进行离散元双轴或者三轴压缩试验计算，根据新计算得到的弹性模量E*和泊松比ν*计算误差函数的大小，并判断误差函数是否小于误差容许值。如小于，结束计算，当前模型的刚度参数即为标定参数；如未满足误差要求，则按照步骤(7)更新参数值并重新计算。本模型迭代至第7次时，误差函数大小为6.65E-05，误差满足要求，标定过程结束，得到最优刚度参数估计值为：kn＝10533175N/m,ks＝3546570N/m。记录标定中的参数更新的过程如表4所示：
Figure BDA0002486424920000125
表4 二维多粒径随机颗粒试样标定过程
具体实施方式三
将本发明的方法应用于三维单粒径随机颗粒试样的标定(如图4)，试验参数如表5：
Figure BDA0002486424920000131
表5 三维多粒径随机排列颗粒模型参数
模拟步骤如下：
步骤(1)：通过室内试验方式，确定待模拟岩土材料的宏观参数。本实施例中岩土材料的宏观参数为：杨氏模量E为10Gpa，泊松比υ为0.2。
步骤(2)：将试验得到的目标宏观参数输入公式，分别确定三维模型中，颗粒接触法向刚度和切向刚度的初始估计值
Figure BDA0002486424920000132
和
Figure BDA0002486424920000133
Figure BDA0002486424920000134
其中，V为三维试样的体积，Nc为试样中的接触数目，r为颗粒体的中值粒径，φ为孔隙率；
Figure BDA0002486424920000135
为颗粒体的平均配位数。
步骤(3)：建立离散元双轴数值试验(如图3)，将步骤(2)计算出的颗粒接触刚度性参数作为初始参数放入模型进行数值计算，数值试验中采用和室内试验相同的围压，为1MPa。
步骤(4)：在加载之前，将试样颗粒间的摩擦系数设置为1.0，实施小应变三轴试验数值仿真，直到轴向应变达到一个较小值0.05％时，停止加载。检查试样加载过程中的应力应变曲线，确认试样在此加载期间仅发生可恢复的弹性变形。计算数值试样模拟得到的弹性模量E*和泊松比ν*。
步骤(5)：构建误差函数L:
Figure BDA0002486424920000136
若只以弹性模量为校核目标，则β1＝1，β2＝0。
步骤(6)：判断误差函数L的大小是否满足误差允许值。得到第一组误差函数的大小为0.0253，不满足容许误差要求。
步骤(7)：更新刚度参数。若颗粒接触刚度参数为第一次更新，按照如下方式更新：
Figure BDA0002486424920000137
其中，“:＝”为赋值符号，αn和αs为系数，分别取值为0.1和-0.5。
若当前颗粒接触刚度参数为第二次或者更多次更新时，则按照改进的梯度下降方法更新，更新公式如下：
Figure BDA0002486424920000138
其中，η为学习速率，
Figure BDA0002486424920000141
和
Figure BDA0002486424920000142
分别为误差函数L对kn和ks的偏导数。若令以kn和ks为自变量的误差函数为f，即令：
f(kn,ks)＝L(E*,ν*) (7)
则两个偏导数可按照差分方式进行估计得到：
Figure BDA0002486424920000143
步骤(8)：采用更新后的刚度参数重新计算数值模型，更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足，继续执行步骤七和八所述过程，并直到误差函数低于误差允许值。具体过程为：将步骤(7)更新得到的刚度参数进行离散元双轴或者三轴压缩试验计算，根据新计算得到的弹性模量E*和泊松比ν*计算误差函数的大小，并判断误差函数是否小于误差容许值。如小于，结束计算，当前模型的刚度参数即为标定参数；如未满足误差要求，则按照步骤(7)更新参数值并重新计算。本模型迭代至第7次时，误差函数大小为6.65E-05，误差满足要求，标定过程结束，得到最优刚度参数估计值为：kn＝6087605N/m,ks＝2056510N/m。记录标定中的参数更新的过程如表6所示：
Figure BDA0002486424920000144
表6 三维单粒径随机颗粒试样标定过程
综上所述，本发明是一种高效、稳定、实施流程清晰且编程方便的方法，可为标定离散元的线性刚度参数提供技术支撑，促进离散元技术对岩土工程问题的分析和研究；不仅适用于单一粒径颗粒散体和多粒径散体，而且对二维和三维试样均有较好的标定能力。在数个迭代步内，参数的预测精度可达到1％以内。
与现有技术相比，本发明的技术效果在于：(1)能自动快速标定离散元模拟过程中颗粒单元中的细观线性刚度参数，极大促进了离散元对岩土工程问题的分析和研究；(2)基于改进的梯度下降方法来迭代标定过程，标定收敛性好，效率高；(3)基于小应变三轴或者双轴加载试验来标定变形参数，模拟时间短且标定结果可靠。
上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (10)

1.一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于包括以下步骤：
一、通过室内真三轴或者双轴试验方式，确定待模拟岩土材料试样的杨氏模量和泊松比；
二、根据小应变条件下推导的均匀散体宏细观变形参数解析公式，计算颗粒接触法向刚度和切向刚度的初始估计值；
三、使用计算出来的刚度参数初始估计值，建立待模拟岩土材料试样离散元双轴或者三轴数值试验；
四、在加载之前，将试验试样颗粒间的摩擦系数设置为1.0或者更大的值，实施小应变三轴试验数值仿真，直到轴向应变达到一个预设阈值ε时，停止加载，得到当前模拟试样的弹性模量与泊松比；
五、构建误差函数L；
六、判断误差函数L的大小是否满足误差允许值；如满足，当前模型所用刚度参数即为标定参数，标定结束；如不满足，继续执行标定；
七、基于改进的梯度下降方法更新刚度参数；
八、采用更新后的刚度参数重新计算数值模型，更新误差函数大小并评估误差函数大小是否满足误差允许值；如不满足，继续执行步骤七和八所述过程，并直到误差函数低于误差允许值。
2.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于，颗粒接触法向初始估计值
Figure FDA0002486424910000014
和切向刚度的初始估计值
Figure FDA0002486424910000015
的计算公式为：
二维离散元模型：
Figure FDA0002486424910000011
三维离散元模型：
Figure FDA0002486424910000012
其中，E为试验得到的小应变杨氏模量，ν为试验得到的泊松比，S为二维试样的面积，V为三维试样的体积，Nc为试样中的接触数目；当试样为单粒径分布时，r为颗粒半径；当试样为多粒径分布时，r为试样的中值粒径。
3.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于，对于待模拟岩土材料散体试样，面积/体积与接触数目的比值计算方法为：
Figure FDA0002486424910000013
其中，2D和3D分别表示2维和3维模型；φ为孔隙率；
Figure FDA0002486424910000027
为颗粒体的平均配位数。
4.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于，采用生成颗粒时的初始摩擦系数fc来估计最终试样的孔隙率和平均配位数经验值，如下：
Figure FDA0002486424910000021
。
5.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于，对于多粒径散体试样，采用的中值粒径来计算模型的初始刚度。
6.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于，真三轴或者双轴模拟结束后，应检查试样加载过程中的应力应变曲线，确保试样在此加载期间仅发生可恢复的弹性变形。
7.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于，误差函数L为：
Figure FDA0002486424910000022
其中，β1和β2是两个值域为[0,1]的加权系数，且满足β1+β2＝1，E*和ν*分别是数值试样计算得到的弹性模量和泊松比。
8.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于：若颗粒接触刚度参数为第一次更新，按照如下方式更新：
Figure FDA0002486424910000023
其中，“:＝”为赋值符号，αn和αs为系数，取值在-0.5-0.5之间。
9.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于，若当前颗粒接触刚度参数为第二次或者更多次更新时，则按照差分方法计算误差函数L对kn和ks的偏导数
Figure FDA0002486424910000024
和
Figure FDA0002486424910000025
若令以kn和ks为自变量的误差函数为f，即令：
f(kn,ks)＝L(E*,ν*) (7)
则两个偏导数按照如下差分方式估计得到：
Figure FDA0002486424910000026
。
10.根据权利要求1所述的一种模拟岩土材料时离散元线性刚度参数的自动标定方法，其特征在于，若当前颗粒接触刚度参数为第二次或者更多次更新时，则按照改进的梯度下降方法更新，更新公式如下：
Figure FDA0002486424910000031
。
CN202010393302.5A 2020-05-11 2020-05-11 一种模拟岩土材料时离散元线性刚度参数的自动标定方法 Pending CN111611695A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010393302.5A CN111611695A (zh) 2020-05-11 2020-05-11 一种模拟岩土材料时离散元线性刚度参数的自动标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010393302.5A CN111611695A (zh) 2020-05-11 2020-05-11 一种模拟岩土材料时离散元线性刚度参数的自动标定方法

Publications (1)

Publication Number Publication Date
CN111611695A true CN111611695A (zh) 2020-09-01

Family

ID=72200225

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010393302.5A Pending CN111611695A (zh) 2020-05-11 2020-05-11 一种模拟岩土材料时离散元线性刚度参数的自动标定方法

Country Status (1)

Country Link
CN (1) CN111611695A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112098209A (zh) * 2020-09-15 2020-12-18 河北工业大学 一种岩土颗粒破坏局部化识别方法
CN112163328A (zh) * 2020-09-18 2021-01-01 武汉大学 一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法
CN112765895A (zh) * 2021-01-28 2021-05-07 南京大学 一种基于机器学习的岩土材料离散元自动建模方法
CN113191004A (zh) * 2021-05-09 2021-07-30 甘肃省地震局 黄土土体内部刚度计算方法
CN115828747A (zh) * 2022-12-01 2023-03-21 山东大学 考虑颗粒互锁效应的裂隙岩体参数智能标定方法与***

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130289953A1 (en) * 2012-01-24 2013-10-31 The University Of Akron Self-optimizing, inverse analysis method for parameter identification of nonlinear material constitutive models
CN109030202A (zh) * 2018-06-19 2018-12-18 湘潭大学 一种快速确定岩石类脆性材料离散元模型参数的方法
CN109543337A (zh) * 2018-12-05 2019-03-29 陕西理工大学 一种离散元模拟中的细观参数标定方法
CN109612885A (zh) * 2019-01-08 2019-04-12 东北大学 一种基于离散元法的矿物颗粒模型参数标定方法
CN109815599A (zh) * 2019-01-28 2019-05-28 南京大学 一种离散元材料自动训练方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130289953A1 (en) * 2012-01-24 2013-10-31 The University Of Akron Self-optimizing, inverse analysis method for parameter identification of nonlinear material constitutive models
CN109030202A (zh) * 2018-06-19 2018-12-18 湘潭大学 一种快速确定岩石类脆性材料离散元模型参数的方法
CN109543337A (zh) * 2018-12-05 2019-03-29 陕西理工大学 一种离散元模拟中的细观参数标定方法
CN109612885A (zh) * 2019-01-08 2019-04-12 东北大学 一种基于离散元法的矿物颗粒模型参数标定方法
CN109815599A (zh) * 2019-01-28 2019-05-28 南京大学 一种离散元材料自动训练方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
TONGMING QU等: "Calibration of linear contact stiffnesses in discrete element models using a hybrid analytical-computational framework", 《POWDER TECHNOLOGY》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112098209A (zh) * 2020-09-15 2020-12-18 河北工业大学 一种岩土颗粒破坏局部化识别方法
CN112163328A (zh) * 2020-09-18 2021-01-01 武汉大学 一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法
CN112765895A (zh) * 2021-01-28 2021-05-07 南京大学 一种基于机器学习的岩土材料离散元自动建模方法
CN112765895B (zh) * 2021-01-28 2023-10-17 南京大学 一种基于机器学习的岩土材料离散元自动建模方法
CN113191004A (zh) * 2021-05-09 2021-07-30 甘肃省地震局 黄土土体内部刚度计算方法
CN113191004B (zh) * 2021-05-09 2022-08-30 甘肃省地震局 黄土土体内部刚度计算方法
CN115828747A (zh) * 2022-12-01 2023-03-21 山东大学 考虑颗粒互锁效应的裂隙岩体参数智能标定方法与***
CN115828747B (zh) * 2022-12-01 2023-09-26 山东大学 考虑颗粒互锁效应的裂隙岩体参数智能标定方法与***

Similar Documents

Publication Publication Date Title
CN111611695A (zh) 一种模拟岩土材料时离散元线性刚度参数的自动标定方法
Shi et al. Calibration of micro-scaled mechanical parameters of granite based on a bonded-particle model with 2D particle flow code
Armaghani et al. Indirect measure of shale shear strength parameters by means of rock index tests through an optimized artificial neural network
Yu et al. An intelligent displacement back-analysis method for earth-rockfill dams
Kahraman et al. The usability of Cerchar abrasivity index for the prediction of UCS and E of Misis Fault Breccia: regression and artificial neural networks analysis
CN105259331B (zh) 一种节理岩体单轴强度预测方法
Kim et al. Prediction of subgrade resilient modulus using artificial neural network
CN111610091A (zh) 一种模拟岩土材料时离散元赫兹接触参数自动标定方法
CN112163328B (zh) 一种基于深度学习和数据驱动的岩土颗粒材料本构建模方法
Groh et al. Damage simulation of brittle heterogeneous materials at the grain size level
Najjar et al. Simulating the stress–strain behavior of Georgia kaolin via recurrent neuronet approach
Kayadelen Estimation of effective stress parameter of unsaturated soils by using artificial neural networks
JP6789274B2 (ja) Ｂｐｍ理論に基づくメゾの力チェーン粒子モデルの構築方法
Sarmadian et al. Modeling of some soil properties using artificial neural network and multivariate regression in Gorgan Province, North of Iran
CN111475876A (zh) 获取粒料动态回弹力学特征参数的方法
Iacono et al. Estimation of model parameters in nonlocal damage theories by inverse analysis techniques
Ma et al. Crack evolution and acoustic emission characteristics of rock specimens containing random joints under uniaxial compression
CN110457746B (zh) 基于bp神经网络的结构面峰值抗剪强度预测模型构建方法
CN110929412B (zh) 一种基于dda理论的节理摩擦系数动态衰减计算方法
CN115508206B (zh) 一种节理岩质边坡岩体强度参数概率反演方法
Dabeet et al. Simulation of cyclic direct simple shear loading response of soils using discrete element modeling
CN112884739B (zh) 一种基于深度学习网络的堆石体填筑密实度快速检测方法
Frery et al. Stochastic particle packing with specified granulometry and porosity
CN110750901B (zh) 基于离散元模型的土体扰动范围判断方法
Li et al. Parameter Estimation of Particle Flow Model for Soils Using Neural Networks.

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
WD01 Invention patent application deemed withdrawn after publication
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20200901