0.   引言

实际工程建设中，包括准静态条件下的隧道与基坑开挖、堤坝修筑等，以及动态条件下的地震、潮汐与交通等荷载作用都会导致工程中土体发生主应力轴旋转。以单剪^[1-2]和空心圆柱环剪试验^[3-5]为代表的模拟主应力轴旋转作用的土单元试验发现土体变形存在显著的非共轴变形特性，即材料当前的主应力方向与塑性主应变率方向不一致。土体在主应力轴旋转作用下的非共轴变形规律有悖于传统塑性流动理论中包含的塑性应变率与应力共轴的假定，这给基于共轴假定的弹塑性模型进行的岩土工程的变形分析的合理性提出了质疑。对于变形控制要求高的工程问题分析，评估土体非共轴变形特性对结果的影响对于提高工程分析准确性具有重要的意义。

为了能够合理描述非共轴行为，基于不同理论框架的非共轴本构模型被陆续提出。大多数非共轴模型都基于一个共同的假设，即假定塑性变形可以划分为共轴项和非共轴项，其中非共轴项用于定义非共轴流动法则。构建非共轴流动法则主要有3类方法，包括Rudnicki^[6]提出的“角点”理论，Yu等^[7]采用Double shearing理论以及Lashikari等^[8]根据应力空间X-Y ($X = {{({\sigma _{11}} - {\sigma _{33}})} \mathord{\left/ {\vphantom {{({\sigma _{11}} - {\sigma _{33}})} 2}} \right. } 2}$，$Y = {\sigma _{13}}$)内应力圆切向直接指定的方法。Rudnicki提出的方法由于其表述形式的一般性，在当前的相关研究中得到广泛应用。钱建固等^[9]在考虑了第三应力不变量的影响后将该方法扩展到了更一般的三维应力状态。Li等^[10]和Tsutsumi等^[11]将这类非共轴流动法则与各向异性本构模型框架相结合，较好地模拟了空心圆柱纯主应轴旋转试验中的非共轴现象。陈洲泉等^[12-13]和Chen等^[14]进一步改进该方法，解决了模型不能模拟单调定轴剪切试验中非共轴现象的缺陷。

从本构模拟的角度研究非共轴变形，学者们已经开展了丰富而深入的研究。然而将非共轴本构模型应用于有限元分析则相对较少。Yang等^[15-17]较早针对角点型非共轴模型开展数值实现方面的研究。他们比较了共轴与非共轴模型在分析了浅基础地基承载力和变形问题上差异，指出不考虑非共轴塑性的计算结果偏于不安全。另外，还分析了非共轴塑性对仓筒问题中主应力轴旋转和土压力分布所产生的影响。Chang等^[18]将“角点”非共轴理论与Cosserat塑性理论结合分析了应变局部化问题。袁冉等^[19]和Yuan等^[20]基于Yu等^[7]提出的非共轴理论构建的各向异性非共轴模型分别开展了浅基础地基承载和变形的分析，以及隧道开挖诱发地表沉降问题的分析。这些研究的一个共同特点就是都采用基于误差控制的自动划分子增量步的显式积分算法将非共轴模型数值实现。就如Chang等^[18]指出的，由于非共轴流动法则构造上的复杂性给构造非共轴模型的隐式积分算法造成了较大困难。尽管该算法具有自动满足屈服面一致性条件，以及平衡迭代二次收敛速度等优势，但当前还未见详实而确凿的相关研究。

考虑到构造非共轴本构模型全隐式积分算法^[21]的复杂性，本文将采用半隐式算法^[22-23]针对角点型非共轴模型开展数值实现工作，这也是当前相关研究中所鲜见的。值得指出的是，半隐式算法避免了对流动法则的偏导，极大简化了推导算法迭代格式的难度。同时，该算法同样自动满足屈服面一致性条件，收敛速度介于隐式与显式之间，特别适合于塑性流动法则形式复杂的本构模型的数值实现。再基于ABAQUS有限元分析软件，对比分析共轴模型和非共轴模型在数值模拟单剪问题与经典的活动门(Trapdoor)问题上的差异，验证非共轴模型与算法的有效性。

1.   角点型非共轴本构模型

