CN111783282B - 基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法 - Google Patents

Abstract

基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法，涉及计量固体的变形领域。为了对压硬性非线性变化和剪缩突变特性材料的振动累积变形进行仿真，立足于循环本构模型理论及数值实现方法，执行相关循环本构模型参数的获得步骤，执行相关材料的振动累积变形的应力驱动的仿真步骤。本发明能够全面反映材料的刚度和强度随周围压力和相对密实度非线性变化的行为；能够反映材料的剪缩趋势随剪应力的增长而发生突变的特性；仿真步骤具有一阶精确性和无条件线性化稳定性；能够准确预测材料的长期累积轴向变形、剪切变形和体积变形。

Description

基于压硬性非线性变化和剪缩突变特性材料的振动累积变形 的仿真方法

技术领域

本发明涉及计量固体的变形领域，特别是基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法。

背景技术

在动力基础的上部结构或上部设备长期的循环荷载作用下，地基土会发生显著的变形累积。一旦地基土的累积变形足够大，上部结构或上部设备会产生安全性和适用性问题。为了应对材料长期累积变形带来的安全性和适用性问题，需根据材料在长期循环荷载下的发生累积变形的规律，结合循环本构模型理论及数值实现方法，对材料振动累积变形进行仿真，为进一步的加固措施提供依据。

材料的循环本构模型理论及数值实现方法立足于广义塑性力学的分量理论、非线性的屈服条件、循环累积变形的塑性硬化模型、描述体积变形的模型和本构模型数值实现方法等。

1、采用了广义塑性力学的分量理论的现有技术有：沈珠江、段建立、郑颖人、孔亮、王硕等建立的本构模型。{沈珠江.土的弹塑性应力—应变关系的合理形式[J].岩土工程学报,1980,2(2):11-19.}、 {段建立.砂土的剪胀性及其数值模拟研究[D].重庆:中国人民解放军后勤工程学院，2000.}、{Zheng Y. R,Yan D.J..Multi-yield surface model forsoil on the basis of test fitting,Computer Method and Advance inGeomechanics,1994,1(1):97-104.}、{冯嵩,郑颖人,孔亮,等.广义塑性力学多重屈服面模型隐式积分算法及其ABAQUS二次开发[J].岩石力学与工程学报,2011,30(10):2019-2025.}、{王硕,段新胜.一种考虑应力洛德角剪切变形的三维双屈服面模型的初步研究[J].土木工程与管理学报,2013,30(3):59-64.}。

2、非线性的屈服条件的现有技术有：Hoek-Brown条件、Desai模型、Lade模型、中国人民解放军后勤工程学院模型、Saniclay模型，等。{Hoek E.,Brown E.T.J..Empiricalstrength criterion for rock masses[J].Journal of the Geotechnical EngineeringDivision,1980,106(15715):1013-1035.}{Desai C.S., Somasundaram S.,Frantziskonis G..A hierarchical approach for constitutive modelling ofgeologic materials[J]. International Journal for Numerical and AnalyticalMethods in Geomechanics,1986,10(3):225-257.}{Lade P. V..Elasto-plasticstress-strain theory for Cohesionless soil with curved yield surfaces[J].International Journal of Solids and Structures,1977,13(11):1019-1035,.}{LadeP.V.,Kim M.K..Single hardening constitutive model for frictional materialsII.Yield critirion and plastic work contours[J].Computers and Geotechnics,1988,6(1), 13-29.}{郑颖人,孔亮.广义塑性力学及其运用[J].中国工程科学,2005,7(11):21-36.}{Dafalias Y.F., Manzari M.T.,Papadimitriou A.G..SANICLAY:simpleanisotropic clay plasticity model[J].International Journal for Numerical andAnalytical Methods in Geomechanics,2006,30(1):1231-1257.}。

Hoek-Brown条件的不足：由于Hoek-Brown条件需要测定完整岩石单轴抗压强度，该类模型不便应用于土壤。Hoek-Brown条件属于破坏准则，而材料通常在长期低水平循环荷载下变形，远未达到破坏，因此工程上更关心的是后继屈服准则，而非破坏准则。该模型不含背应力项，不便于与随动硬化律结合描述材料在循环荷载下的行为。此外Hoek-Brown条件未考虑相对密实度对剪切屈服面的非线性的影响。

Desai模型的不足：Desai系列模型未考虑密度对剪切屈服面的非线性的影响。该模型不含背应力项，不便于与随动硬化律结合描述材料在循环荷载下的行为。此外，子弹头形状并非所有材料的屈服面的形状，如天然的Ottawa砂密实试样的剪切后继屈服面与子弹头形状相去甚远。其剪切后继屈服应力随平均应力的增大而加速提高。

Lade模型的不足：①Lade双屈服面模型不含背应力项，不便于与随动硬化律结合描述材料在循环荷载下的行为。另外，Lade双屈服面模型未考虑密度对剪切屈服面的非线性的影响。②Lade封闭型单屈服面模型未考虑密度对剪切屈服面的非线性的影响。该模型不含背应力项，不便于与随动硬化律结合描述材料在循环荷载下的行为。此外，水滴形并非所有材料的屈服面的形状，如天然的Ottawa砂密实试样的剪切后继屈服面与水滴形相去甚远。其剪切后继屈服应力随平均应力的增大而加速提高。

中国人民解放军后勤工程学院模型的不足：无论是双曲线型剪切屈服面模型还是抛物线型剪切屈服面模型，均未考虑密度对剪切屈服面的非线性的影响。该模型不含背应力项，不便于与随动硬化律结合描述材料在循环荷载下的行为。双曲线和抛物线并非所有材料的屈服面的形状，如天然的Ottawa砂密实试样的剪切后继屈服面与双曲线和抛物线相去甚远。其剪切后继屈服应力随平均应力的增大而加速提高。

Saniclay模型的不足：Saniclay模型的环形屈服面模型未考虑密度对剪切屈服面的非线性的影响。该模型不含背应力项，不便于与随动硬化律结合描述材料在循环荷载下的行为。此外，环形并非所有材料的屈服面的形状，如天然的Ottawa砂密实试样的剪切后继屈服面与环形相去甚远。其剪切后继屈服应力随平均应力的增大而加速提高。

3、循环累积变形的塑性硬化模型的现有技术有：耦合硬化模型、边界面模型，等。{Chaboche J.L.. A review of some plasticity and viscoplasticity constitutivetheories[J].International Journal of Plasticity,2008, 24(10):.1642-1693.}{Taiebat M.,Dafalias Y.F..A Zero Elastic Range Hypoplasticity Model for Sand[J]. Lecture Notes in Applied and Computational Mechanics,2017,1(1):237-256.}。

耦合硬化模型的不足：耦合硬化模型其中的A-F随动硬化模型和Chaboche等向硬化模型未考虑材料的刚度和强度受到周围压力和相对密实度的影响。

边界面模型的不足：边界面模型其中的Sanisand模型的塑性硬化模量和密度的关系仅为线性关系，而大多数材料的塑性硬化模量和密度的关系为非线性关系，如南宁圆砾。边界面模型其中的Sanisand模型的塑性硬化模量和平均应力的关系为0.5次方函数关系，这意味着当平均应力为0 时，塑性硬化模量为0，不符合具有粘聚性的材料的性质。边界面模型的塑性硬化模量在硬化曲线的初始点处是无穷大的，这不符合的部分材料三轴试验的观测，如南宁圆砾。

4、描述体积变形的模型的现有技术有：Terzaghi、Roscoe、Wang的压缩体变模型，Bishop、Newland、 Rowe、Roscoe、张建民的剪胀模型，等。{Terzaghi K.,Peck R.B.,MesriG..Soil Mechanics in Engineering Practice[M].New York:John Wiley and sons,1996.}{Roscoe K.H.,Schofield A.N.,Thurairajah A..Yielding of clays in stateswetter than critical[J].Géotechnique,1963,13(3):211-240.}{Wang Z.L.,DafaliasY.F.,Shen C.K..Bounding surface hypoplasticity model for sand[J].Journal ofEngineering Mechanics,1990,116(5): 983-1001.}{Bishop A.W..Discussion on“Measurement of shear strengths of soils.”[J].Géotechnique,1950, 2(1):113-116.}{Newland P.L.,Allely B.H..Volume Changes in Drained Taixial Tests onGranular Materials[J]. Géotechnique,1957,7(1):17-34.}{Rowe P.W..The Stress-Dilatancy Relation for Static Equilibrium of an Assembly of Particles inContact[A].Proceedings of the Royal Society A:Mathematical,Physical andEngineering Sciences[C].London,JSTOR,1962.500-527.}{Roscoe K.H.,ThurairajahA.,Schofield A.N.. Yielding of Clays in States Wetter than Critical[J].Géotechnique,1963,13(3):211-240.}{张建民,罗刚.考虑可逆与不可逆剪胀的粗粒土动本构模型[J].岩土工程学报,2005,27(2):178-184.}。

Terzaghi、Roscoe、Wang的压缩体变模型未考虑到剪胀剪缩的因素，不足以描述材料的长期累积体变。Bishop、Newland的剪胀模型描述的是土体破坏时的剪胀量。然而材料通常是在低应力水平下累积长期变形的，远没达到破坏应力。Rowe的剪胀模型立足于单调加载的三轴压缩试验，不适合循环加载的材料。Roscoe的剪胀方程描述当剪应力为0时剪缩的趋势最强。然而部分材料在剪应力为0时的剪缩的趋势并不是最强的，如南宁圆砾。直到剪应力到达某个临界点后，剪缩的趋势才突然转强。张建民的剪胀模型需要通过扭剪试验获得参数，不适合振动三轴试验中的材料。

5、本构模型数值实现方法有向前(显式)欧拉差分法、中点欧拉差分法和向后(隐式) 欧拉差分法，等。

向前(显式)欧拉差分法是无条件不稳定的，会造成解的漂移，但计算过程简单。Dafalias的本构模型采用了向前(显式)欧拉差分法。{Dafalias Y.F.,Kourousis K.I.,Saridis G.J..Multiplicative AF kinematic hardening in plasticity[J].International Journal of Solids and Structures,2008,45(1):2861-2880.}。

向后(隐式)欧拉差分法是无条件稳定的，且具有一阶精确性。康国政的本构模型采用了向后(隐式)欧拉差分法。{康国政.循环稳定材料的棘轮行为:II.隐式应力积分算法和有限元实现[J].工程力学,2005,22(3):204-209.}。

中点欧拉差分法是无条件稳定的，且具有二阶精确性，但计算过程比其他方法复杂。周小义的本构模型采用了中点欧拉差分法。{周小义,邓安福.六面体有限覆盖的三维数值流形方法的非线性分析[J].岩土力学,2010,31(7):2276-2282.}。

发明内容

本发明要解决的技术问题是提供基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法。该方法立足于材料的循环本构模型理论及数值实现方法，即立足于广义塑性力学的分量理论、非线性的屈服条件、循环累积变形的塑性硬化模型、描述体积变形的模型和本构模型数值实现方法。该方法能够克服现有技术方案的不足，即：①能够全面反映材料的刚度和强度随周围压力和相对密实度非线性变化的行为；②能够反映材料的剪缩趋势随剪应力的增长而发生突变的特性；③仿真过程由应力驱动。该方法能够对循环荷载作用下的材料的累积变形进行仿真，比如对上部结构或上部设备的长期循环荷载作用下的地基岩土的累积变形进行仿真。为进一步的加固措施提供依据。

本发明通过以下技术方案解决上述技术问题，基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法，包括如下步骤：一、基于压硬性非线性变化和剪缩突变特性的循环本构模型参数的获得步骤；二、基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤。

