CN113536549B - 一种颗粒流微观力学参数反演方法 - Google Patents

一种颗粒流微观力学参数反演方法 Download PDF

Info

Publication number
CN113536549B
CN113536549B CN202110725094.9A CN202110725094A CN113536549B CN 113536549 B CN113536549 B CN 113536549B CN 202110725094 A CN202110725094 A CN 202110725094A CN 113536549 B CN113536549 B CN 113536549B
Authority
CN
China
Prior art keywords
parameters
micromechanics
judge
max
numerical
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.)
Active
Application number
CN202110725094.9A
Other languages
English (en)
Other versions
CN113536549A (zh
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.)
Hunan University of Science and Technology
Original Assignee
Hunan University of Science and 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 Hunan University of Science and Technology filed Critical Hunan University of Science and Technology
Priority to CN202110725094.9A priority Critical patent/CN113536549B/zh
Publication of CN113536549A publication Critical patent/CN113536549A/zh
Application granted granted Critical
Publication of CN113536549B publication Critical patent/CN113536549B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M10/00—Hydrodynamic testing; Arrangements in or on ship-testing tanks or water tunnels
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Chemical & Material Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Dispersion Chemistry (AREA)
  • Computer Hardware Design (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Fluid Mechanics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
  • Developing Agents For Electrophotography (AREA)

Abstract

本发明公开了一种颗粒流微观力学参数反演方法，首先通过物理实验获取岩石宏观力学参数，然后基于物理宏观力学参数反演出颗粒流微观力学参数，接着通过数值模拟计算获取数值模拟宏观力学参数，最后采用模拟退火算法调整颗粒流微观力学参数，使得数值模拟宏观力学参数与物理实验宏观力学参数不断接近，当数值模拟宏观力学参数与物理实验宏观力学参数之间的误差小于10%时，则此时对应的颗粒流微观力学参数则是所需要确定的微观力学参数。本发明操作简便，微观力学参数调整过程中无主观性，收敛速度较好，适用于颗粒流数值模拟计算领域。

Description

一种颗粒流微观力学参数反演方法
技术领域
本发明涉及一种颗粒流微观力学参数反演方法。
背景技术
在开展颗粒流数值模拟计算之前需要确定颗粒流数值模拟模型相应的微观力学参数，该微观力学参数不能随意选取，为了使所选取的微观力学参数对应的数学模型与工程实践更加贴切，一般通过物理试验获取岩石的相关力学参数：单轴抗压强度、弹性模量及泊松比，然后调节颗粒流微观力学参数，使其数值模拟模型对应的宏观力学参数：单轴抗压强度、弹性模型、泊松比与物理试验宏观力学参数靠近，当数值模拟实验得到的宏观力学参数与物理试验宏观力学参数之间误差较小时，则此时对应的微观力学参数即是所求的微观力学参数。目前一般采用试错法不断调整颗粒流微观力学参数使数值模拟计算的宏观力学参数与物理实验宏观力学参数不断接近，当其误差量小于某一值时，则停止调试。采用试错法主要依赖于试验者对颗粒流软件掌握的熟练程度，且调试过程中存在一定的盲目性，同时获取好的微观力学参数具有较大的偶然性。
发明内容
为了解决上述技术问题，本发明提供一种算法简单、收敛速度好的颗粒流微观力学参数反演方法。
本发明解决上述问题的技术方案是：一种颗粒流微观力学参数反演方法，包括以下步骤：
步骤一：通过物理实验获取岩石物理实验宏观力学参数；
步骤二：基于物理宏观力学参数反演出颗粒流微观力学参数；
步骤三：通过数值模拟计算获取数值模拟宏观力学参数；
步骤四：采用模拟退火算法调整颗粒流微观力学参数，使得数值模拟计算宏观力学参数与物理实验宏观力学参数不断接近，当数值模拟宏观力学参数与物理实验宏观力学参数之间的误差小于10％时，则此时对应的颗粒流微观力学参数则是所需要确定的微观力学参数。
上述颗粒流微观力学参数反演方法，所述步骤一中，物理实验宏观力学参数包括单轴抗压强度UCSexperimental、弹性模量Eexperimental、泊松比vexperimental。
上述颗粒流微观力学参数反演方法，所述步骤二中，颗粒流微观力学参数包括颗粒密度ρ、颗粒最大半径与最小半径比值Rmax/Rmin、颗粒最小半径Rmin、颗粒接触刚度Ec、颗粒法向刚度与切向刚度比值kn/ks、平行连接刚度
Figure BDA0003138286360000021
平行连接法向刚度与切向刚度比值
Figure BDA0003138286360000022
颗粒摩擦系数μ、平行连接法向刚度平均值σc-mean、平行连接法向刚度方差σc-std、平行连接剪切刚度平均值τc-mean、平行连接剪切刚度方差τc-std。
上述颗粒流微观力学参数反演方法，所述步骤三中，将步骤二所得到的流微观力学参数开展数值模拟计算，并获取该组微观力学参数所对应的数值模拟宏观力学参数，即单轴抗压强度UCSnumerical、弹性模量Enumerical、泊松比vnumerical。
上述颗粒流微观力学参数反演方法，所述步骤四中，采用模拟退火算法调整颗粒流微观力学参数的具体步骤为：
步骤1：初始化模拟退火算法超参数：温度Temperature、降温系数Decay、马可夫链长度Markov及数据跳动因子Stepfactor，马可夫链长度Markov表示每一温度条件下循环次数，迭代次数Iteration＝0；
初始数值模拟计算时，微观力学参数：颗粒密度、颗粒最大半径与最小半径比值、颗粒最小半径、颗粒接触刚度、颗粒法向刚度与切向刚度比值、平行连接刚度、平行连接法向刚度与切向刚度比值、颗粒摩擦系数、平行连接法向刚度平均值、平行连接法向刚度方差、平行连接剪切刚度平均值、平行连接剪切刚度方差的最大值和最小值分别为ρmax、Rmax-max/Rmin-max、Rmin-max、Ec-max、kn-max/ks-max、
Figure BDA0003138286360000031
μmax、σc-mean-max、σc-std-max、τc-mean-max、τc-std-max和ρmin、Rmax-min/Rmin-min、Rmin-min、Ec-min、kn-min/ks-min、
Figure BDA0003138286360000032
μmin、σc-mean-min、σc-std-min、τc-mean-min、τc-std-min，然后在微观力学参数最大值与最小值之间选取随机值ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000033
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre作为微观力学参数初始值，采用该组微观力学参数开展数值模拟单轴压缩实验，获得数值模拟宏观力学参数单轴抗压强度、弹性模量和泊松比分别为UCSnumerical-pre、Enumerical-pre、vnumerical-pre；
计算获得的数值模拟宏观力学参数与物理试验获得的宏观力学参数之间的相对误差最大值Judgepre，Judgepre即作为判别值表示为：
Figure BDA0003138286360000034
同时将初始的微观力学参数视为最优的微观力学参数：
Figure BDA0003138286360000041
其对应的判别值也是目前最优的相对误差值：
Judgebest＝Judgepre (3)
步骤2：基于先前的微观力学参数集合(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000042
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)确定下一组微观力学参数的数值，其计算公式如下：
Figure BDA0003138286360000043
上式中rand为-1到1之间的随机数，通过式(4)基于先前的微观力学参数集合(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000051
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)，产生当前的微观力学参数集合(ρnext、Rmax-next/Rmin-next、Rmin-next、Ec-next、kn-next/ks-next、
Figure BDA0003138286360000052
μnext、σc-mean-next、σc-std-next、τc-mean-next、τc-std-next)，采用当前的微观力学参数集合开展单轴压缩数值模拟实验获取相应的数值模拟模型宏观力学参数：单轴抗压强度UCSnumerical-next、弹性模量Enumerical-next和泊松比vnumerical-next，并计算判别函数的值Judgenext：
Figure BDA0003138286360000053
步骤3：如果计算得到的Judgenext<Judgebest，则更新最优的微观力学参数集合及最优的目标值：
Figure BDA0003138286360000054
Judgebest＝Judgenext (7)
如果Judgenext>Judgebest，则无需更新最优微观力学参数集合；
步骤4：如果计算得到的Judgenext<Judgepre，则更新先前颗粒流微观力学参数及先前判断函数值：
Figure BDA0003138286360000061
Judgepre＝Judgenext (9)
而当Judgenext≥Judgepre时，先计算概率值p1，p1表示为：
Figure BDA0003138286360000062
然后产生一个0到1的随机数p2，如果p1>p2则更新先前颗粒流微观力学参数及先前判断函数值，即执行公式(8)和公式(9)，否则不执行任何操作；
步骤5：迭代次数更新，Iteration＝Iteration+1，当Iteration>Markov时：Iteration＝0，当前温度下降Temperature＝Temperature×Decay；执行完步骤5之后，再次执行步骤2，依此循环。
上述颗粒流微观力学参数反演方法，所述步骤四中，数值模拟计算终止条件是数值模拟计算得到的宏观力学参数与物理实验得到的宏观力学参数差值相对误差小于10％；数值模拟计算过程中，每迭代一次就对Judgebest进行判定，当Judgebest≤10％时，则停止计算，其中Judgebest表示为：
Figure BDA0003138286360000063
式中，UCSnumerical-best，Enumerical-best，vnumerical-best表示计算过程中最佳的单轴抗压强度、弹性模量及泊松比组合。
本发明的有益效果在于：本发明首先通过物理实验获取岩石宏观力学参数，然后基于物理宏观力学参数反演出颗粒流微观力学参数，接着通过数值模拟计算获取数值模拟宏观力学参数，最后采用模拟退火算法调整颗粒流微观力学参数，使得数值模拟宏观力学参数与物理实验宏观力学参数不断接近，当数值模拟宏观力学参数与物理实验宏观力学参数之间的误差小于10％时，则此时对应的颗粒流微观力学参数则是所需要确定的微观力学参数。本发明操作简便，微观力学参数调整过程中无主观性，收敛速度较好，适用于颗粒流数值模拟计算领域。
附图说明
图1为本发明的流程图。
图2为本发明采用模拟退火算法调整颗粒流微观力学参数的流程图。
具体实施方式
下面结合附图和实施例对本发明做进一步的说明。
如图1所示，一种颗粒流微观力学参数反演方法，包括以下步骤：
步骤一：通过物理实验获取岩石物理实验宏观力学参数，物理实验宏观力学参数包括单轴抗压强度UCSexperimental、弹性模量Eexperimental、泊松比vexperimental。
步骤二：基于物理宏观力学参数反演出颗粒流微观力学参数，颗粒流微观力学参数包括颗粒密度ρ、颗粒最大半径与最小半径比值Rmax/Rmin、颗粒最小半径Rmin、颗粒接触刚度Ec、颗粒法向刚度与切向刚度比值kn/ks、平行连接刚度
Figure BDA0003138286360000072
平行连接法向刚度与切向刚度比值
Figure BDA0003138286360000071
颗粒摩擦系数μ、平行连接法向刚度平均值σc-mean、平行连接法向刚度方差σc-std、平行连接剪切刚度平均值τc-mean、平行连接剪切刚度方差τc-std。
步骤三：将步骤二所得到的流微观力学参数开展数值模拟计算，并获取该组微观力学参数所对应的数值模拟宏观力学参数，即单轴抗压强度UCSnumerical、弹性模量Enumerical、泊松比vnumerical。
步骤四：采用模拟退火算法调整颗粒流微观力学参数，使得数值模拟计算宏观力学参数与物理实验宏观力学参数不断接近，当数值模拟宏观力学参数与物理实验宏观力学参数之间的误差小于10％时，则此时对应的颗粒流微观力学参数则是所需要确定的微观力学参数。
采用模拟退火算法调整颗粒流微观力学参数的具体步骤为：
步骤1：初始化模拟退火算法超参数：温度Temperature、降温系数Decay、马可夫链长度Markov及数据跳动因子Stepfactor，马可夫链长度Markov表示每一温度条件下循环次数，迭代次数Iteration＝0；
初始数值模拟计算时，微观力学参数：颗粒密度、颗粒最大半径与最小半径比值、颗粒最小半径、颗粒接触刚度、颗粒法向刚度与切向刚度比值、平行连接刚度、平行连接法向刚度与切向刚度比值、颗粒摩擦系数、平行连接法向刚度平均值、平行连接法向刚度方差、平行连接剪切刚度平均值、平行连接剪切刚度方差的最大值和最小值分别为ρmax、Rmax-max/Rmin-max、Rmin-max、Ec-max、kn-max/ks-max、
Figure BDA0003138286360000081
μmax、σc-mean-max、σc-std-max、τc-mean-max、τc-std-max和ρmin、Rmax-min/Rmin-min、Rmin-min、Ec-min、kn-min/ks-min、
Figure BDA0003138286360000082
μmin、σc-mean-min、σc-std-min、τc-mean-min、τc-std-min，然后在微观力学参数最大值与最小值之间选取随机值ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000083
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre作为微观力学参数初始值，采用该组微观力学参数开展数值模拟单轴压缩实验，获得数值模拟宏观力学参数单轴抗压强度、弹性模量和泊松比分别为UCSnumerical-pre、Enumerical-pre、vnumerical-pre；
计算获得的数值模拟宏观力学参数与物理试验获得的宏观力学参数之间的相对误差最大值Judgepre，Judgepre即作为判别值表示为：
Figure BDA0003138286360000091
同时将初始的微观力学参数视为最优的微观力学参数：
Figure BDA0003138286360000092
其对应的判别值也是目前最优的相对误差值：
Judgebest＝Judgepre (3)
步骤2：基于先前的微观力学参数集合(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000093
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)确定下一组微观力学参数的数值，其计算公式如下：
Figure BDA0003138286360000101
上式中rand为-1到1之间的随机数，通过式(4)基于先前的微观力学参数集合(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000102
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)，产生当前的微观力学参数集合(ρnext、Rmax-next/Rmin-next、Rmin-next、Ec-next、kn-next/ks-next、
Figure BDA0003138286360000103
μnext、σc-mean-next、σc-std-next、τc-mean-next、τc-std-next)，采用当前的微观力学参数集合开展单轴压缩数值模拟实验获取相应的数值模拟模型宏观力学参数：单轴抗压强度UCSnumerical-next、弹性模量Enumerical-next和泊松比vnumerical-next，并计算判别函数的值Judgenext：
Figure BDA0003138286360000104
步骤3：如果计算得到的Judgenext<Judgebest，则更新最优的微观力学参数集合及最优的目标值：
Figure BDA0003138286360000111
Judgebest＝Judgenext (7)
如果Judgenext>Judgebest，则无需更新最优微观力学参数集合；
步骤4：如果计算得到的Judgenext<Judgepre，则更新先前颗粒流微观力学参数及先前判断函数值：
Figure BDA0003138286360000112
Judgepre＝Judgenext (9)
而当Judgenext≥Judgepre时，先计算概率值p1，p1表示为：
Figure BDA0003138286360000113
然后产生一个0到1的随机数p2，如果p1>p2则更新先前颗粒流微观力学参数及先前判断函数值，即执行公式(8)和公式(9)，否则不执行任何操作；
步骤5：迭代次数更新，Iteration＝Iteration+1，当Iteration>Markov时：Iteration＝0，当前温度下降Temperature＝Temperature×Decay；执行完步骤5之后，再次执行步骤2，依此循环。
数值模拟计算终止条件是数值模拟计算得到的宏观力学参数与物理实验得到的宏观力学参数差值相对误差小于10％；数值模拟计算过程中，每迭代一次就对Judgebest进行判定，当Judgebest≤5％时，则停止计算，其中Judgebest表示为：
Figure BDA0003138286360000121
式中，UCSnumerical-best，Enumerical-best，vnumerical-best表示计算过程中最佳的单轴抗压强度、弹性模量及泊松比组合。
实施例
一种颗粒流微观力学参数反演方法，包括基于颗粒流微观力学参数开展单轴压缩实验获取数值模拟模型宏观力学参数、根据数值模拟宏观力学参数与物理实验宏观力学参数差值采用模拟退火算法调整颗粒流微观力学参数、数值模拟计算终止条件；结合具体实例其步骤如下：
步骤1：通过室内实验获取岩石的单轴抗压强度UCSexprimental、弹性模量Eexperimental、泊松比vexperimental分别为：37MPa、23.3GPa和0.17。同时确定模拟退火算法的超参数温度Temperature为100、降温系数Decay为0.99、马可夫链长度Markov(每一温度条件下循环次数)为10000，数据跳动因子Stepfactor为0.02。并确定微观力学参数(颗粒密度ρ、颗粒最大半径与最小半径比值Rmax/Rmin、颗粒最小半径Rmin、颗粒接触刚度Ec、颗粒法向刚度与切向刚度比值kn/ks、平行连接刚度
Figure BDA0003138286360000134
平行连接法向刚度与切向刚度比值
Figure BDA0003138286360000135
颗粒摩擦系数μ、平行连接法向刚度平均值σc-mean、平行连接法向刚度方差σc-std、平行连接剪切刚度平均值τc-mean、平行连接剪切刚度方差τc-std)的最小值和最大值分别为：(1000,1,0.1e-3,1e9,0.1,1e9,0.1,0.1,1e6,1e6,1e6,1e6)和(5000,10,10e-3,1000e9,10,1000e9,10,10,1000e6,100e6,1000e6,100e6)。
步骤2：初始微观力学参数(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000131
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)是通过在微观力学参数最小值和最大值直接随机选取，然后得到初始微观力学参数组合。然后采用初始微观力学参数组合开展数值单轴压缩实验获得宏观力学参数：单轴抗压强度、弹性模量和泊松比分别为UCSnumerical-pre、Enumerical-pre、vnumerical-pre。计算出目标函数：
Figure BDA0003138286360000132
同时将初始的微观力学参数视为最优的微观力学参数：
Figure BDA0003138286360000133
其对应的判别值也设置为目前最优的判别值：
Judgebest＝Judgepre
步骤3：基于先前的微观力学参数集合(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000141
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)确定下一组微观力学参数的数值，其计算公式如下：
Figure BDA0003138286360000142
式中rand为-1到1的随机数。通过上式可基于先前的微观力学参数集合(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure BDA0003138286360000143
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)，产生当前的微观力学参数集合(ρnext、Rmax-next/Rmin-next、Rmin-next、Ec-next、kn-next/ks-next、
Figure BDA0003138286360000144
μnext、σc-mean-next、σc-std-next、τc-mean-next、τc-std-next)，采用当前的微观力学参数集合开展单轴压缩数值模拟实验获取相应的数值模拟模型宏观力学参数：单轴抗压强度UCSnumerical-next、弹性模量Enumerical-next和泊松比vnumerical-next。并计算判别函数的值Judgenext：
Figure BDA0003138286360000145
步骤4：如果计算得到的Judgenext<Judgebest，则更新最优的微观力学参数集合及最优的目标值：
Figure BDA0003138286360000151
Judgebest＝Judgenext
如果Judgenext>Judgebest，则无需更新最优微观力学参数集合。
步骤5：如果计算得到的Judgenext<Judgepre，则更新先前颗粒流微观力学参数及先前判断函数值：
Figure BDA0003138286360000152
Judgepre＝Judgenext
而当Judgenext≥Judgepre时，先计算概率值p1，p1可表示为：
Figure BDA0003138286360000153
然后产生一个0到1的随机数p2，如果p1>p2则需更新先前颗粒流微观力学参数及先前判断函数值，即执行公式(17)和公式(18)，否则不执行任何操作。
步骤6：迭代次数更新，Iteration＝Iteration+1，当Iteration>10000时：Iteration＝0，当前温度下降Temperature＝Temperature×0.99。执行完步骤6之后，再次执行步骤3，依次循环。
所述数值模拟计算终止条件是数值模拟计算得到的宏观力学参数与物理实验得到的宏观力学参数差值相对误差小于10％时，则终止数值模拟计算。数值模拟计算过程中，每迭代依次就对Judgebest进行判定，当Judgebest≤10％时，则停止计算。其中Judgebest可表示为：
Figure BDA0003138286360000161
式中，UCSnumerical-best，Enumerical-best，vnumerical-best表示计算过程中最佳的单轴抗压强度、弹性模量及泊松比组合。通过计算得到UCSnumerical-best，Enumerical-best，vnumerical-best分别为35.82MPa、23.18GPa和0.16，其中对应的Judgebest为5.8％，小于实验终止条件10％。
同时得到对应最佳的微观力学参数集合(ρbest、Rmax-best/Rmin-best、Rmin-best、Ec-best、kn-best/ks-best、
Figure BDA0003138286360000162
μbest、σc-mean-best、σc-std-best、τc-mean-best、τc-std-best)为(2411，5.54，0.76e-3，1.76e9，0.57，55.96e9，2.94,6.28，77.63e6，77.76e6，64.04e6，33.52e6)。

Claims (2)

1.一种颗粒流微观力学参数反演方法，其特征在于，包括以下步骤：
步骤一：通过物理实验获取岩石物理实验宏观力学参数；
物理实验宏观力学参数包括单轴抗压强度UCSexperimental、弹性模量Eexperimental、泊松比vexperimental；
步骤二：基于物理宏观力学参数反演出颗粒流微观力学参数；
颗粒流微观力学参数包括颗粒密度ρ、颗粒最大半径与最小半径比值Rmax/Rmin、颗粒最小半径Rmin、颗粒接触刚度Ec、颗粒法向刚度与切向刚度比值kn/ks、平行连接刚度
Figure FDA0003722743800000011
平行连接法向刚度与切向刚度比值
Figure FDA0003722743800000012
颗粒摩擦系数μ、平行连接法向刚度平均值σc-mean、平行连接法向刚度方差σc-std、平行连接剪切刚度平均值τc-mean、平行连接剪切刚度方差τc-std；
步骤三：通过数值模拟计算获取数值模拟宏观力学参数；
将步骤二所得到的流微观力学参数开展数值模拟计算，并获取该组微观力学参数所对应的数值模拟宏观力学参数，即单轴抗压强度UCSnumerical、弹性模量Enumerical、泊松比vnumerical；
步骤四：采用模拟退火算法调整颗粒流微观力学参数，使得数值模拟计算宏观力学参数与物理实验宏观力学参数不断接近，当数值模拟宏观力学参数与物理实验宏观力学参数之间的误差小于10％时，则此时对应的颗粒流微观力学参数则是所需要确定的微观力学参数；
采用模拟退火算法调整颗粒流微观力学参数的具体步骤为：
步骤1：初始化模拟退火算法超参数：温度Temperature、降温系数Decay、马可夫链长度Markov及数据跳动因子Stepfactor，马可夫链长度Markov表示每一温度条件下循环次数，迭代次数Iteration＝0；
初始数值模拟计算时，微观力学参数：颗粒密度、颗粒最大半径与最小半径比值、颗粒最小半径、颗粒接触刚度、颗粒法向刚度与切向刚度比值、平行连接刚度、平行连接法向刚度与切向刚度比值、颗粒摩擦系数、平行连接法向刚度平均值、平行连接法向刚度方差、平行连接剪切刚度平均值、平行连接剪切刚度方差的最大值和最小值分别为ρmax、Rmax-max/Rmin-max、Rmin-max、Ec-max、kn-max/ks-max、
Figure FDA0003722743800000021
μmax、σc-mean-max、σc-std-max、τc-mean-max、τc-std-max和ρmin、Rmax-min/Rmin-min、Rmin-min、Ec-min、kn-min/ks-min、
Figure FDA0003722743800000022
μmin、σc-mean-min、σc-std-min、τc-mean-min、τc-std-min，然后在微观力学参数最大值与最小值之间选取随机值ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure FDA0003722743800000023
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre作为微观力学参数初始值，采用该组微观力学参数开展数值模拟单轴压缩实验，获得数值模拟宏观力学参数单轴抗压强度、弹性模量和泊松比分别为UCSnumerical-pre、Enumerical-pre、vnumerical-pre；
计算获得的数值模拟宏观力学参数与物理试验获得的宏观力学参数之间的相对误差最大值Judgepre，Judgepre即作为判别值表示为：
Figure FDA0003722743800000024
同时将初始的微观力学参数视为最优的微观力学参数：
Figure FDA0003722743800000031
其对应的判别值也是目前最优的相对误差值：
Judgebest＝Judgepre (3)
步骤2：基于先前的微观力学参数集合(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure FDA0003722743800000032
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)确定下一组微观力学参数的数值，其计算公式如下：
Figure FDA0003722743800000033
上式中rand为-1到1之间的随机数，通过式(4)基于先前的微观力学参数集合(ρpre、Rmax-pre/Rmin-pre、Rmin-pre、Ec-pre、kn-pre/ks-pre、
Figure FDA0003722743800000041
μpre、σc-mean-pre、σc-std-pre、τc-mean-pre、τc-std-pre)，产生当前的微观力学参数集合(ρnext、Rmax-next/Rmin-next、Rmin-next、Ec-next、kn-next/ks-next、
Figure FDA0003722743800000042
μnext、σc-mean-next、σc-std-next、τc-mean-next、τc-std-next)，采用当前的微观力学参数集合开展单轴压缩数值模拟实验获取相应的数值模拟模型宏观力学参数：单轴抗压强度UCSnumerical-next、弹性模量Enumerical-next和泊松比vnumerical-next，并计算判别函数的值Judgenext：
Figure FDA0003722743800000043
步骤3：如果计算得到的Judgenext<Judgebest，则更新最优的微观力学参数集合及最优的目标值：
Figure FDA0003722743800000044
Judgebest＝Judgenext (7)
如果Judgenext>Judgebest，则无需更新最优微观力学参数集合；
步骤4：如果计算得到的Judgenext<Judgepre，则更新先前颗粒流微观力学参数及先前判断函数值：
Figure FDA0003722743800000051
Judgepre＝Judgenext (9)
而当Judgenext≥Judgepre时，先计算概率值p1，p1表示为：
Figure FDA0003722743800000052
然后产生一个0到1的随机数p2，如果p1>p2则更新先前颗粒流微观力学参数及先前判断函数值，即执行公式(8)和公式(9)，否则不执行任何操作；
步骤5：迭代次数更新，Iteration＝Iteration+1，当Iteration>Markov时：Iteration＝0，当前温度下降Temperature＝Temperature×Decay；执行完步骤5之后，再次执行步骤2，依此循环。
2.根据权利要求1所述的颗粒流微观力学参数反演方法，其特征在于，所述步骤四中，数值模拟计算终止条件是数值模拟计算得到的宏观力学参数与物理实验得到的宏观力学参数差值相对误差小于10％；数值模拟计算过程中，每迭代一次就对Judgebest进行判定，当Judgebest≤10％时，则停止计算，其中Judgebest表示为：
Figure FDA0003722743800000053
式中，UCSnumerical-best，Enumerical-best，vnumerical-best表示计算过程中最佳的单轴抗压强度、弹性模量及泊松比组合。
CN202110725094.9A 2021-06-29 2021-06-29 一种颗粒流微观力学参数反演方法 Active CN113536549B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110725094.9A CN113536549B (zh) 2021-06-29 2021-06-29 一种颗粒流微观力学参数反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110725094.9A CN113536549B (zh) 2021-06-29 2021-06-29 一种颗粒流微观力学参数反演方法

Publications (2)

Publication Number Publication Date
CN113536549A CN113536549A (zh) 2021-10-22
CN113536549B true CN113536549B (zh) 2022-08-26

Family

ID=78097064

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110725094.9A Active CN113536549B (zh) 2021-06-29 2021-06-29 一种颗粒流微观力学参数反演方法

Country Status (1)

Country Link
CN (1) CN113536549B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117113517B (zh) * 2023-10-23 2024-01-23 天津建城基业集团有限公司 一种预制桩基坑支护颗粒流数值模拟方法和系统

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110112394A1 (en) * 2009-11-11 2011-05-12 Mishelevich David J Neuromodulation of deep-brain targets using focused ultrasound
CN103940666A (zh) * 2014-03-18 2014-07-23 中国矿业大学 一种模拟断续裂隙岩石力学特性的细观参数确定方法
US10394976B2 (en) * 2016-12-12 2019-08-27 IFP Energies Nouvelles Method of exploiting hydrocarbons from a sedimentary basin comprising carbonate rocks, by means of stratigraphic simulation
CN108629126B (zh) * 2018-05-09 2020-04-21 中国地质大学（北京） 一种考虑宏细观缺陷耦合的岩体力学数值建模方法
CN109916754A (zh) * 2019-02-26 2019-06-21 成都理工大学 一种基于岩屑微观特征和钻井参数的储层脆性评价方法

Also Published As

Publication number Publication date
CN113536549A (zh) 2021-10-22

Similar Documents

Publication Publication Date Title
CN111310915B (zh) 一种面向强化学习的数据异常检测防御方法
CN109931678B (zh) 基于深度学习lstm的空调故障诊断方法
CN108161934B (zh) 一种利用深度强化学习实现机器人多轴孔装配的方法
JP4334178B2 (ja) 戦略パラメータの最適化方法
CN113536549B (zh) 一种颗粒流微观力学参数反演方法
WO2023019601A1 (zh) 基于结构优化算法的复值神经网络的信号调制识别方法
WO2019184777A1 (zh) 用于竞速类ai模型的ai参数配置方法、装置、设备及存储介质
US7610250B2 (en) Real-time method of determining eye closure state using off-line adaboost-over-genetic programming
CN113947016A (zh) 针对电网紧急控制系统中深度强化学习模型的脆弱性评估方法
Moreno et al. Using prior knowledge to improve reinforcement learning in mobile robotics
Swazinna et al. Comparing model-free and model-based algorithms for offline reinforcement learning
CN107016239B (zh) 一种汽轮机阀门流量特性分析方法
CN109951327B (zh) 一种基于贝叶斯混合模型的网络故障数据合成方法
CN112149884A (zh) 一种面向大规模学员的学业预警监测方法
CN108428256A (zh) 一种基于柔软度的自适应网格细化的软组织形变仿真方法
CN112518742B (zh) 基于动态模型与事后经验回放的多目标机器人控制方法
CN116700003A (zh) 使用流程工业历史数据构建强化学习环境的方法及系统
CN115905861A (zh) 一种基于多级算子变异的强化学习框架安全性检测方法
CN115972197A (zh) 一种基于关节角轨迹编码的机器人动作演示学习方法
CN115510753A (zh) 群智网络中基于矩阵补全与强化学习的数据收集方法
CN116047904A (zh) 面向机器人操作技能学习的人员仿真现实混合训练方法
CN111445005A (zh) 基于强化学习的神经网络控制方法及强化学习系统
Hošovský et al. Improved modeling of pneumatic muscle actuator using recurrent neural network
CN111582495A (zh) 基于优胜劣汰的深度强化学习策略网络存储方法及设备
Pohland et al. Efficient Safe Learning for Robotic Systems in Unstructured Environments

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
GR01 Patent grant
GR01 Patent grant