当前，建立非共轴本构模型的主流思路就是通过向共轴弹塑性模型中增添非共轴流动法则来实现。其最为基础的理论假设就是将塑性应变增量${\text{d}}\varepsilon _{ij}^{\text{p}}$分解为共轴项${\text{d}}\varepsilon _{ij}^{{\text{pc}}}$和非共轴项${\text{d}}\varepsilon _{ij}^{{\text{pn}}}$，形式如下：

式中：共轴塑性应变增量${\text{d}}\varepsilon _{ij}^{{\text{pc}}}$的演化规律完全与共轴塑性模型相关，而非共轴塑性应变增量${\text{d}}\varepsilon _{ij}^{{\text{pn}}}$则需要建立另外的非共轴流动法则。

本文采用的共轴模型为考虑剪切硬化法则的线性Drucker-Prager模型，其屈服函数的形式为

式中：$p = {{{\sigma _{kk}}} \mathord{\left/ {\vphantom {{{\sigma _{kk}}} 3}} \right. } 3}$为平均主应力；$q = \sqrt {3{J_2}}$为广义剪应力；而${J_2} = {{{s_{ij}}{s_{ij}}} \mathord{\left/ {\vphantom {{{s_{ij}}{s_{ij}}} 2}} \right. } 2}$为第二应力不变量；${s_{ij}} = {\sigma _{ij}} - p{\delta _{ij}}$为应力偏张量；$\eta$为硬化参数，反映土体三轴压缩条件应力应变曲线中的应力比；c为土体的黏聚力；${\phi _{\text{c}}}$为土体破坏时的临界摩擦角。如图 0所示，该模型在p-q应力空间为过固定点($- c\cot {\phi _{\text{c}}}$，0)的斜直线，随着应力比变化而绕固定点旋转。

图  0  p-q平面内屈服面和破坏面

Figure  0.  Yield and failure surfaces in p-q plane

下载: 全尺寸图片 幻灯片

模型采用非关联流动法则，塑性势函数为

式中：$\psi$为剪胀角。那么共轴塑性应变增量可表示为

式中：${r_{ij}} = {{\partial Q} \mathord{\left/ {\vphantom {{\partial Q} {\partial {\sigma _{ij}}}}} \right. } {\partial {\sigma _{ij}}}}$为共轴塑性流动方向；${\text{d}}\lambda$为塑性标量因子。

表征应力比$\eta$变化的硬化规律采用双曲线形式

式中：${M_{\text{c}}} = {{6\sin {\phi _{\text{c}}}} \mathord{\left/ {\vphantom {{6\sin {\phi _{\text{c}}}} {(3 - \sin {\phi _{\text{c}}})}}} \right. } {(3 - \sin {\phi _{\text{c}}})}}$为临界应力比，$\varepsilon _{\text{s}}^{\text{p}}$为等效塑性应变($\varepsilon _{\text{s}}^{\text{p}} = \sqrt {{{2e_{ij}^{\text{p}}e_{ij}^{\text{p}}} \mathord{\left/ {\vphantom {{2e_{ij}^{\text{p}}e_{ij}^{\text{p}}} 3}} \right. } 3}}$，$e_{ij}^{\text{p}} = \varepsilon _{ij}^{\text{p}} - {{\varepsilon _{kk}^{\text{p}}{\delta _{ij}}} \mathord{\left/ {\vphantom {{\varepsilon _{kk}^{\text{p}}{\delta _{ij}}} 3}} \right. } 3}$为塑性应变偏张量)；${h_{\text{c}}}$是与硬化曲线相关的材料参数，通过拟合土体的三轴压缩试验曲线来获得。

根据Rudnicki和Rice的理论，非共轴塑性变形被假定为由应力增量在当前应力张量切向方向的分量引起，根据Gram-Schmidt张量正交化理论，该切向分量定义为

然后，非共轴塑性应变增量定义为

或者，

式中：$\boldsymbol{C}_{ijkl}^{} = {\delta _{ik}}{\delta _{jl}} - \frac{{{\delta _{ij}}{\delta _{kl}}}}{3} - \frac{{{s_{ij}}{s_{kl}}}}{{{s_{mn}}{s_{mn}}}}$为非共轴转换张量；${h_n}$为非共轴模量，表征应力增量引发非共轴塑性变形的难易程度。该参数可以在共轴模型参数确定后，通过拟合单剪试验应力应变旋转角变化曲线来获得。

2.   本构模型数值积分算法

