应变空间简化多重本构模型及应用

5.3.1 引言

目前土体的弹塑性模型中，多数是建立于应力空间上的，虽然一些论著从理论上论证了在应变空间上建立弹塑性本构模型的优越性[25，28]，但由于按经典的弹塑性本构理论来建立模型时，需要确定应变空间上的塑性势函数或屈服函数，但要通过试验来确定这两个函数并不是一件容易的事。因此，真正能方便用于实际分析的应变空间弹塑性模型，目前还较少，现有的一些应变空间模型中，其应变空间上的塑性势函数多是从应力空间中的塑性势函数变换而来的。因此，为克服推求应变空间塑性势函数的困难和建立更简便且又更具有一般性的模型，这里采用以上所建立的多重势面理论[114]或张量定律[119]来建立应变空间的弹塑性本构模型，并用于有限元分析，与邓肯—张模型进行比较[122]。

5.3.2 应变空间上的多重势面弹塑性本构方程

作为有限元等数值方法应用的本构模型，其关键是确定增量本构方程：

{dσ}=[dσx，dσy，dσz，dτxy，dτyz，dτzx]T为应力增量，{dε}=[dεx，dεy，dεz，dγxy，dγyz，dγzx]T为应变增量，[Dep]为6×6阶的弹塑性矩阵，故本构方程的确定最终又归结为弹塑性矩阵的确定。

定义塑性应力增量{dσp}与塑性应变增量{dεp}的关系为：

[De]为弹性矩阵。在应变空间εij中，假定已知塑性主应力增量的方向，并把其作为应变空间上的一个矢量，则由多重势面理论可有[114]：

式中：φk（εij）（k=1，2，3）为应变空间上其梯度矢量线性无关的势函数，φk最简单的形式为取三个应变不变量、εv、ψ，ψ为应变Lode角，为简化，不考虑ψ的影响，直接取φ1=εv，φ2=，则由式（5.3.3）得：

式（5.3.4）也可以直接由张量定律而得到[119]，写成矩阵形式为：

应力张量不变量p、q及应变张量不变量εv同以上5.2节中定义的公式，定义dpp、dqp与dp、dq具有相同公式。

则由式（5.3.4）或式（5.3.5）代入dpp、dqp的定义式中可得：

dpp、dqp分别为塑性应力增量的两个不变量，由式（5.3.2）的定义得：

若通过试验求得以下关系：

代入式（5.3.7）、式（5.3.6）则得：

考虑到关系：

则将式（5.3.10）代入式（5.3.9）然后代回式（5.3.5）得：

式中称其为应变空间的塑性矩阵。

应变空间弹塑性矩阵的展开如下：

由于

于是知：

进而可以得到：

其中

假设式（5.2.12）中B=C，则有则式（5.3.11）中表达式的中间两项可以合并，此时为对称矩阵，可展开为：

其中：

由于

又，弹性矩阵

则总应力增量为：

(www.xing528.com)

由式（5.3.11）代入上式得：

即为应变空间的弹塑性矩阵。[De]为弹性矩阵。式（5.3.16）又可写成为下式：

对于平面应变问题，从平面应变的条件dεz=0，dγyz=0，dγzx=0代入，则得平面应变状态时的关系为：

则式（5.3.13）即是可用于有限元等数值分析的应变空间弹塑性本构模型的方程。式中[]是用应变分量表示的，而通常的位移法有限元中，直接得到的是节点的位移，然后通过位移由几何方程得到单元的应变，再根据本构方程由应变求得单元应力，而当塑性矩阵[]是用应变表示时，则可以直接用应变求得塑性矩阵而无需求应力，其比应力空间中的塑性矩阵可以减少求应力的运算，提高精度，程序实现简便。

5.3.3 参数的确定

由以上的推导过程可见，应变空间的弹塑性本构方程只要确定了式（5.3.8）中的系数即可以直接用于有限元分析。另一方面，式（5.3.8）也可以通过解联立方程得另一表达形式：

只要式（5.3.17）中的A′、B′、C′、D′求得，则式（5.3.8）中的系数即可由式（5.3.17）中求解得到用A′、B′、C′、D′来表达，由于目前的试验中多为应力控制，由应力控制的试验资料来确定式（5.3.17）的系数较为方便，由以上应力空间的简化弹塑性模型的研究可知，若对A′、B′、C′、D′4个系数的关系不作规定，则要用多种应力路径试验才能确定这4个系数，若对其关系作一些假定，则试验确定将会简单一些。设塑性应变增量与应力增量的关系为：

