土是一种具有复杂力学特性的历史作用产物，力学性能上表现出强烈的非线性和弹塑性。姚仰平等^[1-3]基于超固结土的一些基本特性提出了适用于超固结土的UH模型，该模型不仅可以合理反映超固结土的剪缩、剪胀、硬化、软化及应力路径相关特性，而且模型的参数少，物理意义明确，可直接通过常规三轴实验测得，便于实际工程的应用。

甄文战等^[4-5]基于UH模型理论及变换应力方法，采用FORTRAN语言首先编写了UH模型的用户材料子程序 (UMAT)，实现了UH模型与有限元软件的结合；冯兴等^[6-7]在UH模型UMAT的基础上，加入罗汀等^[8]将熊文林^[9]提出的非关联塑性切线刚度矩阵对称化方法应用于UH模型内容，并采用与广义SMP准则相似的处理方法编写出了考虑黏聚力c的UH模型UMAT。

原有的UH模型UMAT都是基于不发生破坏的受压应力状态编写而成，事实上，在有限元计算实际工程问题时，很多时候不可避免地会出现局部单元拉裂、剪坏的情况，这种情况的出现不仅会造成计算结果本身的不合理，也有可能会导致有限元计算时出现中途报错退出计算的情况，从而降低了采用UH模型UMAT进行有限元计算的稳定性和实用性。

本文针对采用UH模型UMAT进行有限元计算时出现的单元破坏问题 (包括拉裂和剪坏)，参考平面问题单元破坏修正方法，并结合不同坐标系下应力之间的变换关系，推导得到适用于UH模型的三维问题单元破坏修正公式，对采用UH模型UMAT进行有限元计算时出现的破坏应力状态进行相应的应力修正，以使得采用UH模型UMAT的有限元计算过程更加稳定，计算结果更为合理。

1.   UH模型及其有限元实现

1.1   UH模型

姚仰平等^[10]提出的UH模型是基于当前屈服面 (当前应力点所在的屈服面) 和参考屈服面 (参考应力点所在的剑桥模型屈服面) 的动态变化关系以及与伏斯列夫面变化有关的潜在强度的不同，来描述超固结土的硬化、软化、临界状态、剪缩、剪胀特性和不同应力路径对超固结土变形特性的影响。

如图 1所示，UH模型的当前屈服面与参考屈服面均为椭圆屈服面。

图  1  UH模型的当前屈服面与参考屈服面

Figure  1.  Current yield surface and referential yield surface of UH model

下载: 全尺寸图片 幻灯片

当前屈服面表达式为

参考屈服面表达式为

式中：，e₀为初始孔隙比, λ和κ分别为e-ln p空间中等向压缩线的斜率和等向回弹线的斜率，e为孔隙比，p为平均正应力；M为特征状态或临界状态应力比；p为参考应力点的平均正应力；η为应力比；p_x0和p_x0分别为初始当前屈服面和初始参考屈服面与p轴的交点；ε_v^p为塑性体积应变；H为UH模型的硬化参量，其表达式为

其中：

1.2   变换应力方法

UH模型UMAT中采用姚仰平等^[11-12]提出的变换应力方法进行三维化 (见图 2)，变换应力与真实应力之间的映射关系为

图  2  变换应力空间

Figure  2.  Transformed stress space

下载: 全尺寸图片 幻灯片

式中：分别为变换应力空间的平均正应力和广义剪应力；q_c为三轴压缩应力状态的破坏剪应力；和θ_σ分别为变换应力空间和真实应力空间的应力洛德角。

根据式 (7) 可推得变换应力与真实应力之间的变换关系为

式中：分别为变换应力张量和真实应力张量；q为广义剪应力；δ_ij为克罗内克符号。

1.3   UH模型的UMAT

UH模型的UMAT^[13-14]是在UH模型理论基础上，引入了变换应力三维化方法，并加入黏聚力c和弹塑性刚度矩阵对称化内容，采用FORTRAN语言编写而成的，UH模型UMAT采用半隐式回映应力更新算法，具体的计算流程如图 3所示。

图  3  UH模型UMAT计算流程图

Figure  3.  Computation flowchart of UH model's UMAT

下载: 全尺寸图片 幻灯片

2.   单元拉裂修正

