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

摘要

本发明公开了一种岩土介质弹塑性本构关系隐式自适应应力积分计算方法，步骤是：S1：基于增量步初始时刻的应力张量和应变增量张量，计算试探弹性应力张量；S2：判定增量步内岩土介质加载状态，弹性加载时，将增量步结束时刻的应力张量更新为试探弹性应力张量；S3：当加载状态从弹性过渡至弹塑性，确定过渡点的应力张量；S4：采用2级2阶对角隐式Runge‑Kutta积分法对弹塑性本构关系进行积分步长可控的自适应应力积分计算；S5：更新增量步结束时刻的应力张量、硬化参数向量。相较于显式应力积分法和主流有限元采用的向后欧拉隐式应力积分法，本发明公开的方法操作简便且易行，具有适用性强、计算精度高以及计算效率高的优势。

说明书

技术领域

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

背景技术

随着计算机技术的高速发展，采用数值分析方法对岩土介质进行非线性分析，是岩土工程领域的一个重要研究方向。岩土介质的本构关系是描述其强度与变形特性的重要工具，是岩土工程问题数值分析中的核心内容之一。应力积分计算是借助一定的数值分析方法对本构关系进行积分，进而得到已知应变增量对应的应力增量。它是本构关系数值化过程中最重要的环节之一。

岩土介质具有复杂的力学特性，如各向异性、剪胀性、结构性等，大量学者提出了众多的弹塑性本构关系对其进行描述。目前，弹塑性本构关系的应力积分方法主要可分为显式方法和隐式方法两大类。经典的显式方法主要有具有一阶精度的向前欧拉法和具有二阶精度的改进向前欧拉法。采用这两种方法在进行应力更新计算时，一般不需要求解非线性方程组，具有单个增量步计算量小且易于数值实现的优点。然而，这两种方法是有条件稳定的，在实际应用中，为防止出现较大的累积误差，需设置很小的积分步长，导致其计算效率较低。为了提高显式算法的计算效率，一些学者提出除了显式自适应应力积分方法，如将经典欧拉法和改进向前欧拉法相结合而提出的自带误差控制的子步积分法(Sloan S.W.,Abbo A.J.,&D Sheng.2001.Refined explicit integration of elastoplastic modelswith 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 TIJmodel with stress path dependent hardening.Computers and Geotechnics,36(4):537-548.)等。这类应力更新方法虽然在一定程度上提高了计算效率，但显式算法的特性使得其不能严格满足屈服函数的一致性条件，导致误差积累和精度低的缺点仍然无法完全消除。相较于显式方法，隐式方法强化了在时间步结束时屈服函数的一致性条件，从而避免了屈服面的漂移，因此该法是强健的。目前，主流有限元软件采用的是向后欧拉法，它属于一阶方法，其精度不高。对于一些高度非线性的本构关系，如弹塑性刚度算子与应力增量方向相关的模型(Dafalias Y.F.,Taiebat M.2016.SANISAND-Z:zero elastic range sandplasticity model.Géotechnique,66(12):999-1013)，刚度算子与应力增量方向的相关性导致模型只能采用隐式方法进行应力更新计算。若采用具有一阶精度的向后欧拉法进行应力更新计算，会出现应力更新结果明显依赖于积分步长的现象(Petalas A.L.,&DafaliasY.F.2019.Implicit integration of incrementally non-linear,zero-elastic range,bounding surface plasticity.Computers and Geotechnics,112,386-402)，因而无法有效地兼顾计算效率和计算精度。

发明内容

为了解决现有显式自适应法、一阶向后欧拉隐式法等岩土介质弹塑性本构关系应力积分方法在计算效率和精度等方面的不足，本发明的目的是在于提供了一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法，方法易行，操作简便，适用范围广，计算精度高，效率更高的应力。

为了实现上述的目的，本发明采取如下的技术方案：