在说明步骤前，需预先说明“基准条件”和“应力驱动”的含义。

一些材料，如岩土材料，其剪切硬化曲线具有随周围压力增大而整体提高的特性，且在两个不同周围压力下的后继剪切屈服应力总是保持恒定的比例关系。此外，剪切硬化曲线具有随初始相对密实度增大而整体变化的特性，且两个不同初始相对密实度试样的后继剪切屈服应力总是保持恒定的比例关系。也就是说，材料存在压硬性。根据上述现象，能够以某条剪切硬化曲线上的点为基准，按照几何相似的原则，以某种比例绘制其他周围压力条件下的剪切硬化曲线。且这个比例系数与周围压力有关。因此作为基准的剪切硬化曲线所处的周围压力称为“基准周围压力”。根据上述现象，能够以某条剪切硬化曲线上的点为基准，按照几何相似的原则，以某种比例绘制其他初始相对密实度时的剪切硬化曲线。且这个比例系数与初始相对密实度有关。因此作为基准的剪切硬化曲线的试样的初始相对密实度称为“基准相对密实度”。“基准周围压力”和“基准相对密实度”统称为“基准条件”。

所述“应力驱动”，是指仿真过程中已知应力状态，求解应变状态。“应力驱动”针对的是材料的质点的应力状态受到控制的情况。

一、基于压硬性非线性变化和剪缩突变特性的循环本构模型参数的获得步骤。

A、对材料开展至少三个不同周围压力的三轴压缩试验，记录应力、应变和孔隙水压力的数据，并获得泊松比ν，

B、对材料开展至少三个不同相对密实度的三轴压缩试验，记录应力、应变和孔隙水压力的数据，

C、对材料开展振动三轴试验，记录体变起始点的孔隙比e_ini、应力、应变和孔隙水压力的数据，

D、通过最大孔隙比试验获得最大孔隙比e_max，

E、通过最小孔隙比试验获得最小孔隙比e_min，

F、剪切屈服条件参数C_A、C_B、C_C的获得步骤，

F.a、根据相对密实度相等但处于至少3个不同周围压力的试样的三轴压缩试验应力路径的特征点，编制描述广义剪应力—塑性广义剪应变—周围压力关系的数据表格，即

关系的数据表格，选择

关系表中一个周围压力σ_r作为基准周围压力，基准周围压力为表中最接近具体实际工程的材料的周围压力的中值的周围压力，以减小预测误差，若动力基础底面的土在振密过程中的σ_r在σ_r._min和σ_r.max之间变化，则设置基准周围压力为表中最接近(σ_r.min+σ_r.max)/2的周围压力，

F.b、选择

关系表中一个塑性广义剪应变

作为参考剪切硬化内变量，即作为参考

参考

为表中最接近具体实际工程的材料的

的中值的

以减小预测误差，若动力基础底面的土在振密过程中的

在γ_min和γ_max之间变化，则忽略弹性广义剪应变，设置参考

为表中最接近(γ_min+γ_max)/2的

F.c、将在各周围压力下的试样的参考

对应的广义剪应力q代入式(1)中的q；将基准周围压力下的试样的参考

对应的q代入式(1)中的q^*；将各试验的周围压力σ_r代入式(1)，形成线性方程组，线性方程的数量与

关系表中σ_r的数量相等，

其中：q为等效剪应力，岩土工程常将其称为广义剪应力；q^*为基准条件下的试样的广义剪应力；C_A、C_B、C_C为剪切屈服条件参数，为常数，通过至少3个不同恒定周围压力的三轴压缩试验进行回归确定；σ_r为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力，其值等于试样在整体上所受到的中主应力σ₂，

F.d、用求解矛盾方程组的方法求解线性方程组，得到剪切屈服条件参数C_A、C_B、C_C，

G、剪切屈服条件参数C_D、C_E、C_F的获得步骤，

G.a、根据处于相同周围压力但相对密实度不等的至少3个试样的三轴压缩试验应力路径的特征点，编制描述广义剪应力—塑性广义剪应变—相对密实度关系的数据表格，即

关系的数据表格，选择

关系表中一个相对密实度D_r作为基准相对密实度，基准相对密实度为表中最接近具体实际工程的材料的相对密实度的中值的相对密实度，以减小预测误差，若动力基础底面的土在振密过程中的D_r在D_r.min和D_r.max之间变化，则设置基准相对密实度为表中最接近(D_r.min+D_r.max)/2的相对密实度，

G.b、在

关系表中的参考

与在

关系表中的参考

相同，

G.c、将各相对密实度的试样的参考

对应的广义剪应力q代入式(2)中的q；将基准相对密实度的试样的参考

对应的q代入式(2)中的q^*；将各试样的相对密实度D_r代入式(2)，形成线性方程组，线性方程的数量与

关系表中D_r的数量相等，

其中：C_D、C_E、C_F为剪切屈服条件参数，为常数，通过至少3个不同相对密实度的试样的三轴压缩试验回归确定；D_r为相对密实度，本发明将其取值为从塑性屈服后到弹性卸载前的一次连续塑性过程的初始相对密实度，

G.d、用求解矛盾方程组的方法求解线性方程组，得到剪切屈服条件参数C_D、C_E、C_F，

H、材料在基准条件下单调压缩时的剪切硬化曲线的初值

取值为材料在基准条件下单调压缩时的初始屈服的广义剪应力，对于岩土材料其取值小于剪切强度极限的1/100，

I、材料在基准条件下的剪切硬化曲线的初始斜率

的取值为材料在基准条件下的振动三轴试验获得的广义剪应力—轴应变偏量关系曲线在

点处的斜率，即q-e_a关系曲线在

点处的斜率，

J、材料在基准条件下单调压缩时的剪切硬化曲线的上限

的取值为材料在基准条件下的三轴压缩试验获得的q-e_a关系曲线的q上限，

K、等效的等向压缩线梯度的分段点处的广义剪应力q_seg，观察从三轴压缩试验获得的平均应力—体积应变—广义剪应力关系曲线，即p-ε_v-q关系曲线，若p-ε_v关系曲线出现明显突变，则设定突变点对应的q作为q_seg；若p-ε_v关系曲线无明显突变，则设定振动三轴试验的q的幅值的一半作为q_seg，

L、等效体变模型的参数λ_eq1和λ_eq2，

L.a、根据q_seg点的位置，对从振动三轴试验获得的平均应力—体积应变关系曲线的第1个滞回圈的上升段，即对p-ε_v关系曲线的第1个滞回圈的上升段进行分段，

L.b、采用第1个滞回圈的平均应力p较小的一段p-ε_v关系曲线的数据，对式(3)进行线性回归得到λ_eq1

e_ini-(e_ini+1)ε_v＝Γ-λ_eq1ln(p_abs.ini+p) (3)

其中：e_ini为体变起始点的孔隙比，在振动三轴试验为剪切阶段的初始孔隙比；ε_v为体积应变；λ_eq1为q≤q_seg时的等效的等向压缩线梯度；Γ为等向压缩线参数；p_abs.ini为体变起点的绝对的有效平均应力，在振动三轴试验为剪切阶段的初始有效平均应力；p为有效平均应力， p是相对于p_abs.ini而增加或减少的静水压力，

L.c、采用第1个滞回圈的平均应力p较大的一段p-ε_v关系曲线的数据，对式(4)进行线性回归得到λ_eq2

e_ini-(e_ini+1)ε_v＝Γ-λ_eq2ln(p_abs.ini+p) (4)

其中：λ_eq2为q＞q_seg时的等效的等向压缩线梯度，

M、等效体变模型的参数κ_eq，采用从振动三轴试验获得的p-ε_v关系曲线的第1个滞回圈的下降段的数据，对式(5)进行线性回归得到κ_eq

e_ini-(e_ini+1)ε_v＝Γ_κ-κ_eqln(p_abs.ini+p) (5)

其中：κ_eq为等效的等向膨胀线的梯度；Γ_κ为等向膨胀线参数，

N、剪切硬化权重系数W_sh，W_sh∈[0,1]，待上述其他参数确定后，根据振动三轴试验获得的广义剪应力—轴应变偏量关系曲线，即根据q-e_a关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_sh，

O、体积硬化权重系数W_vh，W_vh∈[0,1]，待上述其他参数确定后，根据振动三轴试验获得的p-ε_v关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_vh，

二、基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤。以下简称“仿真步骤”。

仿真步骤具体为，当循环执行增量步时，在每个增量步内依顺序执行剪切弹塑性仿真步骤和体积弹塑性仿真步骤：

A、剪切弹塑性仿真步骤

A.a、输入常数：C_A,C_B,C_C,C_D,C_E,C_F,W_sh,

ν,e_ini,κ_eq。

输入变量：σ_n,Δσ_n+1,Δγ_s,

u_n,Δu_n+1,σ_r,D_r。

A.b、为判断剪切屈服做准备：

σ′_n＝σ_n-u_n1 (6)

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

σ_n+1＝σ_n+Δσ_n+1 (9)

p_abs.n+1＝tr[σ′_n+1]/3 (10)

s_n+1＝σ′_n+1-p_abs.n+11 (11)

u_n+1＝u_n+Δu_n+1 (12)

σ′_r＝σ_r-u_n+1 (13)

若||α_s.n||＝0，则n_αs＝n_s (20)

否则n_αs＝α_s.n/||α_s.n|| (21)

若

执行步骤A.c；否则，执行步骤A.d

A.c、.发生剪切屈服时的步骤

A.c.a、确定Δγ_s，

A.c.a.a、初始化

k＝0

A.c.a.b、迭代，执行如下牛顿迭代直到

小余预设容许值，k←k+1 计算迭代

A.c.b、更新变量：若Δγ_s＜0，则取Δγ_s＝0

α_s.n+1＝ζ_M(α_s.n+2C_LΔγ_sn_s/3) (36)

K_s.n+1＝ζ_Q(K_s.n+2C_PΔγ_s/3) (38)

执行步骤A.e

A.d、不发生剪切屈服时的步骤：

执行步骤A.e

A.e、.

A.f、输出变量：σ_n+1,Δγ_s,

e_n+1。执行体积弹塑性仿真步骤。

B、体积弹塑性仿真步骤

B.a、.输入常数：λ_eq1,λ_eq2,κ_eq,e_ini,W_vh,p_abs.ini,q_seg,e_max,e_min。

输入变量：σ_n,Δσ_n+1,Δγ_v,α_v.n,K_v.n,

e_n,e_n+1,u_n,Δu_n+1。

B.b、为判断体积屈服做准备：

σ′_n＝σ_n-u_n1 (6)

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

p_n＝tr[σ′_n]/3-p_abs.ini (45)

Δp_n+1＝tr[Δσ′_n+1]/3 (46)

p_n+1＝p_n+Δp_n+1 (47)

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

p_abs.n+1＝tr[σ′_n+1]/3 (10)

s_n+1＝σ′_n+1-p_abs.n+11 (48)

若

执行步骤B.c；否则，执行步骤B.d

B.c、发生体积屈服时的步骤：

Δα_v.n+1＝(1-W_vh)(1+e_ini)(p_abs.ini+p_n+1)Δγ_vn_v/T₂ (58)

α_v.n+1＝α_v.n+Δα_v.n+1 (59)

ΔK_v.n+1＝W_vh(1+e_ini)(p_abs.ini+p_n+1)Δγ_v/T₂ (60)

K_v.n+1＝K_v.n+ΔK_v.n+1 (61)

执行步骤B.e

B.d、不发生体积屈服时的步骤：

D_r＝(e_max-e)/(e_max-e_min) (65)

执行步骤B.e

B.e、