在使用UH模型UMAT进行有限元计算时，有时局部单元的积分点上会出现拉裂应力 (=σ₃+ccot φ < 0，φ为内摩擦角)，这种拉裂应力的存在是不合理的：一方面，无黏性土是没有抗拉强度的，黏性土所受的拉应力也不应大于黏聚力c所提供的抗拉强度，而且拉裂应力的存在也可能导致计算得到的平均正应力出现负值、应力比η>3的情况出现，而这些对于土而言均是不合理的；另一方面，拉裂情况的出现会导致用UH模型UMAT有限元计算时出现中途报错退出的情况，这是由于UH模型UMAT是基于受压应力状态编写的，而实际有限元计算时往往不可避免地会有某些积分点出现拉裂的情况，这种拉裂情况的存在就有可能使得在调用UH模型UMAT进行有限元计算时出现某些数学错误，例如UH模型UMAT中采用的变换应力方法基于SMP准则，在计算三轴压缩状态破坏剪应力(为由计算得到的应力不变量) 时，如果有拉裂应力的存在，就有可能会导致，使得计算报错中止，从而降低了计算的稳定性。

针对拉裂应力引起的两方面问题，需要对出现的不合理拉裂应力进行修正，本文将在已有的平面问题单元拉裂修正方法^[15]基础上，提出适用于UH模型的三维问题单元拉裂修正方法。

2.1   平面问题的单元拉裂修正

对于平面问题单元拉裂修正，可以通过结合Mohr应力圆，并基于一定的假设条件推导出相应的拉裂修正计算公式。如图 4所示，若假定该积分点的最大主应力σ₁保持不变，将σ₃修正为0，并假定主应力方向不发生变化，即可得到平面问题单元拉裂修正公式为

图  4  平面问题的单元拉裂修正

Figure  4.  Element tension failure correction of plane problem

下载: 全尺寸图片 幻灯片

式中：τ′_xz的符号与τ_xz一致。

2.2   适用于UH模型的三维问题单元拉裂修正

虽然平面问题的单元拉裂修正可以解决平面问题的拉裂情况，但是对于不能简化为平面问题的三维问题则并不适用；而UH模型UMAT本身是基于三维单元编写的，因此，平面问题的拉裂修正并不适合直接在UH模型UMAT中的应用，所以需要对平面问题的拉裂修正进行扩展而得到适用于更一般的三维问题的拉裂修正方法。本文所述的三维问题的拉裂修正方法需要用到不同坐标系下的应力变换关系^[16]。

2.2.1   主应力与应力分量之间的变换关系

设新旧2个笛卡儿直角坐标系分别为Ox′y′z′和Oxyz，则2个坐标系下的应力分量之间的变换关系为

式中：σ′_ij为新坐标系Ox′y′z′下的应力张量；σ_kl为旧坐标系Oxyz下的应力张量；l_ik和l_jl为坐标变换张量。

设新旧坐标系Ox′y′z′和Oxyz坐标轴之间的方向余弦如表 1所示。

表  1  新旧坐标轴之间的方向余弦

Table  1.  Direction cosine between new and old coordinate

坐标轴	x	y	z
x′	l1	m1	n1
y′	l2	m2	n2
z′	l3	m3	n3

下载: 导出CSV 

| 显示表格

将式 (10) 展开即可得

假定旧坐标系Oxyz的x、y、z轴对应于σ₁、σ₂、σ₃方向，新坐标系Ox′y′z′的x′、y′、z′轴与σ_x、σ_y、σ_z方向重合，则有

式中：l₁、l₂和l₃为σ₁在Ox′y′z′坐标系中的方向余弦；m₁、m₂和m₃为σ₂在Ox′y′z′坐标系中的方向余弦；n₁、n₂和n₃为σ₃在Ox′y′z′坐标系中的方向余弦。

2.2.2   拉裂修正原理

采用UH模型UMAT进行有限元计算时，若出现拉裂情况，也即，可以参考平面问题单元拉裂修正的假定条件。但是为了得到修正之后的中主应力σ₂值，则需要对假定条件进行相应的补充，本文为了简化修正公式，以及保证修正前后不改变3个主应力之间的大小次序，假定修正前后中主应力参数b不变。具体假定条件如下：

假定1 拉裂之后，认为土不可以再承受拉应力，也即取σ₃=0。

假定2 中主应力参数不变，也即应力洛德角θ_σ保持不变。