对于岩土介质弹塑性本构关系，针对任一增量步[t

S1：计算增量内的试探弹性应力张量，即基于t

式中：σ

S2：判定增量步内岩土介质的加载状态，即根据岩土介质弹塑性本构关系中的屈服面函数f(σ，H)，分别计算增量步初始和结束时刻的屈服面函数大小：f(σ

若

若

S3：确定增量步内岩土介质由弹性加载过渡至弹塑性加载时过渡点的应力状态，即当f(σ

f(σ

式中：A为反映增量步内过渡点位置的标量值；对于非线性方程(2)，采用二分法迭代求取A值；随后，视过渡点为增量步起始点，将将增量步初始时刻t

若|f(σ

S4：采用隐式自适应积分方法对弹塑性加载状态的本构关系进行应力积分计算，即采用2级2阶对角隐式Runge-Kutta法分2级对弹塑性加载状态下的应力应变关系进行积分计算；依据第1级和第2级应力积分计算结果，确定应力积分计算的局部相对误差标量值，并依据这一误差标量值自动控制积分步长，实现自适应应力积分计算；

S5：更新增量步结束时刻的柯西应力张量σ

进一步地，上述步骤S4中的弹塑性加载状态下的应力应变关系积分计算涉及如下约束条件：

f(σ，H)＝0 (5)

式中：

优选地，上述步骤S4中的2级2阶对角隐式Runge-Kutta法中采用的Butcher表为：

式中：Butcher表包含了两种计算方案。主方案计算采用的系数为a、b和c，具有二阶精度，是A-，S-和L-稳定的。次方案计算采用的系数是a、

进一步地，上述步骤S4中的采用2级2阶对角隐式Runge-Kutta法分2级对岩土体弹塑性应力应变关系进行积分计算具体包括：

计算第1级应力：

f(σ

式中：σ

考虑到柯西应力张量的对称性，非线性方程组(7)、(8)、(9)共包含m+7个方程，待定未知量也为m+7个，分别是6个应力分量、m个硬化参数以及增量塑性乘子。采用牛顿迭代法求解非线性方程组(7)、(8)、(9)，获取σ

计算第2级应力：

f(σ

式中：σ

非线性方程组(10)、(11)、(12)共包含m+7个方程，待定未知量也为m+7个。采用牛顿迭代法求解非线性方程组(10)、(11)、(12)，获取σ

进一步地，由于t

进一步地，上述步骤S4中的依据第1级和第2级应力积分计算结果，确定应力积分计算的局部相对误差标量值R

考虑到式(6)中的主方案具有二阶精度，而次方案具有一阶精度，因而，可利用两种方案的应力张量和硬化参数向量计算结果的差值确定应力更新的局部相对误差。由于两种方案中仅系数b和

式中：

R

式中：MAX{x，y}表示取x和y两者的大值；||X||

进一步地，上述步骤S4中的采用2级2阶对角隐式Runge-Kutta法分2级对岩土体弹塑性应力应变关系进行积分计算；依据第1级和第2级应力积分计算结果，确定应力积分计算的局部相对误差标量值，并依据这一误差标量值自动控制积分步长，实现自适应应力积分计算的完整步骤是：

S41：初始化：

设置T＝0，ΔT＝1；

设置Δε

设置误差允许值R

S42：弹塑性本构关系应力积分与局部相对误差计算：

将Δε更新为ΔT·Δε；

采用2级2阶对角隐式Runge-Kutta法对岩土介质弹塑性本构关系进行应力积分计算，确定柯西应力张量σ和硬化参数向量H在t

根据第1级和第2级柯西应力张量和硬化参数向量计算结果，计算σ

S43：自动控制积分步长：

若R

若R

S44：判断应力积分计算是否结束

若T＜1，则应力积分计算未结束，转至S42；

若T＝1，则应力积分计算结束。

所述上述的技术措施，最关键的步骤是S4，通过采用2级2阶对角隐式Runge-Kutta积分方法，不仅可以进行高精度的应力积分计算，还可以通过局部相对误差实现可控步长的自适应应力积分计算，解决了现有岩土介质弹塑性本构关系应力积分算法无法有效兼具计算精度与计算效率的难题，是本发明的核心创新点。

本发明与现有技术相比，具有以下优点和有益效果：

①步骤S42采用2级2阶对角隐式Runge-Kutta法分2级对岩土体弹塑性本构关系进行应力积分计算，保证了本发明提供的应力积分方法具有2阶计算精度，有效解决了现有主流有限元软件中采用的向后欧拉应力积分法计算精度较低(1阶计算精度)的技术问题；

②应力积分方法的隐式特征使其适用性较显式应力积分方法更广，可有效解决岩土介质增量非线性型弹塑性本构关系应力积分这一技术难点；

③步骤S43使得本发明提供的应力积分方法具有自适应特征，可根据局部相对误差自动控制积分的步长，其计算效率高，解决了传统隐式应力积分方法计算效率较低的技术问题。

因此，本发明提供的岩土介质弹塑性本构关系隐式自适应应力积分计算方法，通过步骤S1、S2和S3对岩土介质加载变形过程中屈服状态的判定与划分，以步骤S4为核心，采用2级2阶对角隐式Runge-Kutta法分2级对岩土体弹塑性本构关系进行应力积分计算，利用第1级和第2级应力积分计算结果确定的应力积分局部相对误差，并基于局部相对误差实现自适应应力积分计算，具有计算精度高以及计算高效的技术特点，可达到快速准确地对应力进行更新计算的技术效果。

附图说明

图1为一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法的流程图。

图2本发明实施例的岩土介质屈服状态判断示意图。

图3本发明实施例的2阶2阶对角隐式Runge-Kutta积分法示意图。

图4本发明实施例的Sanisand-Z本构模型在偏平面内的示意图。

图5(a、b)本发明实施例的不排水三轴压缩试验计算结果。

图6(a、b)本发明实施例的不排水循环三轴压缩试验计算结果。

具体实施方式

为了更好地理解上述技术方案，下面结合附图和实施例对本发明的岩土介质弹塑性本构关系隐式自适应应力积分计算方法的性能作进一步的详细说明。以下实施例用于说明本发明，但不能用来限制本发明的适用范围。

实施例1：

一种岩土介质弹塑性本构关系隐式自适应应力积分计算方法，即在任一增量步[t

根据图1可知，一种岩土介质弹塑性本构关系隐式自适应应力更新计算方法，其步骤是：

S1：计算增量内的试探弹性应力张量，即基于t

式中：D

S2：判定增量步内岩土介质的加载状态，即根据岩土介质弹塑性本构关系中的屈服面函数f(σ，H)，分别计算增量步初始和结束时刻的屈服面函数大小：f(σ

若

若

S3：确定增量步内岩土介质由弹性加载过渡至弹塑性加载时过渡点的应力状态，即当f(σ

f(σ

式中：A为反映增量步内过渡点位置的标量值；

对于非线性方程(2)，采用二分法迭代求取A值，即首先假定A＝0.5，若此时f(σ

为便于后续弹塑性加载状态的应力积分计算，视过渡点为增量步起始点，即将增量步初始t

若|f(σ

S4：采用隐式自适应积分方法对弹塑性加载状态的本构关系进行应力积分计算，即采用2级2阶对角隐式Runge-Kutta法分2级对弹塑性加载状态下的应力应变关系进行积分计算，如图3所示，即先基于增量步初始t

S41：初始化，即设置T＝0，ΔT＝1，Δε

S42：弹塑性本构关系应力积分与局部相对误差计算，其中，弹塑性加载阶段的应力积分涉及如下三个约束条件，即

f(σ，H)＝0 (5)

式中：

此外，采用2级2阶对角隐式Runge-Kutta法分2级对岩土体弹塑性应力应变关系进行积分计算时，采用的Butcher表为：

式中：Butcher表包含了两种计算方案。主方案计算采用的系数为a、b和c，具有二阶精度，是A-，S-和L-稳定的；次方案计算采用的系数是a、

对于t

f(σ

式中：σ

考虑到柯西应力张量的对称性，非线性方程组(7)、(8)、(9)共包含m+7个方程，待定未知量也为m+7个，分别是6个应力分量、m个硬化参数以及增量塑性乘子。采用牛顿迭代法求解非线性方程组(7)、(8)、(9)，获取σ

基于第1级应力积分计算获取σ

f(σ

式中：σ

同理，非线性方程组(10)、(11)、(12)也包含m+7个方程和m+7个待定未知量；采用牛顿迭代法求解非线性方程组(10)、(11)、(12)，获取σ

由于t

考虑到式(6)中的主方案具有二阶精度，而次方案具有一阶精度，因而，可利用两种方案的应力张量和硬化参数向量计算结果的差值确定应力更新的局部相对误差。由于两种方案中仅系数b和

式中：

R

式中：MAX{x，y}表示取x和y两者的大值；||X||

S43：自动控制积分步长，即通过对比应力积分的局部相对误差R

若R

若R

S44：判断应力积分计算是否结束，即根据T值的大小，判断弹塑性加载状态下的应力积分计算是否结束，具体为：

若T＜1，则应力积分计算未结束，转至步骤S42；

若T＝1，则应力积分计算结束；

S5：更新增量步结束t

通过逐级执行上述岩土介质弹塑性本构关系隐式自适应应力积分计算方法的步骤，对岩土工程数值分析中广泛使用的Sanisand系列模型中的Sanisand-Z边界面弹塑性本构模型进行了应力积分计算。Sanisand-Z边界面弹塑性本构模型(Dafalias Y.F.，TaiebatM.2016.SANISAND-Z：zero elastic range sand plasticity model.Géotechnique，66(12)：999-1013)包括边界面、剪胀面和临界状态面，如图4所示(图中r和

值得注意的是，Sanisand-Z边界面弹塑性本构模型中，描述土体弹塑性应力应变关系的四阶弹塑性刚度张量D

式中：D

综上所述，本发明提供的一种岩土介质弹塑性本构关系应力积分计算方法，无论在计算精度或计算效率方面，还是在适用性方面，均较隐式自适应应力积分法或主流有限元软件采用的隐式向后欧拉应力积分法更具优势，是一种极具推广价值的岩土介质弹塑性本构关系应力积分计算方法。

以上所述仅为本发明的优选实施例而已，并不用于限制发明，凡是根据本发明实质对以上实施例做任何简单修改、变更以及等效结构变换，均仍属本发明技术方案的保护范围内。