本文采用Moran等^[24]提出的半隐式积分算法。这种算法在应力更新过程中，塑性标量因子作为未知量，采用隐式表达；而塑性流动方向已知，为显式表达。这意味着在一个增量步中，塑性因子要在增量步结束时求得，而塑性流动方向则是在增量步开始时就计算好的。同时，为了防止显式算法存在屈服面偏移问题，该算法限定计算结束时应力状态满足屈服方程。这种算法构建方式有效回避了隐式算法推导塑性流动函数偏导存在的困难，特别是对数学形式复杂的本构模型，同时也能克服显式算法屈服面偏移问题，是一种具有广泛适用性和准确性的本构模型数值积分算法。

在ABAQUS/Explicit计算平台下，应力更新过程为根据在${t_n}$~${t_{n + 1}}$时间步中给定的应变增量${\text{d}}{\varepsilon _{ij}}$来计算应力增量${\text{d}}{\sigma _{ij}}$。按照弹塑性理论，将应变增量分解为弹性部分${\text{d}}\varepsilon _{ij}^{\text{e}}$和塑性部分${\text{d}}\varepsilon _{ij}^{\text{p}}$：

那么，增量步${t_n}$到${t_{n + 1}}$的应力更新格式即为

式中：$\boldsymbol{D}_{ijkl}^{\text{e}} = G({\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}}) + (K - {{2G} \mathord{\left/ {\vphantom {{2G} 3}} \right. } 3}){\delta _{ij}}{\delta _{kl}}$为弹性刚度矩阵；$\boldsymbol{C}_{ijkl}^n$代表矩阵由${t_n}$的应力状态定义。

考虑到${\text{d}}{\sigma _{kl}} = {}^{n + 1}{\sigma _{ij}} - {}^n{\sigma _{ij}}$，$\boldsymbol{C}_{ijkl}^n{}^n{\sigma _{kl}} = 0$（公式相关推导参考附录1），则式（9）可转换为

式中：${}^n\sigma _{ij}^{{\text{trial}}} = {}^n{\sigma _{ij}} + \boldsymbol{D}_{ijkl}^{\text{e}}{\text{d}}{\varepsilon _{kl}}$为弹性预测应力。再结合硬化规律和屈服准则，完整的应力更新过程可由下列方程组表示：

值得指出的是方程组（12）包含了11个式子，式中与${t_{n + 1}}$对应的未知量构成方程组解向量X，表示为

指定方程组（12）对应的残差向量为R(X)，在采用Newton-Raphson迭代法求解该非线性方程组时，该残差向量的Jacobi矩阵可表示为$\partial R/\partial X$，是一个11×11的方阵，其形式为

那么，该模型的半隐式应力积分算法的迭代求解过程如下：

（1）初始化：k=0，赋值初始解向量${X^0}$

（2）判断是否满足屈服条件：

如果$f({}^{n + 1}\sigma _{ij}^{(k)}, {}^{n + 1}{\eta ^{(k)}}) \leqslant 0$，那么不进行牛顿迭代，更新应力，${ }^{n+1} \sigma_{i j}={ }^n \sigma_{i j}$，结束计算。

否则进入下步。

（3）判断残差向量的范数$\left\| {R({X^k})} \right\|$是否满足收敛条件：

如果$\left\| {R({X^k})} \right\|$ > 收敛容差，进行牛顿迭代：

设置$k \leftarrow k + 1$，返回第3步。

否则进入下步。

（4）更新应力及相关变量：

迭代结束。

3.   数值算例

3.1   单剪试验模拟

将在ABAQUS/Explicit分析模块中，基于用户材料子程序开发接口VUMAT编写文中介绍的非共轴模型的应力积分算法，并通过模拟单剪试验的过程来验证算法的有效性及模型描述非共轴变形的能力。单剪有限元分析模型采用单个减缩积分平面应变单元CPE4R。本构模型参数选取参考Yang等^[15]等工作，如表 1所示。试验加载方式如图 1所示，竖向压力保持不变，${\sigma _{22}}$=500 kPa，初始静止侧压力系数${K_0} = 0.5$，在单元上边缘施加水平位移实现剪切过程。

表  1  模型参数

Table  1.  Model parameters

弹性参数	塑性参数
	共轴参数		非共轴参数
