岩土介质弹塑性本构关系隐式自适应应力积分计算方法技术

本发明专利技术公开了一种岩土介质弹塑性本构关系隐式自适应应力积分计算方法，步骤是：S1：基于增量步初始时刻的应力张量和应变增量张量，计算试探弹性应力张量；S2：判定增量步内岩土介质加载状态，弹性加载时，将增量步结束时刻的应力张量更新为试探弹性应力张量；S3：当加载状态从弹性过渡至弹塑性，确定过渡点的应力张量；S4：采用2级2阶对角隐式Runge

【技术实现步骤摘要】
岩土介质弹塑性本构关系隐式自适应应力积分计算方法


[0001]本专利技术属于计算岩土力学
，更具体是涉及一种岩土介质弹塑性本构关系隐式自适应应力积分计算方法。

技术介绍

[0002]随着计算机技术的高速发展，采用数值分析方法对岩土介质进行非线性分析，是岩土工程领域的一个重要研究方向。岩土介质的本构关系是描述其强度与变形特性的重要工具，是岩土工程问题数值分析中的核心内容之一。应力积分计算是借助一定的数值分析方法对本构关系进行积分，进而得到已知应变增量对应的应力增量。它是本构关系数值化过程中最重要的环节之一。
[0003]岩土介质具有复杂的力学特性，如各向异性、剪胀性、结构性等，大量学者提出了众多的弹塑性本构关系对其进行描述。目前，弹塑性本构关系的应力积分方法主要可分为显式方法和隐式方法两大类。经典的显式方法主要有具有一阶精度的向前欧拉法和具有二阶精度的改进向前欧拉法。采用这两种方法在进行应力更新计算时，一般不需要求解非线性方程组，具有单个增量步计算量小且易于数值实现的优点。然而，这两种方法是有条件稳定的，在实际应用中，为防止出现较大的累积误差，需设置很小的积分步长，导致其计算效率较低。为了提高显式算法的计算效率，一些学者提出除了显式自适应应力积分方法，如将经典欧拉法和改进向前欧拉法相结合而提出的自带误差控制的子步积分法(Sloan S.W.,Abbo A.J.,&D Sheng.2001.Refined explicit integration of elastoplastic models with automatic error control.Engineering Computations,18(1/2):121
‑
194)，以及将经典向前欧拉法和Runge
‑
Kutta法相结合而提出的子步积分法(Farias M.M.,Pedroso,D.M.,&Nakai T.2009.Automatic substepping integration of the subloading TIJ model with stress path dependent hardening.Computers and Geotechnics,36(4):537
‑
548.)等。这类应力更新方法虽然在一定程度上提高了计算效率，但显式算法的特性使得其不能严格满足屈服函数的一致性条件，导致误差积累和精度低的缺点仍然无法完全消除。相较于显式方法，隐式方法强化了在时间步结束时屈服函数的一致性条件，从而避免了屈服面的漂移，因此该法是强健的。目前，主流有限元软件采用的是向后欧拉法，它属于一阶方法，其精度不高。对于一些高度非线性的本构关系，如弹塑性刚度算子与应力增量方向相关的模型(Dafalias Y.F.,Taiebat M.2016.SANISAND
‑
Z:zero elastic range sand plasticity model.G
é
otechnique,66(12):999
‑
1013)，刚度算子与应力增量方向的相关性导致模型只能采用隐式方法进行应力更新计算。若采用具有一阶精度的向后欧拉法进行应力更新计算，会出现应力更新结果明显依赖于积分步长的现象(Petalas A.L.,&Dafalias Y.F.2019.Implicit integration of incrementally non
‑
linear,zero
‑
elastic range,bounding surface plasticity.Computers and Geotechnics,112,386
‑
402)，因而无法有效地兼顾计算效率和计算精度。

技术实现思路