ε_n+1＝e_n+1+ε_v.n+11/3 (67)

B.f、输出变量：ε_n+1,Δγ_v,α_v.n+1,K_v.n+1,

D_r。结束当前增量步。

以上仿真步骤中符号的含义：变量右下标_n指上一增量步；变量右下标_n+1指当前增量步；变量右上标trial指该变量是采用上一增量步的硬化参数试计算得到的；变量右上标^*指该变量处于基准条件；变量前的Δ指该变量为增量；符号:为内积符号，是对张量缩并；变量右上标^(k)指第(k)次牛顿迭代；变量右上标'指该变量为有效应力；||·||指二范数；tr[·]指对张量求迹；sign(·)为符号函数。

以上仿真步骤中的变量，粗体符号为张量，非粗体符号为标量。变量的含义：α_s为试样在实际条件下的背应力偏张量；

为试样在基准条件下的背应力偏张量；α_v为体积背应力；B_s为与周围压力有关的比例系数；C_A、C_B、C_C、C_D、C_E、C_F为剪切屈服条件参数；C_L、C_M为A-F模型的随动硬化参数；C_P、C_Q为Chaboche模型的等向硬化参数；D_r为相对密实度；D_s为与相对密实度有关的比例系数；e为应变偏张量；e^e为弹性应变偏张量；e^p为塑性应变偏张量；e为孔隙比；e_ini为体变起始点的孔隙比；e_max为最大孔隙比；e_min为最小孔隙比；ε为应变张量；ε^p为塑性应变张量；ε_v为体积应变；

为弹性体积应变；

为塑性体积应变；f_s为剪切屈服函数；f_v为体积屈服函数；G为剪切弹性模量；γ_s为剪切塑性滑移率；γ_v为体积塑性滑移率；k为迭代次数指示变量；K_s为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力的各向同性硬化部分；

为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力的各向同性硬化部分；K_v为体积等向塑性流动应力；κ_eq为等效的等向膨胀线的梯度；ξ_s为相对应力的偏张量；ξ_v为相对球应力；λ_eq1为q≤q_seg时的等效的等向压缩线梯度；λ_eq2为q＞q_seg时的等效的等向压缩线梯度；

为体积弹性模量；n_v为体积塑性流动方向；n_s为剪切塑性流动方向；n_αs为α_s方向的单位向量；ν为泊松比；o(κ_eq)为远小于κ_eq的非零正数，o(κ_eq)∈(0,κ_eq×10^-4]；p_abs为绝对的有效平均应力；p_abs.ini为体变起点的绝对的有效平均应力；p为有效平均应力，p是相对于p_abs.ini而增加或减少的静水压力；q为等效剪应力，岩土工程常将其称为广义剪应力；q_seg为等效的等向压缩线梯度的分段点处的广义剪应力；

为材料在基准条件下单调压缩时的剪切硬化曲线的初值；

为材料在基准条件下的剪切硬化曲线的初始斜率；

为材料在基准条件下单调压缩时的剪切硬化曲线的上限；s为应力偏张量；σ为应力张量；σ_r为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力；T₂为由式(50)定义的函数；u为孔隙水压力；W_sh为剪切硬化权重系数，W_sh∈[0,1]；W_vh为体积硬化权重系数，W_vh∈[0,1]；ζ_M为式(28)定义的函数；ζ_Q为式(26)定义的函数；1为二阶单位张量，1的矩阵形式表示为[1 1 1 0 0 0]^T。

符号和变量的补充说明，仿真步骤中大部分变量是由上述符号和上述变量复合而成，其含义也由各部分的含义复合而成。如

是由变量f_v、符号_n+1、符号

复合而成，因此其含义为：体积屈服函数，该变量处于当前增量步，该变量为弹性试探值。其余变量依此类推。

本发明的附带功能及其特征：

一、通过调整“基于压硬性非线性变化和剪缩突变特性的循环本构模型参数的获得步骤”，能够对金属材料的振动累积变形进行仿真：

A、对金属材料开展单轴拉伸试验，记录应力、应变的数据，并获得泊松比ν，

B、无步骤B，步骤A紧接步骤C，

C、对金属材料开展循环加载试验，记录应力、应变的数据，设定体变起始点的孔隙比e_ini＝0，

D、设定最大孔隙比e_max＝0，

E、设定最小孔隙比e_min＝0，

F、剪切屈服条件参数C_A、C_B、C_C的获得步骤，设定C_A＝0；C_B＝0；C_C＝1，

G、剪切屈服条件参数C_D、C_E、C_F的获得步骤，设定C_D＝0；C_E＝0；C_F＝1，

H、材料在基准条件下单调压缩时的剪切硬化曲线的初值

取值为金属材料的初始剪切屈服强度，

I、材料在基准条件下的剪切硬化曲线的初始斜率

的取值为金属材料在循环加载试验获得的广义剪应力—轴应变偏量关系曲线在

点处的斜率，即q-e_a关系曲线在

点处的斜率，

J、材料在基准条件下单调压缩时的剪切硬化曲线的上限

取值为金属的剪切强度极限，

K、等效的等向压缩线梯度的分段点处的广义剪应力q_seg，q_seg取值为0至剪切强度极限，

L、等效体变模型的参数λ_eq1和λ_eq2，λ_eq1和λ_eq2取值小于1×10-¹⁵且大于0，

M、等效体变模型的参数κ_eq，κ_eq取值小于1×10-¹⁵且大于0，

N、剪切硬化权重系数W_sh，W_sh∈[0,1]，待上述其他参数确定后，根据金属循环加载试验获得的 q-ε_a关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_sh，

O、体积硬化权重系数W_vh，W_vh∈[0,1]，

二、“基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤”不做调整。

本发明的原理：

一、广义塑性力学的分量理论

由于应力张量和应变张量能够分解为线性无关的球张量和偏张量，本发明应用广义塑性势理论，对塑性应变进行分解

其中：ε^p为塑性应变张量；e^p为塑性应变偏张量；

为塑性体积应变；1为二阶单位张量；γ_s为剪切塑性滑移率；γ_v为体积塑性滑移率；Q_s为剪切塑性势；Q_v为体积塑性势；s为应力偏张量；p为平均应力。本发明基于这个分解，分别在剪切方向和体积方向建立屈服面、硬化律、塑性流动向量。

二、基于压硬性非线性变化的材料循环本构模型的剪切分量

1.线弹性本构关系

本发明采用广义胡克定律描述材料的剪切弹性。广义胡克定律的应力张偏量—弹性应变偏张量关系，即s-e^e关系表述为：

e^e＝0.5s/G (69)

其中：e^e为弹性应变偏张量；s为应力偏张量；G为剪切弹性模量，根据弹性理论，表示为

其中：

为体积弹性模量；ν为泊松比。

2.非线性剪切屈服条件

包含背应力项和等向塑性流动应力项的非线性剪切屈服条件的表达式为：

其中：f_s为剪切屈服函数；s为应力偏张量；α_s为背应力偏张量；K_s为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力q的各向同性硬化部分；

为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力q的各向同性硬化部分；H_s为剪切硬化内变量，

为塑性等效剪应变，岩土工程常将其称为塑性广义剪应变；σ_r为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力，其值等于试样在整体上所受到的中主应力σ₂；D_r为相对密实度，本发明将其取值为从塑性屈服后到弹性卸载前的一次连续塑性过程的初始相对密实度；B_s为与周围压力有关的比例系数；D_s为与相对密实度有关的比例系数；C_A、C_B、 C_C为剪切屈服条件参数，为常数，通过至少3个不同恒定周围压力的三轴压缩试验进行回归确定。C_D、 C_E、C_F为剪切屈服条件参数，为常数，通过至少3个不同相对密实度的试样的三轴压缩试验回归确定。根据具体材料，比例系数B_s和D_s表示为直线函数、双曲线函数、指数函数、幂函数或对数函数。当C_A为0时，B_s退化为直线函数，能够描述剪切屈服应力随周围压力线性增长。当C_D为0时，D_s退化为直线函数，能够描述剪切屈服应力随初始相对密实度线性增长。

3.关联相对应力偏量的剪切塑性流动法则

建立循环加载本构模型时采用与相对应力偏量(s-α_s)同向的单位向量作为塑性流动方向。

4.基于耦合硬化模型的剪切硬化律

采用组合随动/等向硬化律描述材料的剪切塑性硬化，并按权重分配每个无穷小增量q的随动/等向硬化的占比，即：

其中：q为剪切硬化曲线上的等效剪应力，岩土工程常将其称为广义剪应力；W_sh为剪切硬化权重系数，W_sh∈[0,1]，根据振动三轴试验获得的广义剪应力—轴应变偏量关系，即q-e_a关系曲线确定。

采用A-F随动硬化模型和Chaboche等向硬化模型进一步描述材料的剪切硬化曲线：

其中：γ_s为剪切塑性滑移率，

C_L和C_M为A-F随动硬化参数；n_s为塑性流动方向；C_P和C_Q为Chaboche等向硬化参数。

一些材料，如岩土材料，其剪切硬化曲线的形状受到周围压力和相对密实度影响，因此A-F 随动硬化模型和Chaboche等向硬化模型的参数也会随着这些条件的变化而变化。结合材料的压硬性的非线性变化性质拓展这两个模型，得到等向硬化参数和随动硬化参数为

其中：

为材料在基准条件下单调压缩时的剪切硬化曲线的初值，对于岩土材料其取值小于剪切强度极限的1/100；

为材料在基准条件下单调压缩时的剪切硬化曲线的上限，通过相应的三轴压缩试验获得；

为材料在基准条件下的剪切硬化曲线的初始斜率，通过相应的振动三轴试验第1个滞回圈的初始上升段获得。

三、基于剪缩突变特性的圆砾循环本构模型的体积分量

1.等效体变模型

将材料的体变ε_v分解为弹性体变

和塑性体变

并采用Roscoe等提出的公式描述

和

即

其中：

为弹性体变；

为塑性体变；p_abs.ini为体变起点的绝对的有效平均应力，在振动三轴试验为剪切阶段的初始有效平均应力；e_ini为体变起始点的孔隙比，在振动三轴试验为剪切阶段的初始孔隙比；p为有效平均应力，p是相对于p_abs.ini而增加或减少的静水压力的量；κ_eq为等效的等向膨胀线的梯度；λ_eq为等效的等向压缩线的梯度。

2.体积屈服条件和体积塑性流动法则

由于p和

为标量，本发明采用一维的屈服条件描述体积屈服，即

f_v＝|p-α_v(H_v)-K_v(H_v) (80)

其中：f_v为体积屈服函数；α_v为体积背应力；K_v为体积等向塑性流动应力；H_v为体积硬化内变量，

采用关联相对平均应力(p-α_v)方向的塑性流动法则描述体积塑性流动，即

n_v＝sign(p-α_v) (81)

其中：n_v为体积塑性流动方向；sign(·)为符号函数。

3.分段梯度的体积硬化律

采用组合随动/等向硬化律描述体积塑性硬化。其中体积的组合随动/等向硬化律按权重分配每个无穷小增量p的随动/等向硬化的占比。

其中：W_vh为体积硬化权重系数，W_vh∈[0,1]。W_vh通过振动三轴试验获得的p-ε_v关系曲线确定。针对材料的剪缩趋势不连续变化的现象，本发明提出了分段梯度的体积硬化律。联立式(78)、式(79)、式(82)、式(83)、式(84)得到体积随动硬化模型和体积等向硬化模型