G=16 MPa	${\phi _{\text{c}}}$=30°	$\psi$=0°	hn=2G, 1G,
v=0.25	c =5 kPa	${h_{\text{c}}}$=0.001	0.5G, 0.2G

下载: 导出CSV 

| 显示表格

图  1  单剪试验和张量主轴旋转

Figure  1.  Simple shear tests and rotation of principal axis of tensors

下载: 全尺寸图片 幻灯片

图 1（b）说明了平面应变条件下力学张量主轴方向构成的局部坐标系(${x'_1}$, ${x'_2}$)与固定不变的整体坐标(${x_1}$, ${x_2}$)的夹角关系。对于本问题中考虑的整体坐标系下力学张量，主要为应力张量${\sigma _{ij}}$和应变张量${\varepsilon _{ij}}$，求解对应的主分量可以表述为将与力学张量相关的实对称矩阵转换为对角阵的过程，其数学关系为

式中：${\boldsymbol{T}_{kl}}$可以代表为应力张量或者应变张量，其矩阵形式为

${\boldsymbol{\bar T}_{ij}}$为${\boldsymbol{T}_{kl}}$对应的主分量，矩阵形式为

${\boldsymbol{Q}_{ik}}$为正交矩阵，表示坐标转换过程，矩阵形式为

$\boldsymbol{Q}_{lj}^{\text{T}}$为${\boldsymbol{Q}_{ik}}$的转置。值得指出的是下文中主应力方向角取为$\alpha$，主应变率方向角取为$\beta$。并且规定逆时针旋转整体坐标系得到主轴坐标系时，主轴方向角取值为正，顺时针为负，旋转角取值范围为[-180°，180°]。

图 2比较了共轴模型与取不同非共轴模量的非共轴模型计算的剪应力与剪应变的关系。从图 2可以看出非共轴模量越小，剪切应力增长越慢，说明非共轴塑性流动机制引起了更大的塑性剪切应变。随着非共轴模量的增大，应力应变曲线越来越接近共轴模型的计算结果。

图  2  应力应变关系

Figure  2.  Relationship between shear stress and shear strain

下载: 全尺寸图片 幻灯片

图 3比较了共轴模型与非共轴模型分析得到的主应力旋转角和主应变率旋转角在剪切过程中的变化。图 4则进一步总结了图 3中5种情况中主应变率方向角与主应力方向角之差（即非共轴角$\beta - \alpha$）的变化情况。这些结果都表明随着非共轴模量的增大，非共轴效应越小，非共轴模型计算的结果与共轴模型的计算结果越接近一致。这些分析结果有效验证了模型和算法的合理性和可靠性。

图  3  主应力偏角和主应变偏角

Figure  3.  Rotation angles of principal stress and principal strain rate

下载: 全尺寸图片 幻灯片

图  4  非共轴角的变化

Figure  4.  Variation of the non-coaxial angle

下载: 全尺寸图片 幻灯片

3.2   活动门问题有限元分析

活动门问题是与隧道开挖等地下工程密切相关的力学分析模型。采用有限元方法分析该问题能够为理解隧道工程引起的应力和位移变化规律提供了更为深刻的认识。本文在平面应变条件下建立该问题的1/2平面对称有限元分析模型，其有限元网格如图 5所示，单元数目总共为6000个，单元类型为CPE4R。模型中活动门宽度B=1 m，上覆土体长度为10B，深度为6B。材料模型的参数沿用表 1中的取值。并在本节中根据上一节获得的结论，就非共轴变形特性对活动门问题计算结果的影响开展参数分析。

图  5  活动门问题计算模型

Figure  5.  FEM model for trapdoor problem

下载: 全尺寸图片 幻灯片

分析过程划分为两步：第一步进行地应力平衡。首先，约束计算模型左右竖向边界的水平位移。在土体的底部布置两块刚性挡板，通过建立挡板与土体的接触条件并挡板施加固定约束来限制土体的竖向位移。整个土体施加重力，重力加速度为g=9.8 m/s²。然后模型内部定义沿竖直方向线性分布的地应力，该分布力可表示为${\sigma _{\text{v}}} = \rho gH$，其中H表示土体内某位置距离模型上表面的距离, $\rho$=1.63 t/m³。水平地应力定义为${\sigma _{\text{h}}} = {K_0}{\sigma _{\text{v}}}$，式中${K_0}$为静止土压力系数，取为0.5。地应力平衡后，在第二步进行位移加载，解除模型下部活动门的竖向位移约束，并施加0.02 m的竖向位移。

