CN111610091A - 一种模拟岩土材料时离散元赫兹接触参数自动标定方法 - Google Patents

一种模拟岩土材料时离散元赫兹接触参数自动标定方法 Download PDF

Info

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
Application number
CN202010394118.2A
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 CN202010394118.2A priority Critical patent/CN111610091A/zh
Publication of CN111610091A publication Critical patent/CN111610091A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N3/00Investigating strength properties of solid materials by application of mechanical stress
    • G01N3/08Investigating strength properties of solid materials by application of mechanical stress by applying steady tensile or compressive forces
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/0058Kind of property studied
    • G01N2203/0069Fatigue, creep, strain-stress relations or elastic constants
    • G01N2203/0075Strain-stress relations or elastic constants
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/02Details not specific for a particular testing method
    • G01N2203/025Geometry of the test
    • G01N2203/0254Biaxial, the forces being applied along two normal axes of the specimen
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/02Details not specific for a particular testing method
    • G01N2203/025Geometry of the test
    • G01N2203/0256Triaxial, i.e. the forces being applied along three normal axes of the specimen
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force 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关于颗粒剪切模量和泊松比变化的梯度值;
九、计算未经偏差修正的一阶和二阶矩估计值;
十、修正一阶和二阶矩估计;
十一、基于修正的一阶和二阶矩估计值,计算第二次及更多次迭代步中颗粒剪切模量和泊松比的更新值;
十二、采用更新后的颗粒剪切模量和泊松比重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值;如不满足,继续执行步骤七到步骤十一所述过程,直到误差函数低于误差允许值。
其中,颗粒接触剪切模量和泊松比的初始估计值
Figure BDA0002486734950000029
Figure BDA00024867349500000210
按如下方式更新:
二维离散元模型:
Figure BDA0002486734950000021
Figure BDA0002486734950000022
三维离散元模型:
Figure BDA0002486734950000023
Figure BDA0002486734950000024
其中,
Figure BDA0002486734950000025
Figure BDA0002486734950000026
分别为试验得到的小应变杨氏模量和泊松比,φ为试样中的孔隙率;
Figure BDA0002486734950000027
为颗粒试样的平均配位数;当试样单粒径分布时,r为颗粒半径;当试样为多粒径分布时,r为试样的中值粒径。
其中,采用生成颗粒时的初始摩擦系数fc来估计最终试样的孔隙率和平均配位数经验值,如下:
Figure BDA0002486734950000028
其中,对于多粒径散体试样,采用试样的中值粒径来计算模型的初始接触剪切模量和泊松比。
其中,真三轴或者双轴模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
其中,误差函数L为:
Figure BDA0002486734950000031
其中,E*是将颗粒间摩擦系数设置为1.0时,使试样在整体发生较小轴向变形时,颗粒间不发生滑动变形,计算数值试样模拟得到的弹性模量。
其中,若颗粒剪切模量G和泊松比ν为第一次迭代更新,按照如下方式更新:
Figure BDA0002486734950000032
其中,“:=”为赋值符号,η为系数,取值在-0.5-0.5之间。
其中,若当前颗粒接触刚度参数为第二次或者更多次更新时,需按照差分方法计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值gG和gυ,令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
Figure BDA0002486734950000033
其中,若当前颗粒接触刚度参数为第二次或者更多次更新时,基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
Figure BDA0002486734950000034
其中,
Figure BDA0002486734950000035
Figure BDA0002486734950000036
分别为第t次迭代时,经偏差修正后的一阶和二阶矩估计值;ε是为了保持数值稳定性而设置的非常小的常数,一般取值e-7;α为学习率。
区别于现有技术,本发明的模拟岩土材料时离散元赫兹接触参数自动标定方法首先将离散元标定的目标宏观参数代入解析公式计算出颗粒剪切模量和泊松比。以此解析结果作为颗粒参数初始估计值,建立数值三轴或双轴试验,得到数值模型的宏观杨氏模量和泊松比。计算误差函数大小,判断是否更新颗粒参数。根据迭代次数的不同,分别采取不同的策略。初次迭代时,通过人工给予一个与初始估计值(解析结果)成正比的较小扰动;二次及更多次迭代时则采用自适应矩估计策略来更新参数。通过本发明,能够自动快速标定离散元模拟过程中颗粒剪切模量和泊松比,极大促进了离散元对岩土工程问题的分析和研究;本发明基于自适应的矩估计策略来迭代标定过程,标定收敛性好,效率高;基于小应变三轴或者双轴加载试验来标定变形参数,模拟时间短且标定结果可靠。本发明提供了一种高效、稳定、实施流程清晰且编程方便的方法,可为标定离散元的非线性赫兹变形参数提供技术支撑,促进离散元技术对岩土工程问题的分析和研究。
附图说明
附图1是所述的离散元非线性赫兹接触变形参数自动标定技术路线图;
附图2是所述的二维单粒径随机颗粒试样模型图;
附图3是所述的三维单粒径随机颗粒试样模型图;
附图4是所述的三维多粒径随机颗粒试样模型图;
具体实施方式
为了对本发明的技术特征、目的和效果有更加清楚的理解,现对照附图详细说明本发明的具体实施方式。
参阅图1,本发明提供了一种模拟岩土材料时离散元赫兹接触参数自动标定方法,包括步骤:
步骤(1):在使用离散元模拟岩土材料前,需采用室内真三轴或者双轴试验等方式,确定待模拟岩土材料的宏观参数。所用宏观参数为小应变杨氏模量
Figure BDA0002486734950000041
与泊松比
Figure BDA0002486734950000042
步骤(2):将试验得到的目标宏观参数输入以下公式,分别确定二维和三维模拟情况下,颗粒接触剪切模量和泊松比的初始估计值
Figure BDA0002486734950000043
Figure BDA0002486734950000044
二维离散元模型:
Figure BDA0002486734950000045
Figure BDA0002486734950000046
三维离散元模型:
Figure BDA0002486734950000047
Figure BDA0002486734950000048
其中,φ为试样中的孔隙率;
Figure BDA0002486734950000049
为颗粒试样的平均配位数;当试样单粒径分布时,r为颗粒半径;当试样为多粒径分布时,r为试样的中值粒径。
对于离散元散体试样,生成颗粒时的初始摩擦系数往往对颗粒体的最终密实程度有很大影响。当初始摩擦系数为fc时,最终试样的孔隙率和平均配位数的值可按下列经验公式估计:
Figure BDA00024867349500000410
其中,2D和3D分别表示2维和3维模型;
步骤(3):建立离散元双轴或者真三轴数值试验,将步骤(2)计算出的颗粒的剪切模量和泊松比作为初始参数投入模型进行数值计算,数值试验中采用和室内试验相同的围压。
步骤(4):将试样颗粒间的摩擦系数设置为1.0,实施小应变三轴试验数值仿真,直到轴向应变达到一个较小值ε时,停止加载。(例如:轴向应变ε=0.05%)。模拟结束后,应检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。将颗粒间摩擦系数设置为1.0的目的是使试样在整体发生较小轴向变形时,颗粒间不发生滑动变形(塑性变形)。计算数值试样模拟得到的弹性模量E*和泊松比ν*。计算方法为:
Figure BDA0002486734950000051
Figure BDA0002486734950000052
其中,Δσ为偏应力,Δε为偏应变,下标x、y、z分别表示笛卡尔坐标系下的三个坐标方向。
步骤(5):构建以弹性模量为目标宏观参数的误差函数L:
Figure BDA0002486734950000053
步骤(6):判断误差函数L的大小是否满足误差允许值。如满足,当前模型所用刚度参数即为标定参数,标定结束;如不满足,继续执行标定步骤。具体过程为:将步骤(4)模拟得到的弹性模量E*代入步骤(5)构建的误差函数,判断误差函数是否小于误差允许值(例如:0.01%)。如满足,则使用的细观参数为模拟该试验的合理估计参数;如不满足,则按照步骤(7)的描述更新参数。
步骤(7):首次更新颗粒剪切模量和泊松比并再次计算误差函数大小。若是为第一次更新颗粒剪切模量和泊松比,按照如下方式更新:
Figure BDA0002486734950000054
其中,“:=”为赋值符号,η为系数,一般取值为-0.5-0.5之间。若当前颗粒赫兹接触参数为第二次或者更多次更新时,则按照步骤(8)到步骤(11)进行更新。
步骤(8):计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值。若令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
Figure BDA0002486734950000055
步骤(9):计算未经偏差修正的一阶和二阶矩估计值。令Mt和Gt分别为第t次迭代时,两个偏导数关于未经偏差修正的一阶和二阶矩估计值按如下方式赋值得到:
Mt=β1Mt-1+(1-β1)gt (13)
Figure BDA0002486734950000056
其中:gt表示在第t次迭代时误差函数的梯度值,
Figure BDA0002486734950000057
代表各个梯度分量依次相乘,β1和β2一阶矩估计Mt和二阶矩估计Gt的衰减率,通常取值为0.9和0.99。
步骤(10):修正一阶和二阶矩估计。
Figure BDA0002486734950000061
Figure BDA0002486734950000062
分别为第t次迭代时,经偏差修正后的一阶和二阶矩估计值,分别按如下方式计算:
Figure BDA0002486734950000063
Figure BDA0002486734950000064
步骤(11):基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
Figure BDA0002486734950000065
其中,ε是为了保持数值稳定性而设置的非常小的常数,一般取值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之间的关系:
Figure BDA0002486734950000066
其中,Nc为在体积V内的参与荷载传递的接触数目;c表示该范围内的第c个接触。c=1,2,…Nc
(2)基于均匀应力假设的宏细观参数关系
假设颗粒试样内部应力均匀分布,则颗粒间的接触力可根据如下公式通过颗粒试样中的平均应力得出:
Figure BDA0002486734950000067
其中,Ajk是一个与颗粒内部结构相关的张量,
Figure BDA0002486734950000068
为第c个接触平面的单位外法线方向矢量在第k个基向量方向的分量。
将公式(17)代入公式(18),得到:
Figure BDA0002486734950000071
其中,δij为克罗内克函数。
根据公式(19),可求出张量Ajk的表达式:
Figure BDA0002486734950000072
(3)颗粒试样的应变描述
根据能量守恒原理,在不考虑颗粒体力的情况下,外力对试样的做功等于颗粒接触变形积累的应变能,即为:
Figure BDA0002486734950000073
通过将公式(18)代入公式(21),可得:
Figure BDA0002486734950000074
由于应力张量Δσij可以为物理材料条件允许下的任意非零值,因此可以得到颗粒试样应变的表达式为:
Figure BDA0002486734950000075
(4)本构方程与柔度张量
Figure BDA0002486734950000076
是第c个接触的柔度,则颗粒间接触的重叠量
Figure BDA0002486734950000077
与接触力fi c的关系可表示为:
Figure BDA0002486734950000078
根据弹性力学原理,固体的应力应变关系可用柔度张量Sijkl的形式表达:
εij=Sijklσkl (25)
结合公式(18),(23),(24)和(25),可计算得到该颗粒材料的应力应变柔度张量的大小为:
Figure BDA0002486734950000079
(5)柔度张量的积分形式与应力应变关系
在颗粒材料中颗粒的大小相同时,颗粒接触支向量Lj c可表达为:
Figure BDA00024867349500000710
其中,r为颗粒的半径。
将公式(24)和公式(27)代入公式(26)。另外,对于一个包含大量颗粒数目的试样,接触法线的方向在统计上可被近似为各向均匀的。根据接触法向方向均匀分布的假设,柔度张量Sijkl的积分形式可计算得:
Figure BDA0002486734950000081
通过对公式(28)积分,可得到散体的应力应变关系为:
2D模型:
Figure BDA0002486734950000082
其中:
Figure BDA0002486734950000083
3D模型:
Figure BDA0002486734950000084
其中
Figure BDA0002486734950000085
对于一般各向同性的弹性固体,其应力应变关系为:
2D平面应力模型:
Figure BDA0002486734950000086
3D模型:
Figure BDA0002486734950000087
其中,
Figure BDA0002486734950000088
Figure BDA0002486734950000089
为一般各向同性弹性材料的宏观弹性模量与泊松比。
通过对比公式(29)和公式(31),对比公式(30)和公式(32),可得到一般均匀散体的宏观弹性常数与细观参数的关系为:
2D平面应力模型:
Figure BDA00024867349500000810
3D模型:
Figure BDA00024867349500000811
其中,
Figure BDA00024867349500000812
为一般各向同性弹性材料的宏观剪切模量。
(6)基于赫兹接触模型的散体试样小应变模量
对于承受均匀正压力的各向同性材料,可认为在水平和竖直方向,材料中只有主应力,没有剪切应力。基于颗粒间无剪切力,且各个接触法向作用力的大小相等的假设,公式(18)的全量形式可写为:
Figure BDA0002486734950000091
其中,Nc为颗粒试样中参与荷载传递的有效接触数目;σc为颗粒试样承受的正应力大小。
考虑到赫兹接触模型的接触力与接触位移的作用关系为:
Figure BDA0002486734950000092
Figure BDA0002486734950000093
其中,G和ν是细观尺度上的颗粒剪切模量与泊松比。
结合公式(35),(36)和(37),可得到接触刚度的表达形式为:
Figure BDA0002486734950000094
将公式(38)代入公式(33)和(34),可得:
2D模型:
Figure BDA0002486734950000095
Figure BDA0002486734950000096
3D模型:
Figure BDA0002486734950000097
Figure BDA0002486734950000098
联立公式(39)和(40),当变形参数的标定目标为离散元试样的小应变杨氏模量
Figure BDA0002486734950000099
与泊松比
Figure BDA00024867349500000910
时,可得:
2D模型:
Figure BDA00024867349500000911
Figure BDA00024867349500000912
3D模型:
Figure BDA0002486734950000101
Figure BDA0002486734950000102
此公式基于各向同性均匀单粒径颗粒系统推导得出,对于具有多粒径颗粒系统,其内部组构不可避免存在一定程度的各向异性特征。因此,本解析公式在本发明中仅用于对真实颗粒赫兹变形参数的估计值。精确的结果将由本发明提出的迭代技术方案自动训练得出。
实施方式一
将本发明的方法应用于二维单粒径随机颗粒试样(如图2)的标定,试验参数如下:
Figure BDA0002486734950000103
表1二维单粒径随机排列颗粒模型参数
标定步骤如下:
步骤(1):通过室内试验确定岩土材料宏观参数。本实施例将以小应变杨氏模量
Figure BDA0002486734950000104
为10Gpa,泊松比
Figure BDA0002486734950000105
为0.2为标定目标。
步骤(2):将目标宏观参数输入以下公式,在二维模拟情况下,颗粒接触剪切模量和泊松比的初始估计值
Figure BDA0002486734950000106
Figure BDA0002486734950000107
二维离散元模型:
Figure BDA0002486734950000108
Figure BDA0002486734950000109
其中,φ为试样中的孔隙率;
Figure BDA00024867349500001010
为颗粒试样的平均配位数;当试样单粒径分布时,r为颗粒半径;当试样为多粒径分布时,r为试样的中值粒径。
步骤(3):建立离散元双轴数值试验,将步骤(2)计算出的颗粒接触参数作为初始参数投入模型进行数值计算,数值试验中采用和室内试验相同的围压。
步骤(4):将试样颗粒间的摩擦系数设置为1.0,实施小应变三轴试验数值仿真,直到轴向应变达到0.05%时,停止加载。模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
步骤(5):构建以弹性模量为目标宏观参数的误差函数L:
Figure BDA00024867349500001011
步骤(6):判断误差函数L的大小是否满足误差允许值。如满足,当前模型所用刚度参数即为标定参数,标定结束;如不满足,继续执行标定步骤。具体过程为:将步骤(4)模拟得到的弹性模量E*代入步骤(5)构建的误差函数,判断误差函数是否小于误差允许值(例如:0.01%)。如满足,则使用的细观参数为模拟该试验的合理估计参数;如不满足,则按照步骤(7)的描述更新参数。
步骤(7):首次更新颗粒剪切模量和泊松比并再次计算误差函数大小。若是为第一次更新颗粒剪切模量和泊松比,按照如下方式更新:
Figure BDA0002486734950000111
其中,“:=”为赋值符号,η为系数,取值为-0.03。若当前颗粒赫兹接触参数为第二次或者更多次更新时,则按照步骤(8)到步骤(11)进行更新。
步骤(8):计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值。若令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
Figure BDA0002486734950000112
步骤(9):计算未经偏差修正的一阶和二阶矩估计值。令Mt和Gt分别为第t次迭代时,两个偏导数关于未经偏差修正的一阶和二阶矩估计值按如下方式赋值得到:
Mt=β1Mt-1+(1-β1)gt (43)
Figure BDA0002486734950000113
其中:gt表示在第t次迭代时误差函数的梯度值,
Figure BDA0002486734950000114
代表各个梯度分量依次相乘,β1和β2一阶矩估计Mt和二阶矩估计Gt的衰减率,分别取值为0.9和0.99。
步骤(10):修正一阶和二阶矩估计。
Figure BDA0002486734950000115
Figure BDA0002486734950000116
分别为第t次迭代时,经偏差修正后的一阶和二阶矩估计值,分别按如下方式计算:
Figure BDA0002486734950000117
Figure BDA0002486734950000118
步骤(11):基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
Figure BDA0002486734950000119
其中,ε是为了保持数值稳定性而设置的非常小的常数,取值e-7;α为学习率,取0.01。
步骤(12):采用更新后的刚度参数重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足,继续执行步骤(7)到步骤(12)所述过程,并直到误差函数低于误差允许值。具体过程为:将步骤(7)更新得到的刚度参数进行离散元双轴或者真三轴压缩试验计算,根据新计算得到的弹性模量E*计算误差函数的大小,并判断误差函数是否小于误差容许值。如小于,结束计算,当前模型的颗粒剪切模量和泊松比即为标定参数;如未满足误差要求,则按照步骤(7)到步骤(12)更新参数值并重新计算,直至误差满足要求。
该标定案列的参数迭代过程如表2所示,在第五次迭代后,数值试验的杨氏模量与目标标定模量的误差仅为0.2%左右,此时对应的颗粒剪切模量和数值的颗粒泊松比为:59270227Pa和0.71。注意到:此颗粒泊松比为数值泊松比,而非物理泊松比;数值泊松比包含了对实际颗粒接触作用和颗粒形状等因素在数值上被简化的补偿作用,因此其取值范围与常规物理泊松比并不相同。
Figure BDA0002486734950000121
表2二维单粒径随机颗粒试样标定过程
实施方式二
将本发明的方法应用于三维单粒径随机颗粒试样(如图3)的标定,试验参数如下:
Figure BDA0002486734950000122
表3三维单粒径随机排列颗粒模型参数
标定步骤如下:
步骤(1):通过室内试验确定岩土材料宏观参数。本实施例将以小应变杨氏模量
Figure BDA0002486734950000123
为10Gpa,泊松比
Figure BDA0002486734950000124
为0.2为标定目标。
步骤(2):将目标宏观参数输入以下公式,在二维模拟情况下,颗粒接触剪切模量和泊松比的初始估计值
Figure BDA0002486734950000125
Figure BDA0002486734950000126
三维离散元模型:
Figure BDA0002486734950000127
Figure BDA0002486734950000128
其中,φ为试样中的孔隙率;
Figure BDA0002486734950000129
为颗粒试样的平均配位数;当试样单粒径分布时,r为颗粒半径;当试样为多粒径分布时,r为试样的中值粒径。
步骤(3):建立离散元双轴数值试验,将步骤(2)计算出的颗粒接触参数作为初始参数投入模型进行数值计算,数值试验中采用和室内试验相同的围压。
步骤(4):将试样颗粒间的摩擦系数设置为1.0,实施小应变三轴试验数值仿真,直到轴向应变达到0.05%时,停止加载。模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
步骤(5):构建以弹性模量为目标宏观参数的误差函数L:
Figure BDA0002486734950000131
步骤(6):判断误差函数L的大小是否满足误差允许值。当前误差函数大小为1.68E-05,已满足误差容许要求,当前模型所用刚度参数即为标定参数,标定结束;
该标定案列的参数迭代过程如表4所示,在使用估计解析值进行计算时,数值试验的杨氏模量与目标标定模量的误差仅为0.4%左右,此时对应的颗粒剪切模量和数值的颗粒泊松比为:18151797Pa和0.77。注意到:此颗粒泊松比为数值泊松比,而非物理泊松比;数值泊松比包含了对实际颗粒接触作用和颗粒形状等因素在数值上被简化的补偿作用,因此其取值范围与常规物理泊松比并不相同。
Figure BDA0002486734950000132
表4三维单粒径随机颗粒试样标定过程
实施方式三
将本发明的方法应用于三维多粒径随机颗粒试样(如图4)的标定,试验参数如下:
Figure BDA0002486734950000133
表5三维多粒径随机排列颗粒模型参数
标定步骤如下:
步骤(1):通过室内试验确定岩土材料宏观参数。本实施例将以小应变杨氏模量
Figure BDA0002486734950000134
为10Gpa,泊松比
Figure BDA0002486734950000135
为0.2为标定目标。
步骤(2):将目标宏观参数输入以下公式,在二维模拟情况下,颗粒接触剪切模量和泊松比的初始估计值
Figure BDA0002486734950000136
Figure BDA0002486734950000137
三维离散元模型:
Figure BDA0002486734950000138
Figure BDA0002486734950000139
其中,φ为试样中的孔隙率;
Figure BDA00024867349500001310
为颗粒试样的平均配位数;当试样单粒径分布时,r为颗粒半径;当试样为多粒径分布时,r为试样的中值粒径。
步骤(3):建立离散元双轴数值试验,将步骤(2)计算出的颗粒接触参数作为初始参数投入模型进行数值计算,数值试验中采用和室内试验相同的围压。
步骤(4):将试样颗粒间的摩擦系数设置为1.0,实施小应变三轴试验数值仿真,直到轴向应变达到0.05%时,停止加载。模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
步骤(5):构建以弹性模量为目标宏观参数的误差函数L:
Figure BDA0002486734950000141
步骤(6):判断误差函数L的大小是否满足误差允许值。如满足,当前模型所用刚度参数即为标定参数,标定结束;如不满足,继续执行标定步骤。具体过程为:将步骤(4)模拟得到的弹性模量E*代入步骤(5)构建的误差函数,判断误差函数是否小于误差允许值(例如:0.01%)。如满足,则使用的细观参数为模拟该试验的合理估计参数;如不满足,则按照步骤(7)的描述更新参数。
步骤(7):首次更新颗粒剪切模量和泊松比并再次计算误差函数大小。若是为第一次更新颗粒剪切模量和泊松比,按照如下方式更新:
Figure BDA0002486734950000142
其中,“:=”为赋值符号,η为系数,取值为0.1。若当前颗粒赫兹接触参数为第二次或者更多次更新时,则按照步骤(8)到步骤(11)进行更新。
步骤(8):计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值。若令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
Figure BDA0002486734950000143
步骤(9):计算未经偏差修正的一阶和二阶矩估计值。令Mt和Gt分别为第t次迭代时,两个偏导数关于未经偏差修正的一阶和二阶矩估计值按如下方式赋值得到:
Mt=β1Mt-1+(1-β1)gt (37)
Figure BDA0002486734950000144
其中:gt表示在第t次迭代时误差函数的梯度值,
Figure BDA0002486734950000145
代表各个梯度分量依次相乘,β1和β2一阶矩估计Mt和二阶矩估计Gt的衰减率,分别取值为0.9和0.99。
步骤(10):修正一阶和二阶矩估计。
Figure BDA0002486734950000146
Figure BDA0002486734950000147
分别为第t次迭代时,经偏差修正后的一阶和二阶矩估计值,分别按如下方式计算:
Figure BDA0002486734950000148
Figure BDA0002486734950000149
步骤(11):基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
Figure BDA0002486734950000151
其中,ε是为了保持数值稳定性而设置的非常小的常数,取值e-7;α为学习率,取0.01。
步骤(12):采用更新后的刚度参数重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值。如不满足,继续执行步骤(7)到步骤(12)所述过程,并直到误差函数低于误差允许值。具体过程为:将步骤(7)更新得到的刚度参数进行离散元双轴或者真三轴压缩试验计算,根据新计算得到的弹性模量E*计算误差函数的大小,并判断误差函数是否小于误差容许值。如小于,结束计算,当前模型的颗粒剪切模量和泊松比即为标定参数;如未满足误差要求,则按照步骤(7)到步骤(12)更新参数值并重新计算,直至误差满足要求。
该标定案列的参数迭代过程如表2所示,在第五次迭代后,数值试验的杨氏模量与目标标定模量的误差在1%以内,此时对应的颗粒剪切模量和数值的颗粒泊松比为:15792063Pa和0.67。注意到:此颗粒泊松比为数值泊松比,而非物理泊松比;数值泊松比包含了对实际颗粒接触作用和颗粒形状等因素在数值上被简化的补偿作用,因此其取值范围与常规物理泊松比并不相同。
Figure BDA0002486734950000152
表6三维多粒径随机颗粒试样标定过程
区别于现有技术,本发明的模拟岩土材料时离散元赫兹接触参数自动标定方法首先将离散元标定的目标宏观参数代入解析公式计算出颗粒剪切模量和泊松比。以此解析结果作为颗粒参数初始估计值,建立数值三轴或双轴试验,得到数值模型的宏观杨氏模量和泊松比。计算误差函数大小,判断是否更新颗粒参数。根据迭代次数的不同,分别采取不同的策略。初次迭代时,通过人工给予一个与初始估计值(解析结果)成正比的较小扰动;二次及更多次迭代时则采用自适应矩估计策略来更新参数。通过本发明,能够自动快速标定离散元模拟过程中颗粒剪切模量和泊松比,极大促进了离散元对岩土工程问题的分析和研究;本发明基于自适应的矩估计策略来迭代标定过程,标定收敛性好,效率高;基于小应变三轴或者双轴加载试验来标定变形参数,模拟时间短且标定结果可靠。本发明提供了一种高效、稳定、实施流程清晰且编程方便的方法,可为标定离散元的非线性赫兹变形参数提供技术支撑,促进离散元技术对岩土工程问题的分析和研究。
上面结合附图对本发明的实施例进行了描述,但是本发明并不局限于上述的具体实施方式,上述的具体实施方式仅仅是示意性的,而不是限制性的,本领域的普通技术人员在本发明的启示下,在不脱离本发明宗旨和权利要求所保护的范围情况下,还可做出很多形式,这些均属于本发明的保护之内。