其中：o(κ_eq)为远小于κ_eq的非零正数，o(κ_eq)∈(0,κ_eq×10^-4]；q_seg为等效的等向压缩线梯度的分段点处的广义剪应力；λ_eq1为q≤q_seg时的等效的等向压缩线梯度；λ_eq2为q＞q_seg时的等效的等向压缩线梯度。由于剪缩突变，在较高的q水平时和较低的q水平时圆砾的p-ε_v曲线斜率存在明显差异，本发明将λ_eq分成2段，通过式(87)的前两式表示这种差异。为了描述膨胀时的包辛格效应，其中式 (87)的第三式控制着膨胀且屈服时几乎不发生塑性变形。

四、基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤中部分公式的说明

1.有效应力增量在当前增量步的形式，即式(7)

证明：由向后欧拉差分法有

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

σ_n+1＝σ_n+Δσ_n+1 (9)

u_n+1＝u_n+Δu_n+1 (12)

其中：变量右下标_n指上一增量步；变量右下标_n+1指当前增量步；变量前的Δ指该变量为增量；变量右上标'指该变量为有效应力。整理式(8)得到

Δσ′_n+1＝σ′_n+1-σ′_n (88)

由有效应力原理有

σ′_n＝σ_n-u_n1 (6)

σ′_n+1＝σ_n+1-u_n+11 (89)

其中：1为二阶单位张量，1的矩阵形式表示为[1 1 1 0 0 0]^T。

将式(6)和式(89)代入式(88)得到

Δσ′_n+1＝(σ_n+1-u_n+11)-(σ_n-u_n1)＝(σ_n+1-σ_n)-(u_n+1-u_n)1 (90)

整理式(9)和式(12)得到

σ_n+1-σ_n＝Δσ_n+1 (91)

u_n+1-u_n＝Δu_n+1 (92)

将式(91)和式(92)代入式(90)得到

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

证毕。

2.剪切背应力的向后欧拉差分形式，即式(36)和式(28)

证明：对A-F随动硬化模型式(75)的等号两边乘以时间增量Δt，有

式(93)联立

得到

式(95)联立向后欧拉差分

α_s.n+1＝α_s.n+Δα_s.n+1 (96)

得到

解方程式(97)得到

即：α_s.n+1＝ζ_M(α_s.n+2C_LΔγ_sn_s/3) (36)

其中：

证毕。

3.剪切等向塑性流动应力的向后欧拉差分形式，即式(38)和式(26)

证明：将Chaboche等向硬化模型式(76)等号两边乘以时间增量Δt，有

式(99)联立式(94)得到

式(100)联立向后欧拉差分

K_s.n+1＝K_s.n+ΔK_s.n+1 (101)

得到

解方程式(102)得到

即：K_s.n+1＝ζ_Q(K_s.n+2C_PΔγ_s/3) (38) 其中：

证毕。

4.用于判断屈服的剪切屈服条件的差分形式，即式(27)和式(26)

证明：剪切屈服条件式(70)在当前增量步的形式为

其中：||·||指二范数。式(103)代入式(104)得到

其中：

变量右上标^trial指该变量是采用上一增量步的硬化参数试计算得到的。由于其中的Δγ_s是从上一增量步传来的，该函数为试剪切屈服函数。证毕。

5.用于求解塑性滑移的剪切屈服条件的差分形式，即式(32)、式(28)、式(26)和式(18) 证明：剪切相对应力的定义为：ξ_s.n+1＝s_n+1-α_s.n+1 (105)

联立式(105)和式(98)得到

整理得到

式(107)等号两边内积径向流动向量n_s得到

其中：符号:为内积符号，是对张量缩并。这里假设了ξ_s.n+1的方向和(s_n+1-α_s.n)的方向同向。材料在三轴压缩、三轴拉伸和三轴卸载时满足π平面上的s_n+1、α_s.n+1和α_s.n在一条直线上，即这些变量都在Lode角θ＝-π/6或θ＝π/6的位置上，因此ξ_s.n+1的方向和(s_n+1-α_s.n)的方向同向，式(108) 在本发明的范围内成立。式(108)代入剪切屈服条件式(104)得到

将式(103)代入式(109)得到

即

其中：

证毕。

6.体积弹性模量在当前增量步的形式，即式(41)

证明：对式(78)的等号两边同时乘时间微分dt得到

整理(111)得到

其中：

为体积弹性模量。在当前增量步有

证毕。

7.塑性体积应变增量在当前增量步的形式，即式(55)

证明：对式(79)的等号两边乘以时间增量Δt，有

由式(113)得到

在当前增量步有

联立式(87)和式(115)得到

证毕。

8.体积背应力增量在当前增量步的形式，即式(58)

证明：将体积塑性滑移率

代入式(85)，得到

对式(117)的等号两边乘以时间增量Δt，有

由式(118)得到

在当前增量步有Δα_v.n+1＝(1-W_vh)(1+e_ini)(p_abs.ini+p_n+1)Δγ_vn_v/T₂ (58)

证毕。

9.体积等向塑性流动应力增量在当前增量步的形式，即式(60)

证明：将体积塑性滑移率

代入式(86)，得到

对式(120)的等号两边乘以时间增量Δt，有

由式(121)得到

在当前增量步有ΔK_v.n+1＝W_vh(1+e_ini)(p_abs.ini+p_n+1)Δγ_v/T₂ (60)

证毕。

10.弹性体积应变在当前增量步的形式，即式(62)

证明：对等效体变模型式(78)的等号两边积分，并应用边界条件

得到

证毕。

11.仿真步骤的其他补充说明

需要补充说明的是，仿真步骤里的式(22)的C_L比式(77)的C_L多了一个绝对值符号。这是因为，在三轴卸载时，塑性流动方向n_s发生了逆转。这时n_s的方向与n_αs的方向正好相反。这时n_s:n_αs＝-1。为了避免材料参数的正负符号的剧烈变化导致数值实现的困难，仿真步骤对C_L取绝对值。

本发明的有益效果是：

(1)非线性剪切屈服条件，即式(14)～式(27)能够在仿真过程中反映周围压力和相对密实度对剪切屈服面的非线性的影响；

(2)拓展的A-F随动硬化模型和Chaboche等向硬化模型，即式(36)～式(39)能够在仿真过程中反映塑性硬化模量随着周围压力和相对密实度的变化而非线性地变化；

(3)分段梯度的体积硬化律，即式(50)、式(58)和式(60)能够在仿真过程中准确反映剪缩趋势不连续变化的特性；

(4)基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤立足于向后欧拉差分法，具有一阶精确性和无条件(线性化)稳定性；

(5)采用本发明能够准确预测材料的长期累积轴向变形、剪切变形和体积变形。

附图说明

图1为南宁圆砾D_r＝0.5试样在三轴压缩时的广义剪应力—塑性广义剪应变

关系曲线。

图2为南宁圆砾在σ_r＝0.2MPa下三轴压缩时的广义剪应力—塑性广义剪应变

关系曲线。

图3为南宁圆砾在σ_r＝0.2MPa下三轴压缩时的平均应力—体积应变(p-ε_v)关系曲线。

图4为南宁圆砾D_r＝0.3试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数 (e_a-ε_v-ε_a-N)关系仿真曲线与试验曲线对比。

图5为南宁圆砾D_r＝0.5试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数 (e_a-ε_v-ε_a-N)关系仿真曲线与试验曲线对比。

图6为南宁圆砾D_r＝0.7试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数 (e_a-ε_v-ε_a-N)关系仿真曲线与试验曲线对比。

图7为南宁圆砾D_r＝0.5试样在q-σ_r平面上的剪切后继屈服面。

图8为南宁圆砾在σ_r＝0.2MPa时q-D_r平面上的剪切后继屈服面。

图9为SS304钢在单轴拉伸时的广义剪应力—广义剪应变

关系仿真曲线与试验曲线对比。

图10为SS304钢在循环加载时的轴应变—振动次数(ε_a-N)关系仿真曲线与试验曲线对比。

具体实施方式

以下通过实施例对本发明的技术方案作进一步说明。

实施例1

本发明所述的基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法的具体应用实例，是对南宁圆砾的振动三轴试验所测累积变形进行仿真，在每个增量步内，按顺序执行如下步骤：

一、基于压硬性非线性变化和剪缩突变特性的循环本构模型参数的获得步骤。

A、执行《土工试验规程》SL237-1999及《土工试验方法标准》GB/T 50123-1999，对材料开展至少三个不同周围压力的三轴压缩试验，记录应力、应变和孔隙水压力的数据，并获得泊松比ν，从南宁圆砾的4个不同周围压力的三轴压缩试验获得的广义剪应力—塑性广义剪应变

关系曲线见图1。南宁圆砾ν＝0.15。

B、执行《土工试验规程》SL237-1999及《土工试验方法标准》GB/T 50123-1999，对材料开展至少三个不同相对密实度的三轴压缩试验，记录应力、应变和孔隙水压力的数据，

从南宁圆砾的3个不同相对密实度试样的三轴压缩试验获得的广义剪应力—塑性广义剪应变

关系曲线和平均应力—体积应变(p-ε_v)关系曲线分别见图2和图3。

C、执行《土工试验规程》SL237-1999及《土工试验方法标准》GB/T 50123-1999，对材料开展振动三轴试验，记录体变起始点的孔隙比e_ini、应力、应变和孔隙水压力的数据。南宁圆砾D_r＝0.3的试样，e_ini＝0.6076；南宁圆砾D_r＝0.5的试样，e_ini＝0.5558；南宁圆砾D_r＝0.5的试样，e_ini＝0.5290。

从南宁圆砾D_r＝0.3试样的振动三轴试验获得的轴应变偏量—体积应变—轴应变—振动次数 (e_a-ε_v-ε_a-N)关系曲线见图4；从南宁圆砾D_r＝0.5试样的振动三轴试验获得的轴应变偏量—体积应变—轴应变—振动次数(e_a-ε_v-ε_a-N)关系曲线见图5；从南宁圆砾D_r＝0.7试样的振动三轴试验获得的轴应变偏量—体积应变—轴应变—振动次数(e_a-ε_v-ε_a-N)关系曲线见图6。

D、执行《土工试验规程》SL237-1999及《土工试验方法标准》GB/T 50123-1999，通过最大孔隙比试验获得最大孔隙比e_max。南宁圆砾e_max＝0.684。

E、执行《土工试验规程》SL237-1999及《土工试验方法标准》GB/T 50123-1999，通过最小孔隙比试验获得最小孔隙比e_min。南宁圆砾e_min＝0.411。

F、剪切屈服条件参数C_A、C_B、C_C的获得步骤。

F.a、根据相对密实度相等但处于至少3个不同周围压力的试样的三轴压缩试验应力路径的特征点，编制描述广义剪应力—塑性广义剪应变—周围压力关系的数据表格，即

关系的数据表格，选择

关系表中一个周围压力σ_r作为基准周围压力，基准周围压力为表中最接近具体实际工程的材料的周围压力的中值的周围压力，以减小预测误差，若动力基础底面的土在振密过程中的σ_r在σ_r.min和σ_r.max之间变化，则设置基准周围压力为表中最接近(σ_r.min+σ_r.max)/2的周围压力，以南宁圆砾D_r＝0.5的试样为例，在恒定周围压力下的固结排水的三轴压缩试验采样得到特征点的

关系见表1。

表1南宁圆砾D_r＝0.5试样的

关系表/MPa

表1中的数据在在q-σ_r平面上显示南宁圆砾D_r＝0.5试样剪切后继屈服面，见图7。剪切后继屈服面即图中的趋势线。

设置σ_r＝0.2MPa为基准周围压力。

F.b、选择

关系表中一个塑性广义剪应变

作为参考剪切硬化内变量，即作为参考

参考

为表中最接近具体实际工程的材料的

