一种岩土材料本构模型建立方法

1.本发明属于岩土工程领域，具体涉及一种岩土材料本构模型建立方法。


背景技术：

2.本构模型是评价岩土材料力学性能和变形性能的重要工具之一，对于合理设计和精确评价岩石结构的稳定性和可靠性具有重要的指导意义。本构模型简称为sdc模型。现有本构模型大致可分为线性模型和非线性模型两类，如图1 为现有本构模型示意图；如图1a为弹脆性模型，图1b为三线性模型，图1c为四线性模型，图1d为五线性模型，图1e为duncan模型，图1f为mazars模型，图1g为loland模型，图1h为统计损伤。其中，线性线性模型包括弹脆性损伤模型、三线性模型、四线性模型(即双线性应变软化-残余理想塑性模型)和五线性模型；非线性本构模型包括duncan模型、loland模型、mazars模型和scd 模型。这些模型均对轴向应力-应变的压密阶段作了简化假设，这限制了这些模型只能定性分析而不能定量描述完整的轴向应力-应变曲线。压密阶段的非线性特征可以归结为原始微裂纹的变形。对于花岗岩、辉绿岩、大理岩等致密岩石，原始微裂纹较少，非线性变形较小，忽略压密阶段的简化假设与对致密岩石的实验观测结果基本一致。而对于多孔岩石，如碳质页岩、砂质板岩、碎屑砂岩等则不适用该假设，因为它们轴向应力-应变曲线的压密阶段明显，并表现出明显的非线性特征。由于压密阶段是岩石变形的关键阶段，在建立本构模型时应考虑压密阶段。
3.一般而言，岩土材料在压缩作用下的应力-应变曲线可分为压密阶段、弹性阶段、塑性屈服阶段和峰后阶段四个阶段，如图2为典型压缩试验轴向应力-应变示意图。如图2所示压密点a和屈服点b的应变分别为ε
cc
和εy，通过识别3 个特殊的应力阈值(裂纹闭合应力σ
cc
、屈服应力σy和峰值应力σ
p
)，可以确定4 个阶段的范围：(1)裂纹闭合应力代表原始微裂纹的完全闭合和弹性变形的开始；(2)弹性阶段结束时出现屈服应力，岩石开始发生非线性变形；(3)峰值应力为整个本构曲线的最大应力。值得注意的是，在三轴压缩试验中，轴向应力σi指的是偏应力σ
1-σ3，而在单轴压缩试验中，σi即为最大主应力σ1。
4.目前关于考虑压密阶段的岩土材料本构模型的一些代表性的研究如下：
5.cao,et al.[k.cao,l.ma,y.wu,a.spearing,n.khan,and y.xie,"thedetermination of a damage model for mudstone under uniaxial loading in acidicconditions,"geofluids,vol.2020,pp.1-12,10/08 2020.]将压密阶段的本构关系唯像地简化为基于裂纹闭合应力的二次函数。liu,et al.[h.h.liu,j.rutqvist,and j. g.berryman,"on the relationship between stress and elastic strain for porous andfractured rock,"int j rock mech min sci,vol.46,no.2,pp.289-296,2009.]开发了双应变胡克模型来模拟岩石的应力-应变关系。cao,et al.[w.g.cao,h.zhao,x. li,and y.j.zhang,"statistical damage model with strain softening and hardening forrocks under the influence of voids and volume changes,"(in english),can.geotech. j.,article vol.47,no.8,pp.857-871,aug 2010.]将岩石材料抽象为孔隙
和岩石基体两部分(包括损伤部分和未损伤部分)，并对经典lemaitre损伤模型进行修正，假设损伤部分也能承受载荷。这在一定程度上提高了岩石整体变形的计算精度。 peng,et al.[j.peng,g.rong,c.b.zhou,and k.peng,"a study of crack closureeffect of rocks and its quantitative model,"(in chinese),rock and soil mechanics, article vol.37,no.1,pp.126-132,jan 2016],[j.peng,g.rong,m.cai,and c.b. zhou,"a model for characterizing crack closure effect of rocks,"eng geol,vol.189, pp.48-57,2015/04/30/2015.]基于北非辉长岩和卡拉拉大理岩的压缩试验结果，针对岩石压密阶段提出了一个经验本构模型。
[0006]
但是以上研究尚未明确一个公认的具有明确物理数学含义的模型，每一方法均带有较强经验性与主观性，模型描述结果的精确度低，误差大，偏离实际结果较大。