Claims (9)

1.一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,包括以下步骤:
一、在使用离散元模拟岩土材料前,采用室内三轴或者双轴试验,确定待模拟岩土材料试样的杨氏模量和泊松比;
二、以试验得到的宏观参数为目标参数,根据解析公式计算颗粒剪切模量和泊松比的初始估计值;
三、使用计算出来的颗粒剪切模量和泊松比作为输入参数,建立离散元双轴或者三轴数值试验;
四、在加载之前,将试样颗粒间的摩擦系数设置为1.0或者更大的值,实施小应变三轴/双轴试验数值仿真,直到轴向应变达到预设阈值,停止加载;
五、构建误差函数L;
六、判断误差函数L的大小是否满足误差允许值;如满足,当前模型所用颗粒剪切模量和泊松比即为标定参数,标定结束;如不满足,继续执行标定;
七、首次更新颗粒剪切模量和泊松比并再次评估误差函数大小;
八、计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值;
九、计算未经偏差修正的一阶和二阶矩估计值;
十、修正一阶和二阶矩估计;
十一、基于修正的一阶和二阶矩估计值,计算第二次及更多次迭代步中颗粒剪切模量和泊松比的更新值;
十二、采用更新后的颗粒剪切模量和泊松比重新计算数值模型,更新误差函数大小并评估误差函数大小是否满足误差允许值;如不满足,继续执行步骤七到步骤十一所述过程,直到误差函数低于误差允许值。
2.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,颗粒接触剪切模量和泊松比的初始估计值
Figure FDA0002486734940000011
Figure FDA0002486734940000012
按如下方式更新:
二维离散元模型:
Figure FDA0002486734940000013
Figure FDA0002486734940000014
三维离散元模型:
Figure FDA0002486734940000015
Figure FDA0002486734940000016
其中,
Figure FDA0002486734940000017
Figure FDA0002486734940000018
分别为试验得到的小应变杨氏模量和泊松比,φ为试样中的孔隙率;
Figure FDA0002486734940000019
为颗粒试样的平均配位数;当试样单粒径分布时,r为颗粒半径;当试样为多粒径分布时,r为试样的中值粒径。
3.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,采用生成颗粒时的初始摩擦系数fc来估计最终试样的孔隙率和平均配位数经验值,如下:
Figure FDA0002486734940000021
4.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,对于多粒径散体试样,采用试样的中值粒径来计算模型的初始刚度。
5.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,真三轴或者双轴模拟结束后,检查试样加载过程中的应力应变曲线,确保试样在此加载期间仅发生可恢复的弹性变形。
6.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,误差函数L为:
Figure FDA0002486734940000022
其中,E*是将颗粒间摩擦系数设置为1.0时,使试样在整体发生较小轴向变形时,颗粒间不发生滑动变形,计算数值试样模拟得到的弹性模量。
7.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,若颗粒剪切模量G和泊松比ν为第一次迭代更新,按照如下方式更新:
Figure FDA0002486734940000023
其中,“:=”为赋值符号,η为系数,取值在-0.5-0.5之间。
8.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,若当前颗粒接触刚度参数为第二次或者更多次更新时,需按照差分方法计算误差函数L关于颗粒剪切模量和泊松比变化的梯度值gG和gυ,令以G和υ为自变量的误差函数为f,即令:
f(G,ν)=L(E*) (8)
则两个偏导数可按照差分方式进行估计得到:
Figure FDA0002486734940000024
9.根据权利要求1所述的一种模拟岩土材料时离散元赫兹接触参数自动标定方法,其特征在于,若当前颗粒接触刚度参数为第二次或者更多次更新时,基于修正的一阶和二阶矩估计值,计算每次迭代步中颗粒剪切模量和泊松比的更新值,更新公式如下:
Figure FDA0002486734940000025
其中,
Figure FDA0002486734940000031
Figure FDA0002486734940000032
分别为第t次迭代时,经偏差修正后的一阶和二阶矩估计值;ε是为了保持数值稳定性而设置的非常小的常数,一般取值e-7;α为学习率。
CN202010394118.2A 2020-05-11 2020-05-11 一种模拟岩土材料时离散元赫兹接触参数自动标定方法 Pending CN111610091A (zh)

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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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 南京大学 一种离散元材料自动训练方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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