的中值的

以减小预测误差，若动力基础底面的土在振密过程中的

在γ_min和γ_max之间变化，则忽略弹性广义剪应变，设置参考

为表中最接近 (γ_min+γ_max)/2的

以南宁圆砾D_r＝0.5的试样为例，设置参考

F.c、将在各周围压力下的试样的参考

对应的广义剪应力q代入式(1)中的q；将基准周围压力下的试样的参考

对应的q代入式(1)中的q^*；将各试验的周围压力σ_r代入式(1)。形成线性方程组，线性方程的数量与

关系表中σ_r的数量相等。

其中：q为等效剪应力，岩土工程常将其称为广义剪应力；q^*为基准条件下的试样的广义剪应力；C_A、C_B、C_C为剪切屈服条件参数，为常数，通过至少3个不同恒定周围压力的三轴压缩试验进行回归确定；σ_r为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力，其值等于试样在整体上所受到的中主应力σ₂。

F.d、用求解矛盾方程组的方法求解线性方程组，得到剪切屈服条件参数C_A、C_B、C_C。

以南宁圆砾D_r＝0.5的试样为例，C_A＝-2.3455、C_B＝5.3433、C_C＝0.0252。

G、剪切屈服条件参数C_D、C_E、C_F的获得步骤。

G.a、根据处于相同周围压力但相对密实度不等的至少3个试样的三轴压缩试验应力路径的特征点，编制描述广义剪应力—塑性广义剪应变—相对密实度关系的数据表格，即

关系的数据表格，选择

关系表中一个相对密实度D_r作为基准相对密实度，基准相对密实度为表中最接近具体实际工程的材料的相对密实度的中值的相对密实度，以减小预测误差，若动力基础底面的土在振密过程中的D_r在D_r.min和D_r.max之间变化，则设置基准相对密实度为表中最接近(D_r.min+D_r.max)/2的相对密实度，从南宁圆砾不同相对密实度试样在0.2MPa周围压力下的固结排水的三轴压缩试验采样得到

关系见表2。

表2南宁圆砾在σ_r＝0.2MPa周围压力下的

关系表/MPa

表2中的数据在在q-D_r平面上显示南宁圆砾在σ_r＝0.2MPa时剪切后继屈服面，见图8。剪切后继屈服面即图中的趋势线。

设置D_r＝0.5为基准相对密实度。

G.b、在

关系表中的参考

与在

关系表中的参考

相同。以南宁圆砾D_r＝0.5的试样为例，设置参考

G.c、将各相对密实度的试样的参考

对应的广义剪应力q代入式(2)中的q；将基准相对密实度的试样的参考

对应的q代入式(2)中的q^*；将各试样的相对密实度D_r代入式(2)。形成线性方程组，线性方程的数量与

关系表中D_r的数量相等。

其中：C_D、C_E、C_F为剪切屈服条件参数，为常数，通过至少3个不同相对密实度的试样的三轴压缩试验回归确定；D_r为相对密实度，本发明将其取值为从塑性屈服后到弹性卸载前的一次连续塑性过程的初始相对密实度。

G.d、用求解矛盾方程组的方法求解线性方程组，得到剪切屈服条件参数C_D、C_E、C_F。

C_D＝-0.3571、C_E＝0.7143、C_F＝0.7321。

H、材料在基准条件下单调压缩时的剪切硬化曲线的初值

取值为材料在基准条件下单调压缩时的初始屈服的广义剪应力，对于岩土材料其取值小于剪切强度极限的1/100，

以南宁圆砾D_r＝0.5的试样为例，

I、材料在基准条件下的剪切硬化曲线的初始斜率

的取值为材料在基准条件下的振动三轴试验获得的广义剪应力—轴应变偏量关系曲线在

点处的斜率，即q-e_a关系曲线在

点处的斜率。以南宁圆砾D_r＝0.5的试样为例，

J、材料在基准条件下单调压缩时的剪切硬化曲线的上限

的取值为材料在基准条件下的三轴压缩试验获得的q-e_a关系曲线的q上限。以南宁圆砾D_r＝0.5的试样为例，

K、等效的等向压缩线梯度的分段点处的广义剪应力q_seg。观察从三轴压缩试验获得的平均应力—体积应变—广义剪应力关系曲线，即p-ε_v-q关系曲线。若p-ε_v关系曲线出现明显突变，则设定突变点对应的q作为q_seg；若p-ε_v关系曲线无明显突变，则设定振动三轴试验的q的幅值的一半作为 q_seg。以南宁圆砾D_r＝0.5的试样为例，q_seg＝0.055MPa。

L、等效体变模型的参数λ_eq1和λ_eq2

L.a、根据q_seg点的位置，对从振动三轴试验获得的平均应力—体积应变关系曲线的第1个滞回圈的上升段，即对p-ε_v关系曲线的第1个滞回圈的上升段进行分段。

L.b、采用第1个滞回圈的平均应力p较小的一段p-ε_v关系曲线的数据，对式(3)进行线性回归得到λ_eq1

e_ini-(e_ini+1)ε_v＝Γ-λ_eq1ln(p_abs.ini+p) (3)

其中：e_ini为体变起始点的孔隙比，在振动三轴试验为剪切阶段的初始孔隙比；ε_v为体积应变；λ_eq1为q≤q_seg时的等效的等向压缩线梯度；Γ为等向压缩线参数；p_abs.ini为体变起点的绝对的有效平均应力，在振动三轴试验为剪切阶段的初始有效平均应力；p为有效平均应力，p是相对于p_abs.ini而增加或减少的静水压力。以南宁圆砾D_r＝0.5的试样为例，λ_eq1＝9.2226×10^-4。

L.c、采用第1个滞回圈的平均应力p较大的一段p-ε_v关系曲线的数据，对式(4)进行线性回归得到λ_eq2

e_ini-(e_ini+1)ε_v＝Γ-λ_eq2ln(p_abs.ini+p) (4)

其中：λ_eq2为q＞q_seg时的等效的等向压缩线梯度。

以南宁圆砾D_r＝0.5的试样为例，λ_eq2＝1.7154×10^-3。

M、等效体变模型的参数κ_eq。采用从振动三轴试验获得的p-ε_v关系曲线的第1个滞回圈的下降段的数据，对式(5)进行线性回归得到κ_eq

e_ini-(e_ini+1)ε_v＝Γ_κ-κ_eqln(p_abs.ini+p) (5)

其中：κ_eq为等效的等向膨胀线的梯度；Γ_κ为等向膨胀线参数。

以南宁圆砾D_r＝0.5的试样为例，κ_eq＝7.6730×10^-4。

N、剪切硬化权重系数W_sh，W_sh∈[0,1]。待上述其他参数确定后，根据振动三轴试验获得的广义剪应力—轴应变偏量关系曲线，即根据q-e_a关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_sh。以南宁圆砾D_r＝0.5的试样为例，W_sh＝1.02×10^-5。

O、体积硬化权重系数W_vh，W_vh∈[0,1]。待上述其他参数确定后，根据振动三轴试验获得的p-ε_v关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_vh。以南宁圆砾D_r＝0.5的试样为例，W_vh＝0.00077。汇总南宁圆砾的循环本构模型参数，见表3。

表3南宁圆砾的循环本构模型参数

二、基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤。以下简称“仿真步骤”。

仿真步骤具体为，当循环执行增量步时，在每个增量步内依顺序执行剪切弹塑性仿真步骤和体积弹塑性仿真步骤：以南宁圆砾D_r＝0.5的试样的第8个增量步为例，

A、剪切弹塑性仿真步骤

A.a、输入常数：C_A,C_B,C_C,C_D,C_E,C_F,W_sh,

ν,e_ini,κ_eq。

输入变量：σ_n,Δσ_n+1,Δγ_s,

u_n,Δu_n+1,σ_r,D_r。

A.b、为判断剪切屈服做准备：

σ′_n＝σ_n-u_n1 (6)

σ′_n＝[0.2261 0.1926 0.1926 0 0 0]

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

Δσ′_n+1＝[0.0034 -0.0003 -0.0003 0 0 0]

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

σ′_n+1＝[0.2295 0.1923 0.1923 0 0 0]

σ_n+1＝σ_n+Δσ_n+1 (9)

σ_n+1＝[0.2365 0.1993 0.1993 0 0 0]

p_abs.n+1＝tr[σ′_n+1]/3 (10)

p_abs.n+1＝0.2047

s_n+1＝σ′_n+1-p_abs.n+11 (11)

s_n+1＝[0.0248 -0.0124 -0.0124 0 0 0]

u_n+1＝u_n+Δu_n+1 (12)

u_n+1＝0.0070

σ′_r＝σ_r-u_n+1 (13)

σ′_r＝0.1923

B_s＝0.9661

D_s＝1

α_s.n＝[0.0216 -0.0108 -0.0108 0 0 0]

K_s.n＝9.6633×10-⁴

n_s＝[0.8165 -0.4082 -0.4082 0 0 0]

若||α_s.n||＝0，则n_αs＝n_s (20)

否则n_αs＝α_s.n/||α_s.n|| (21)

n_αs＝[0.8165 -0.4082 -0.4082 0 0 0]

C_L＝212.9678

C_M＝375

C_P＝-1.0688×10-⁵

C_Q＝-1.8450

ζ_Q＝1.0000

若

执行步骤A.c；否则，执行步骤A.d

A.c、.发生剪切屈服时的步骤

A.c.a、确定Δγ_s，

A.c.a.a、初始化

k＝0

A.c.a.b、迭代，执行如下牛顿迭代直到

小余预设容许值，k←k+1 计算迭代

Δγ_s＝2.2997×10^-5

A.c.b、更新变量：若Δγ_s＜0，则取Δγ_s＝0

α_s.n+1＝ζ_M(α_s.n+2C_LΔγ_sn_s/3) (36)

α_s.n+1＝[0.0241 -0.0121 -0.0121 0 0 0]

K_s.n+1＝ζ_Q(K_s.n+2C_PΔγ_s/3) (38)

K_s.n+1＝9.6636×10^-4

执行步骤A.e

A.d、不发生剪切屈服时的步骤：

执行步骤A.e

A.e、.

G＝378.9748

e_n+1＝[0.2090 -0.1045 -0.1045 0 0 0]×10^-3

A.f、输出变量：σ_n+1,Δγ_s,

e_n+1。执行体积弹塑性仿真步骤。

B、体积弹塑性仿真步骤

B.a、.输入常数：λ_eq1,λ_eq2,κ_eq,e_ini,W_vh,p_abs.ini,q_seg,e_max,e_min。

输入变量：σ_n,Δσ_n+1,Δγ_v,α_v.n,K_v.n,

e_n,e_n+1,u_n,Δu_n+1。

B.b、为判断体积屈服做准备：

σ′_n＝σ_n-u_n1 (6)

σ′_n＝[0.2261 0.1926 0.1926 0 0 0]

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

Δσ′_n+1＝[0.0034 -0.0003 -0.0003 0 0 0]

p_n＝tr[σ′_n]|/3-p_abs.ini (45)

p_n＝0.0038

Δp_n+1＝tr[Δσ′_n+1]/3 (46)

Δp_n+1＝9.5084×10^-4

p_n+1＝p_n+Δp_n+1 (47)

p_n+1＝0.0047

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

σ′_n+1＝[0.2295 0.1923 0.1923 0 0 0]

p_abs.n+1＝tr[σ′_n+1]|/3 (10)

p_abs.n+1＝0.2047

s_n+1＝σ′_n+1-p_abs.n+11 (48)

s_n+1＝[0.0248 -0.0124 -0.0124 0 0 0]

q_n+1＝0.0372

