CN111783332B - 压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法 - Google Patents

压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法 Download PDF

Info

Publication number
CN111783332B
CN111783332B CN202010538076.5A CN202010538076A CN111783332B CN 111783332 B CN111783332 B CN 111783332B CN 202010538076 A CN202010538076 A CN 202010538076A CN 111783332 B CN111783332 B CN 111783332B
Authority
CN
China
Prior art keywords
shear
stress
simulation
finite element
tensor
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
CN202010538076.5A
Other languages
English (en)
Other versions
CN111783332A (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.)
Nanning Urban And Rural Planning And Design Institute
Guangxi University
Original Assignee
Nanning Urban And Rural Planning And Design Institute
Guangxi University
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 Nanning Urban And Rural Planning And Design Institute, Guangxi University filed Critical Nanning Urban And Rural Planning And Design Institute
Priority to CN202010538076.5A priority Critical patent/CN111783332B/zh
Publication of CN111783332A publication Critical patent/CN111783332A/zh
Application granted granted Critical
Publication of CN111783332B publication Critical patent/CN111783332B/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
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • 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

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 hierarchicalapproach 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 Equilibriumof an Assembly of ParticlesinContact[A].Proceedings of the Royal Society A:Mathematical,PhysicalandEngineering Sciences[C].London,JSTOR,1962.500-527.}{Roscoe K.H.,ThurairajahA.,SchofieldA.N.. Yielding of Clays inStates Wetter thanCritical[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.}。
通常要对具体结构各部分的累积变形进行仿真，需要借助有限元法、有限差分法、样条函数法、边界元法等数值实现方法。其中有限元法是应用最广泛的方法。{王勋成.有限单元法[M].北京:清华大学出版社,2006.}。有限元法的积分点的应力的计算通常是应变驱动的，即已知应变状态，求解应力状态。有限元法计算单元刚度时，需要应用到具体材料的一致弹塑性模量。
发明内容
本发明要解决的技术问题是提供压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法。该方法立足于材料的循环本构模型理论及数值实现方法，即立足于广义塑性力学的分量理论、非线性的屈服条件、循环累积变形的塑性硬化模型、描述体积变形的模型和本构模型数值实现方法。该方法能够克服现有技术方案的不足，即：①能够全面反映材料的刚度和强度随周围压力和相对密实度非线性变化的行为；②能够反映材料的剪缩趋势随剪应力的增长而发生突变的特性；③仿真过程由应变驱动；④仿真过程能够向有限元程序提供一致弹塑性模量。因此仿真步骤能够与有限元法结合，对含有压硬性非线性变化和剪缩突变特性材料的具体结构各部位在循环荷载下产生的累积变形进行仿真。比如对上部结构或上部设备的长期循环荷载作用下的不同深度的地基岩土的累积变形进行仿真。为进一步的加固措施提供依据。
所述“应变驱动”，是指仿真过程中已知应变状态，求解应力状态。“应变驱动”针对的是有限元模型的积分点的应变状态受到控制的情况。
本发明通过以下技术方案解决上述技术问题，压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法，包括如下步骤：
基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应变驱动的仿真步骤。以下简称“仿真步骤”。
在执行仿真步骤前，需要建立有限元模型。
仿真步骤具体为，当循环执行增量步时，在每个增量步内依顺序执行剪切弹塑性仿真步骤、体积弹塑性仿真步骤和向有限元法传递变量的仿真步骤：
A、剪切弹塑性仿真步骤
A.a、输入常数：CA,CB,CC,CD,CE,CF,Wsh,
Figure GDA0003778943730000041
ν,eini,κeq。
输入变量：
Figure GDA0003778943730000042
Δγs,
Figure GDA0003778943730000043
Dr。
A.b、为判断剪切屈服做准备：
Figure GDA0003778943730000051
Figure GDA0003778943730000052
Figure GDA0003778943730000053
Figure GDA0003778943730000054
Figure GDA0003778943730000055
Figure GDA0003778943730000056
εn+1＝εn+Δεn+1 (7)
εv.n+1＝tr[εn+1] (8)
en+1＝εn+1-(εv.n+1/3)1 (9)
σ′n＝σn-un1 (10)
pabs.n＝tr[σ′n]/3 (11)
Figure GDA0003778943730000057
Figure GDA0003778943730000058
Figure GDA0003778943730000059
un+1＝un+Δun+1 (15)
σ′r＝σr-un+1 (16)
Figure GDA00037789437300000510
Figure GDA00037789437300000511
Figure GDA00037789437300000512
Figure GDA00037789437300000513
Figure GDA00037789437300000514
Figure GDA00037789437300000515
若||αs.n||＝0，则nαs＝ns (23)
否则nαs＝αs.n/||αs.n|| (24)
Figure GDA00037789437300000516
Figure GDA00037789437300000517
Figure GDA00037789437300000518
Figure GDA00037789437300000519
Figure GDA00037789437300000520
Figure GDA00037789437300000521
Figure GDA00037789437300000522
若
Figure GDA00037789437300000523
执行步骤3；否则，执行步骤4
A.c、发生剪切屈服时的步骤
A.c.a、确定Δγs，
A.c.a.a、初始化
Figure GDA0003778943730000061
k＝0
A.c.a.b、迭代，执行如下牛顿迭代直到
Figure GDA0003778943730000062
小于预设容许值，k←k+1
计算迭代
Figure GDA0003778943730000063
Figure GDA0003778943730000064
Figure GDA0003778943730000065
Figure GDA0003778943730000066
Figure GDA0003778943730000067
Figure GDA0003778943730000068
Figure GDA0003778943730000069
Figure GDA00037789437300000610
A.c.b、更新变量：若Δγs<0，则取Δγs＝0
Figure GDA00037789437300000611
Figure GDA00037789437300000612
αs.n+1＝ζM(αs.n+2CLΔγsns/3) (41)
Figure GDA00037789437300000613
Ks.n+1＝ζQ(Ks.n+2CPΔγs/3) (43)
Figure GDA00037789437300000614
A.c.c、计算一致弹塑性切线模量的偏张量
Figure GDA00037789437300000615
Figure GDA00037789437300000616
Figure GDA00037789437300000617
Figure GDA0003778943730000071
Figure GDA0003778943730000072
执行步骤A.e
A.d、不发生剪切屈服时的步骤：
Figure GDA0003778943730000073
执行步骤A.e
A.e、输出变量：sn+1、Δγs、
Figure GDA0003778943730000074
执行体积弹塑性仿真步骤。
B、体积弹塑性仿真步骤
B.a、输入常数：λeq1,λeq2,κeq,eini,Wvh,pabs.ini,qseg,emax,emin。
输入变量：sn+1,
Figure GDA0003778943730000075
Δγv,αv.n,Kv.n,
Figure GDA0003778943730000076
Dr,Wr。
B.b、为判断体积屈服做准备：
Figure GDA0003778943730000077
Figure GDA0003778943730000078
Figure GDA0003778943730000079
Δεv.n+1＝tr[Δεn+1] （50）
Figure GDA00037789437300000711
εn+1＝εn+Δεn+1 (7)
εv.n+1＝tr[εn+1] (8)
Figure GDA00037789437300000712
Figure GDA00037789437300000713
Figure GDA00037789437300000714
Figure GDA00037789437300000715
Figure GDA00037789437300000716
若
Figure GDA00037789437300000717
执行步骤3；否则，执行步骤4
B.c、发生体积屈服时的步骤：
B.c.a、确定Δγv，
B.c.a.a、初始化
Figure GDA00037789437300000718
k＝0
B.c.a.b、迭代，执行如下不动点迭代直到
Figure GDA00037789437300000719
小于预设容许值，k←k+1计算迭代
Figure GDA00037789437300000819
Figure GDA0003778943730000082
Figure GDA0003778943730000083
Figure GDA0003778943730000084
B.c.b、更新变量：若Δγv<0，则取Δγv＝0
Figure GDA0003778943730000085
p′n+1＝pabs.ini(T1-1) (58)
Δαv.n+1＝(1-Wvh)(1+eini)(pabs.ini+p′n+1)Δγvnv/T2 (60)
ΔKv.n+1＝Wvh(1+eini)(pabs.ini+p′n+1)Δγv/T2 (61)
αv.n+1＝αv.n+Δαv.n+1 (62)
Kv.n+1＝Kv.n+ΔKv.n+1 (63)
Figure GDA00037789437300000820
执行步骤B.e
B.d、不发生体积屈服时的步骤：
Figure GDA0003778943730000086
Figure GDA0003778943730000087
Figure GDA0003778943730000088
Dr＝(emax-e)/(emax-emin) (68)
执行步骤B.e
B.e、计算一致弹塑性切线模量的球分量
Figure GDA0003778943730000089
Figure GDA00037789437300000810
Figure GDA00037789437300000811
B.f、输出变量：p′n+1,Δγv,αv.n+1,Kv.n+1,
Figure GDA00037789437300000812
Dr。执行向有限元法传递变量的仿真步骤。
C、向有限元法传递变量的仿真步骤
C.a、输入常数：pabs.ini
输入变量：
Figure GDA00037789437300000813
p′n+1,sn+1,
Figure GDA00037789437300000814
Wr。
Figure GDA00037789437300000815
Figure GDA00037789437300000816
Figure GDA00037789437300000817
Figure GDA00037789437300000818
un+1＝un+Δun+1 (15)
p′abs.n+1＝pabs.ini+p′n+1 (72)
σ′n+1＝sn+1+p′abs.n+11 (73)
σn+1＝σ′n+1+un+11 (74)
εn+1＝εn+Δεn+1 (7)
Figure GDA0003778943730000091
Figure GDA0003778943730000092
Figure GDA0003778943730000093
Figure GDA0003778943730000094
Figure GDA0003778943730000095
Figure GDA0003778943730000096
C.c、输出变量：
Figure GDA0003778943730000097
结束当前增量步。
为了说明仿真步骤中符号和变量的含义，需预先说明“基准条件”的含义。
一些材料，如岩土材料，其剪切硬化曲线具有随周围压力增大而整体提高的特性，且在两个不同周围压力下的后继剪切屈服应力总是保持恒定的比例关系。此外，剪切硬化曲线具有随初始相对密实度增大而整体变化的特性，且两个不同初始相对密实度试样的后继剪切屈服应力总是保持恒定的比例关系。也就是说，材料存在压硬性。根据上述现象，能够以某条剪切硬化曲线上的点为基准，按照几何相似的原则，以某种比例绘制其他周围压力条件下的剪切硬化曲线。且这个比例系数与周围压力有关。因此作为基准的剪切硬化曲线所处的周围压力称为“基准周围压力”。根据上述现象，能够以某条剪切硬化曲线上的点为基准，按照几何相似的原则，以某种比例绘制其他初始相对密实度时的剪切硬化曲线。且这个比例系数与初始相对密实度有关。因此作为基准的剪切硬化曲线的试样的初始相对密实度称为“基准相对密实度”。“基准周围压力”和“基准相对密实度”统称为“基准条件”。
仿真步骤中符号的含义：变量右下标n指上一增量步；变量右下标n+1指当前增量步；变量右上标trial指该变量为弹性试探值；变量右上标*指该变量处于基准条件；变量前的Δ指该变量为增量；符号:为内积符号，是对张量缩并；变量右上标(k)指第(k)次牛顿迭代；变量右上标'指该变量为有效应力；||·||指二范数；tr[·]指对张量求迹；sign(·)为符号函数；符号
Figure GDA0003778943730000098
为张量积；exp(·)为以自然常数e为底的指数函数。
仿真步骤中的变量，粗体符号为张量，非粗体符号为标量。变量的含义：αs为试样在实际条件下的背应力偏张量；
Figure GDA0003778943730000099
为试样在基准条件下的背应力偏张量；αv为体积背应力；Bs为与周围压力有关的比例系数；CA、CB、CC、CD、CE、CF为剪切屈服条件参数；CL、CM为A-F模型的随动硬化参数；CP、CQ为Chaboche模型的等向硬化参数；Dr为相对密实度；Ds为与相对密实度有关的比例系数；e为应变偏张量；ee为弹性应变偏张量；ep为塑性应变偏张量；e为孔隙比；eini为体变起始点的孔隙比；emax为最大孔隙比；emin为最小孔隙比；ε为应变张量；εp为塑性应变张量；εv为体积应变；
Figure GDA0003778943730000101
为弹性体积应变；
Figure GDA0003778943730000102
为塑性体积应变；
Figure GDA0003778943730000103
为从有限元程序传递到仿真步骤的应变张量；
Figure GDA0003778943730000104
为从有限元程序传递到仿真步骤的应变张量的增量；
Figure GDA0003778943730000105
为从仿真步骤传递到有限元程序的应变张量； fs为剪切屈服函数；fv为体积屈服函数；G为剪切弹性模量；γs为剪切塑性滑移率；γv为体积塑性滑移率；I为四阶单位张量，I的矩阵形式表示为6×6的对角矩阵，对角线上的元素均为1；k为迭代次数指示变量；Ks为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力的各向同性硬化部分；
Figure GDA0003778943730000106
为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力的各向同性硬化部分；Kv为体积等向塑性流动应力；κeq为等效的等向膨胀线的梯度；ξs为相对应力的偏张量；ξv为相对球应力；λeq1为q≤qseg时的等效的等向压缩线梯度；λeq2为 q>qseg时的等效的等向压缩线梯度；Me为弹性刚度张量；Mep为一致弹塑性模量张量，该模量被传递给有限元程序；Mtan为一致弹塑性切线模量张量；
Figure GDA0003778943730000107
为一致弹塑性切线模量的偏张量；
Figure GDA0003778943730000108
为一致弹塑性切线模量的球分量；
Figure GDA0003778943730000109
为体积弹性模量；nv为体积塑性流动方向；ns为剪切塑性流动方向；nαs为αs方向的单位向量；ν为泊松比；o(κeq)为远小于κeq的非零正数，o(κeq)∈(0,κeq×10-4]； pabs为绝对的有效平均应力；pabs.ini为体变起点的绝对的有效平均应力；p是相对于pabs.ini而增加或减少的静水压力；Pdev为求偏张量的投影算子；q为等效剪应力，岩土工程常将其称为广义剪应力；qseg为等效的等向压缩线梯度的分段点处的广义剪应力；
Figure GDA00037789437300001010
为材料在基准条件下单调压缩时的剪切硬化曲线的初值；
Figure GDA00037789437300001011
为材料在基准条件下的剪切硬化曲线的初始斜率；
Figure GDA00037789437300001012
为材料在基准条件下单调压缩时的剪切硬化曲线的上限；s为应力偏张量；σ为应力张量；σr为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力；
Figure GDA00037789437300001013
为从有限元程序传递到仿真步骤的周围压力；
Figure GDA00037789437300001014
为从有限元程序传递到仿真步骤的应力张量；
Figure GDA00037789437300001015
为从仿真步骤传递到有限元程序的应力张量；T1为由式(57)定义的函数；T2为由式(51)定义的函数；T7为由式(69)定义的函数；u为孔隙水压力；
Figure GDA00037789437300001016
为从有限元程序传递到仿真步骤的孔隙水压力；
Figure GDA00037789437300001017
为从有限元程序传递到仿真步骤的孔隙水压力的增量；Wsh为剪切硬化权重系数，Wsh∈[0,1]；Wvh为体积硬化权重系数，Wvh∈[0,1]；Wr为权重系数，Wr∈[0,1]，使得Mep的大小介于Mtan和Me之间；ζM为式(32)定义的函数；ζQ为式(29)定义的函数；1为二阶单位张量，1的矩阵形式表示为[1 1 1 0 0 0]T。
符号和变量的补充说明，仿真步骤中大部分变量是由上述符号和上述变量复合而成，其含义也由各部分的含义复合而成。如
Figure GDA00037789437300001018
是由变量
Figure GDA00037789437300001019
符号n+1、符号trial复合而成，因此其含义为：一致弹塑性切线模量的偏张量，该变量处于当前增量步，该变量为弹性试探值。其余变量依此类推。
本发明的原理：
一、广义塑性力学的分量理论
由于应力张量和应变张量能够分解为线性无关的球张量和偏张量，本发明应用广义塑性势理论，对塑性应变进行分解
Figure GDA0003778943730000111
其中：εp为塑性应变张量；ep为塑性应变偏张量；
Figure GDA0003778943730000112
为塑性体积应变；1为二阶单位张量；γs为剪切塑性滑移率；γv为体积塑性滑移率；Qs为剪切塑性势；Qv为体积塑性势；s为应力偏张量；p为平均应力。本发明基于这个分解，分别在剪切方向和体积方向建立屈服面、硬化律、塑性流动向量。
二、基于压硬性非线性变化的材料循环本构模型的剪切分量
1.线弹性本构关系
本发明采用广义胡克定律描述材料的剪切弹性。广义胡克定律的应力张偏量—弹性应变偏张量关系，即s-ee关系表述为：
ee＝0.5s/G (82)
其中：ee为弹性应变偏张量；s为应力偏张量；G为剪切弹性模量，根据弹性理论，表示为
Figure GDA0003778943730000113
其中：
Figure GDA0003778943730000114
为体积弹性模量；ν为泊松比。
2.非线性剪切屈服条件
包含背应力项和等向塑性流动应力项的非线性剪切屈服条件的表达式为：
Figure GDA0003778943730000115
Figure GDA0003778943730000116
Figure GDA0003778943730000117
Figure GDA0003778943730000118
其中：fs为剪切屈服函数；s为应力偏张量；αs为背应力偏张量；Ks为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力q的各向同性硬化部分；
Figure GDA0003778943730000119
为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力q的各向同性硬化部分；Hs为剪切硬化内变量，
Figure GDA00037789437300001110
Figure GDA00037789437300001111
为塑性等效剪应变，岩土工程常将其称为塑性广义剪应变；σr为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力，其值等于试样在整体上所受到的中主应力σ2；Dr为相对密实度，本发明将其取值为从塑性屈服后到弹性卸载前的一次连续塑性过程的初始相对密实度；Bs为与周围压力有关的比例系数；Ds为与相对密实度有关的比例系数；CA、CB、 CC为剪切屈服条件参数，为常数，通过至少3个不同恒定周围压力的三轴压缩试验进行回归确定。CD、 CE、CF为剪切屈服条件参数，为常数，通过至少3个不同相对密实度的试样的三轴压缩试验回归确定。根据具体材料，比例系数Bs和Ds表示为直线函数、双曲线函数、指数函数、幂函数或对数函数。当CA为0时，Bs退化为直线函数，能够描述剪切屈服应力随周围压力线性增长。当CD为0时，Ds退化为直线函数，能够描述剪切屈服应力随初始相对密实度线性增长。
3.关联相对应力偏量的剪切塑性流动法则
建立循环加载本构模型时采用与相对应力偏量(s-αs)同向的单位向量作为塑性流动方向。
4.基于耦合硬化模型的剪切硬化律
采用组合随动/等向硬化律描述材料的剪切塑性硬化，并按权重分配每个无穷小增量q的随动/等向硬化的占比，即：
Figure GDA0003778943730000121
Figure GDA0003778943730000122
Figure GDA0003778943730000123
其中：q为剪切硬化曲线上的等效剪应力，岩土工程常将其称为广义剪应力；Wsh为剪切硬化权重系数，Wsh∈[0,1]，根据振动三轴试验获得的广义剪应力—轴应变偏量关系，即q-ea关系曲线确定。
采用A-F随动硬化模型和Chaboche等向硬化模型进一步描述材料的剪切硬化曲线：
Figure GDA0003778943730000124
Figure GDA0003778943730000125
其中：γs为剪切塑性滑移率，
Figure GDA0003778943730000126
CL和CM为A-F随动硬化参数；ns为塑性流动方向；CP和CQ为Chaboche等向硬化参数。
一些材料，如岩土材料，其剪切硬化曲线的形状受到周围压力和相对密实度影响，因此A-F 随动硬化模型和Chaboche等向硬化模型的参数也会随着这些条件的变化而变化。结合材料的压硬性的非线性变化性质拓展这两个模型，得到等向硬化参数和随动硬化参数为
Figure GDA0003778943730000127
其中：
Figure GDA0003778943730000128
为材料在基准条件下单调压缩时的剪切硬化曲线的初值，对于岩土材料其取值小于剪切强度极限的1/100；
Figure GDA0003778943730000129
为材料在基准条件下单调压缩时的剪切硬化曲线的上限，通过相应的三轴压缩试验获得；
Figure GDA00037789437300001210
为材料在基准条件下的剪切硬化曲线的初始斜率，通过相应的振动三轴试验第1个滞回圈的初始上升段获得。
三、基于剪缩突变特性的圆砾循环本构模型的体积分量
1.等效体变模型
将材料的体变εv分解为弹性体变
Figure GDA00037789437300001211
和塑性体变
Figure GDA00037789437300001212
并采用Roscoe等提出的公式描述
Figure GDA00037789437300001213
和
Figure GDA0003778943730000131
即
Figure GDA0003778943730000132
Figure GDA0003778943730000133
其中：
Figure GDA0003778943730000134
为弹性体变；
Figure GDA0003778943730000135
为塑性体变；pabs.ini为体变起点的绝对的有效平均应力，在振动三轴试验为剪切阶段的初始有效平均应力；eini为体变起始点的孔隙比，在振动三轴试验为剪切阶段的初始孔隙比；p是相对于pabs.ini而增加或减少的静水压力的量；κeq为等效的等向膨胀线的梯度；λeq为等效的等向压缩线的梯度。
2.体积屈服条件和体积塑性流动法则
由于p和
Figure GDA0003778943730000136
为标量，本发明采用一维的屈服条件描述体积屈服，即
fv＝|p-αv(Hv)-Kv(Hv) (93)
其中：fv为体积屈服函数；αv为体积背应力；Kv为体积等向塑性流动应力；Hv为体积硬化内变量，
Figure GDA0003778943730000137
采用关联相对平均应力(p-αv)方向的塑性流动法则描述体积塑性流动，即
nv＝sign(p-αv) (94)
其中：nv为体积塑性流动方向；sign(·)为符号函数。
3.分段梯度的体积硬化律
采用组合随动/等向硬化律描述体积塑性硬化。其中体积的组合随动/等向硬化律按权重分配每个无穷小增量p的随动/等向硬化的占比。
Figure GDA0003778943730000138
Figure GDA0003778943730000139
Figure GDA00037789437300001310
其中：Wvh为体积硬化权重系数，Wvh∈[0,1]。Wvh通过振动三轴试验获得的p-εv关系曲线确定。针对材料的剪缩趋势不连续变化的现象，本发明提出了分段梯度的体积硬化律。联立式(91)、式(92)、式(95)、式(96)、式(97)得到体积随动硬化模型和体积等向硬化模型
Figure GDA00037789437300001311
Figure GDA00037789437300001312
Figure GDA00037789437300001313
其中：o(κeq)为远小于κeq的非零正数，o(κeq)∈(0,κeq×10-4]；qseg为等效的等向压缩线梯度的分段点处的广义剪应力；λeq1为q≤qseg时的等效的等向压缩线梯度；λeq2为q>qseg时的等效的等向压缩线梯度。由于剪缩突变，在较高的q水平时和较低的q水平时圆砾的p-ev曲线斜率存在明显差异，本发明将λeq分成2段，通过式(100)的前两式表示这种差异。为了描述膨胀时的包辛格效应，其中式(100)的第三式控制着膨胀且屈服时几乎不发生塑性变形。
四、基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应变驱动的仿真步骤中部分公式的说明
1.剪切背应力的向后欧拉差分形式，即式(41)和式(32)
证明：对A-F随动硬化模型式(88)的等号两边乘以时间增量Δt，有
Figure GDA0003778943730000141
式(101)联立
Figure GDA0003778943730000142
得到
Figure GDA0003778943730000143
式(103)联立向后欧拉差分
αs.n+1＝αs.n+Δαs.n+1 (104)
得到
Figure GDA0003778943730000144
解方程式(105)得到
Figure GDA0003778943730000145
即：αs.n+1＝ζM(αs.n+2CLΔγsns/3) (41)
其中：
Figure GDA0003778943730000146
证毕。
2.剪切等向塑性流动应力的向后欧拉差分形式，即式(43)和式(29)
证明：将Chaboche等向硬化模型式(89)等号两边乘以时间增量Δt，有
Figure GDA0003778943730000147
式(107)联立式(102)得到
Figure GDA0003778943730000148
式(108)联立向后欧拉差分
Ks.n+1＝Ks.n+ΔKs.n+1 (109)
得到
Figure GDA0003778943730000149
解方程式(110)得到
Figure GDA00037789437300001410
即：Ks.n+1＝ζQ(Ks.n+2CPΔγs/3) (43)
其中：
Figure GDA00037789437300001411
证毕。
3.体积背应力增量在当前增量步的形式，即式(60)
证明：将体积塑性滑移率
Figure GDA00037789437300001412
代入式(98)，得到
Figure GDA00037789437300001413
对式(113)的等号两边乘以时间增量Δt，有
Figure GDA0003778943730000151
由式(114)得到
Figure GDA0003778943730000152
在当前增量步有Δαv.n+1＝(1-Wvh)(1+eini)(pabs.ini+pn+1)Δγvnv/T2 (60)
证毕。
4.体积等向塑性流动应力增量在当前增量步的形式，即式(61)
证明：将体积塑性滑移率
Figure GDA0003778943730000153
代入式(99)，得到
Figure GDA0003778943730000154
对式(116)的等号两边乘以时间增量Δt，有
Figure GDA0003778943730000155
由式(117)得到
Figure GDA0003778943730000156
在当前增量步有ΔKv.n+1＝Wvh(1+eini)(pabs.ini+pn+1)Δγv/T2 (61)
证毕。
5.用于求解塑性滑移的剪切屈服条件的差分形式，即式(36)、式(32)、式(33)和式(21)。证明：应力偏张量的“最近点投影”的形式可通过如下方式推导
Figure GDA0003778943730000157
其中：变量右上标trial是指弹性试探值。联立式(39)以及剪切相对应力的定义
ξs.n+1＝sn+1-αs.n+1 (119)
以及αs.n+1＝ζM(αs.n+2CLΔγsns/3) (41)
以及
Figure GDA0003778943730000158
得到
Figure GDA0003778943730000159
整理得
Figure GDA00037789437300001510
式(121)等号两边内积径向流动向量ns得到
Figure GDA00037789437300001511
这里假设了ξs.n+1的方向和
Figure GDA00037789437300001512
的方向同向。圆砾在三轴压缩、三轴拉伸和三轴卸载时满足π平面上的
Figure GDA00037789437300001513
sn+1、αs.n+1和αs.n在一条直线上，即这些变量都在Lode角θ＝-π/6或θ＝π/6的位置上，因此ξs.n+1的方向和
Figure GDA0003778943730000161
的方向同向，式(122)在本发明范围内成立。
联立式(119)、式(122)以及剪切屈服条件式(83)在当前增量步的形式
Figure GDA0003778943730000162
得到
Figure GDA0003778943730000163
其中：||·||指二范数。将式(43)和式(29)代入式(124)得到用于求解塑性滑移的剪切屈服条件的差分形式为
Figure GDA0003778943730000164
即
Figure GDA0003778943730000165
其中：
Figure GDA0003778943730000166
Figure GDA0003778943730000167
Figure GDA0003778943730000168
证毕。
6.剪切塑性滑移对应变的偏导的差分形式，即式(46)。
证明：对式(36)求相对于εn+1的偏导数
Figure GDA0003778943730000169
根据弹性理论，有
Figure GDA00037789437300001610
Figure GDA00037789437300001611
Figure GDA00037789437300001612
以及引理：
Figure GDA00037789437300001613
联立式(126)、式(127)、式(128)、式(129)、式(130)得到
Figure GDA0003778943730000171
其中：Pdev为求偏张量的投影算子。I为四阶单位张量，I的矩阵形式表示为6×6的对角矩阵，对角线上的元素均为1。符号
Figure GDA0003778943730000172
为张量积。根据一致性条件
当fs＝0时，
Figure GDA0003778943730000173
有
Figure GDA0003778943730000174
解方程得到剪切塑性滑移对应变的偏导的差分形式为
Figure GDA0003778943730000175
证毕。
7.一致弹塑性切线模量的偏张量在当前增量步的形式，即式(47)。
证明：结合式(128)、式(129)、式(130)，对式(39)等号两边求微分
Figure GDA0003778943730000176
因此一致弹塑性切线模量的偏张量在当前增量步的形式为
Figure GDA0003778943730000177
证毕。
8.平均应力的弹性试探值，即式(52)。
证明：对Roscoe等提出的弹性体变模型式(91)的等号两边同时乘时间微分dt得到
Figure GDA0003778943730000178
对式(134)的等号两边同时积分得到：
Figure GDA0003778943730000179
其中：C为待定常数。应用初始条件
Figure GDA00037789437300001710
得到C＝1/pabs.ini，则：
Figure GDA0003778943730000181
其中：exp(·)为以自然常数e为底的指数函数。
由于本文将圆砾的体变εv分解为弹性体变
Figure GDA0003778943730000182
和塑性体变
Figure GDA0003778943730000183
有
Figure GDA0003778943730000184
联立式(136)和式(137)可得到有效平均应力的弹性试探值为
Figure GDA0003778943730000185
证毕。
9.有效平均应力的向后欧拉差分形式，即式(57)和式(58)。
证明：联立式(136)和式(137)得到
Figure GDA0003778943730000186
由弹塑性理论有
Figure GDA0003778943730000187
则
Figure GDA0003778943730000188
的向后欧拉差分形式为
Figure GDA0003778943730000189
联立式(138)和式(64)得到有效平均应力的向后欧拉差分形式为
Figure GDA00037789437300001810
p′n+1＝pabs.ini(T1-1) (58)
证毕。
10.体积塑性滑移增量，即式(59)。
证明：将式(52)写成如下形式：
Figure GDA00037789437300001811
其中：
Figure GDA00037789437300001812
将式(57)和式(58)写成如下形式：
Figure GDA00037789437300001813
其中：ΔX＝-Δγvnv。对式(141)进行一阶泰勒展开得到：
Figure GDA00037789437300001814
将式(142)写成向后欧拉差分形式得到：
Figure GDA00037789437300001815
相对平均应力的定义为：ξv.n+1＝p′n+1-αv.n+1 (144)
联立式(144)和式(62)得到：
ξv.n+1＝p′n+1-αv.n-Δαv.n+1 (145)
联立式(145)和式(60)得到：
ξv.n+1＝p′n+1-αv.n-(1-Wvh)(1+eini)(pabs.ini+p′n+1)Δγvnv/T2 (146)
相对平均应力的弹性试探值的定义为：
Figure GDA0003778943730000191
联立式(146)、式(53)和式(144)得到：
Figure GDA0003778943730000192
式(147)成立的前提是ξv.n+1和
Figure GDA0003778943730000193
的方向同向。由于体积变形只有膨胀和压缩两种状态，体积弹塑变形是个一维问题。
Figure GDA0003778943730000194
p′n+1、αv.n+1、αv.n自然在一条直线上，因此ξv.n+1和
Figure GDA0003778943730000195
的方向均为nv，式(147)成立。由于ξv.n+1＝|ξv.n+1|nv
Figure GDA0003778943730000196
Figure GDA0003778943730000197
根据一致性条件，在材料屈服时，体积屈服函数fv＝0。
因此联立式(149)和式(61)得到：
Figure GDA0003778943730000201
解方程式(150)得到体积塑性滑移增量为：
Figure GDA0003778943730000202
证毕。
11.一致弹塑性切线模量的球分量在当前增量步的形式，即式(69)和式(70)。证明：(1)膨胀时的一致弹塑性切线模量的球分量
整理式(134)得到
Figure GDA0003778943730000203
在当前增量步有
Figure GDA0003778943730000204
(2)压缩时的一致弹塑性切线模量的球分量
对Roscoe等提出的塑性体变模型式(92)的等号两边同时乘时间微分dt得到
Figure GDA0003778943730000205
联立式(153)和式(134)得到
Figure GDA0003778943730000206
整理式(153)得到
Figure GDA0003778943730000207
在当前增量步有
Figure GDA0003778943730000208
(3)一致弹塑性切线模量的球分量在当前增量步的形式
对压缩时的
Figure GDA0003778943730000209
按剪应力大小分段，并且联立膨胀时的
Figure GDA00037789437300002010
有
Figure GDA00037789437300002011
Figure GDA00037789437300002012
证毕。
12.体积弹性模量在当前增量步的形式，即式(71)。
证明：整理式(134)得到
Figure GDA0003778943730000211
在当前增量步，体积弹性模量在当前增量步的形式为
Figure GDA0003778943730000212
证毕。
13.不动点存在性证明。“B、体积弹塑性仿真步骤”的中的式(57)、式(58)和式(59)的迭代收敛的前提是不动点存在。以下将证明该迭代的不动点存在。
证明：设
Figure GDA0003778943730000213
是映射F的定义域D的闭子集，
Figure GDA0003778943730000214
令y＝Δγv，将式(57)、式(58)和式(59)写成y＝F(Δγv)的形式：
Figure GDA0003778943730000215
其中：
Figure GDA0003778943730000216
则：
Figure GDA0003778943730000217
其中：
Figure GDA0003778943730000218
∵eini≥0，κeq>0，
Figure GDA0003778943730000219
∴0≤|exp(T3)|≤1
∵κeq≥0；
根据式(51)，T2≥0
∴0≤|T2/(κeq+T2)|≤1
∵在正向加载时nv＝1；在反向加载时nv＝-1
∴|nv|＝1
∵通常在工程实践中，总量远大于增量，即：
Figure GDA00037789437300002110
∴
Figure GDA00037789437300002111
∴
Figure GDA00037789437300002112
根据柯西中值定理
Figure GDA00037789437300002113
其中：
Figure GDA00037789437300002114
根据压缩映射原理，从式(162)和式(163)可知：
采用式(57)、式(58)和式(59)进行迭代时，存在不动点。或者表述为，式(57)、式(58)和式(59)为不动点迭代。证毕。
14.仿真步骤的其他补充说明
式(12)其中绝对平均应力pabs采用的是上一增量步pabs.n而不是当前增量步pabs.n+1，这是因为pabs.n+1并不能事先知道，只要Δp′n+1趋于无穷小，
Figure GDA0003778943730000221
的误差也趋于无穷小，因此
Figure GDA0003778943730000222
的误差为Δp′n+1的高阶无穷小o(Δp′n+1)。式(31)为式(83)的剪切屈服函数的差分形式，由于其中的Δγs是从上一增量步传来的，该函数为试剪切屈服函数。式(51)是在式(100)基础上调整了分段条件：因为本仿真步骤是应变驱动的，平均应力增量Δp并不能事先知道，所以式(51)采用Δεv.n+1判断拉压趋势。式(56)为式(93)的体积屈服函数的差分形式，由于其中的Δγv是从上一增量步传来的，该函数为试体积屈服函数。式(76)为式(13)在当前增量步的形式。
由于通常有限元法规定的力和变形的符号与部分工程学科规定的符号相反，在有限元程序向仿真步骤传递应力、应变、孔隙水压力和周围压力的变量时，将这些变量取相反数，即式(1)、式(2)、式(3)、式(4)、式(5)、式(6)；在仿真步骤向有限元程序传递应力和应变的变量时，将这些变量取相反数，即式(79)和式(80)。
本仿真步骤在向有限元程序传递材料的弹塑性模量矩阵时，传递的并不仅限于切线模量
Figure GDA0003778943730000223
(式 (75))，而是介于切线模量
Figure GDA0003778943730000224
与弹性刚度
Figure GDA0003778943730000225
(式(77))之间的模量
Figure GDA0003778943730000226
(式(78))。通常不宜采用切线模量，而是比切线模量稍大的模量，以避免可能的迭代发散。其具体原因是，通常在计算切线模量时，由于存在数值误差，所得到的模量可能比切线模量稍小，这可能造成迭代收敛困难。根据有限元法，只要计算能够收敛，传递给有限元程序的一致弹塑性模量的大小并不影响计算结果，只影响收敛速度。通常采用切线模量计算的收敛速度最快。
本发明的有益效果是：
(1)非线性剪切屈服条件，即式(17)～式(31)能够在仿真过程中反映周围压力和相对密实度对剪切屈服面的非线性的影响；
(2)拓展的A-F随动硬化模型和Chaboche等向硬化模型，即式(41)～式(44)能够在仿真过程中反映塑性硬化模量随着周围压力和相对密实度的变化而非线性地变化；
(3)分段梯度的体积硬化律，即式(51)、式(60)和式(61)能够在仿真过程中准确反映剪缩趋势不连续变化的特性；
(4)仿真过程为应变驱动，即通过式(1)和式(2)输入应变状态，通过式(80)输出应力状态，使得仿真步骤能够与有限元程序良好对接；
(5)仿真过程提供了一致弹塑性模量
Figure GDA0003778943730000227
(式(78))，用于有限元程序计算单元刚度，这使得仿真步骤能够与有限元程序良好对接；
(6)仿真步骤是基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真步骤，而仿真步骤又能够与有限元法结合，因此本发明能够对含有压硬性非线性变化和剪缩突变特性材料的具体结构的各部位在循环荷载下产生的累积变形进行仿真。
(7)仿真步骤立足于向后欧拉差分法，具有一阶精确性和无条件(线性化)稳定性；
(8)采用本发明能够准确预测结构各部位的长期累积变形。
附图说明
图1为圆砾试样的有限元模型。
图2为有限元模型的单元与结点编号。
图3为南宁圆砾Dr＝0.3试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数 (ea-εv-εa-N)关系的试验曲线与仿真曲线对比图。
图4为南宁圆砾Dr＝0.5试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数 (ea-εv-εa-N)关系的试验曲线与仿真曲线对比图。
图5为南宁圆砾Dr＝0.7试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数 (ea-εv-εa-N)关系的试验曲线与仿真曲线对比图。
图6为南宁圆砾Dr＝0.3试样在循环荷载下结构各单元的轴应变偏量—振动次数(ea-N) 关系的仿真曲线图。
图7为南宁圆砾Dr＝0.5试样在循环荷载下结构各单元的轴应变偏量—振动次数(ea-N) 关系的仿真曲线图。
图8为南宁圆砾Dr＝0.7试样在循环荷载下结构各单元的轴应变偏量—振动次数(ea-N) 关系的仿真曲线图。
图9为南宁圆砾Dr＝0.3试样在循环荷载下结构各单元的侧应变偏量—振动次数(er-N) 关系的仿真曲线图。
图10为南宁圆砾Dr＝0.5试样在循环荷载下结构各单元的侧应变偏量—振动次数(er-N) 关系的仿真曲线图。
图11为南宁圆砾Dr＝0.7试样在循环荷载下结构各单元的侧应变偏量—振动次数(er-N) 关系的仿真曲线图。
图12为南宁圆砾Dr＝0.3试样在循环荷载下结构各单元的轴应变—振动次数(εa-N)关系的仿真曲线图。
图13为南宁圆砾Dr＝0.5试样在循环荷载下结构各单元的轴应变—振动次数(εa-N)关系的仿真曲线图。
图14为南宁圆砾Dr＝0.7试样在循环荷载下结构各单元的轴应变—振动次数(εa-N)关系的仿真曲线图。
图15为南宁圆砾Dr＝0.3试样在循环荷载下结构各单元的侧应变—振动次数(εr-N)关系的仿真曲线图。
图16为南宁圆砾Dr＝0.5试样在循环荷载下结构各单元的侧应变—振动次数(εr-N)关系的仿真曲线图。
图17为南宁圆砾Dr＝0.7试样在循环荷载下结构各单元的侧应变—振动次数(εr-N)关系的仿真曲线图。
图18为南宁圆砾Dr＝0.3试样在循环荷载下结构各单元的体积应变—振动次数(εv-N) 关系的仿真曲线图。
图19为南宁圆砾Dr＝0.5试样在循环荷载下结构各单元的体积应变—振动次数(εv-N) 关系的仿真曲线图。
图20为南宁圆砾Dr＝0.7试样在循环荷载下结构各单元的体积应变—振动次数(εv-N) 关系的仿真曲线图。
具体实施方式
以下通过实施例对本发明的技术方案作进一步说明。
实施例
本发明所述的压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法的具体应用实例，是对南宁圆砾的振动三轴试验所测累积变形进行仿真，在每个增量步内，按顺序执行如下步骤：
基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应变驱动的仿真步骤。以下简称“仿真步骤”。
在执行仿真步骤前，需要建立有限元模型。采用三节点三角形环状单元建立有限元模型。有限元模型的形状和尺寸与圆砾振动三轴试验的试样相同，即Ф300mm×H750mm的圆柱形，如图1 所示。整个有限元模型为4个单元，单元与结点编号如图2所示。圆柱形的底面有竖向约束，顶面作用有均布变化轴压力σa，中轴线有水平约束，外侧面作用有均布变化周围压力σr。采用对角元素改1法强制位移边界条件。数值积分方法为单元中点近似法。采用考虑平衡校正的隐式欧拉增量解法求解非线性代数方程，并且为了减少数值误差，当轴应力到达极值时进行收敛迭代。
仿真步骤具体为，当循环执行增量步时，在每个增量步内依顺序执行剪切弹塑性仿真步骤、体积弹塑性仿真步骤和向有限元法传递变量的仿真步骤：
以南宁圆砾Dr＝0.5的试样的第4个单元的第60个增量步为例，
A、剪切弹塑性仿真步骤
A.a、输入常数：CA,CB,CC,CD,CE,CF,Wsh,
Figure GDA0003778943730000251
ν,eini,κeq。
CA＝-2.3455；CB＝5.3433；CC＝0.0252；CD＝-0.3571；CE＝0.7143；CF＝0.7321；Wsh＝1.07×10-5；
Figure GDA0003778943730000252
ν＝0.15；eini＝0.5150；κeq＝7.6730×10-4
输入变量：
Figure GDA0003778943730000253
Δγs,
Figure GDA0003778943730000254
Dr。
A.b、为判断剪切屈服做准备：
Figure GDA0003778943730000255
εn＝[0.1221 -0.3418 -0.0102 -0.1059 0 0]×10-5
Figure GDA0003778943730000256
Δεn+1＝[-0.0950 0.2174 -0.0949 0.0001 0 0]×10-4
Figure GDA0003778943730000257
σn＝[0.1996 0.1994 0.1995-0.0001 0 0]
Figure GDA0003778943730000258
un＝0.0078
Figure GDA0003778943730000259
Δun+1＝-0.0033
Figure GDA00037789437300002510
σr＝0.1993
εn+1＝εn+Δεn+1 (7)
εn+1＝[-0.0828 0.1832 -0.0959 -0.0105 0 0]×10-4
εv.n+1＝tr[εn+1] (8)
εv.n+1＝4.4400×10-7
en+1＝εn+1-(εv.n+1/3)1 (9)
en+1＝[-0.0843 0.1817 -0.0974 -0.0105 0 0]×10-4
σ′n＝σn-un1 (10)
σ′n＝[0.1918 0.1916 0.1917 -0.0001 0 0]
pabs.n＝tr[σ′n]/3 (11)
pabs.n＝0.1917
Figure GDA00037789437300002511
Figure GDA00037789437300002512
Figure GDA00037789437300002513
G＝354.8929
Figure GDA00037789437300002514
Figure GDA00037789437300002515
un+1＝un+Δun+1 (15)
un+1＝0.0045
σ′r＝σr-un+1 (16)
σ′r＝0.1949
Figure GDA0003778943730000261
Bs＝0.9774
Figure GDA0003778943730000262
Ds＝0.9887
Figure GDA0003778943730000263
αs.n＝[0.4311 -0.7377 0.3066 -0.0996 0 0]×10-3
Figure GDA0003778943730000264
Ks.n＝9.6799×10-4
Figure GDA0003778943730000265
Figure GDA0003778943730000266
Figure GDA0003778943730000267
ns＝[-0.4084 0.8165 -0.4081 0.0002 0 0]
若||αs.n||＝0，则nαs＝ns (23)
否则nαs＝αs.n/||αs.n|| (24)
nαs＝[0.4721 -0.8078 0.3357 -0.1091 0 0]
Figure GDA0003778943730000268
CL＝213.0316
Figure GDA0003778943730000269
CM＝375
Figure GDA00037789437300002610
CP＝-1.1768×10-5
Figure GDA00037789437300002611
CQ＝-1.9359
Figure GDA00037789437300002612
ζQ＝1.0000
Figure GDA00037789437300002613
Figure GDA00037789437300002614
Figure GDA00037789437300002615
Figure GDA00037789437300002616
若
Figure GDA0003778943730000271
执行步骤3；否则，执行步骤4
A.c、发生剪切屈服时的步骤
A.c.a、确定Δγs，
A.c.a.a、初始化
Figure GDA0003778943730000272
k＝0
A.c.a.b、迭代，执行如下牛顿迭代直到
Figure GDA0003778943730000273
小于预设容许值，k←k+1
计算迭代
Figure GDA0003778943730000274
Figure GDA0003778943730000275
Figure GDA0003778943730000276
Figure GDA0003778943730000277
Figure GDA0003778943730000278
Figure GDA0003778943730000279
Figure GDA00037789437300002710
Figure GDA00037789437300002711
Δγs＝2.1265×10-5
A.c.b、更新变量：若Δγs<0，则取Δγs＝0
Figure GDA00037789437300002712
sn+1＝[-0.0011 0.0024 -0.0012 -0.0001 0 0]
Figure GDA00037789437300002713
Figure GDA00037789437300002714
αs.n+1＝ζM(αs.n+2CLΔγsns/3) (41)
αs.n+1＝[-0.0008 0.0017 -0.0009 -0.0001 0 0]
Figure GDA00037789437300002715
Figure GDA00037789437300002716
Ks.n+1＝ζQ(Ks.n+2CPΔγs/3) (43)
Ks.n+1＝9.6802×10-4
Figure GDA00037789437300002717
Figure GDA00037789437300002718
A.c.c、计算一致弹塑性切线模量的偏张量
Figure GDA0003778943730000281
Figure GDA0003778943730000282
Figure GDA0003778943730000283
Figure GDA0003778943730000284
Figure GDA0003778943730000285
Figure GDA0003778943730000286
Figure GDA0003778943730000287
Figure GDA0003778943730000288
Figure GDA0003778943730000289
Figure GDA00037789437300002810
执行步骤A.e
A.d、不发生剪切屈服时的步骤：
Figure GDA00037789437300002811
执行步骤A.e
A.e、输出变量：sn+1、Δγs、
Figure GDA00037789437300002812
执行体积弹塑性仿真步骤。
B、体积弹塑性仿真步骤
B.a、输入常数：λeq1,λeq2,κeq,eini,Wvh,pabs.ini,qseg,emax,emin。
λeq1＝9.2226×10-4；λeq2＝1.7154×10-3；κeq＝7.6730×10-4；eini＝0.5150；Wvh＝0.00088； pabs.ini＝0.2MPa；qseg＝0.0550MPa；emax＝0.684；emin＝0.411。
输入变量：sn+1,
Figure GDA0003778943730000291
Δγv,αv.n,Kv.n,
Figure GDA0003778943730000292
Dr,Wr。
B.b、为判断体积屈服做准备：
Figure GDA0003778943730000293
εn＝[0.1221 -0.3418 -0.0102 -0.1059 0 0]×10-5
Figure GDA0003778943730000294
Δεn+1＝[-0.0950 0.2174 -0.0949 0.0001 0 0]×10-4
Figure GDA0003778943730000295
qn+1＝0.0035
Δεv.n+1＝tr[Δεn+1] (50)
Δεv.n+1＝2.7432×10-6
Figure GDA0003778943730000296
T2＝1.5495×10-4
εn+1＝εn+Δεn+1 (7)
εn+1＝[-0.0828 0.1832 -0.0959 -0.0105 0 0]×10-4
εv.n+1＝tr[εn+1] (8)
εvn+1＝4.4400×10-7
Figure GDA0003778943730000297
Figure GDA0003778943730000298
Figure GDA0003778943730000299
Figure GDA00037789437300002910
Figure GDA00037789437300002911
nv＝1
Figure GDA00037789437300002912
Figure GDA00037789437300002913
Figure GDA00037789437300002914
Figure GDA00037789437300002915
若
Figure GDA00037789437300002916
执行步骤3；否则，执行步骤4
B.c、发生体积屈服时的步骤：
B.c.a、确定Δγv，
B.c.a.a、初始化
Figure GDA00037789437300002917
k＝0
B.c.a.b、迭代，执行如下不动点迭代直到
Figure GDA00037789437300002918
小于预设容许值，k←k+1计算迭代
Figure GDA00037789437300003013
Figure GDA0003778943730000302
Figure GDA0003778943730000303
Figure GDA0003778943730000304
Δγv＝4.6000×10-7
B.c.b、更新变量：若Δγv<0，则取Δγv＝0
Figure GDA0003778943730000305
T1＝0.9813
p′n+1＝pabs.ini(T1-1) (58)
p′n+1＝-0.0037
Δαv.n+1＝(1-Wvh)(1+eini)(pabs.ini+p′n+1)Δγvnv/T2 (60)
Δαv.n+1＝9.0566×10-4
ΔKv.n+1＝Wvh(1+eini)(pabs.ini+p′n+1)Δγv/T2 (61)
ΔKv.n+1＝7.9768×10-7
αv.n+1＝αv.n+Δαv.n+1 (62)
αv.n+1＝-0.0038
Kv.n+1＝Kv.n+ΔKv.n+1 (63)
Kv.n+1＝6.6963×10-5
Figure GDA0003778943730000306
Figure GDA0003778943730000307
执行步骤B.e
B.d、不发生体积屈服时的步骤：
Figure GDA0003778943730000308
Figure GDA0003778943730000309
Figure GDA00037789437300003010
Dr＝(emax-e)/(emax-emin) (68)
执行步骤B.e
B.e、计算一致弹塑性切线模量的球分量
Figure GDA00037789437300003011
T7＝9.2226×10-4
Figure GDA00037789437300003012
Figure GDA0003778943730000311
Figure GDA0003778943730000312
Figure GDA0003778943730000313
B.f、输出变量：p′n+1,Δγv,αv.n+1,Kv.n+1,
Figure GDA0003778943730000314
Dr。执行向有限元法传递变量的仿真步骤。
C、向有限元法传递变量的仿真步骤
C.a、输入常数：pabs.ini
输入变量：
Figure GDA0003778943730000315
p′n+1,sn+1,
Figure GDA0003778943730000316
Wr。
Figure GDA0003778943730000317
εn＝[0.1221 -0.3418 -0.0102 -0.1059 0 0]×10-5
Figure GDA0003778943730000318
Δεn+1＝[-0.0950 0.2174 -0.0949 0.0001 0 0]×10-4
Figure GDA0003778943730000319
un＝0.0078
Figure GDA00037789437300003110
Δun+1＝-0.0033
un+1＝un+Δun+1 (15)
un+1＝0.0045
p′abs.n+1＝pabs.ini+p′n+1 (72)
p′abs.n+1＝0.1963
σ′n+1＝sn+1+p′abs.n+11 (73)
σ′n+1＝[0.1951 0.1986 0.1950 -0.0001 0 0]
σn+1＝σ′n+1+un+11 (74)
σn+1＝[0.1996 0.2031 0.1995 -0.0001 0 0]
εn+1＝εn+Δεn+1 (7)
εn+1＝[-0.0828 0.1832 -0.0959 -0.0105 0 0]
Figure GDA00037789437300003111
Figure GDA00037789437300003112
Figure GDA00037789437300003113
G＝363.3349
Figure GDA00037789437300003114
Figure GDA0003778943730000321
Figure GDA0003778943730000322
Figure GDA0003778943730000323
Figure GDA0003778943730000324
Figure GDA0003778943730000325
Figure GDA0003778943730000326
Figure GDA0003778943730000327
C.c、输出变量：
Figure GDA0003778943730000328
结束当前增量步。
以上仿真步骤中符号的含义：变量右下标n指上一增量步；变量右下标n+1指当前增量步；变量右上标trial指该变量为弹性试探值；变量右上标*指该变量处于基准条件；变量前的Δ指该变量为增量；符号:为内积符号，是对张量缩并；变量右上标(k)指第(k)次牛顿迭代；变量右上标'指该变量为有效应力；||·||指二范数；tr[·]指对张量求迹；sign(·)为符号函数；符号
Figure GDA0003778943730000329
为张量积；exp(·)为以自然常数e为底的指数函数。
以上仿真步骤中的变量，粗体符号为张量，非粗体符号为标量。变量的含义：αs为试样在实际条件下的背应力偏张量；
Figure GDA00037789437300003210
为试样在基准条件下的背应力偏张量；αv为体积背应力；Bs为与周围压力有关的比例系数；CA、CB、CC、CD、CE、CF为剪切屈服条件参数；CL、CM为A-F模型的随动硬化参数；CP、CQ为Chaboche模型的等向硬化参数；Dr为相对密实度；Ds为与相对密实度有关的比例系数；e为应变偏张量；ee为弹性应变偏张量；ep为塑性应变偏张量；e为孔隙比；eini为体变起始点的孔隙比；emax为最大孔隙比；emin为最小孔隙比；ε为应变张量；εp为塑性应变张量；εv为体积应变；
Figure GDA00037789437300003211
为弹性体积应变；
Figure GDA00037789437300003212
为塑性体积应变；
Figure GDA00037789437300003213
为从有限元程序传递到仿真步骤的应变张量；
Figure GDA00037789437300003214
为从有限元程序传递到仿真步骤的应变张量的增量；
Figure GDA00037789437300003215
为从仿真步骤传递到有限元程序的应变张量；fs为剪切屈服函数；fv为体积屈服函数；G为剪切弹性模量；γs为剪切塑性滑移率；γv为体积塑性滑移率；I为四阶单位张量，I的矩阵形式表示为6×6的对角矩阵，对角线上的元素均为1； k为迭代次数指示变量；Ks为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力的各向同性硬化部分；
Figure GDA0003778943730000331
为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力的各向同性硬化部分；Kv为体积等向塑性流动应力；κeq为等效的等向膨胀线的梯度；ξs为相对应力的偏张量；ξv为相对球应力；λeq1为q≤qseg时的等效的等向压缩线梯度；λeq2为q>qseg时的等效的等向压缩线梯度；Me为弹性刚度张量；Mep为一致弹塑性模量张量，该模量被传递给有限元程序；Mtan为一致弹塑性切线模量张量；
Figure GDA0003778943730000332
为一致弹塑性切线模量的偏张量；
Figure GDA0003778943730000333
为一致弹塑性切线模量的球分量；
Figure GDA0003778943730000334
为体积弹性模量；nv为体积塑性流动方向；ns为剪切塑性流动方向；nαs为αs方向的单位向量；ν为泊松比；o(κeq)为远小于κeq的非零正数， o(κeq)∈(0,κeq×10-4]；pabs为绝对的有效平均应力；pabs.ini为体变起点的绝对的有效平均应力；p是相对于pabs.ini而增加或减少的静水压力；Pdev为求偏张量的投影算子；q为等效剪应力，岩土工程常将其称为广义剪应力；qseg为等效的等向压缩线梯度的分段点处的广义剪应力；
Figure GDA0003778943730000335
为材料在基准条件下单调压缩时的剪切硬化曲线的初值；
Figure GDA0003778943730000336
为材料在基准条件下的剪切硬化曲线的初始斜率；
Figure GDA0003778943730000337
为材料在基准条件下单调压缩时的剪切硬化曲线的上限；s为应力偏张量；σ为应力张量；σr为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力；
Figure GDA0003778943730000338
为从有限元程序传递到仿真步骤的周围压力；
Figure GDA0003778943730000339
为从有限元程序传递到仿真步骤的应力张量；
Figure GDA00037789437300003310
为从仿真步骤传递到有限元程序的应力张量； T1为由式(57)定义的函数；T2为由式(51)定义的函数；T7为由式(69)定义的函数；u为孔隙水压力；
Figure GDA00037789437300003311
为从有限元程序传递到仿真步骤的孔隙水压力；
Figure GDA00037789437300003312
为从有限元程序传递到仿真步骤的孔隙水压力的增量；Wsh为剪切硬化权重系数，Wsh∈[0,1]；Wvh为体积硬化权重系数，Wvh∈[0,1]；Wr为权重系数，Wr∈[0,1]，使得Mep的大小介于Mtan和Me之间；ζM为式(32)定义的函数；ζQ为式 (29)定义的函数；1为二阶单位张量，1的矩阵形式表示为[1 1 1 0 0 0]T。
符号和变量的补充说明，仿真步骤中大部分变量是由上述符号和上述变量复合而成，其含义也由各部分的含义复合而成。如
Figure GDA00037789437300003313
是由变量
Figure GDA00037789437300003314
符号n+1、符号trial复合而成，因此其含义为：一致弹塑性切线模量的偏张量，该变量处于当前增量步，该变量为弹性试探值。其余变量依此类推。
南宁圆砾Dr＝0.3试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数 (ea-εv-εa-N)关系的试验曲线与仿真曲线对比见图3；南宁圆砾Dr＝0.5试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数(ea-εv-εa-N)关系的试验曲线与仿真曲线对比见图4；南宁圆砾Dr＝0.7试样在循环荷载下的轴应变偏量—体积应变—轴应变—振动次数(ea-εv-εa-N)关系的试验曲线与仿真曲线对比见图5。因此采用本发明能够准确预测材料的长期累积轴向变形、剪切变形和体积变形。
南宁圆砾Dr＝0.3、Dr＝0.5、Dr＝0.7试样在循环荷载下结构各单元的轴应变偏量—振动次数 (ea-N)关系的仿真曲线分别见图6、图7、图8；南宁圆砾Dr＝0.3、Dr＝0.5、Dr＝0.7试样在循环荷载下结构各单元的侧应变偏量—振动次数(er-N)关系的仿真曲线分别见图9、图10、图11；南宁圆砾Dr＝0.3、Dr＝0.5、Dr＝0.7试样在循环荷载下结构各单元的轴应变—振动次数 (εa-N)关系的仿真曲线分别见图12、图13、图14；南宁圆砾Dr＝0.3、Dr＝0.5、Dr＝0.7试样在循环荷载下结构各单元的侧应变—振动次数(εr-N)关系的仿真曲线分别见图15、图16、图 17；南宁圆砾Dr＝0.3、Dr＝0.5、Dr＝0.7试样在循环荷载下结构各单元的体积应变—振动次数 (εv-N)关系的仿真曲线分别见图18、图19、图20。因此本发明能够对含有压硬性非线性变化和剪缩突变特性材料的具体结构的各部位在循环荷载下产生的累积变形进行仿真。
本发明所述“压硬性非线性变化”并不局限于抛物线变化，还指双曲线函数、指数函数、幂函数或对数函数变化。本发明也涵盖压硬性线性变化。本发明所述“剪缩突变”并不局限于体积应变曲线的一处突变，也指多处突变。只要是体积应变曲线不连续变化，都是本发明涵盖的范围。本发明也涵盖体积应变曲线无突变。

Claims (1)

1.压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法，包括基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应变驱动的仿真步骤，其特征在于：
基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的应变驱动的仿真步骤，以下简称“仿真步骤”，
仿真步骤具体为，当循环执行增量步时，在每个增量步内依顺序执行剪切弹塑性仿真步骤、体积弹塑性仿真步骤和向有限元法传递变量的仿真步骤：
A、剪切弹塑性仿真步骤
A.a、输入常数：CA,CB,CC,CD,CE,CF,Wsh,
Figure FDA0003778943720000011
v,eini,κeq，
输入变量：
Figure FDA0003778943720000012
Δγs,
Figure FDA0003778943720000013
Dr，
A.b、为判断剪切屈服做准备：
Figure FDA0003778943720000014
Figure FDA0003778943720000015
Figure FDA0003778943720000016
Figure FDA0003778943720000017
Figure FDA0003778943720000018
Figure FDA0003778943720000019
εn+1＝εn+Δεn+1 (7)
εv.n+1＝tr[εn+1] (8)
en+1＝εn+1-(εv.n+1/3)1 (9)
σ′n＝σn-un1 (10)
pabs.n＝tr[σ′n]|/3 (11)
Figure FDA00037789437200000111
Figure FDA00037789437200000112
Figure FDA00037789437200000113
un+1＝un+Δun+1 (15)
σ′r＝σr-un+1 (16)
Figure FDA00037789437200000114
Figure FDA00037789437200000115
Figure FDA00037789437200000116
Figure FDA00037789437200000117
Figure FDA00037789437200000118
Figure FDA00037789437200000119
若||αs.n||＝0，则nαs＝ns (23)
否则nαs＝αs.n/||αs.n|| (24)
Figure FDA0003778943720000021
Figure FDA0003778943720000022
Figure FDA0003778943720000023
Figure FDA0003778943720000024
Figure FDA0003778943720000025
Figure FDA0003778943720000026
Figure FDA0003778943720000027
若
Figure FDA0003778943720000028
执行步骤3；否则，执行步骤4
A.c、发生剪切屈服时的步骤
A.c.a、确定Δγs，
A.c.a.a、初始化
Figure FDA0003778943720000029
k＝0
A.c.a.b、迭代，执行如下牛顿迭代直到
Figure FDA00037789437200000210
小于预设容许值，k←k+1计算迭代
Figure FDA00037789437200000211
Figure FDA00037789437200000212
Figure FDA00037789437200000213
Figure FDA00037789437200000214
Figure FDA00037789437200000215
Figure FDA00037789437200000216
Figure FDA00037789437200000217
Figure FDA00037789437200000218
A.c.b、更新变量：若Δγs<0，则取Δγs＝0
Figure FDA00037789437200000219
Figure FDA00037789437200000220
αs.n+1＝ζM(αs.n+2CLΔγsns/3) (41)
Figure FDA0003778943720000031
Ks.n+1＝ζQ(Ks.n+2CPΔγs/3) (43)
Figure FDA0003778943720000032
A.c.c、计算一致弹塑性切线模量的偏张量
Figure FDA0003778943720000033
Figure FDA0003778943720000034
Figure FDA0003778943720000035
Figure FDA0003778943720000036
Figure FDA0003778943720000037
执行步骤A.e
A.d、不发生剪切屈服时的步骤：
Figure FDA0003778943720000038
执行步骤A.e
A.e、输出变量：sn+1、Δγs、
Figure FDA0003778943720000039
Figure FDA00037789437200000310
执行体积弹塑性仿真步骤，
B、体积弹塑性仿真步骤
B.a、输入常数：λeq1,λeq2,κeq,eini,Wvh,pabs.ini,qseg,emax,emin，
输入变量：sn+1,
Figure FDA00037789437200000311
Δγv,αv.n,Kv.n,
Figure FDA00037789437200000312
Dr,Wr，
B.b、为判断体积屈服做准备：
Figure FDA00037789437200000313
Figure FDA00037789437200000314
Figure FDA00037789437200000315
Δεv.n+1＝tr[Δεn+1] （50）
Figure FDA00037789437200000317
εn+1＝εn+Δεn+1 (7)
εv.n+1＝tr[εn+1] (8)
Figure FDA00037789437200000319
Figure FDA0003778943720000041
Figure FDA0003778943720000042
Figure FDA0003778943720000043
Figure FDA0003778943720000044
若
Figure FDA0003778943720000045
执行步骤3；否则，执行步骤4
B.c、发生体积屈服时的步骤：
B.c.a、确定Δγv，
B.c.a.a、初始化
Figure FDA0003778943720000046
k＝0
B.c.a.b、迭代，执行如下不动点迭代直到
Figure FDA0003778943720000047
小于预设容许值，k←k+1计算迭代
Figure FDA00037789437200000419
Figure FDA0003778943720000049
Figure FDA00037789437200000410
Figure FDA00037789437200000411
B.c.b、更新变量：若Δγv<0，则取Δγv＝0
Figure FDA00037789437200000412
p′n+1＝pabs.ini(T1-1) (58)
Δαv.n+1＝(1-Wvh)(1+eini)(pabs.ini+p′n+1)Δγvnv/T2 (60)
ΔKv.n+1＝Wvh(1+eini)(pabs.ini+p′n+1)Δγv/T2 (61)
αv.n+1＝αv.n+Δαv.n+1 (62)
Kv.n+1＝Kv.n+ΔKv.n+1 (63)
Figure FDA00037789437200000413
执行步骤B.e
B.d、不发生体积屈服时的步骤：
Figure FDA00037789437200000414
Figure FDA00037789437200000415
Figure FDA00037789437200000416
Dr＝(emax-e)/(emax-emin) (68)
执行步骤B.e
B.e、计算一致弹塑性切线模量的球分量
Figure FDA00037789437200000417
Figure FDA00037789437200000418
Figure FDA0003778943720000051
B.f、输出变量：p′n+1,Δγv,αv.n+1,Kv.n+1,
Figure FDA0003778943720000052
Dr，执行向有限元法传递变量的仿真步骤，
C、向有限元法传递变量的仿真步骤
C.a、输入常数：pabs.ini
输入变量：
Figure FDA0003778943720000053
p′n+1,sn+1,
Figure FDA0003778943720000054
Wr，
C.b、
Figure FDA0003778943720000055
Figure FDA0003778943720000056
Figure FDA0003778943720000057
Figure FDA0003778943720000058
un+1＝un+Δun+1 (15)
p′abs.n+1＝pabs.ini+p′n+1 (72)
σ′n+1＝sn+1+p′abs.n+11 (73)
σn+1＝σ′n+1+un+11 (74)
εn+1＝εn+Δεn+1 (7)
Figure FDA0003778943720000059
Figure FDA00037789437200000510
Figure FDA00037789437200000511
Figure FDA00037789437200000512
Figure FDA00037789437200000513
Figure FDA00037789437200000514
C.c、输出变量：
Figure FDA00037789437200000515
结束当前增量步，
以上仿真步骤中符号的含义：变量右下标n指上一增量步；变量右下标n+1指当前增量步；变量右上标trial指该变量为弹性试探值；变量右上标*指该变量处于基准条件；变量前的Δ指该变量为增量；符号:为内积符号，是对张量缩并；变量右上标(k)指第(k)次牛顿迭代；变量右上标'指该变量为有效应力；||·||指二范数；tr[·]指对张量求迹；sign(·)为符号函数；符号
Figure FDA00037789437200000516
为张量积；exp(·)为以自然常数e为底的指数函数，
以上仿真步骤中的变量，粗体符号为张量，非粗体符号为标量，变量的含义：αs为试样在实际条件下的背应力偏张量；
Figure FDA00037789437200000517
为试样在基准条件下的背应力偏张量；αv为体积背应力；Bs为与周围压力有关的比例系数；CA、CB、CC、CD、CE、CF为剪切屈服条件参数；CL、CM为A-F模型的随动硬化参数；CP、CQ为Chaboche模型的等向硬化参数；Dr为相对密实度；Ds为与相对密实度有关的比例系数；e为应变偏张量；ee为弹性应变偏张量；ep为塑性应变偏张量；e为孔隙比；eini为体变起始点的孔隙比；emax为最大孔隙比；emin为最小孔隙比；ε为应变张量；εp为塑性应变张量；εv为体积应变；
Figure FDA0003778943720000061
为弹性体积应变；
Figure FDA0003778943720000062
为塑性体积应变；
Figure FDA0003778943720000063
为从有限元程序传递到仿真步骤的应变张量；
Figure FDA0003778943720000064
为从有限元程序传递到仿真步骤的应变张量的增量；
Figure FDA0003778943720000065
为从仿真步骤传递到有限元程序的应变张量；fs为剪切屈服函数；fv为体积屈服函数；G为剪切弹性模量；γs为剪切塑性滑移率；γv为体积塑性滑移率；I为四阶单位张量，I的矩阵形式表示为6×6的对角矩阵，对角线上的元素均为1；k为迭代次数指示变量；Ks为试样在实际条件下的剪切等向塑性流动应力，即试样在实际条件下屈服时的广义剪应力的各向同性硬化部分；
Figure FDA0003778943720000066
为试样在基准条件下的剪切等向塑性流动应力，即试样在基准条件下屈服时的广义剪应力的各向同性硬化部分；Kv为体积等向塑性流动应力；κeq为等效的等向膨胀线的梯度；ξs为相对应力的偏张量；ξv为相对球应力；λeq1为q≤qseg时的等效的等向压缩线梯度；λeq2为q>qseg时的等效的等向压缩线梯度；Me为弹性刚度张量；Mep为一致弹塑性模量张量，该模量被传递给有限元程序；Mtan为一致弹塑性切线模量张量；
Figure FDA0003778943720000067
为一致弹塑性切线模量的偏张量；
Figure FDA0003778943720000068
为一致弹塑性切线模量的球分量；
Figure FDA0003778943720000069
为体积弹性模量；nv为体积塑性流动方向；ns为剪切塑性流动方向；nαs为αs方向的单位向量；v为泊松比；o(κeq)为远小于κeq的非零正数，o(κeq)∈(0,κeq×10-4]；pabs为绝对的有效平均应力；pabs.ini为体变起点的绝对的有效平均应力；p是相对于pabs.ini而增加或减少的静水压力；Pdev为求偏张量的投影算子；q为等效剪应力，岩土工程常将其称为广义剪应力；qseg为等效的等向压缩线梯度的分段点处的广义剪应力；
Figure FDA00037789437200000610
为材料在基准条件下单调压缩时的剪切硬化曲线的初值；
Figure FDA00037789437200000611
为材料在基准条件下的剪切硬化曲线的初始斜率；
Figure FDA00037789437200000612
为材料在基准条件下单调压缩时的剪切硬化曲线的上限；s为应力偏张量；σ为应力张量；σr为三轴压缩试验和振动三轴试验中试样受到的有效的周围压力；
Figure FDA00037789437200000613
为从有限元程序传递到仿真步骤的周围压力
Figure FDA00037789437200000614
为从有限元程序传递到仿真步骤的应力张量；
Figure FDA00037789437200000615
为从仿真步骤传递到有限元程序的应力张量；T1为由式(57)定义的函数；T2为由式(51)定义的函数；T7为由式(69)定义的函数；u为孔隙水压力；
Figure FDA00037789437200000616
为从有限元程序传递到仿真步骤的孔隙水压力；
Figure FDA00037789437200000617
为从有限元程序传递到仿真步骤的孔隙水压力的增量；Wsh为剪切硬化权重系数，Wsh∈[0,1]；Wvh为体积硬化权重系数，Wvh∈[0,1]；Wr为权重系数，Wr∈[0,1]，使得Mep的大小介于Mtan和Me之间；ζM为式(32)定义的函数；ζQ为式(29)定义的函数；1为二阶单位张量，1的矩阵形式表示为[1 1 1 0 0 0]T，
符号和变量的补充说明，仿真步骤中大部分变量是由上述符号和上述变量复合而成，其含义也由各部分的含义复合而成，如
Figure FDA00037789437200000618
是由变量
Figure FDA00037789437200000619
符号n+1、符号trial复合而成，因此其含义为：一致弹塑性切线模量的偏张量，该变量处于当前增量步，该变量为弹性试探值，其余变量依此类推。
CN202010538076.5A 2020-06-12 2020-06-12 压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法 Active CN111783332B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010538076.5A CN111783332B (zh) 2020-06-12 2020-06-12 压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010538076.5A CN111783332B (zh) 2020-06-12 2020-06-12 压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法

Publications (2)

Publication Number Publication Date
CN111783332A CN111783332A (zh) 2020-10-16
CN111783332B true CN111783332B (zh) 2022-10-21

Family

ID=72756247

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010538076.5A Active CN111783332B (zh) 2020-06-12 2020-06-12 压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法

Country Status (1)

Country Link
CN (1) CN111783332B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113155612B (zh) * 2021-04-16 2022-09-30 浙江科技学院 微纤维混合硅溶胶固化钙质砂的变形预报方法
CN114912314B (zh) * 2022-04-21 2024-04-02 中国科学院武汉岩土力学研究所 岩土介质弹塑性本构关系隐式自适应应力积分计算方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN203594126U (zh) * 2013-10-10 2014-05-14 上海大学 厚度以二次抛物线变化的钢剪切板阻尼器
CN109580388A (zh) * 2019-01-21 2019-04-05 广西大学 一种岩土材料剪切屈服面与体积屈服面的测定方法
CN111062162A (zh) * 2019-12-12 2020-04-24 王靖涛 一种岩土材料精确本构模型的数值建模与应用方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN203594126U (zh) * 2013-10-10 2014-05-14 上海大学 厚度以二次抛物线变化的钢剪切板阻尼器
CN109580388A (zh) * 2019-01-21 2019-04-05 广西大学 一种岩土材料剪切屈服面与体积屈服面的测定方法
CN111062162A (zh) * 2019-12-12 2020-04-24 王靖涛 一种岩土材料精确本构模型的数值建模与应用方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
In situ study of static and dynamic strain energy density at notch roots and fatigue cracks using digital image correlation;Holycross, C.M等;《2016 | 57th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference》;20160225;1-19 *
土本构关系及数值建模;张光永;《中国博士学位论文全文数据库 (工程科技Ⅱ辑)》;20090515(第5期);C038-35 *

Also Published As

Publication number Publication date
CN111783332A (zh) 2020-10-16

Similar Documents

Publication Publication Date Title
Zheng et al. Reformulation of dynamic crack propagation using the numerical manifold method
CN111783282B (zh) 基于压硬性非线性变化和剪缩突变特性材料的振动累积变形的仿真方法
Onate et al. A local constitutive model for the discrete element method. Application to geomaterials and concrete
Wang et al. A constitutive model for rock interfaces and joints
Li et al. Nonlinear vibration analysis of fiber metal laminated plates with multiple viscoelastic layers
Zhao et al. A generic approach to modelling flexible confined boundary conditions in SPH and its application
Goodman et al. Finite element analysis for discontinuous rocks
CN111783332B (zh) 压硬性非线性变化和剪缩突变特性材料振动累积变形的有限元仿真方法
Uzuoka et al. Three-dimensional numerical simulation of earthquake damage to group-piles in a liquefied ground
Tamagnini et al. Plasticity with generalized hardening: constitutive modeling and computational aspects
Cheng et al. Cavity expansion in unsaturated soils of finite radial extent
Fernandes et al. Multi-scale modelling for bending analysis of heterogeneous plates by coupling BEM and FEM
Loukidis Advanced constitutive modeling of sands and applications to foundation engineering
Jalali et al. Using finite element method for pile-soil interface (through PLAXIS and ANSYS)
Lackner et al. Constitutive modeling of cementitious materials in the framework of chemoplasticity
Yuan et al. Stabilized smoothed particle finite element method for coupled large deformation problems in geotechnics
Ebrahimian et al. A numerical study on interface shearing of granular Cosserat materials
Ma et al. Unified elastoplastic finite difference and its application
Li et al. Limit deformation analysis of unsaturated expansive soils during wetting and drying cycles
Yao et al. An SBFEM-Based model for hydraulic fracturing in quasi-brittle materials
Lagioia et al. The ‘I 3’generalization of the Galileo–Rankine tension criterion
Zhang et al. Generalized plasticity model considering grain crushing and anisotropy for rockfill materials
Dong Performance of explicit substepping integration scheme for complex constitutive models in finite element analysis
Kolmayer et al. Numerical implementation of a new rheological law for argilites
Molenkamp Application of non‐linear elastic model

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