图 6显示了活动门反力与竖向位移的关系曲线，即地层特征曲线^[25]。从图 6可以发现随着活动门向下移动，作用在其上的反力逐渐减小，并趋向于一个稳定值，即极限状态。将图中红框范围内的曲线放大可以更清晰地看出，相对于共轴模型计算得到反力位移曲线，非共轴模型计算得到的结果需要经历更大的位移才能达到稳定值。非共轴模量越小，曲线下降的趋势越平缓，达到极限状态所需的位移更大。值得注意的是当非共轴模量小到一定程度时，反力曲线不再与共轴模型的结果趋于相同的稳定值，而是在经历了更大的位移后达到一个更大的极限反力值。

图  6  活动门反力位移曲线

Figure  6.  Reaction curves of tropdoor

下载: 全尺寸图片 幻灯片

图 7为共轴和非共轴模型条件下分析活动门上方土体塑性区发展分布情况。图中采用应力比$\eta$表征土的塑性发展程度，并且每幅图的右边展示了从图左侧提取的应力比介于1与1.2的分布情况，这些塑性发展显著的红色区域可以认为是土体极限破坏时产生的塑性滑动面。图 7（a）显示了在共轴情况下，土体的塑性变形展现出径直向上发展，然后向两侧延伸的发展趋势。右边的塑性集中区则显示滑动面在向上发展时会更倾向于向对称轴一侧延伸出分支，从而形成拱形滑动面，并且在滑动面继续竖向发展的过程形成拱形分支。当考虑非共轴变形后，土体中滑动面竖直向上延伸并向对称轴一侧成拱的发展趋势随着非共轴效应的增强（或非共轴模量的减小）而减弱，但是沿着背离对称轴斜向上方向发展的趋势却在增强。如图 7（e）所示，非共轴塑性变形足够强（或非共轴模量足够小）时，滑动面已经变成两条交叉的直线。活动门竖向范围内形成三角滑动面，同时另一条滑动面向外不断延伸形成一个斜坡状。图 7（a），（e）两种破坏模式的显著差异在一定程度上决定了图 6中两种模型对应的反力曲线极限值的差异。图 7（a）中的拱形滑动面主要反映的是活动门竖向范围内松动土体的影响，而图 7（e）中的向外延伸的滑动面则反映了活动门竖向范围以外的土体向活动门滑动而产生的影响。

图  7  不同模型条件下土体破坏时的塑性区分布

Figure  7.  Distribution of plastic parameter at failure with different model parameters

下载: 全尺寸图片 幻灯片

图 8展示了图 5中对称轴位置处不同模型参数条件下计算得到的土体竖向位移沿深度的分布情况。每条曲线中两个红色实心符号与图 5中轴线上塑性集中区域的起始位置相对应，并且举例说明了与红色实心符号标记图 7（a）中共轴情况下塑性区范围的情况。图 8中红框的放大图更清晰地展示了各曲线的相对位置关系以及各曲线中塑性区的起始位置。总体上，图中各曲线都表现为沉降随深度开始变化平缓，然后出现明显的变化，再又出现变化放缓的趋势。

图  8  沿深度分布的沉降值

Figure  8.  Distribution of settlements along depth

下载: 全尺寸图片 幻灯片

结合图 7进行对比分析可以发现曲线变化放缓的区域都是弹性区，曲线陡变的区域正是夹在两个红色实心符号标记之间的塑性区域。并且，随着非共轴模量的减小，塑性区的起始位置越深，而终止的位置则相差不大，即竖向的塑性区变得范围越小。

基于图 7，8确定不同模型条件下活动门上方塑性区高度和剪切带高度展示在图 9中。由该图可知，在共轴条件下，活动门上方土体的塑性区高度和剪切带高度均大于等于非共轴情况，并且随着非共轴模量由2G减小到0.2G，塑性区高度减小了1.3 m（由3.5 m减小到2.2 m），剪切带高度减小了0.7 m（由2.6 m减小到1.9 m）。以上结果表明，非共轴特性显著影响了土体内部塑性区和剪切带的发展。

图  9  不同模型条件下活动门上方塑性区高度和剪切带高度