假定3 最大主应力σ₁不变。

假定4 主应力方向不变。

事实上，假定3可以有2种情况：竖向应力σ_z不变和最大主应力σ₁不变。竖向应力σ_z保持不变的假定对于以竖向载荷为主的情况更为合理，但不适用于大正应力非竖向的情况，例如出现被动土压力情况；而假定最大主应力σ₁不变的拉裂修正则不会存在这个问题，其具有更大的适用性；而且对于以竖向载荷为主的情况，假定最大主应力σ₁不变的拉裂修正和假定竖向应力σ_z不变的拉裂修正得到结果是很接近的；假定σ₁不变的拉裂修正公式相对而言也更为简单，所以本文所述的拉裂修正采用的是基于最大主应力σ₁不变的假定。通过假定1~假定3即可推得修正后的主应力与修正前的主应力之间的关系为

得到拉裂修正之后的主应力，再根据假定4和式 (13) 即可得到修正改之后的六维应力分量为

式中：为拉裂修正之后的3个主应力分量；。

3.   单元剪坏修正

经过第2节所述的单元拉裂修正虽然消除了不合理的拉裂应力状态，但该应力还有可能超出抗剪强度包线，这种应力状态理论上也是不合理的，因此需要再对其进行剪坏修正。

3.1   平面问题的单元剪坏修正

与平面问题的单元拉裂修正方法相似，可以结合一定的假设条件和Mohr应力圆，推导出基于Mohr-Coulomb强度理论的平面问题单元剪坏修正公式^[15]。如图 5所示，若假定该积分点的最大主应力σ₁不变，修正之后Mohr圆与强度破坏包线相切，且假定修正前后的主应力方向不发生改变，则可以推得平面问题单元剪坏修正的计算公式为

图  5  平面问题的单元剪坏修正

Figure  5.  Element shear failure correction of plane problem

下载: 全尺寸图片 幻灯片

式中：τ′_xz的符号与τ_xz一致。

3.2   适用于UH模型的三维问题单元剪坏修正

对于三维问题的单元剪坏修正，可以参考平面问题单元剪坏修正的假设条件和三维问题单元拉裂修正的处理方法得到。在采用UH模型UMAT进行有限元计算时，若出现(为由计算得到的应力比和潜在强度) 的情况，则采取如下假定：

假定1 。

假定2 中主应力参数不变，也即应力洛德角θ_σ不变。

假定3 最大主应力σ₁不变。

假定4 主应力方向不变。

事实上，假定3可以有3种情况：平均正应力p不变、竖向应力σ_z不变和最大主应力σ₁不变。若假定修正前后平均正应力p不变，则剪坏修正会使得大正应力值变化很大，而一般来说，外载荷主要由与之方向相一致的大主应力承担，因此这种假定的剪坏修正对于计算收敛性会有较大的影响；假定竖向应力σ_z不变的剪坏修正对于以竖向载荷为主的情况来说是更为合理的，但实际问题的大正应力方向不一定都是竖向的，对于大正应力非竖向则该假定也就不再合适了；最大主应力σ₁不变的假定对大正应力值改动很小，所以可以降低修正对计算收敛性的影响，而且该假定不仅适用于大正应力竖向的情况，也适用于大正应力非竖向的情况，具有更大的适用性，所以本文所述的剪坏修正是基于最大主应力σ₁不变假定的。通过假定1~假定3，有

式中：

联合式 (17)、式 (18)，可推得更新后的主应力分量为

式中：

得到修正之后的主应力分量，再根据式 (15) 即可推导得到经过剪坏修正之后的六维应力分量。

4.   破坏修正的有限元验证

UH模型UMAT的单元验证可以参见参考文献[6]，本文不再赘述。本文将采用一个有限元计算基坑开挖的例子来说明拉裂修正和剪坏修正的有效性。UH模型UMAT是基于三维单元编写的，这里通过对边界条件的特殊设置来实现用三维单元描述平面问题。基坑的各个部分几何尺寸如图 6所示。土材料的本构模型选用UH模型，模型参数取值如表 2所示。表中：ν为泊松比；N为e-lnp空间中正常压缩线的截距 (即当p=1时的孔隙比)。支护桩采用线弹性模型，弹性模量E=30 GPa，泊松比ν=0.1。土体与支护桩之间采用简单的绑定约束。图 7为计算模型的网格划分，土单元与墙单元均采用8节点六面体实体单元，单元总数为1 380，开挖区土单元数为150，墙单元数为20。