考虑弹性应变为：

Ke、Ge分别为弹性体积模量和剪切模量，则由式（5.3.18）、式（5.3.19）可得：

称A、B、C、D为塑性系数，由以上的研究结果表明，通常的关联或非模型的塑性系数只是在多重势面理论基础上，当规定了A、B、C、D4个系数满足一定的关系时的特殊情况，为便于参数的确定和程序的方便，可对A、B、C、D4个系数的关系规定与关联模型时的关系一致，即为：

此时塑性矩阵]由式（5.3.11）可见将是对称的。作了这一假设后，则塑性系数A、B、C、D即可以像邓肯—张模型那样简单地确定，也可以直接用邓肯—张模型的Et、μt参数来表示了。如在常规三轴试验中，由于σ3=const.，dσ2=dσ3=0，dpdq=dσ1，dε2=dε3=-μtdε1，d（1+μt）dε1，dεv=（1-2μt）dε1，代入式（5.3.20）得：

考虑到Et=dσ1/dε1，则以上两式成为：

由式（5.3.21）、式（5.3.22）可解得A、B、C、D4个系数为：

E，μ为弹性部分的弹性模量和弹性泊松比。

A、B、C、D求得后代回式（5.3.20），再由式（5.3.20）可解得：

比较式（5.3.24）与式（5.3.8），则可得：

由此，应变空间弹塑性本构模型中的塑性矩阵中的系数即可以由邓肯—张模型所定义的切线模量Et、切线泊松比μt及弹性体积模量Ke和剪切模量Ge而确定，则应变空间弹塑性模型的增量本构方程即完全确定了。

5.3.4 计算实例

为检验模型及程序的可靠性，以图5.3.1所示的一个软土路堤的固结沉降为例，应用了以上的应变空间弹塑性模型和邓肯—张非线性模型进行了计算，邓肯—张模型方程及相应的软土参数为：

K=4.25MPa，n=0.308，Rf=0.897，c=3kPa，φ=16.9°，G=0.358，F=0.138，D=0.48。

对于应变空间的弹塑性模型，由邓肯—张模型公式（5.3.26）求得Et、μt，对弹性模量Ke、Ge，应通过卸荷试验来确定，本文在此假定取弹性模量Ee为：

同时，取弹性泊松比μ=0.3，则由弹性关系有：

图5.3.1 有限元网格

把Et、μt、Ke、Ge代入式（5.3.23），即求得A、B、C、D，再代入式（5.3.25）即得到应变空间弹塑性模型所需的参数

分别用以上的参数及邓肯—张模型和本文的应变空间弹塑性模型，用有限元计算了图5.3.1所示的软土路堤的填筑过程的固结沉降情况，采用渗透系数为K=10-7cm/s。两个模型计算所得的原地面沉降情况如图5.3.2所示，路堤中点原地面处的固结沉降过程如图5.3.3所示，由图5.3.3可见，最大沉降点的沉降过程及大小基本相同，由图5.3.2可见，路面沉降也基本一致，但路堤脚两侧的隆起数值，则弹塑性模型比邓肯—张非线性弹性模型要大一些，这主要是弹塑性模型的弹塑性矩阵是满阵的，其对侧向变形或剪缩剪胀等特性的反映要比非线性弹性模型较强，故隆起量稍大一些是正常的。

图5.3.2 原地面路堤中线点两种模型计算的沉降过程

图5.3.3 路堤下原地面两种模型计算的最终沉降

由于路堤填筑是一个加载过程，应力路径也较简单，两种模型计算结果较一致是合理的，这里的简单情况下的计算比较的主要目的只是检验应变空间弹塑性模型的可靠性和程序的正确性，通过以上的计算比较，可以认为以上所得出的应变空间弹塑性本构模型及计算程序是可靠的和正确的。

5.3.5 小结

本文在多重势面理论或张量定律基础上建立的应变空间弹塑性本构模型，具有参数确定容易，无需推求塑性势函数的优点，在程序实现上，与应力空间模型比较，具有减少运算，提高精度和简便的特点，更易于应用。