T₂＝1.5495×10-⁴

n_v＝1

若

执行步骤B.c；否则，执行步骤B.d

B.c、发生体积屈服时的步骤：

Δγ_v＝4.6261×10^-7

Δα_v.n+1＝(1-W_vh)(1+e_ini)(p_abs.ini+p_n+1)Δγ_vn_v/T₂ (58)

Δα_v.n+1＝9.5010×10^-4

α_v.n+1＝α_v.n+Δα_v.n+1 (59)

α_v.n+1＝0.0056

ΔK_v.n+1＝W_vh(1+e_ini)(p_abs.ini+p_n+1)Δγ_v/T₂ (60)

ΔK_v.n+1＝7.3214×10^-7

K_v.n+1＝K_v.n+ΔK_v.n+1 (61)

K_v.n+1＝4.2944×10^-6

执行步骤B.e

B.d、不发生体积屈服时的步骤：

D_r＝(e_max-e)/(e_max-e_min) (65)

执行步骤B.e

B.e、

ε_v.n+1＝1.6215×10^-5

ε_n+1＝e_n+1+ε_v.n+11/3 (67)

ε_n+1＝[0.2144 -0.0991 -0.0991 0 0 0]×10^-3

B.f、输出变量：ε_n+1,Δγ_v,α_v.n+1,K_v.n+1,

D_r。结束当前增量步。

以上仿真步骤中符号的含义：变量右下标_n指上一增量步；变量右下标_n+1指当前增量步；变量右上标trial指该变量是采用上一增量步的硬化参数试计算得到的；变量右上标^*指该变量处于基准条件；变量前的Δ指该变量为增量；符号:为内积符号，是对张量缩并；变量右上标^(k)指第(k)次牛顿迭代；变量右上标'指该变量为有效应力；||·||指二范数；tr[·]指对张量求迹；sign(·)为符号函数。

以上仿真步骤中的变量，粗体符号为张量，非粗体符号为标量。变量的含义：α_s为试样在实际条件下的背应力偏张量；

为试样在基准条件下的背应力偏张量；α_v为体积背应力；B_s为与周围压力有关的比例系数；C_A、C_B、C_C、C_D、C_E、C_F为剪切屈服条件参数；C_L、C_M为A-F模型的随动硬化参数；C_P、C_Q为Chaboche模型的等向硬化参数；D_r为相对密实度；D_s为与相对密实度有关的比例系数；e为应变偏张量；e^e为弹性应变偏张量；e^p为塑性应变偏张量；e为孔隙比；e_ini为体变起始点的孔隙比；e_max为最大孔隙比；e_min为最小孔隙比；ε为应变张量；ε^p为塑性应变张量；ε_v为体积应变；

为弹性体积应变；

为塑性体积应变；f_s为剪切屈服函数；f_v为体积屈服函数；G为剪切弹性模量；γ_s为剪切塑性滑移率；γ_v为体积塑性滑移率；k为迭代次数指示变量；K_s为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力的各向同性硬化部分；

为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力的各向同性硬化部分；K_v为体积等向塑性流动应力；κ_eq为等效的等向膨胀线的梯度；ξ_s为相对应力的偏张量；ξ_v为相对球应力；λ_eq1为q≤q_seg时的等效的等向压缩线梯度；λ_eq2为q＞q_seg时的等效的等向压缩线梯度；

为体积弹性模量；n_v为体积塑性流动方向；n_s为剪切塑性流动方向；n_αs为α_s方向的单位向量；ν为泊松比；o(κ_eq)为远小于κ_eq的非零正数，o(κ_eq)∈(0,κ_eq×10^-4]；p_abs为绝对的有效平均应力；p_abs.ini为体变起点的绝对的有效平均应力；p为有效平均应力，p是相对于p_abs.ini而增加或减少的静水压力；q为等效剪应力，岩土工程常将其称为广义剪应力；q_seg为等效的等向压缩线梯度的分段点处的广义剪应力；

为材料在基准条件下单调压缩时的剪切硬化曲线的初值；

为材料在基准条件下的剪切硬化曲线的初始斜率；

为材料在基准条件下单调压缩时的剪切硬化曲线的上限；s为应力偏张量；σ为应力张量；σ_r为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力；T₂为由式(50)定义的函数；u为孔隙水压力；W_sh为剪切硬化权重系数，W_sh∈[0,1]；W_vh为体积硬化权重系数，W_vh∈[0,1]；ζ_M为式(28)定义的函数；ζ_Q为式(26)定义的函数；1为二阶单位张量，1的矩阵形式表示为[1 1 1 0 0 0]^T。

符号和变量的补充说明，仿真步骤中大部分变量是由上述符号和上述变量复合而成，其含义也由各部分的含义复合而成。如

是由变量f_v、符号_n+1、符号

复合而成，因此其含义为：体积屈服函数，该变量处于当前增量步，该变量为弹性试探值。其余变量依此类推。

对南宁圆砾在4个不同周围压力下三轴压缩时的广义剪应力—塑性广义剪应变

关系曲线的仿真见图1；对南宁圆砾的3个不同相对密实度试样在三轴压缩时的广义剪应力—塑性广义剪应变

关系曲线和平均应力—体积应变(p-ε_v)关系曲线的仿真分别见图2和图3；对南宁圆砾D_r＝0.3试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数(e_a-ε_v-ε_a-N)关系曲线的仿真见图4；对南宁圆砾D_r＝0.5试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数(e_a-ε_v-ε_a-N)关系曲线的仿真见图5；对南宁圆砾D_r＝0.7 试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数(e_a-ε_v-ε_a-N)关系曲线的仿真见图6。因此采用本发明能够准确预测材料的长期累积轴向变形、剪切变形和体积变形。

从图7和图8可见，本发明能够全面反映材料的强度随周围压力和相对密实度非线性变化的行为，即能够在仿真过程中反映圆砾的剪切屈服面随周围压力和相对密实度的变化而非线性变化；从图1和图2可见，本发明能够全面反映材料的刚度随周围压力和相对密实度非线性变化的行为，即能够在仿真过程中反映圆砾的塑性硬化模量随着周围压力和相对密实度的变化而非线性地变化；从图3可见，本发明能够反映材料的剪缩趋势随剪应力的增长而发生突变的特性，即能够在仿真过程中准确反映圆砾的体积应变曲线的斜率不连续变化的特性。

实施例2

本发明所述的基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法的具体应用实例，是对SS304钢的循环加载试验所测累积变形进行仿真，在每个增量步内，按顺序执行如下步骤：

一、基于压硬性非线性变化和剪缩突变特性的循环本构模型参数的获得步骤。

A、执行《金属材料拉伸试验第1部分：室温试验方法》GB/T 228.1-2010及《实验力学》，对金属材料开展单轴拉伸试验，记录应力、应变的数据，并获得泊松比ν，SS304钢ν＝0.3，

从SS304钢的单轴拉伸试验获得的广义剪应力—广义剪应变

关系曲线见图9，

B、无步骤B，步骤A紧接步骤C，

C、执行《金属材料拉伸试验第1部分：室温试验方法》GB/T 228.1-2010及《实验力学》，对金属材料开展循环加载试验，记录应力、应变的数据，设定体变起始点的孔隙比e_ini＝0，

从SS304钢的循环加载试验获得的轴应变—振动次数(ε_a-N)关系曲线见图10，

D、设定最大孔隙比e_max＝0，

E、设定最小孔隙比e_min＝0，

F、剪切屈服条件参数C_A、C_B、C_C的获得步骤，设定C_A＝0；C_B＝0；C_C＝1，

G、剪切屈服条件参数C_D、C_E、C_F的获得步骤，设定C_D＝0；C_E＝0；C_F＝1，

H、材料在基准条件下单调压缩时的剪切硬化曲线的初值

取值为金属材料的初始剪切屈服强度，以SS304钢为例，

I、材料在基准条件下的剪切硬化曲线的初始斜率

的取值为金属材料在循环加载试验获得的广义剪应力—轴应变偏量关系曲线在

点处的斜率，即q-e_a关系曲线在

点处的斜率，以SS304钢为例，

J、材料在基准条件下单调压缩时的剪切硬化曲线的上限

取值为金属的剪切强度极限，以SS304钢为例，

K、等效的等向压缩线梯度的分段点处的广义剪应力q_seg，q_seg取值为0至剪切强度极限，以SS304钢为例，q_seg＝0。

L、等效体变模型的参数λ_eq1和λ_eq2，λ_eq1和λ_eq2取值小于1×10^-15且大于0，

以SS304钢为例，λ_eq1＝λ_eq2＝2.2204×10^-16，

M、等效体变模型的参数κ_eq，κ_eq取值小于1×10^-15且大于0，以SS304钢为例，κ_eq＝2.2204×10^-16，

N、剪切硬化权重系数W_sh，W_sh∈[0,1]，待上述其他参数确定后，根据金属循环加载试验获得的q-ε_a关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_sh，以SS304钢为例，W_sh＝0.009，

O、体积硬化权重系数W_vh，W_vh∈[0,1]，以SS304钢为例，W_vh＝1。

汇总SS304钢的循环本构模型参数，见表4。

表4 SS304钢的循环本构模型参数

二、基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤。以下简称“仿真步骤”。

仿真步骤具体为，当循环执行增量步时，在每个增量步内依顺序执行剪切弹塑性仿真步骤和体积弹塑性仿真步骤：以SS304钢的第8个增量步为例，

A、剪切弹塑性仿真步骤

A.a、输入常数：C_A,C_B,C_C,C_D,C_E,C_F,W_sh,

ν,e_ini,κ_eq。

输入变量：σ_n,Δσ_n+1,Δγ_s,

u_n,Δu_n+1,σ_r,D_r。

A.b、为判断剪切屈服做准备：

σ′_n＝σ_n-u_n1 (6)

σ′_n＝[78 0 0 0 0 0]

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

Δσ′_n+1＝[24.4597 0 0 0 0 0]

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

σ′_n+1＝[102.4597 0 0 0 0 0]

σ_n+1＝σ_n+Δσ_n+1 (9)

σ_n+1＝[102.4597 0 0 0 0 0]

p_abs.n+1＝tr[σ′_n+1]/3 (10)

p_abs.n+1＝34.1532

s_n+1＝σ′_n+1-p_abs.n+11 (11)

s_n+1＝[68.3064 0 0 0 0 0]

u_n+1＝u_n+Δu_n+1 (12)

u_n+1＝0

σ′_r＝σ_r-u_n+1 (13)

σ′_r＝0

B_s＝1

D_s＝1

α_s.n＝[79.3794 0 0 0 0 0]

K_s.n＝343.8649×10^-4

n_s＝[1.0000 0 0 0 0 0]

若||α_s.n||＝0，则n_αs＝n_s (20)

否则n_αs＝α_s.n/||α_s.n|| (21)

n_αs＝[1.0000 0 0 0 0 0]

C_L＝9.5760×10⁴

C_M＝150.3615

C_P＝0

C_Q＝0

ζ_Q＝1.0000

若

执行步骤A.c；否则，执行步骤A.d

A.c、.发生剪切屈服时的步骤

A.c.a、确定Δγ_s，

A.c.a.a、初始化

A.c.a.b、迭代，执行如下牛顿迭代直到

小余预设容许值，k←k+1

计算迭代

Δγ_s＝6.6548×10^-9

A.c.b、更新变量：若Δγ_s＜0，则取Δγ_s＝0

α_s.n+1＝ζ_M(α_s.n+2C_LΔγ_sn_s/3) (36)

α_s.n+1＝[84.0745 0 0 0 0 0]

K_s.n+1＝ζ_Q(K_s.n+2C_PΔγ_s/3) (38)

