CN111610091A - 一种模拟岩土材料时离散元赫兹接触参数自动标定方法 - Google Patents
一种模拟岩土材料时离散元赫兹接触参数自动标定方法 Download PDFInfo
- Publication number
- CN111610091A CN111610091A CN202010394118.2A CN202010394118A CN111610091A CN 111610091 A CN111610091 A CN 111610091A CN 202010394118 A CN202010394118 A CN 202010394118A CN 111610091 A CN111610091 A CN 111610091A
- Authority
- CN
- China
- Prior art keywords
- particle
- sample
- discrete element
- poisson
- error function
- 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
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N3/00—Investigating strength properties of solid materials by application of mechanical stress
- G01N3/08—Investigating strength properties of solid materials by application of mechanical stress by applying steady tensile or compressive forces
-
- 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
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
- G01N2203/0058—Kind of property studied
- G01N2203/0069—Fatigue, creep, strain-stress relations or elastic constants
- G01N2203/0075—Strain-stress relations or elastic constants
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
- G01N2203/02—Details not specific for a particular testing method
- G01N2203/025—Geometry of the test
- G01N2203/0254—Biaxial, the forces being applied along two normal axes of the specimen
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
- G01N2203/02—Details not specific for a particular testing method
- G01N2203/025—Geometry of the test
- G01N2203/0256—Triaxial, i.e. the forces being applied along three normal axes of the specimen
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- 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)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Analytical Chemistry (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- General Engineering & Computer Science (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
Abstract
本发明公开了一种模拟岩土材料时离散元赫兹接触参数自动标定方法,将离散元标定的目标宏观参数代入解析公式计算出颗粒剪切模量和泊松比,以此作为颗粒参数初始估计值,建立数值三轴或双轴试验,得到宏观杨氏模量和泊松比;计算误差函数大小,判断是否更新颗粒参数;根据迭代次数的不同,分别采取不同的策略;初次迭代时,通过人工给予一个与初始估计值成正比的较小扰动;更多次迭代时采用自适应矩估计策略更新参数。通过实施案例表明,单粒径或多粒径离散元试样,在有限次迭代后,标定参数的误差被有效控制在一定范围内。该标定方法具备自动标定能力,收敛能力强,标定效率高,能有效解决离散元模拟岩土体材料时赫兹变形参数的标定问题。
Description
技术领域
本发明属于岩土工程领域,具体的说,本发明涉及一种模拟岩土材料时离散元赫兹接触参数自动标定方法。
背景技术
离散元方法广泛应用于地基加固、隧道开挖等岩土工程问题的分析和研究中。然而,由于离散元算法采用的细观尺度参数不易直接通过物理试验测得,当前大多数针对岩土工程的离散元模拟,首先都要通过模拟某类常规岩土试验(如三轴试验),通过不断调试离散元细观参数,直到模拟对象在目标宏观指标上和物理试验基本一致时,所用的模拟参数才被当作是一组可靠的输入参数。传统的人工标定方法不仅费时费力,而且标定精度相当有限。
线性接触模型和非线性的赫兹模型是广泛应用于离散元模拟的两个基本接触模型。各个模型所采用的输入参数不一。其中,线性接触模型采用法向接触刚度和切向接触刚度作为控制变形的参数;非线性的赫兹模型采用颗粒剪切模量和泊松比来作为控制变形的参数。由于细观变形机理的不同,两者在标定策略上存在差异。本发明主要解决非线性的赫兹模型中细观变形参数(即:颗粒剪切模量和泊松比)的参数标定问题。
近年来,离散元参数的标定问题得到了广泛关注。
中国发明专利(申请号:201910080559.2,专利名称:一种离散元材料自动训练方法)提供了一种离散元材料参数自动训练方法。该方法将基于规则排列的颗粒体宏细观参数解析关系作为材料训练初始值,采用实际宏观参数和目标宏观参数的比值乘以当前细观参数来得到下一步更新的细观参数,并不断训练直到达到目标值。该方法上实现了自动训练,相对传统方法取得了明显进步。然而也存在以下不足:(1)所有细观参数同时训练,过程中各个参数相互干扰,影响标定效率;(2)更新策略来源于经验性的假设,缺乏严格数学基础,无法保证足够的收敛效率和收敛精度;(3)输入细观初始值来源于单一粒径且规则排列颗粒的变形关系,这与实际岩土材料中颗粒组构和粒径大小都随机分布的现实情况差异较大。过于偏离实际散体的宏观响应的初始估计,不仅使得训练时间过长,而且材料在训练过程中面对局部最优陷阱的风险明显增加;(4)参数更新过程忽略了散体的宏细观参数关系,例如:法向刚度和切向刚度同时对宏观参数中的弹性模量和泊松比存在影响,根据该方法使用的参数更新策略,单独只用弹性模量或者泊松比来更新法向刚度或者切向刚度,实际上都存在不妥;(5)不区分二维模型和三维模型在参数训练策略和过程上的差异。例如:对二维和三维模型采用同样的初值估计,模型不能达到最佳收敛效率。
发明内容
本发明的目的在于,为离散元模拟砂土和堆石料等散粒体时,提供一套适用于二维和三维模型的、兼具高效和稳定性能的非线性赫兹接触参数自动标定方法。
本发明的目的通过以下步骤实现:本发明提供了一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,包括以下步骤:
一、在使用离散元模拟岩土材料前,需采用室内三轴或者双轴试验,确定待模拟岩土材料的杨氏模量和泊松比;
二、以试验得到的宏观参数为目标参数,根据解析公式计算颗粒剪切模量和泊松比的初始估计值;
三、使用计算出来的颗粒剪切模量和泊松比作为输入参数,建立离散元双轴或者三轴数值试验;
四、在加载之前,将试样颗粒间的摩擦系数设置为1.0或者更大的值,实施小应变三轴/双轴试验数值仿真,直到轴向应变达到预设阈值,停止加载;
五、构建误差函数L;
六、判断误差函数L的大小是否满足误差允许值;如满足,当前模型所用颗粒剪切模量和泊松比即为标定参数,标定结束;如不满足,继续执行标定;
七、首次更新颗粒剪切模量和泊松比并再次评估误差函数大小;
八、计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值;
九、计算未经偏差修正的一阶和二阶矩估计值;
十、修正一阶和二阶矩估计;
十一、基于修正的一阶和二阶矩估计值,计算第二次及更多次迭代步中颗粒剪切模量和泊松比的更新值;
十二、采用更新后的颗粒剪切模量和泊松比重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值;如不满足,继续执行步骤七到步骤十一所述过程,直到误差函数低于误差允许值。
二维离散元模型:
三维离散元模型:
其中,采用生成颗粒时的初始摩擦系数fc来估计最终试样的孔隙率和平均配位数经验值,如下:
其中,对于多粒径散体试样,采用试样的中值粒径来计算模型的初始接触剪切模量和泊松比。
其中,真三轴或者双轴模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
其中,误差函数L为:
其中,E*是将颗粒间摩擦系数设置为1.0时,使试样在整体发生较小轴向变形时,颗粒间不发生滑动变形,计算数值试样模拟得到的弹性模量。
其中,若颗粒剪切模量G和泊松比ν为第一次迭代更新,按照如下方式更新:
其中,“:=”为赋值符号,η为系数,取值在-0.5-0.5之间。
其中,若当前颗粒接触刚度参数为第二次或者更多次更新时,需按照差分方法计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值gG和gυ,令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
其中,若当前颗粒接触刚度参数为第二次或者更多次更新时,基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
区别于现有技术,本发明的模拟岩土材料时离散元赫兹接触参数自动标定方法首先将离散元标定的目标宏观参数代入解析公式计算出颗粒剪切模量和泊松比。以此解析结果作为颗粒参数初始估计值,建立数值三轴或双轴试验,得到数值模型的宏观杨氏模量和泊松比。计算误差函数大小,判断是否更新颗粒参数。根据迭代次数的不同,分别采取不同的策略。初次迭代时,通过人工给予一个与初始估计值(解析结果)成正比的较小扰动;二次及更多次迭代时则采用自适应矩估计策略来更新参数。通过本发明,能够自动快速标定离散元模拟过程中颗粒剪切模量和泊松比,极大促进了离散元对岩土工程问题的分析和研究;本发明基于自适应的矩估计策略来迭代标定过程,标定收敛性好,效率高;基于小应变三轴或者双轴加载试验来标定变形参数,模拟时间短且标定结果可靠。本发明提供了一种高效、稳定、实施流程清晰且编程方便的方法,可为标定离散元的非线性赫兹变形参数提供技术支撑,促进离散元技术对岩土工程问题的分析和研究。
附图说明
附图1是所述的离散元非线性赫兹接触变形参数自动标定技术路线图;
附图2是所述的二维单粒径随机颗粒试样模型图;
附图3是所述的三维单粒径随机颗粒试样模型图;
附图4是所述的三维多粒径随机颗粒试样模型图;
具体实施方式
为了对本发明的技术特征、目的和效果有更加清楚的理解,现对照附图详细说明本发明的具体实施方式。
参阅图1,本发明提供了一种模拟岩土材料时离散元赫兹接触参数自动标定方法,包括步骤:
二维离散元模型:
三维离散元模型:
对于离散元散体试样,生成颗粒时的初始摩擦系数往往对颗粒体的最终密实程度有很大影响。当初始摩擦系数为fc时,最终试样的孔隙率和平均配位数的值可按下列经验公式估计:
其中,2D和3D分别表示2维和3维模型;
步骤(3):建立离散元双轴或者真三轴数值试验,将步骤(2)计算出的颗粒的剪切模量和泊松比作为初始参数投入模型进行数值计算,数值试验中采用和室内试验相同的围压。
步骤(4):将试样颗粒间的摩擦系数设置为1.0,实施小应变三轴试验数值仿真,直到轴向应变达到一个较小值ε时,停止加载。(例如:轴向应变ε=0.05%)。模拟结束后,应检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。将颗粒间摩擦系数设置为1.0的目的是使试样在整体发生较小轴向变形时,颗粒间不发生滑动变形(塑性变形)。计算数值试样模拟得到的弹性模量E*和泊松比ν*。计算方法为:
其中,Δσ为偏应力,Δε为偏应变,下标x、y、z分别表示笛卡尔坐标系下的三个坐标方向。
步骤(5):构建以弹性模量为目标宏观参数的误差函数L:
步骤(6):判断误差函数L的大小是否满足误差允许值。如满足,当前模型所用刚度参数即为标定参数,标定结束;如不满足,继续执行标定步骤。具体过程为:将步骤(4)模拟得到的弹性模量E*代入步骤(5)构建的误差函数,判断误差函数是否小于误差允许值(例如:0.01%)。如满足,则使用的细观参数为模拟该试验的合理估计参数;如不满足,则按照步骤(7)的描述更新参数。
步骤(7):首次更新颗粒剪切模量和泊松比并再次计算误差函数大小。若是为第一次更新颗粒剪切模量和泊松比,按照如下方式更新:
其中,“:=”为赋值符号,η为系数,一般取值为-0.5-0.5之间。若当前颗粒赫兹接触参数为第二次或者更多次更新时,则按照步骤(8)到步骤(11)进行更新。
步骤(8):计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值。若令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
步骤(9):计算未经偏差修正的一阶和二阶矩估计值。令Mt和Gt分别为第t次迭代时,两个偏导数关于未经偏差修正的一阶和二阶矩估计值按如下方式赋值得到:
Mt=β1Mt-1+(1-β1)gt (13)
步骤(11):基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
其中,ε是为了保持数值稳定性而设置的非常小的常数,一般取值e-7;α为学习率,一般通过尝试0.001,0.01,0.1等几个学习率后确定。
步骤(12):采用更新后的刚度参数重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足,继续执行步骤(7)到步骤(12)所述过程,并直到误差函数低于误差允许值。具体过程为:将步骤(7)更新得到的刚度参数进行离散元双轴或者真三轴压缩试验计算,根据新计算得到的弹性模量E*计算误差函数的大小,并判断误差函数是否小于误差容许值。如小于,结束计算,当前模型的颗粒剪切模量和泊松比即为标定参数;如未满足误差要求,则按照步骤(7)到步骤(12)更新参数值并重新计算,直至误差满足要求。
以下推导过程是本发明采用的以赫兹模型为接触法则的颗粒试样宏细观变形解析关系推导过程
以下推导过程基于张量运算完成。张量计算中,笛卡尔坐标系下的x轴,y轴和z轴,通常采用下标1,2,3表示;张量x的分量记法为xi(一阶)或者xij(二阶),其中i和j可能是1,2,3中任意一项。
(1)颗粒物质力学中的应力表达
应力属于连续介质力学范畴,为连接连续介质力学中的应力概念与细观力学中的颗粒接触力概念,不可避免需要对颗粒模型采取均质化或者平均化的手段。通常情况下,在一定体积V内,等效平均应力增量Δσij可表示为在该体积内颗粒间的接触力增量Δfi c和接触支向量Lj c之间的关系:
其中,Nc为在体积V内的参与荷载传递的接触数目;c表示该范围内的第c个接触。c=1,2,…Nc;
(2)基于均匀应力假设的宏细观参数关系
假设颗粒试样内部应力均匀分布,则颗粒间的接触力可根据如下公式通过颗粒试样中的平均应力得出:
将公式(17)代入公式(18),得到:
其中,δij为克罗内克函数。
根据公式(19),可求出张量Ajk的表达式:
(3)颗粒试样的应变描述
根据能量守恒原理,在不考虑颗粒体力的情况下,外力对试样的做功等于颗粒接触变形积累的应变能,即为:
通过将公式(18)代入公式(21),可得:
由于应力张量Δσij可以为物理材料条件允许下的任意非零值,因此可以得到颗粒试样应变的表达式为:
(4)本构方程与柔度张量
根据弹性力学原理,固体的应力应变关系可用柔度张量Sijkl的形式表达:
εij=Sijklσkl (25)
结合公式(18),(23),(24)和(25),可计算得到该颗粒材料的应力应变柔度张量的大小为:
(5)柔度张量的积分形式与应力应变关系
在颗粒材料中颗粒的大小相同时,颗粒接触支向量Lj c可表达为:
其中,r为颗粒的半径。
将公式(24)和公式(27)代入公式(26)。另外,对于一个包含大量颗粒数目的试样,接触法线的方向在统计上可被近似为各向均匀的。根据接触法向方向均匀分布的假设,柔度张量Sijkl的积分形式可计算得:
通过对公式(28)积分,可得到散体的应力应变关系为:
2D模型:
3D模型:
对于一般各向同性的弹性固体,其应力应变关系为:
2D平面应力模型:
3D模型:
通过对比公式(29)和公式(31),对比公式(30)和公式(32),可得到一般均匀散体的宏观弹性常数与细观参数的关系为:
2D平面应力模型:
3D模型:
(6)基于赫兹接触模型的散体试样小应变模量
对于承受均匀正压力的各向同性材料,可认为在水平和竖直方向,材料中只有主应力,没有剪切应力。基于颗粒间无剪切力,且各个接触法向作用力的大小相等的假设,公式(18)的全量形式可写为:
其中,Nc为颗粒试样中参与荷载传递的有效接触数目;σc为颗粒试样承受的正应力大小。
考虑到赫兹接触模型的接触力与接触位移的作用关系为:
其中,G和ν是细观尺度上的颗粒剪切模量与泊松比。
结合公式(35),(36)和(37),可得到接触刚度的表达形式为:
将公式(38)代入公式(33)和(34),可得:
2D模型:
3D模型:
2D模型:
3D模型:
此公式基于各向同性均匀单粒径颗粒系统推导得出,对于具有多粒径颗粒系统,其内部组构不可避免存在一定程度的各向异性特征。因此,本解析公式在本发明中仅用于对真实颗粒赫兹变形参数的估计值。精确的结果将由本发明提出的迭代技术方案自动训练得出。
实施方式一
将本发明的方法应用于二维单粒径随机颗粒试样(如图2)的标定,试验参数如下:
表1二维单粒径随机排列颗粒模型参数
标定步骤如下:
二维离散元模型:
步骤(3):建立离散元双轴数值试验,将步骤(2)计算出的颗粒接触参数作为初始参数投入模型进行数值计算,数值试验中采用和室内试验相同的围压。
步骤(4):将试样颗粒间的摩擦系数设置为1.0,实施小应变三轴试验数值仿真,直到轴向应变达到0.05%时,停止加载。模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
步骤(5):构建以弹性模量为目标宏观参数的误差函数L:
步骤(6):判断误差函数L的大小是否满足误差允许值。如满足,当前模型所用刚度参数即为标定参数,标定结束;如不满足,继续执行标定步骤。具体过程为:将步骤(4)模拟得到的弹性模量E*代入步骤(5)构建的误差函数,判断误差函数是否小于误差允许值(例如:0.01%)。如满足,则使用的细观参数为模拟该试验的合理估计参数;如不满足,则按照步骤(7)的描述更新参数。
步骤(7):首次更新颗粒剪切模量和泊松比并再次计算误差函数大小。若是为第一次更新颗粒剪切模量和泊松比,按照如下方式更新:
其中,“:=”为赋值符号,η为系数,取值为-0.03。若当前颗粒赫兹接触参数为第二次或者更多次更新时,则按照步骤(8)到步骤(11)进行更新。
步骤(8):计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值。若令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
步骤(9):计算未经偏差修正的一阶和二阶矩估计值。令Mt和Gt分别为第t次迭代时,两个偏导数关于未经偏差修正的一阶和二阶矩估计值按如下方式赋值得到:
Mt=β1Mt-1+(1-β1)gt (43)
步骤(11):基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
其中,ε是为了保持数值稳定性而设置的非常小的常数,取值e-7;α为学习率,取0.01。
步骤(12):采用更新后的刚度参数重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足,继续执行步骤(7)到步骤(12)所述过程,并直到误差函数低于误差允许值。具体过程为:将步骤(7)更新得到的刚度参数进行离散元双轴或者真三轴压缩试验计算,根据新计算得到的弹性模量E*计算误差函数的大小,并判断误差函数是否小于误差容许值。如小于,结束计算,当前模型的颗粒剪切模量和泊松比即为标定参数;如未满足误差要求,则按照步骤(7)到步骤(12)更新参数值并重新计算,直至误差满足要求。
该标定案列的参数迭代过程如表2所示,在第五次迭代后,数值试验的杨氏模量与目标标定模量的误差仅为0.2%左右,此时对应的颗粒剪切模量和数值的颗粒泊松比为:59270227Pa和0.71。注意到:此颗粒泊松比为数值泊松比,而非物理泊松比;数值泊松比包含了对实际颗粒接触作用和颗粒形状等因素在数值上被简化的补偿作用,因此其取值范围与常规物理泊松比并不相同。
表2二维单粒径随机颗粒试样标定过程
实施方式二
将本发明的方法应用于三维单粒径随机颗粒试样(如图3)的标定,试验参数如下:
表3三维单粒径随机排列颗粒模型参数
标定步骤如下:
三维离散元模型:
步骤(3):建立离散元双轴数值试验,将步骤(2)计算出的颗粒接触参数作为初始参数投入模型进行数值计算,数值试验中采用和室内试验相同的围压。
步骤(4):将试样颗粒间的摩擦系数设置为1.0,实施小应变三轴试验数值仿真,直到轴向应变达到0.05%时,停止加载。模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
步骤(5):构建以弹性模量为目标宏观参数的误差函数L:
步骤(6):判断误差函数L的大小是否满足误差允许值。当前误差函数大小为1.68E-05,已满足误差容许要求,当前模型所用刚度参数即为标定参数,标定结束;
该标定案列的参数迭代过程如表4所示,在使用估计解析值进行计算时,数值试验的杨氏模量与目标标定模量的误差仅为0.4%左右,此时对应的颗粒剪切模量和数值的颗粒泊松比为:18151797Pa和0.77。注意到:此颗粒泊松比为数值泊松比,而非物理泊松比;数值泊松比包含了对实际颗粒接触作用和颗粒形状等因素在数值上被简化的补偿作用,因此其取值范围与常规物理泊松比并不相同。
表4三维单粒径随机颗粒试样标定过程
实施方式三
将本发明的方法应用于三维多粒径随机颗粒试样(如图4)的标定,试验参数如下:
表5三维多粒径随机排列颗粒模型参数
标定步骤如下:
三维离散元模型:
步骤(3):建立离散元双轴数值试验,将步骤(2)计算出的颗粒接触参数作为初始参数投入模型进行数值计算,数值试验中采用和室内试验相同的围压。
步骤(4):将试样颗粒间的摩擦系数设置为1.0,实施小应变三轴试验数值仿真,直到轴向应变达到0.05%时,停止加载。模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
步骤(5):构建以弹性模量为目标宏观参数的误差函数L:
步骤(6):判断误差函数L的大小是否满足误差允许值。如满足,当前模型所用刚度参数即为标定参数,标定结束;如不满足,继续执行标定步骤。具体过程为:将步骤(4)模拟得到的弹性模量E*代入步骤(5)构建的误差函数,判断误差函数是否小于误差允许值(例如:0.01%)。如满足,则使用的细观参数为模拟该试验的合理估计参数;如不满足,则按照步骤(7)的描述更新参数。
步骤(7):首次更新颗粒剪切模量和泊松比并再次计算误差函数大小。若是为第一次更新颗粒剪切模量和泊松比,按照如下方式更新:
其中,“:=”为赋值符号,η为系数,取值为0.1。若当前颗粒赫兹接触参数为第二次或者更多次更新时,则按照步骤(8)到步骤(11)进行更新。
步骤(8):计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值。若令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
步骤(9):计算未经偏差修正的一阶和二阶矩估计值。令Mt和Gt分别为第t次迭代时,两个偏导数关于未经偏差修正的一阶和二阶矩估计值按如下方式赋值得到:
Mt=β1Mt-1+(1-β1)gt (37)
步骤(11):基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
其中,ε是为了保持数值稳定性而设置的非常小的常数,取值e-7;α为学习率,取0.01。
步骤(12):采用更新后的刚度参数重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足,继续执行步骤(7)到步骤(12)所述过程,并直到误差函数低于误差允许值。具体过程为:将步骤(7)更新得到的刚度参数进行离散元双轴或者真三轴压缩试验计算,根据新计算得到的弹性模量E*计算误差函数的大小,并判断误差函数是否小于误差容许值。如小于,结束计算,当前模型的颗粒剪切模量和泊松比即为标定参数;如未满足误差要求,则按照步骤(7)到步骤(12)更新参数值并重新计算,直至误差满足要求。
该标定案列的参数迭代过程如表2所示,在第五次迭代后,数值试验的杨氏模量与目标标定模量的误差在1%以内,此时对应的颗粒剪切模量和数值的颗粒泊松比为:15792063Pa和0.67。注意到:此颗粒泊松比为数值泊松比,而非物理泊松比;数值泊松比包含了对实际颗粒接触作用和颗粒形状等因素在数值上被简化的补偿作用,因此其取值范围与常规物理泊松比并不相同。
表6三维多粒径随机颗粒试样标定过程
区别于现有技术,本发明的模拟岩土材料时离散元赫兹接触参数自动标定方法首先将离散元标定的目标宏观参数代入解析公式计算出颗粒剪切模量和泊松比。以此解析结果作为颗粒参数初始估计值,建立数值三轴或双轴试验,得到数值模型的宏观杨氏模量和泊松比。计算误差函数大小,判断是否更新颗粒参数。根据迭代次数的不同,分别采取不同的策略。初次迭代时,通过人工给予一个与初始估计值(解析结果)成正比的较小扰动;二次及更多次迭代时则采用自适应矩估计策略来更新参数。通过本发明,能够自动快速标定离散元模拟过程中颗粒剪切模量和泊松比,极大促进了离散元对岩土工程问题的分析和研究;本发明基于自适应的矩估计策略来迭代标定过程,标定收敛性好,效率高;基于小应变三轴或者双轴加载试验来标定变形参数,模拟时间短且标定结果可靠。本发明提供了一种高效、稳定、实施流程清晰且编程方便的方法,可为标定离散元的非线性赫兹变形参数提供技术支撑,促进离散元技术对岩土工程问题的分析和研究。
上面结合附图对本发明的实施例进行了描述,但是本发明并不局限于上述的具体实施方式,上述的具体实施方式仅仅是示意性的,而不是限制性的,本领域的普通技术人员在本发明的启示下,在不脱离本发明宗旨和权利要求所保护的范围情况下,还可做出很多形式,这些均属于本发明的保护之内。
Claims (9)
1.一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,包括以下步骤:
一、在使用离散元模拟岩土材料前,采用室内三轴或者双轴试验,确定待模拟岩土材料试样的杨氏模量和泊松比;
二、以试验得到的宏观参数为目标参数,根据解析公式计算颗粒剪切模量和泊松比的初始估计值;
三、使用计算出来的颗粒剪切模量和泊松比作为输入参数,建立离散元双轴或者三轴数值试验;
四、在加载之前,将试样颗粒间的摩擦系数设置为1.0或者更大的值,实施小应变三轴/双轴试验数值仿真,直到轴向应变达到预设阈值,停止加载;
五、构建误差函数L;
六、判断误差函数L的大小是否满足误差允许值;如满足,当前模型所用颗粒剪切模量和泊松比即为标定参数,标定结束;如不满足,继续执行标定;
七、首次更新颗粒剪切模量和泊松比并再次评估误差函数大小;
八、计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值;
九、计算未经偏差修正的一阶和二阶矩估计值;
十、修正一阶和二阶矩估计;
十一、基于修正的一阶和二阶矩估计值,计算第二次及更多次迭代步中颗粒剪切模量和泊松比的更新值;
十二、采用更新后的颗粒剪切模量和泊松比重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值;如不满足,继续执行步骤七到步骤十一所述过程,直到误差函数低于误差允许值。
4.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,对于多粒径散体试样,采用试样的中值粒径来计算模型的初始刚度。
5.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,真三轴或者双轴模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010394118.2A CN111610091A (zh) | 2020-05-11 | 2020-05-11 | 一种模拟岩土材料时离散元赫兹接触参数自动标定方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010394118.2A CN111610091A (zh) | 2020-05-11 | 2020-05-11 | 一种模拟岩土材料时离散元赫兹接触参数自动标定方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN111610091A true CN111610091A (zh) | 2020-09-01 |
Family
ID=72201886
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010394118.2A Pending CN111610091A (zh) | 2020-05-11 | 2020-05-11 | 一种模拟岩土材料时离散元赫兹接触参数自动标定方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111610091A (zh) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112098209A (zh) * | 2020-09-15 | 2020-12-18 | 河北工业大学 | 一种岩土颗粒破坏局部化识别方法 |
CN113408124A (zh) * | 2021-06-16 | 2021-09-17 | 河海大学 | 一种不改变边界形状的颗粒体系伺服控制方法 |
CN113887092A (zh) * | 2021-08-31 | 2022-01-04 | 浙江工业大学 | 土工织物包裹泥浆加排水板抽水的离散元模型构建方法 |
CN114542043A (zh) * | 2022-04-28 | 2022-05-27 | 太原理工大学 | 基于压裂液粘度优化改进岩层压裂增渗的方法及装置 |
CN115931548A (zh) * | 2022-10-13 | 2023-04-07 | 荣梓华 | 杨氏模量确定方法、装置和电子设备 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109030202A (zh) * | 2018-06-19 | 2018-12-18 | 湘潭大学 | 一种快速确定岩石类脆性材料离散元模型参数的方法 |
US20190057168A1 (en) * | 2016-03-28 | 2019-02-21 | Baker Hughes, A Ge Company, Llc | Obtaining micro- and macro-rock properties with a calibrated rock deformation simulation |
CN109815599A (zh) * | 2019-01-28 | 2019-05-28 | 南京大学 | 一种离散元材料自动训练方法 |
-
2020
- 2020-05-11 CN CN202010394118.2A patent/CN111610091A/zh active Pending
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20190057168A1 (en) * | 2016-03-28 | 2019-02-21 | Baker Hughes, A Ge Company, Llc | Obtaining micro- and macro-rock properties with a calibrated rock deformation simulation |
CN109030202A (zh) * | 2018-06-19 | 2018-12-18 | 湘潭大学 | 一种快速确定岩石类脆性材料离散元模型参数的方法 |
CN109815599A (zh) * | 2019-01-28 | 2019-05-28 | 南京大学 | 一种离散元材料自动训练方法 |
Non-Patent Citations (1)
Title |
---|
TONGMING QU 等: "A hybrid calibration approach to Hertz-type contact parameters for discrete element models", 《INTERNATIONAL JOURNAL FOR NUMERICAL AND ANALYTICAL METHODS IN GEOMECHANICS》 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112098209A (zh) * | 2020-09-15 | 2020-12-18 | 河北工业大学 | 一种岩土颗粒破坏局部化识别方法 |
CN113408124A (zh) * | 2021-06-16 | 2021-09-17 | 河海大学 | 一种不改变边界形状的颗粒体系伺服控制方法 |
CN113408124B (zh) * | 2021-06-16 | 2023-08-22 | 河海大学 | 一种不改变边界形状的颗粒体系伺服控制方法 |
CN113887092A (zh) * | 2021-08-31 | 2022-01-04 | 浙江工业大学 | 土工织物包裹泥浆加排水板抽水的离散元模型构建方法 |
CN113887092B (zh) * | 2021-08-31 | 2023-03-31 | 浙江工业大学 | 土工织物包裹泥浆加排水板抽水的离散元模型构建方法 |
CN114542043A (zh) * | 2022-04-28 | 2022-05-27 | 太原理工大学 | 基于压裂液粘度优化改进岩层压裂增渗的方法及装置 |
CN114542043B (zh) * | 2022-04-28 | 2022-08-12 | 太原理工大学 | 基于压裂液粘度优化改进岩层压裂增渗的方法及装置 |
CN115931548A (zh) * | 2022-10-13 | 2023-04-07 | 荣梓华 | 杨氏模量确定方法、装置和电子设备 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111610091A (zh) | 一种模拟岩土材料时离散元赫兹接触参数自动标定方法 | |
CN111611695A (zh) | 一种模拟岩土材料时离散元线性刚度参数的自动标定方法 | |
CN107729706B (zh) | 一种非线性机械系统的动力学模型构建方法 | |
Leblicq et al. | A discrete element approach for modelling the compression of crop stems | |
CN108959686A (zh) | 一种基于灵敏度分析的有限元模型修正方法 | |
CN110631792B (zh) | 基于卷积神经网络的抗震混合试验模型更新方法 | |
CN108984829B (zh) | 堆石混凝土堆石体堆积过程的计算方法和系统 | |
CN111666621B (zh) | 粘土地层隧道开挖面安全支护压力区间确定方法 | |
CN113536623B (zh) | 一种材料不确定性结构稳健性拓扑优化设计方法 | |
Zhang et al. | A continuous contact force model for impact analysis | |
Goodarzi et al. | Modelling slope failure using a quasi-static MPM with a non-local strain softening approach | |
Mousavi et al. | Imperfection sensitivity analysis of double domes free form space structures | |
CN104899642B (zh) | 基于混合神经网络模型的大变形柔性体动态应力补偿方法 | |
CN113468466A (zh) | 基于神经网络的多工况一维波动方程求解方法 | |
CN111158059B (zh) | 基于三次b样条函数的重力反演方法 | |
CN108446507A (zh) | 基于网格质量反馈优化的弹性体网格变形方法 | |
CN108920737B (zh) | 一种水动力模型的粒子滤波同化方法、装置和计算设备 | |
CN109885896B (zh) | 一种基于复变差分灵敏度的非线性结构有限元模型修正方法 | |
CN108491636B (zh) | 基于几何约束的弹性体网格变形方法 | |
Moosavi et al. | Isogeometric meshless finite volume method in nonlinear elasticity | |
Jovanović | Identification of dynamic system using neural network | |
Dabeet et al. | Simulation of cyclic direct simple shear loading response of soils using discrete element modeling | |
CN111159946B (zh) | 一种基于最小势能原理的非连续性问题分区求解方法 | |
CN101634618B (zh) | 一种三维弹性模量成像方法 | |
CN114372387A (zh) | 镍基单晶合金微观组织演化相场模型的有限元法 |
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 |
Application publication date: 20200901 |
|
WD01 | Invention patent application deemed withdrawn after publication |