图  6  计算模型的几何尺寸

Figure  6.  Geometric dimension of computation model

下载: 全尺寸图片 幻灯片

图  7  计算模型网格划分

Figure  7.  Computation model meshing

下载: 全尺寸图片 幻灯片

表  2  UH模型参数取值

Table  2.  Parameter value of UH model

参数	M	ν	κ	λ	e0	N	c/kPa
数值	0.984	0.3	0.016	0.09	0.6	1.15	10

下载: 导出CSV 

| 显示表格

图 8为有无破坏修正的土体等值云图。

图  8  有无破坏修正的土体等值云图

Figure  8.  contour of soils with and without failure correction

下载: 全尺寸图片 幻灯片

表 3为图 8(a)中1 244号单元第1积分点和18号单元第4积分点破坏修正前后的主应力分量、应力分量、应力比、潜在强度在某个增量步的计算值。由破坏修正前的结果可以看出，该积分点是存在拉裂情况的，该应力状态也会导致，如果不进行处理，计算就会出现中途中断的情况，而且由该应力状态计算得到的应力比，也不满足的要求；进行破坏修正后，最大主应力σ₁的值没有改变，最小主应力σ₃被修正为0(=21.44 kPa)，应力比，破坏修正之后的值也均大于0。所以，破坏修正既消除了不合理的应力状态，也消除了由于破坏应力导致的子程序计算中途报错退出的问题。

表  3  破坏修正前后计算结果对比

Table  3.  Comparison of calculation results before and after failure correction

变量	1 244号单元第1积分点			18号单元第4积分点
	修正前	修正后		修正前	修正后
	32.53	32.53		42.11	42.11
	7.39	23.97		5.07	23.93
	-14.16	21.44		-7.47	21.44
σx/kPa	-11.16	3.51		-15.61	2.86
σy/kPa	-35.600	0.007		-28.910	0
σz/kPa	8.19	10.10		19.92	20.30
τxz/kPa	-8.03	-2.74		-5.24	-2.57
	4.714	0.387		3.373	0.670
	1.908	1.908		1.801	1.801
	-0.223	29.290		-0.444	10.830

下载: 导出CSV 

| 显示表格

由第2节和第3节的描述可知，拉裂和剪坏的判定均是基于考虑了黏聚力之后的应力，所以采用基于计算得到的变量的等值云图，而不是基于真实应力计算得到的变量的等值云图。

图 8(a)和图 8(b)分别为基坑开挖至3.75 m时没有破坏修正 (为了避免由于破坏应力而导致的中途退出计算情况，需对由于破坏应力导致的数学错误做相应的处理) 和有破坏修正计算得到的等值云图对比 (应力以压为正、拉为负)。由图 8(a)可以看出，如果没有破坏修正，基坑内部最表层的土体确实存在拉裂应力，这对于土而言是不合理的，需要对其进行修正；而由图 8(b)可以看出，有破坏修正的有限元计算提高了σ₃的值，使得基坑内部表层不合理的拉裂应力状态得以消除。由图 8(a)和图 8(b)对比还可以发现，破坏修正使得σ₃的值整体提高较大，没有破坏修正的约为1/6~1/5，有破坏修正的约为2/5，这主要是由于剪坏修正引起的，参考图 5，如果没有剪坏修正，应力比的值没有限制，可以达到很大，也即相当于提高了土的抗剪强度，根据土压力理论可知，土的侧向应力也就相应的会变小，而剪坏修正使得土的应力比不能过大，相当于限制住了土的抗剪强度，所以侧向应力会相应提高。

图 9和图 10分别为有无破坏修正的应力比和潜在强度等值云图。图 9(a)和图 10(a)分别为没有破坏修正的应力比和潜在强度等值云图，由两者的云图对比可以看出，计算结果中存在应力比大于潜在强度的情况，也即存在剪坏的应力状态。图 9(b)和图 10(b)分别为有破坏修正的应力比和潜在强度等值云图，可以看出，经过修正之后，消除了计算结果中的不合理剪坏应力状态，使得土单元上的应力比均符合要求。

图  9  有无破坏修正的土体的等值云图

Figure  9.  contour of soils with and without failure correction