K_s.n+1＝351.7456×10^-4

执行步骤A.e

A.d、不发生剪切屈服时的步骤：

执行步骤A.e

A.e、.

G＝78188

e_n+1＝[0.8103 -0.4051 -0.4051 0 0 0]×10^-9

A.f、输出变量：σ_n+1,Δγ_s,

e_n+1。执行体积弹塑性仿真步骤。

B、体积弹塑性仿真步骤

B.a、.输入常数：λ_eq1,λ_eq2,κ_eq,e_ini,W_vh,p_abs.ini,q_seg,e_max,e_min。

输入变量：σ_n,Δσ_n+1,Δγ_v,α_v.n,K_v.n,

e_n,e_n+1,u_n,Δu_n+1。

B.b、为判断体积屈服做准备：

σ′_n＝σ_n-u_n1 (6)

σ′_n＝[234.5766 0 0 0 0 0]

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

Δσ′_n+1＝[17.3193 0 0 0 0 0]

p_n＝tr[σ′_n]/3-p_abs.ini (45)

p_n＝78.1912

Δp_n+1＝tr[Δσ′_n+1]/3 (46)

Δp_n+1＝5.7731

p_n+1＝p_n+Δp_n+1 (47)

p_n+1＝83.9643

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

σ′_n+1＝[251.8959 0 0 0 0 0]

p_abs.n+1＝tr[σ′_n+1]/3 (10)

p_abs.n+1＝83.9653

s_n+1＝σ′_n+1-p_abs.n+11 (48)

s_n+1＝[167.9306 -83.9653 -83.9653 0 0 0]

q_n+1＝251.8959

T₂＝9.4809×10^-10

n_v＝1

若

执行步骤B.c；否则，执行步骤B.d B.c、发生体积屈服时的步骤：

Δα_v.n+1＝(1-W_vh)(1+e_ini)(p_abs.ini+p_n+1)Δγ_vn_v/T₂ (58)

α_v.n+1＝α_v.n+Δα_v.n+1 (59)

ΔK_v.n+1＝W_vh(1+e_ini)(p_abs.ini+p_n+1)Δγ_v/T₂ (60)

K_v.n+1＝K_v.n+ΔK_v.n+1 (61)

执行步骤B.e

B.d、不发生体积屈服时的步骤：

Δγ_v＝0,α_v.n+1＝0,K_v.n+1＝∞

e＝0

D_r＝(e_max-e)/(e_max-e_min) (65)

D_r＝1

执行步骤B.e

B.e、

ε_v.n+1＝5.5919×10^-9

ε_n+1＝e_n+1+ε_v.n+11/3 (67)

ε_n+1＝[0.2674 0.1459 0.1459 0 0 0]×10^-8

B.f、输出变量：ε_n+1,Δγ_v,α_v.n+1,K_v.n+1,

D_r。结束当前增量步。

以上仿真步骤中符号的含义：变量右下标_n指上一增量步；变量右下标_n+1指当前增量步；变量右上标trial指该变量是采用上一增量步的硬化参数试计算得到的；变量右上标^*指该变量处于基准条件；变量前的Δ指该变量为增量；符号:为内积符号，是对张量缩并；变量右上标^(k)指第(k)次牛顿迭代；变量右上标'指该变量为有效应力；||·||指二范数；tr[·]指对张量求迹；sign(·)为符号函数。

以上仿真步骤中的变量，粗体符号为张量，非粗体符号为标量。变量的含义：α_s为试样在实际条件下的背应力偏张量；

为试样在基准条件下的背应力偏张量；α_v为体积背应力；B_s为与周围压力有关的比例系数；C_A、C_B、C_C、C_D、C_E、C_F为剪切屈服条件参数；C_L、C_M为A-F模型的随动硬化参数；C_P、C_Q为Chaboche模型的等向硬化参数；D_r为相对密实度；D_s为与相对密实度有关的比例系数；e为应变偏张量；e^e为弹性应变偏张量；e^p为塑性应变偏张量；e为孔隙比；e_ini为体变起始点的孔隙比；e_max为最大孔隙比；e_min为最小孔隙比；ε为应变张量；ε^p为塑性应变张量；ε_v为体积应变；

为弹性体积应变；

为塑性体积应变；f_s为剪切屈服函数；f_v为体积屈服函数；G为剪切弹性模量；γ_s为剪切塑性滑移率；γ_v为体积塑性滑移率；k为迭代次数指示变量；K_s为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力的各向同性硬化部分；

为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力的各向同性硬化部分；K_v为体积等向塑性流动应力；κ_eq为等效的等向膨胀线的梯度；ξ_s为相对应力的偏张量；ξ_v为相对球应力；λ_eq1为q≤q_seg时的等效的等向压缩线梯度；λ_eq2为q＞q_seg时的等效的等向压缩线梯度；

为体积弹性模量；n_v为体积塑性流动方向；n_s为剪切塑性流动方向；n_αs为α_s方向的单位向量；ν为泊松比；o(κ_eq)为远小于κ_eq的非零正数，o(κ_eq)∈(0,κ_eq×10^-4]；p_abs为绝对的有效平均应力；p_abs.ini为体变起点的绝对的有效平均应力；p为有效平均应力，p是相对于p_abs.ini而增加或减少的静水压力；q为等效剪应力，岩土工程常将其称为广义剪应力；q_seg为等效的等向压缩线梯度的分段点处的广义剪应力；

为材料在基准条件下单调压缩时的剪切硬化曲线的初值；

为材料在基准条件下的剪切硬化曲线的初始斜率；

为材料在基准条件下单调压缩时的剪切硬化曲线的上限；s为应力偏张量；σ为应力张量；σ_r为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力；T₂为由式(50)定义的函数；u为孔隙水压力；W_sh为剪切硬化权重系数，W_sh∈[0,1]；W_vh为体积硬化权重系数，W_vh∈[0,1]；ζ_M为式(28)定义的函数；ζ_Q为式(26)定义的函数；1为二阶单位张量，1的矩阵形式表示为[1 1 1 0 0 0]^T。

符号和变量的补充说明，仿真步骤中大部分变量是由上述符号和上述变量复合而成，其含义也由各部分的含义复合而成。如

是由变量f_v、符号_n+1、符号

复合而成，因此其含义为：体积屈服函数，该变量处于当前增量步，该变量为弹性试探值。其余变量依此类推。

对SS304钢在单轴拉伸时的广义剪应力—广义剪应变

关系曲线的仿真见图9；对 SS304钢在循环加载时的轴应变—振动次数(ε_a-N)关系曲线的仿真见图10。因此采用本发明能够准确预测材料的长期累积变形。

本发明所述“压硬性非线性变化”并不局限于抛物线变化，还指双曲线函数、指数函数、幂函数或对数函数变化。本发明也涵盖压硬性线性变化。本发明所述“剪缩突变”并不局限于体积应变曲线的一处突变，也指多处突变。只要是体积应变曲线不连续变化，都是本发明涵盖的范围。本发明也涵盖体积应变曲线无突变。

Claims (2)

1.基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法，包括基于压硬性非线性变化和剪缩突变特性的循环本构模型参数的获得步骤、基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤，其特征在于：

基于压硬性非线性变化和剪缩突变特性的循环本构模型参数的获得步骤，

A、对材料开展至少三个不同周围压力的三轴压缩试验，记录应力、应变和孔隙水压力的数据，并获得泊松比ν，

B、对材料开展至少三个不同相对密实度的三轴压缩试验，记录应力、应变和孔隙水压力的数据，

C、对材料开展振动三轴试验，记录体变起始点的孔隙比e_ini、应力、应变和孔隙水压力的数据，

D、通过最大孔隙比试验获得最大孔隙比e_max，

E、通过最小孔隙比试验获得最小孔隙比e_min，

F、剪切屈服条件参数C_A、C_B、C_C的获得步骤，

F.a、根据相对密实度相等但处于至少3个不同周围压力的试样的三轴压缩试验应力路径的特征点，编制描述广义剪应力—塑性广义剪应变—周围压力关系的数据表格，即

关系的数据表格，选择

关系表中一个周围压力σ_r作为基准周围压力，基准周围压力为表中最接近具体实际工程的材料的周围压力的中值的周围压力，以减小预测误差，若动力基础底面的土在振密过程中的σ_r在σ_r.min和σ_r.max之间变化，则设置基准周围压力为表中最接近(σ_r.min+σ_r.max)/2的周围压力，

F.b、选择

关系表中一个塑性广义剪应变

作为参考剪切硬化内变量，即作为参考

参考

为表中最接近具体实际工程的材料的

的中值的

以减小预测误差，若动力基础底面的土在振密过程中的

在γ_min和γ_max之间变化，则忽略弹性广义剪应变，设置参考

为表中最接近(γ_min+γ_max)/2的

F.c、将在各周围压力下的试样的参考

对应的广义剪应力q代入式(1)中的q；将基准周围压力下的试样的参考

对应的q代入式(1)中的q^*；将各试验的周围压力σ_r代入式(1)，形成线性方程组，线性方程的数量与

关系表中σ_r的数量相等，

其中：q为等效剪应力，岩土工程常将其称为广义剪应力；q^*为基准条件下的试样的广义剪应力；C_A、C_B、C_C为剪切屈服条件参数，为常数，通过至少3个不同恒定周围压力的三轴压缩试验进行回归确定；σ_r为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力，其值等于试样在整体上所受到的中主应力σ₂，

F.d、用求解矛盾方程组的方法求解线性方程组，得到剪切屈服条件参数C_A、C_B、C_C，

G、剪切屈服条件参数C_D、C_E、C_F的获得步骤，

G.a、根据处于相同周围压力但相对密实度不等的至少3个试样的三轴压缩试验应力路径的特征点，编制描述广义剪应力—塑性广义剪应变—相对密实度关系的数据表格，即

关系的数据表格，选择

关系表中一个相对密实度D_r作为基准相对密实度，基准相对密实度为表中最接近具体实际工程的材料的相对密实度的中值的相对密实度，以减小预测误差，若动力基础底面的土在振密过程中的D_r在D_r.min和D_r.max之间变化，则设置基准相对密实度为表中最接近(D_r.min+D_r.max)/2的相对密实度，

G.b、在

关系表中的参考

与在

关系表中的参考

相同，

G.c、将各相对密实度的试样的参考

对应的广义剪应力q代入式(2)中的q；将基准相对密实度的试样的参考

对应的q代入式(2)中的q^*；将各试样的相对密实度D_r代入式(2)，形成线性方程组，线性方程的数量与

关系表中D_r的数量相等，

其中：C_D、C_E、C_F为剪切屈服条件参数，为常数，通过至少3个不同相对密实度的试样的三轴压缩试验回归确定；D_r为相对密实度，将其取值为从塑性屈服后到弹性卸载前的一次连续塑性过程的初始相对密实度，

G.d、用求解矛盾方程组的方法求解线性方程组，得到剪切屈服条件参数C_D、C_E、C_F，

H、材料在基准条件下单调压缩时的剪切硬化曲线的初值

取值为材料在基准条件下单调压缩时的初始屈服的广义剪应力，对于岩土材料其取值小于剪切强度极限的1/100，

I、材料在基准条件下的剪切硬化曲线的初始斜率

的取值为材料在基准条件下的振动三轴试验获得的广义剪应力—轴应变偏量关系曲线在

点处的斜率，即q-e_a关系曲线在

点处的斜率，

J、材料在基准条件下单调压缩时的剪切硬化曲线的上限

的取值为材料在基准条件下的三轴压缩试验获得的q-e_a关系曲线的q上限，