[0004]为了解决现有显式自适应法、一阶向后欧拉隐式法等岩土介质弹塑性本构关系应力积分方法在计算效率和精度等方面的不足，本专利技术的目的是在于提供了一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法，方法易行，操作简便，适用范围广，计算精度高，效率更高的应力。
[0005]为了实现上述的目的，本专利技术采取如下的技术方案：
[0006]对于岩土介质弹塑性本构关系，针对任一增量步[t
n
，t
n+1
]，已知增量步初始t
n
时刻的柯西应力张量σ
n
和包含m个分量的硬化参数向量以及增量步内的应变增量张量Δε，采用隐式自适应应力积分方法计算增量步结束t
n+1
时刻的柯西应力张量σ
n+1
和硬化参数向量H
n+1
，一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法，其步骤是：
[0007]S1：计算增量内的试探弹性应力张量，即基于t
n
时刻的柯西应力张量σ
n
和增量步内的应变增量张量Δε，计算t
n+1
时刻的试探弹性应力张量
[0008][0009]式中：σ
n
为柯西应力张量、Δs为应变增量张量、为试探弹性应力张量、D
e
为岩土介质弹塑性本构关系中描述弹性应力应变关系的四阶刚度张量；“：”表示张量内积；
[0010]S2：判定增量步内岩土介质的加载状态，即根据岩土介质弹塑性本构关系中的屈服面函数f(σ，H)，分别计算增量步初始和结束时刻的屈服面函数大小：f(σ
n
，H
n
)和
[0011]若则增量步内岩土介质处于完全弹性状态，增量步结束时刻t
n+1
的柯西应力张量σ
n+1
和硬化张量向量H
n+1
分别为：H
n+1
＝H
n
；跳转至步骤S5；
[0012]若进入下一步；
[0013]S3：确定增量步内岩土介质由弹性加载过渡至弹塑性加载时过渡点的应力状态，即当f(σ
n
，H
n
)＜
‑
FTOL，表明增量步内岩土介质从完全弹性加载状态过渡至弹塑性加载状态，过渡点处柯西应力张量σ
int
满足如下条件：
[0014]f(σ
int
，H
n
)＝f(σ
n
+A
·
D
e
：Δε，H
n
)＝0
ꢀꢀꢀꢀ
(2)
[0015]式中：A为反映增量步内过渡点位置的标量值；对于非线性方程(2)，采用二分法迭代求取A值；随后，视过渡点为增量步起始点，将将增量步初始时刻t
n
更新为t
n
+A(t
n+1
‑
t
n
)，柯西应力张量σ
n
更新为本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法，其步骤是：S1：计算增量内的试探弹性应力张量，即基于增量步初始t
n
时刻的应力张量σ
n
和增量步应变增量张量Δε，计算增量步结束t
n+1
时刻的试探弹性应力张量时刻的试探弹性应力张量式中：D
e
为岩土介质弹塑性本构关系中描述弹性应力应变关系的四阶刚度张量、符号“：”为张量内积；S2：判定增量步内岩土介质的加载状态，即根据弹塑性本构关系中的屈服面函数，分别计算增量步初始和结束时刻的屈服面函数大小：f(σ
n
，H
n
)和其中H
n
表示增量步初始t
n
时刻的硬化参数向量，包含m个分量，即即则增量步内岩土介质处于完全弹性状态，增量步结束时刻t
n+1
的柯西应力张量σ
n+1
和硬化张量向量H
n+1
分别为：H
n+1
＝H
n
；跳转至步骤S5；则进入下一步；S3：确定增量步内岩土介质由弹性加载过渡至弹塑性加载时过渡点的应力状态，即当f(σ
n
，H
n
)＜
‑
FTOL时，表明增量步内岩土介质从完全弹性加载状态过渡至弹塑性加载状态，过渡点处柯西应力张量σ
int
如下：f(σ
int
，H
n
)＝f(σ
n
+A
·
D
e
：Δε，H
n
)＝0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)式中：A为反映增量步内过渡点位置的标量值；对于非线性方程(2)，采用二分法迭代求取A值；随后，视过渡点为增量步起始点，将将增量步初始时刻t
n
更新为t
n
+A(t
n+1
‑
t
n
)，柯西应力张量σ
n
更新为σ
int
＝σ
n
+A
·
D
e
：Δε，硬化参数向量仍为H
n
，增量步内的应变增量张量Δε更新为(1
‑
A)
·
Δε；进入下一步；若|f(σ
n
，H
n
)|≤FTOL，则增量步内岩土介质处于弹塑性加载状态，进入下一步；S4：采用隐式自适应积分方法对弹塑性加载状态的本构关系进行应力积分计算，即采用2级2阶对角隐式Runge
‑
Kutta法分2级对弹塑性加载状态下的应力应变关系进行积分计算；依据第1级和第2级应力积分计算结果，确定应力积分计算的局部相对误差标量值，依据这一误差标量值自动控制积分步长，实现自适应应力积分计算；S5：更新增量步结束时刻的柯西应力张量σ
n+1
、硬化参数向量H
n+1
。2.根据权利要求1所述的一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法，其特征在于：所述步骤(S4)中的弹塑性加载状态下的应力应变关系积分计算涉及如下约束条件：约束条件：f(σ，H)＝0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)式中：和分别为柯西应力张量和硬化参数向量随时间的变化率；和分别为总应变张量和塑性应变张量随时间变化率；表示塑性势函数g(σ，H)对柯西应力张量的偏导数，为一描述塑性应变流动方向的2阶张量；为塑性乘子随时间变化率；为塑性乘子随时间变化率；其中，为与弹塑性本构关系中第i个硬化参数相关且形式已知的硬化函数；f(σ，H)为屈服面函数。3.根据权利要求1所述的一种岩土介质弹塑性本构关系隐式自适应应力更新计算方
法，其特征在于：所述步骤(S4)中的2级2阶对角隐式Runge
‑
Kutta法中采用的Butcher为：式中：Butcher表包含两种计算，主方案计算采用的系数为a、b和c，次方案计算采用的系数是a、和c；利用主方案和次方案计算精度的不同，确定出应力积分计算的局部相对误差。4.根据权利要求1所述的一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法，其特征在于：所述步骤(S4)中的采用2级2阶对角隐式Runge
‑
Kutta法分2级对岩土体弹塑性应力应变关系进行积分计算中，第1级应力的计算为求解下述非线性方程组：第1级应力的计算为求解下述非线性方程组：f(σ
n1
，H
n1
)＝0
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)式中：σ
n1
、H
n1
分别为柯西应力张量σ和硬化参数向量H在t
n1
＝t
n
+α
·
(t
n+1
‑
t
n
)时刻值；为t
n1
时刻塑性势函数g(σ，H)对柯西应力张量的偏导数；为t
n1
时刻向量函数的值；Δλ
I
为与第1级应力积分求解相关的增量塑性乘子；α为Butcher表中参数。5.根据权利要求4所述的一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法，其特征在于：所述的柯西应力张量具有对称性，故非线性方程组(7)、(8)、(9)共包含m+7个方程，待定未知量为m+7个，分别是6个应力分量、m个硬化参数以及增量塑性乘子，采用牛顿迭代法求解非线性方程组(7)、(8)、(9)，获取σ
n1
、H
n1
和Δλ
I

【专利技术属性】
技术研发人员：陈成，罗玛诗艺，吴勋，张先伟，刘晓敏，
申请(专利权)人：中国科学院武汉岩土力学研究所，
类型：发明
国别省市：