Figure  9.  Heights of plastic zone and shear band above trapdoor with different model parameters

下载: 全尺寸图片 幻灯片

图 10显示了不同模型参数条件下计算得到的土体表面竖向位移沿水平方向的分布情况。显然，共轴情况下对称中心处的土体沉降最小。随着考虑非共轴塑性，非共轴模量越小，对称中心处的土体则沉降越来越大。同时，接近对称中心处的沉降曲线变化越陡。而在水平距离4 m以外，非共轴效应越强的沉降曲线沉降值则要小于弱非共轴效应的沉降曲线。

图  10  土体表面竖向位移分布曲线

Figure  10.  Curves of vertical displacement on soil surface

下载: 全尺寸图片 幻灯片

图 11显示了活动门下移时共轴模型计算得到土体应力主轴旋转分布情况。在初始时刻，土体内所有应力的主轴方向都与整体坐标方向平行。在比较绝对值的情况下，大主应力为竖向压应力。从图 11中可以发现，左侧对称轴附近存在“长颈”漏斗型的红色区域。该区域内的土体应力主轴基本上都逆时针旋转了90°，意味着土体的大主应力方向都是沿水平方向，在距漏斗型区域以外一定距离的绿色和蓝色区域，主应力偏转逐渐减小，表明在活动门上方形成了拱形大主应力迹线，即产生了土拱效应。在远离对称轴的区域，主应力轴的旋转程度则较小，甚至出现逆时针旋转的情况（黑色区域）。由此可见，拱形主应力迹线只能在一定范围内形成。

图  11  共轴情况下土体主应力轴旋转角分布

Figure  11.  Distribution of rotation angle of principal stress axis for coaxial model

下载: 全尺寸图片 幻灯片

模型考虑非共轴变形后，获得的土体内部应力主轴旋转的整体分布特征与图 11是相似的。但为了体现共轴与非共轴情况的局部差异，笔者将不同参数非共轴模型的结果与图 11中的结果作差来进行比较。图 12则说明了不同参数非共轴模型结果与共轴模型结果的差异。非共轴模型会在特定区域引起额外的应力轴旋转。随着非共轴模量的减小，应力主轴旋转的程度加剧。值得注意的是，参照图 7和图 11的结果，非共轴变形引起额外的应力主轴旋转主要发生在塑性发展较为集中的区域。

图  12  其他非共轴模型与共轴模型的旋转角分布差异

Figure  12.  Differences of rotation angle between coaxial model and other non-coaxial models

下载: 全尺寸图片 幻灯片

4.   结语

本文在ABAQUS/Explicit分析模块下，基于用户材料子程序VUMAT开发了基于双曲剪切硬化Drucker- Prager模型建立的非共轴模型的半隐式应力更新算法，推导了该算法的牛顿迭代格式。基于该Drucker-Prager模型和相应的非共轴模型，分析了土体单剪试验中的非共轴现象和活动门下移问题。分析结果表明该应力积分算法合理可行、强健有效。同时，采用共轴模型和非共轴模型对活动门问题的分析揭示出：引入非共轴塑性变形机制会使活动门的反力位移曲线在经历更大的位移时，才能达到极限状态。同时，在存在显著主应力轴旋转的区域产生更大的变形，如引起更大土表沉降。然而，当非共轴变形过于强烈时则会改变土体极限破坏形式，并造成反力位移曲线极限值的改变。由此可见，土体非共轴变形特性是工程变形预测以及极限破坏问题中不可忽视的一个因素。

附录A

关于第二节公式$\boldsymbol{C}_{ijkl}^n{}^n{\sigma _{kl}} = 0$的证明：

因为$\boldsymbol{C}_{i j k l}=\delta_{i k} \delta_{j l}-\frac{\delta_{i j} \delta_{k l}}{3}-\frac{S_{i j} S_{k l}}{S_{m n} S_{m n}}$，

以及$\sigma_{k l}=p \delta_{k l}+s_{k l}$，那么

考虑到$\delta_{i k} \delta_{j l} \delta_{k l}=\delta_{i j}, \delta_{i j} \delta_{i j}=3$，$\delta_{i k} \delta_{j l} s_{k l}=s_{i j}, \quad \delta_{i j} s_{i j}=0$，上式（A1）可化简为

命题得证。