下载: 全尺寸图片 幻灯片

图  10  有无破坏修正的土体的等值云图

Figure  10.  contour of soils with and without failure correction

下载: 全尺寸图片 幻灯片

图 11为有无破坏修正和有破坏修正的等值云图 (应力以压为正、拉为负)。可以看出，破坏修正对于的分布大小影响很小，这是因为破坏修正 (包括拉裂修正和剪坏修正) 都是基于最大主应力σ₁不变假定进行的，而本例中最大主应力σ₁与竖向应力σ_z的大小方向又是很接近的，在有限元计算过程中，竖向的载荷主要由σ_z承担，假定最大主应力σ₁不变，也即意味着σ_z基本不变，因此根据式 (24) 可知，由于破坏修正而产生的竖向不平衡力也是很小的，所以破坏修正对竖向应力σ_z的大小分布影响很小，也即对最大主应力σ₁的影响也很小。

图  11  有无破坏修正的土体的等值云图

Figure  11.  contour of soils with and without failure correction

下载: 全尺寸图片 幻灯片

图 12为有无破坏修正与有破坏修正的竖向位移等值云图。可以看出，两者的位移分布形式是一致的，但破坏修正使得整体变形值变小，这是因为UH模型中，弹塑性刚度矩阵是与应力状态相关的，一般来说，平均正应力p越大，刚度也越大，破坏修正会提高平均正应力p的值 (参见图 4和图 5)，也就提高了破坏单元的刚度矩阵，而其他单元的刚度矩阵是基本不变的，因此破坏修正会提高整体的刚度矩阵值，从而导致整体变形值的下降。

图  12  有无破坏修正的土体的竖向位移等值云图

Figure  12.  Vertical displacement contour of soils with and without failure correction

下载: 全尺寸图片 幻灯片

对于采用有破坏修正UH模型UMAT的有限元计算，一个迭代步下实际的计算流程为：先调用UH模型UMAT进行应力更新得到应力σ；再调用破坏修正子程序进行破坏修正得到应力σ′；由破坏修正的应力σ′可以计算得到不平衡力ΔR′，再用不平衡力ΔR′进行有限元的平衡迭代。不平衡力ΔR′具体表达式为

式 (21) 又可以化为

式中：

其中：ΔR为没有破坏修正的不平衡力；Δr为破坏修正引起的附加不平衡力。因此，如果破坏修正准则较为合理，也有可能使得破坏修正引起的附加不平衡力Δr降低总的不平衡力，提高迭代计算的收敛性。如本文所述的计算基坑开挖的例子，有破坏修正的计算收敛性就好于没有破坏修正的计算收敛性，在计算第1、2、3层开挖时，没有破坏修正分别用了4、7、92个增量步，而有破坏修正则分别用了1、7、8个增量步，这也说明了本文所采用的破坏修正假定是较为合理的。

综上所述，采用有破坏修正的UH模型UMAT的有限元计算不仅避免了由于破坏应力状态引起的计算中断问题，也消除了计算结果中的不合理破坏应力状态，使得计算过程更加稳定，计算结果更为合理。

5.   结论

在采用UH模型UMAT进行有限元计算时，经常会遇到局部单元出现破坏应力的情况 (主要包括拉裂和剪坏)，这种破坏应力的存在是不合理的：一方面，土单元所受的应力不应该超出理论上的强度包线；另一方面，破坏应力的存在会导致有限元计算上出现相应的稳定性问题。

1) 针对在有限元计算过程中出现的局部单元破坏问题，本文介绍了一种适用于UH模型的三维问题单元破坏修正方法，该方法不仅可以消除有限元计算中出现的不合理的破坏应力状态，以及由于这种应力状态而导致的计算稳定性问题，而且该方法并没有改变原有UH模型UMAT的算法流程，不会由于这种修正而导致原有UH模型的UMAT产生新的理论上的问题。

2) 本文所述的单元破坏修正方法对于有限元计算收敛性是存在影响的，由于修正过程强制修改了应力分量，会产生了新的不平衡力，理论上会增加有限元迭代的步数。但一方面这种修正一般只是对局部少数单元的修正，因此产生的不平衡力一般不会很大；另一方面，修正之后的应力如果变得更为合理，也有可能降低了迭代的不平衡力，降低有限元计算整体的迭代步数，提高计算的收敛性。