技术实现要素：

[0007]
本发明的目的在于提供一种岩土材料本构模型建立方法，该方法考虑了孔隙压实特性，建立了更加高效准确的岩土材料本构模型。
[0008]
本发明提供的这种岩土材料本构模型建立方法，包括如下步骤：
[0009]
s1.进行单轴压缩试验，基于轴向应力差确定压密点和屈服点；
[0010]
s2.基于孔隙的压实特性，计算岩石基质的轴向应变和孔隙的轴向应变，建立多级加载后，屈服点之前的本构模型；
[0011]
s3.通过基于屈服点修正的宏观应变的概率密度函数，建立屈服点之后的本构模型；
[0012]
s4.结合屈服点之前的本构模型和屈服点之后的本构模型，完成岩土材料本构模型的建立。
[0013]
所述的步骤s1，包括对岩石试样进行单轴压缩试验，基于试验数据建立岩石试样的轴向应变-轴向应力曲线图；绘制基于弹性模量的参考直线；绘制应力差-应变曲线，包括计算轴向应力-应变曲线上的每个试验数据点与参考直线之间的轴向应力差；求得的应力差-应变曲线水平端的一端为压密点，另一端为屈服点。
[0014]
所述的步骤s1，包括对岩石试样进行单轴压缩试验，基于试验数据建立岩石试样的轴向应变-轴向应力曲线图；绘制基于弹性模量e的参考直线，参考直线的公式为σi*＝eεi，其中，σi*为参考应力，εi表示应变；绘制应力差-应变曲线，包括计算轴向应力-应变曲线上的每个试验数据点(εi,σi)与参考直线σi*＝eεi之间的轴向应力差(σi*-σi)；σi表示应力；求得的应力差-应变曲线水平端的一端为压密点a，另一端为屈服点b。
[0015]
所述的步骤s2，包括如下步骤：
[0016]
a1.建立柱状岩石单元的初始模型；
[0017]
a2.基于孔隙的压实特性，计算岩石基质的轴向应变和孔隙的轴向应变；
[0018]
a3.将轴向应力下的变形代入柱状岩石模型中，获取岩石的基础本构模型；
[0019]
a4.将步骤a2获取的岩石基质的轴向应变和孔隙的轴向应变代入步骤a3 岩石的基础本构模型，获取屈服点前的岩石的最终本构模型。
[0020]
所述的步骤a1，包括沿轴向应力方向取一个代表性的柱状岩石单元，设置柱状岩石单元的轴向长度为l，岩石基质长度为lm，孔隙长度为lv，γ为孔隙在柱状岩石单元的比
例；所述的步骤a2，包括求得加载后的岩石基质的轴向应变εm：
[0021][0022]
其中，em表示岩石基质的弹性模量；μm表示岩石基质的泊松比；σ1表示最大主应力；σ2表示中间主应力；σ3表示最小主应力；并设压缩后孔隙的长度为则孔隙的轴向应变εv为：
[0023][0024]
其中，lv为孔隙长度；
[0025]
应力σi在压缩试验中为多级加载，第i级的应力增量为趋于无穷小；孔隙的轴向应变εv为每一级应力增量作用下的应变δε
v,i
之和：
[0026][0027]
其中，e
vi
表示第i级应力增量对应的孔隙变形模量；根据推导得到压缩后孔隙的长度
[0028][0029]
其中，lv为孔隙长度；e
vi
表示第i级应力增量对应的孔隙变形模量；从而求得压缩后的孔隙的轴向应变εv：
[0030][0031]
其中，σ1表示最大主应力；ev为孔隙的弹性模量。
[0032]
所述的步骤a3，包括设岩石基质的变形量为δlm，孔隙的变形量为δlv，计算获取岩石基质的轴向应变εm与岩石基质的变形量δlm的关系和孔隙的轴向应变εv与孔隙的变形量δlv的关系：
[0033][0034]
其中，l为柱状岩石单元的轴向长度；δlm为岩石基质变形量；δlv为孔隙变形；γ为孔隙在柱状岩石单元的比例；设柱状岩石单元在轴向应力作用下的宏观变形δl＝δlm δlv，求得岩石的基础本构模型：
[0035]
[0036]
其中，ε为宏观应变，l为柱状岩石单元的轴向长度；δlm为岩石基质变形量；δlv为孔隙变形；γ为孔隙在柱状岩石单元的比例。
[0037]
所述的步骤a4，包括将步骤a2获取的岩石基质的轴向应变εm和孔隙的轴向应变εv代入步骤a3岩石的基础本构模型，获取屈服点前的岩石的最终本构模型：
[0038][0039]
其中，e表示岩石材料的弹性模量；μ表示屈服点前岩石材料的泊松比；γ为孔隙在柱状岩石单元的比例；ev为孔隙的弹性模量；σ1表示最大主应力；σ2表示中间主应力；σ3表示最小主应力；
[0040]
计算孔隙的弹性模量ev和孔隙在柱状岩石单元的比例γ，包括：对于单轴压缩下的应力-应变曲线，在屈服点b(εy,σy)处存在轴向应力-应变曲线的导数等于屈服点前岩石材料的岩石材料的弹性模量e：
[0041][0042]
其中，e表示岩石材料的弹性模量；(εy,σy)为屈服点b的应力-应变坐标；γ为孔隙在柱状岩石单元的比例；ev为孔隙的弹性模量；σ1表示最大主应力；ε表示宏观应变；
[0043]
求解得到：设参考方程参考方程的交点的横坐标为孔隙的弹性模量ev。
[0044]
所述的步骤s3，包括如下步骤：
[0045]
b1.计算改进后的损伤模型：
[0046]
b2.采用修正weibull分布，计算宏观应变的概率密度函数；
[0047]
b3.通过宏观应变的概率密度函数计算损伤变量，并将损伤变量代入步骤 b1的改进后的损伤模型，得到屈服点之后的本构模型。
[0048]
所述的步骤b1，包括轴向应力由岩石基质和损伤岩石基质共同承担：
[0049]
σia＝e(ε-ε
cc
)a
*
σr(a-a
*
)
[0050]
其中，σi表示应力；σr表示损伤岩石基质的应力；a为轴向应力作用区域；ε表示宏观应变；e表示岩石材料的弹性模量；ε
cc
表示压密点a的应变；a*表示未损伤岩石基质的分布区域；
[0051]
设置损伤变量d为：d＝(a-a*)/a，则把公式化简为
[0052]
σi＝e(ε-ε
cc
)(1-d) σ
rd[0053]
其中，σi表示应力；σr表示损伤岩石基质的应力；ε表示宏观应变；e表示岩石材料的弹性模量；ε
cc
表示压密点a的应变；结合裂纹闭合应力σ
cc
，得到改进后的损伤模型：
[0054]
σ
i-σ
cc
＝e(ε-ε
cc
)(1-d) (σ
r-σ
cc
)d
[0055]
其中，σi表示应力；σr表示损伤岩石基质的应力；ε表示宏观应变；e表示岩石材料的弹性模量；ε
cc
表示压密点a的应变；σ
cc
为裂纹闭合应力；d表示损伤变量。
[0056]
所述的步骤b2，包括将屈服点b(εy,σy)作为岩石损伤的起始点；因此，采用屈服应力σy修正weibull分布：
[0057][0058]
其中，ω(ε)表示宏观应变ε的概率密度函数；m为模型形状参数；ε0为模型尺度参数；(εy,σy)为屈服点b的应力-应变坐标；
[0059]
所述的步骤b3，包括对宏观应变ε的概率密度函数ω(ε)进行积分，得到损伤变量d：
[0060][0061]
其中，ω(ε)表示宏观应变ε的概率密度函数；m为模型形状参数；ε0为模型尺度参数；(εy,σy)为屈服点b的应力-应变坐标；
[0062]
将损伤变量d代入步骤b1的改进后的损伤模型，得到屈服点b(εy,σy)之后的本构模型：
[0063][0064]
其中，σi表示应力；ε表示宏观应变；σr表示损伤岩石基质的应力；e表示岩石材料的弹性模量；ε
cc
表示压密点a的应变；σ
cc
为裂纹闭合应力；(εy,σy)为屈服点b的应力-应变坐标；m为模型形状参数；ε0为模型尺度参数；拟合模型形状参数和模型尺度参数ε0，包括：
[0065]
将屈服点b(εy,σy)代入屈服点b(εy,σy)之后的本构模型：
[0066]
σy＝e(ε
y-ε
cc
) σ
cc
[0067]
其中，e表示岩石材料的弹性模量；ε
cc
表示压密点a的应变；σ
cc
为裂纹闭合应力；(εy,σy)为屈服点b的应力-应变坐标；采用数学运算得到模型形状参数 m和模型尺度参数ε0的关系：
[0068][0069]
将两边看作是新变量x,y和z：
[0070][0071]
则屈服点b(εy,σy)之后的本构模型可简化为：
[0072]
y＝mx z
[0073]
将应力-应变数据简化后的屈服点b(εy,σy)之后的本构模型，拟合线的斜率和截距分别为参数m和-mlnε0，并通过线性回归得到m和-mlnε0，求得模型形状参数m和模型尺度参数ε0的值。
[0074]
本发明提供的这种岩土材料本构模型建立方法，利用对岩石压密阶段的分析，并采用了屈服点对现有应变的概率密度函数进行修正，保证了对岩石压缩变形行为描述的精确性，尤其适用于岩土工程领域的理论建模以及工程稳定性分析。
附图说明
[0075]
图1为现有本构模型示意图；
[0076]
图2为典型压缩试验轴向应力-应变示意图。
[0077]
图3为本发明方法的流程示意图。
[0078]
图4为本发明方法的压密点和屈服点位置示意图。
[0079]
图5为本发明方法的柱状岩石单元的结构示意图。
[0080]
图6为本发明方法的参考方程交点示意图。
[0081]
图7为本发明方法的损伤岩石基质部分示意图。
[0082]
图8为本发明实施例的现有理论本构模型与试验结果比较示意图。
[0083]
图9为本发明实施例的a1试样的单轴压缩应力-应变曲线示意图。
[0084]
图10为本发明实施例的确定a1试样相关参数示意图。
[0085]
图11为本发明实施例的参考方程交点示意图。
[0086]
图12为本发明实施例的模型参数拟合示意图。
[0087]
图13为本发明实施例与试验结果的对比示意图。
具体实施方式
[0088]
如图3为本发明方法的流程示意图：本发明提供的这种岩土材料本构模型建立方法，包括如下步骤：
[0089]
s1.进行单轴压缩试验，基于轴向应力差确定压密点和屈服点；
[0090]
s2.基于孔隙的压实特性，计算岩石基质的轴向应变和孔隙的轴向应变，建立多级加载后，屈服点之前的本构模型；
[0091]
s3.通过基于屈服点修正的宏观应变的概率密度函数，建立屈服点之后的本构模型；
[0092]
s4.结合屈服点之前的本构模型和屈服点之后的本构模型，完成岩土材料本构模
型的建立。
[0093]
如图4为本发明方法的压密点和屈服点位置示意图。所述的步骤s1，包括对岩石试样进行单轴压缩试验，基于试验数据建立岩石试样的轴向应变-轴向应力曲线图；绘制基于弹性模量e的参考直线，参考直线的公式为σi*＝eεi，其中，σi*为参考应力，εi表示应变；绘制应力差-应变曲线，包括计算轴向应力
‑ꢀ
应变曲线上的每个试验数据点(εi,σi)与参考直线σi*＝eεi之间的轴向应力差 (σi*-σi)；σi表示应力；求得的应力差-应变曲线水平端的一端为压密点a，另一端为屈服点b。
[0094]
所述的步骤s2，包括如下步骤：
[0095]
如图5为本发明方法的柱状岩石单元的结构示意图。a1.建立柱状岩石单元的初始模型；为便于分析，将岩石基体视为由孔隙和岩石基质组成的复合地质材料；如图3所示，沿轴向应力方向取一个代表性的柱状岩石单元，设置柱状岩石单元的轴向长度为l，岩石基质长度为lm，孔隙长度为lv，γ为孔隙在柱状岩石单元的比例，因此：
[0096][0097]
a2.基于孔隙的压实特性，计算岩石基质的轴向应变和孔隙的轴向应变：由于岩石基体在加载后呈线性变形，卸载后可以恢复到屈服点之前的原始状态，因此，根据广义胡克定律，得到加载后的岩石基质的轴向应变εm：
[0098][0099]
其中，em表示岩石基质的弹性模量；μm表示岩石基质的泊松比；σ1表示最大主应力；σ2表示中间主应力；σ3表示最小主应力；在常规三轴压缩试验中，σ2＝σ3；当变形达到弹性阶段时，孔隙已接近闭合，因此，无需考虑孔隙的影响，这一阶段的变形可以近似于基岩基质的变形。因此，em和μm可以近似地被屈服点前岩石材料的弹性模量e和屈服点前岩石材料的泊松比μ所代替。
[0100]
从应变本身的定义出发，准确描述压缩后孔隙的变形：假设压缩后孔隙的长度为则孔隙的轴向应变εv为：
[0101][0102]
其中，lv为孔隙长度；
[0103]
假设应力σi在压缩试验中视为多级加载，第i级的应力增量为趋于无穷小；则孔隙的轴向应变εv可以转换为关于应力σi的函数；不考虑侧向应力的影响，孔隙的轴向应变εv为每一级应力增量作用下的应变δε
v,i
之和：
[0104][0105]
其中，e
vi
表示第i级应力增量对应的孔隙变形模量；则得到等式：
[0106][0107]
求解得到压缩后孔隙的长度
[0108][0109]
其中，lv为孔隙长度；e
vi
表示第i级应力增量对应的孔隙变形模量；从而求得压缩后的孔隙的轴向应变εv[0110][0111]
其中，e
vi
表示第i级应力增量对应的孔隙变形模量；
[0112]
若每一级应力增量作用时对应的孔隙变形模量e
vi
都保持不变，将化简为其中，ev为孔隙的弹性模量。
[0113]
a3.将轴向应力下的变形代入柱状岩石模型中，获取柱状岩石单元在轴向应力作用下的宏观应变：设岩石基质的变形量为δlm，孔隙的变形量为δlv，计算获取岩石基质的轴向应变εm与岩石基质的变形量δlm的关系和孔隙的轴向应变εv与孔隙的变形量δlv的关系：
[0114][0115]
其中，l为柱状岩石单元的轴向长度；δlm为岩石基质变形量；δlv为孔隙变形；γ为孔隙在柱状岩石单元的比例；设柱状岩石单元在轴向应力作用下的宏观变形δl＝δlm δlv，求得柱状岩石单元在轴向应力作用下的宏观应变ε，即岩石的基础本构模型：
[0116][0117]
其中，l为柱状岩石单元的轴向长度；δlm为岩石基质变形量；δlv为孔隙变形；γ为孔隙在柱状岩石单元的比例；在现有技术中，宏观应变ε为ε＝εm εv，本方法加入了孔隙在柱状岩石单元的比例对柱状岩石单元的影响，预测更加准确，更接近于岩石材料的实际变形。
[0118]
a4.将步骤a2获取的岩石基质的轴向应变εm和孔隙的轴向应变εv代入步骤a3岩石
的基础本构模型，获取屈服点前的岩石的最终本构模型：
[0119][0120]
其中，e表示岩石材料的弹性模量；μ表示屈服点前岩石材料的泊松比；γ为孔隙在柱状岩石单元的比例；ev为孔隙的弹性模量；σ1表示最大主应力；σ2表示中间主应力；σ3表示最小主应力；
[0121]
其中，孔隙的弹性模量ev和孔隙在柱状岩石单元的比例γ的求解方法如下：如图2所示屈服点b的应变为εy；对于单轴压缩下的应力-应变曲线，在屈服点 b(εy,σy)处存在特殊边界条件，即轴向应力-应变曲线的导数等于屈服点前岩石材料的岩石材料的弹性模量e：
[0122][0123]
其中，e表示岩石材料的弹性模量；(εy,σy)为屈服点b的应力-应变坐标；γ为孔隙在柱状岩石单元的比例；ev为孔隙的弹性模量；σ1表示最大主应力；ε表示宏观应变；
[0124]
求解得到：如图6为本发明方法的参考方程交点示意图，设参考方程参考方程的交点的横坐标为孔隙的弹性模量ev。
[0125]
如图7为本发明方法的损伤岩石基质部分示意图。因为岩石变形在屈服点之后进入塑性屈服阶段，因此岩石基质在弹性阶段未发生损伤，但在屈服点之后发生损伤。考虑到在三轴压缩试验中，如图2所示，峰值后阶段的岩样在摩擦和围压的影响下仍然具有承载能力，所述的步骤s3，包括如下步骤：
[0126]
b1.计算改进后的损伤模型：如图7所示，轴向应力由岩石基质和损伤岩石基质共同承担：
[0127]
σia＝e(ε-ε
cc
)a
*
σr(a-a
*
)
[0128]
其中，σi表示应力；σr表示损伤岩石基质的应力；a为轴向应力作用区域；ε表示宏观应变；e表示岩石材料的弹性模量；如图2所示ε
cc
表示压密点a的应变；a*表示未损伤岩石基质的分布区域；
[0129]
设置损伤变量d为：d＝(a-a*)/a，则把公式化简为
[0130]
σi＝e(ε-ε
cc
)(1-d) σ
rd[0131]
其中，σi表示应力；σr表示损伤岩石基质的应力；ε表示宏观应变；e表示岩石材料的弹性模量；ε
cc
表示压密点a的应变；结合裂纹闭合应力σ
cc
，可得到改进后的损伤模型：
[0132]
σ
i-σ
cc
＝e(ε-ε
cc
)(1-d) (σ
r-σ
cc
)d
[0133]
其中，σi表示应力；σr表示损伤岩石基质的应力；ε表示宏观应变；e表示岩石材料的弹性模量；ε
cc
表示压密点a的应变；σ
cc
为裂纹闭合应力；d表示损伤变量；
[0134]
b2.采用修正weibull分布，计算宏观应变ε的概率密度函数；根据损伤统计理论，损伤变量d采用概率密度函数求解，现有技术中基于weibull分布的统计本构模型适合于岩石材料的损伤变量d，即：
[0135][0136]
其中，m为模型形状参数；ε0为模型尺度参数。如图8为本发明实施例的现有理论本构模型与试验结果比较示意图。但根据现有技术建立的理论本构模型在峰前阶段与实验结果存在较大偏差，原因为现有技术在忽略压密阶段的基础上，将损伤起始点认为是(0,0)，即岩石材料一发生宏观压缩变形就会发生损伤，这与岩石试样在弹性阶段卸载后可以恢复应力-应变曲线的事实相冲突。
[0137]
因此，由于屈服点的应力σy为岩石塑性屈服变形的应力阈值，将屈服点 b(εy,σy)作为岩石损伤的起始点；因此，采用屈服点的应力σy修正weibull分布：
[0138][0139]
其中，ω(ε)表示宏观应变ε的概率密度函数；m为模型形状参数；ε0为模型尺度参数；(εy,σy)为屈服点b的应力-应变坐标；
[0140]
b3.通过宏观应变ε的概率密度函数计算损伤变量d，包括对宏观应变ε的概率密度函数ω(ε)进行积分，得到损伤变量d：
[0141][0142]
其中，ω(ε)表示宏观应变ε的概率密度函数；m为模型形状参数；ε0为模型尺度参数；(εy,σy)为屈服点b的应力-应变坐标；
[0143]
将损伤变量d代入步骤b1的改进后的损伤模型，得到屈服点b(εy,σy)之后的本构模型：
[0144][0145]
其中，σi表示应力；ε表示宏观应变；σr表示损伤岩石基质的应力；e表示岩石材料的
弹性模量；ε
cc
表示压密点a的应变；σ
cc
为裂纹闭合应力；m为模型形状参数；ε0为模型尺度参数；(εy,σy)为屈服点b的应力-应变坐标。
[0146]
由于本方法采用屈服点b(εy,σy)作为分界点分别建立本构模型，因此需要证明屈服点之前的本构模型和屈服点之后的本构模型能否在分界点衔接，在步骤a4已将屈服点的坐标代入了屈服点之前的本构模型，因此，仅需证明屈服点 b(εy,σy)符合屈服点之后的本构模型，将屈服点b(εy,σy)代入屈服点b(εy,σy)之后的本构模型：
[0147]
σy＝e(ε
y-ε
cc
) σ
cc
[0148]
其中，e表示岩石材料的弹性模量；ε
cc
表示压密点a的应变；σ
cc
为裂纹闭合应力；(εy,σy)为屈服点b的应力-应变坐标；如图2所示，在弹性阶段轴向应力的增量与轴向应变的增量呈线性关系，应力-应变比为弹性模量e，即：σ
i-σ
cc
＝e(ε-ε
cc
)，由于屈服点b(εy,σy)为弹性阶段的终点，代入得σ
y-σ
cc
＝e(ε
y-ε
cc
)，因此，屈服点b(εy,σy)代入屈服点b(εy,σy)之后的本构模型，即σy＝e(ε
y-ε
cc
) σ
cc
成立，因此，屈服点b(εy,σy)符合屈服点之后的本构模型。
[0149]
拟合模型形状参数m和模型尺度参数ε0，将屈服点b(εy,σy)代入屈服点 b(εy,σy)之后的本构模型：σy＝e(ε
y-ε
cc
) σ
cc
改写为：
[0150][0151]
并两边同时除以e(ε-ε
cc
)-(σ
r-σ
cc
)：
[0152][0153]
两边同时取对数：
[0154][0155]
两边同时再次取对数：
[0156][0157]
将两边看作是新变量x,y和z：
[0158][0159]
则屈服点b(εy,σy)之后的本构模型可简化为：
[0160]
y＝mx z
[0161]
然后将应力-应变数据简化后的屈服点b(εy,σy)之后的本构模型，拟合线的斜率和
截距分别为参数m和-mlnε0。因此，通过线性回归得到m和-mlnε0可以进一步得到模型形状参数m和ε0为模型尺度参数ε0的值。
[0162]
如图9为本发明实施例的a1试样的单轴压缩应力-应变曲线示意图。在本实施例中，采用bai,et al.[y.bai,r.l.shan,y.ju,y.x.wu,p.f.sun,and z.e. wang,"study on the mechanical properties and damage constitutive model of frozenweakly cemented red sandstone,"cold reg sci technol,vol.171,p.102980, 2020/03/01/2020.]单轴压缩下红砂岩a1试样的应力-应变数据对新模型进行验证。
[0163]
如图10为本发明实施例的确定a1试样相关参数示意图。首先，按照步骤 s1的方法确定压密点a和屈服点b。由弹性模量e＝1.302gpa，残余强度σr＝7.083mpa。由弹性模量可以确定参考线σ＝1302
·
ε。由此计算出参考线与轴向应力-应变曲线之间的应力差，并绘制在图10中。从图10可以直接确定水平段的两个端点：压密点a和屈服点b，由此可以确定σ
cc
＝4.254mpa，ε
cc
＝0.01058，σy＝12.61mpa、εy＝0.017；如图11为本发明实施例的参考方程交点示意图。其中参考方程计算求得 ev＝1.936mpa，γ＝0.0074。
[0164]
求得屈服点之前的模型表达式为:
[0165][0166]
如图12为本发明实施例的模型参数拟合示意图。拟合得到m＝2.147，ε0＝0.0067。因此，可以得到屈服点后的模型表达式：
[0167][0168]
得到屈服点前后的模型表达式。如图12为本发明实施例与试验结果的对比示意图。图13的对比结果表明，模型曲线与试验结果吻合较好。