K、等效的等向压缩线梯度的分段点处的广义剪应力q_seg，观察从三轴压缩试验获得的平均应力—体积应变—广义剪应力关系曲线，即p-ε_v-q关系曲线，若p-ε_v关系曲线出现明显突变，则设定突变点对应的q作为q_seg；若p-ε_v关系曲线无明显突变，则设定振动三轴试验的q的幅值的一半作为q_seg，

L、等效体变模型的参数λ_eq1和λ_eq2，

L.a、根据q_seg点的位置，对从振动三轴试验获得的平均应力—体积应变关系曲线的第1个滞回圈的上升段，即对p-ε_v关系曲线的第1个滞回圈的上升段进行分段，

L.b、采用第1个滞回圈的平均应力p较小的一段p-ε_v关系曲线的数据，对式(3)进行线性回归得到λ_eq1

e_ini-(e_ini+1)ε_v＝Γ-λ_eq1ln(p_abs.ini+p) (3)

其中：e_ini为体变起始点的孔隙比，在振动三轴试验为剪切阶段的初始孔隙比；ε_v为体积应变；λ_eq1为q≤q_seg时的等效的等向压缩线梯度；Γ为等向压缩线参数；p_abs.ini为体变起点的绝对的有效平均应力，在振动三轴试验为剪切阶段的初始有效平均应力；p是相对于p_abs.ini而增加或减少的静水压力，

L.c、采用第1个滞回圈的平均应力p较大的一段p-ε_v关系曲线的数据，对式(4)进行线性回归得到λ_eq2

e_ini-(e_ini+1)ε_v＝Γ-λ_eq2ln(p_abs.ini+p) (4)

其中：λ_eq2为q＞q_seg时的等效的等向压缩线梯度，

M、等效体变模型的参数κ_eq，采用从振动三轴试验获得的p-ε_v关系曲线的第1个滞回圈的下降段的数据，对式(5)进行线性回归得到κ_eq

e_ini-(e_ini+1)ε_v＝Γ_κ-κ_eqln(p_abs.ini+p) (5)

其中：κ_eq为等效的等向膨胀线的梯度；Γ_κ为等向膨胀线参数，

N、剪切硬化权重系数W_sh，W_sh∈[0,1]，待上述其他参数确定后，根据振动三轴试验获得的广义剪应力—轴应变偏量关系曲线，即根据q-e_a关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_sh，

O、体积硬化权重系数W_vh，W_vh∈[0,1]，待上述其他参数确定后，根据振动三轴试验获得的p-ε_v关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_vh，

基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤，以下简称“仿真步骤”，

仿真步骤具体为，当循环执行增量步时，在每个增量步内依顺序执行剪切弹塑性仿真步骤和体积弹塑性仿真步骤：

A、剪切弹塑性仿真步骤

A.a、输入常数：C_A,C_B,C_C,C_D,C_E,C_F,W_sh,

v,e_ini,κ_eq，

输入变量：σ_n,Δσ_n+1,Δγ_s,

u_n,Δu_n+1,σ_r,D_r，

A.b、为判断剪切屈服做准备：

σ′_n＝σ_n-u_n1 (6)

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

σ_n+1＝σ_n+Δσ_n+1 (9)

p_abs.n+1＝tr[σ′_n+1]/3 (10)

s_n+1＝σ′_n+1-p_abs.n+11 (11)

u_n+1＝u_n+Δu_n+1 (12)

σ_r′＝σ_r-u_n+1 (13)

若||α_s.n||＝0，则n_αs＝n_s (20)

否则n_αs＝α_s.n/||α_s.n|| (21)

若

执行步骤A.c；否则，执行步骤A.d

A.c、发生剪切屈服时的步骤

A.c.a、确定Δγ_s，

A.c.a.a、初始化

k＝0

A.c.a.b、迭代，执行如下牛顿迭代直到

小于预设容许值，k←k+1计算迭代

A.c.b、更新变量：若Δγ_s＜0，则取Δγ_s＝0

α_s.n+1＝ζ_M(α_s.n+2C_LΔγ_sn_s/3) (36)

K_s.n+1＝ζ_Q(K_s.n+2C_PΔγ_s/3) (38)

执行步骤A.e

A.d、不发生剪切屈服时的步骤：

执行步骤A.e

A.e、

A.f、输出变量：σ_n+1,Δγ_s,

e_n+1，执行体积弹塑性仿真步骤，

B、体积弹塑性仿真步骤

B.a、输入常数：λ_eq1,λ_eq2,k_eq,e_ini,W_vh,p_abs.ini,q_seg,e_max,e_min，

输入变量：σ_n,Δσ_n+1,Δγ_v,a_v.n,K_v.n,

e_n,e_n+1,u_n,Δu_n+1，

B.b、为判断体积屈服做准备：

σ′_n＝σ_n-u_n1 (6)

Δσ′_n+1＝Δσ_n+1-Δu_n+11 (7)

p_n＝tr[σ′_n]/3-p_abs.ini (45)

Δp_n+1＝tr[Δσ′_n+1]/3 (46)

p_n+1＝p_n+Δp_n+1 (47)

σ′_n+1＝σ′_n+Δσ′_n+1 (8)

p_abs.n+1＝tr[σ′_n+1]/3 (10)

s_n+1＝σ′_n+1-p_abs.n+11 (48)

若

执行步骤B.c；否则，执行步骤B.d

B.c、发生体积屈服时的步骤：

Δα_v.n+1＝(1-W_vh)(1+e_ini)(p_abs.ini+p_n+1)Δγ_vn_v/T₂ (58)

α_v.n+1＝α_v.n+Δα_v.n+1 (59)

ΔK_v.n+1＝W_vh(1+e_ini)(p_abs.ini+p_n+1)Δγ_v/T₂ (60)

K_v.n+1＝K_v.n+ΔK_v.n+1 (61)

执行步骤B.e

B.d、不发生体积屈服时的步骤：

D_r＝(e_max-e)/(e_max-e_min) (65)

执行步骤B.e

B.e、

ε_n+1＝e_n+1+ε_v.n+11/3 (67)

B.f、输出变量：ε_n+1,Δγ_v,α_v.n+1,K_v.n+1,

D_r，结束当前增量步，

以上仿真步骤中符号的含义：变量右下标_n指上一增量步；变量右下标_n+1指当前增量步；变量右上标^trial指该变量是采用上一增量步的硬化参数试计算得到的；变量右上标^*指该变量处于基准条件；变量前的Δ指该变量为增量；符号:为内积符号，是对张量缩并；变量右上标^(k)指第(k)次牛顿迭代；变量右上标'指该变量为有效应力；||·||指二范数；tr[·]指对张量求迹；sign(·)为符号函数，

以上仿真步骤中的变量，粗体符号为张量，非粗体符号为标量，变量的含义：α_s为试样在实际条件下的背应力偏张量；

为试样在基准条件下的背应力偏张量；α_v为体积背应力；B_s为与周围压力有关的比例系数；C_A、C_B、C_C、C_D、C_E、C_F为剪切屈服条件参数；C_L、C_M为A-F模型的随动硬化参数；C_P、C_Q为Chaboche模型的等向硬化参数；D_r为相对密实度；D_s为与相对密实度有关的比例系数；e为应变偏张量；e^e为弹性应变偏张量；e^p为塑性应变偏张量；e为孔隙比；e_ini为体变起始点的孔隙比；e_max为最大孔隙比；e_min为最小孔隙比；ε为应变张量；ε^p为塑性应变张量；ε_v为体积应变；

为弹性体积应变；

为塑性体积应变；f_s为剪切屈服函数；f_v为体积屈服函数；G为剪切弹性模量；γ_s为剪切塑性滑移率；γ_v为体积塑性滑移率；k为迭代次数指示变量；K_s为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力的各向同性硬化部分；

为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力的各向同性硬化部分；K_v为体积等向塑性流动应力；k_eq为等效的等向膨胀线的梯度；ξ_s为相对应力的偏张量；ξ_v为相对球应力；λ_eq1为q≤q_seg时的等效的等向压缩线梯度；λ_eq2为q＞q_seg时的等效的等向压缩线梯度；

为体积弹性模量；n_v为体积塑性流动方向；n_s为剪切塑性流动方向；n_αs为α_s方向的单位向量；v为泊松比；o(k_eq)为远小于k_eq的非零正数，o(k_eq)∈(0,k_eq×10^-4]；p_abs为绝对的有效平均应力；p_abs.ini为体变起点的绝对的有效平均应力；p是相对于p_abs.ini而增加或减少的静水压力；q为等效剪应力，岩土工程常将其称为广义剪应力；q_seg为等效的等向压缩线梯度的分段点处的广义剪应力；

为材料在基准条件下单调压缩时的剪切硬化曲线的初值；

为材料在基准条件下的剪切硬化曲线的初始斜率；

为材料在基准条件下单调压缩时的剪切硬化曲线的上限；s为应力偏张量；σ为应力张量；σ_r为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力；T₂为由式(50)定义的函数；u为孔隙水压力；W_sh为剪切硬化权重系数，W_sh∈[0,1]；W_vh为体积硬化权重系数，W_vh∈[0,1]；ζ_M为式(28)定义的函数；ζ_Q为式(26)定义的函数；1为二阶单位张量，1的矩阵形式表示为[1 1 1 0 0 0]^T，

符号和变量的补充说明，仿真步骤中大部分变量是由上述符号和上述变量复合而成，其含义也由各部分的含义复合而成，如

是由变量f_v、符号_n+1、符号^trial复合而成，因此其含义为：体积屈服函数，该变量处于当前增量步，该变量为弹性试探值，其余变量依此类推。

2.根据权利要求1所述的基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法，其特征在于：

通过调整“基于压硬性非线性变化和剪缩突变特性的循环本构模型参数的获得步骤”，能够对金属材料的振动累积变形进行仿真：

A、对金属材料开展单轴拉伸试验，记录应力、应变的数据，并获得泊松比v，

B、无步骤B，步骤A紧接步骤C，

C、对金属材料开展循环加载试验，记录应力、应变的数据，设定体变起始点的孔隙比e_ini＝0，

D、设定最大孔隙比e_max＝0，

E、设定最小孔隙比e_min＝0，

F、剪切屈服条件参数C_A、C_B、C_C的获得步骤，设定C_A＝0；C_B＝0；C_C＝1，

G、剪切屈服条件参数C_D、C_E、C_F的获得步骤，设定C_D＝0；C_E＝0；C_F＝1，

H、材料在基准条件下单调压缩时的剪切硬化曲线的初值

取值为金属材料的初始剪切屈服强度，

I、材料在基准条件下的剪切硬化曲线的初始斜率

的取值为金属材料在循环加载试验获得的广义剪应力—轴应变偏量关系曲线在

点处的斜率，即q-e_a关系曲线在

点处的斜率，

J、材料在基准条件下单调压缩时的剪切硬化曲线的上限

取值为金属的剪切强度极限，

K、等效的等向压缩线梯度的分段点处的广义剪应力q_seg，q_seg取值为0至剪切强度极限，

L、等效体变模型的参数λ_eq1和λ_eq2，λ_eq1和λ_eq2取值小于1×10^-15且大于0，

M、等效体变模型的参数k_eq，k_eq取值小于1×10^-15且大于0，

N、剪切硬化权重系数W_sh，W_sh∈[0,1]，待上述其他参数确定后，根据金属循环加载试验获得的q-ε_a关系曲线的弹性区域的扩张速度和塑性区域的收敛速度，对比模型计算结果，用试错法确定W_sh，

O、体积硬化权重系数W_vh，W_vh∈[0,1]，

“基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应力驱动的仿真步骤”不做调整。

Images