针对砂土应变软化强非线性问题的动态松弛有限元法研究

彭芳乐，李福林，白晓宇

（1.同济大学 地下建筑与工程系，上海 200092；2.同济大学 岩土及地下工程教育部重点实验室，上海 200092；3. 中国矿业大学 力学与建筑工程学院，江苏 徐州 221116）

1 引 言

砂土材料，尤其是密实砂土在峰值后具有显著的应变软化特性。软化引起的材料强非线性问题在现有的有限元分析中仍然非常难以处理。目前，大部分有限元计算程序仍使用增量加载的步骤和修正Newton-Raphson迭代来消除不平衡荷载[1]。该法与Newton-Raphson法的区别在于刚度矩阵不再需要在每个增量步进行重组或重新分解，但修正Newton-Raphson带来的高效却因为收敛性的降低而抵消，特别是在后期迭代时。Newton型求解方法虽然取得了很多重要的成果（如共扼牛顿法，割线牛顿法，准牛顿法等），但由于这种方法需要大量的存储空间以及计算复杂，不是非常适合材料的强非线性问题的求解。而动态松弛法由于每次求解过程中不需要形成整体刚度矩阵，且要求内存空间少，迭代计算简单，从而为材料的强非线性问题的求解提供了一条有效的数值计算途径。动态松弛首先出现在20世纪60年代中期Otter[2]和Day[3]的论文中。Welsh[4]、Cassel等[5]引入虚拟质量对动态松弛法做了改进。Rushton[6]将这种方法用于处理非线性问题。Underwood[7]提出了可适应的动态松弛方法，Papadrakakis[8]提出如何自动选取动态松弛迭代参数的方法。Tanaka和 Kawamoto[9－10]将动态松弛法应用到三维弹塑性有限元中，成功地预测了土体结构的破坏荷载。Siddiquee[11]首先提出动态松弛法的相同路径追踪算法，用来求解材料的非线性问题。

为解决砂土应变软化引起的强非线性问题，本文提出一种动态松弛有限元法。文中对动态松弛有限元技术进行详细阐述，给出了动态松弛有限元法的基本控制方程、直接位移控制的平衡路径追踪以及应力更新的回归映射算法等。最后将动态松弛有限元法应用于砂土平面应变压缩试验的数值计算，验证了动态松弛有限元法在求解砂土应变软化特性引起的材料强非线性问题方面的优势。在应力更新算法和数值分析中，材料模型采用砂土的3要素弹黏塑性本构模型[12－15]。

2 动态松弛有限元求解技术

2.1 动态松弛法

非线性有限元的平衡控制方程为

式中：P为内力矢量；Pinit为初始应力引起的节点力矢量；F为外力矢量；B为几何矩阵；N为有限元中离散的单元个数；σ为单元高斯点上的应力矢量；Ve为每个单元的体积。

动态松弛法的核心思想就是将上述静力学问题看作动力学问题来处理，即通过求解动力学运动方程的稳态响应来获得上述平衡控制方程的解。动态松弛法在式（1）中引入假想力，即惯性力和阻尼力

动态松弛法是一种显式算法。将中心差分法应用于式（2），可得t时刻的速度矢量和加速度矢量，分别为

式中：Δt为时间增量；ut-Δt、ut和ut+Δt分别为t-Δt、t和t+Δt时刻的位移矢量。

显式算法的每一步求解均为矩阵乘法，系数矩阵均可以转换成对角阵。为了保证动态松弛法的可靠性，质量矩阵m选用对角阵，则由如下关系可得阻尼矩阵为

其中，α为阻尼因子，它是确定动态松弛法求解时能否快速达到稳态解的最关键因素。如果任意给定一个阻尼因子，就可能需要大量的计算时间来获得临界阻尼值。因此，许多研究者建议了多种确定临界阻尼因子的方法[6,8,16－18]。考虑到动态松弛法中阻尼力是假想力的特点，只要能够使系统振荡趋于收敛，任意估计的临界阻尼因子均可接受。在本研究中，采用Rayleigh商确定临界阻尼因子[19]

式中：x为特征向量；K为整体刚度矩阵；M为整体质量矩阵。这里，x可用t时刻的速度矢量t近似。M用对角阵m来近似可达到理想效果，如集中质量。K可用对角切线刚度矩阵Klt来近似取值[7]，具体形式如下

式中：Pt和Pt-Δt分别为t和t-Δt时刻的内力矢量。将上述这些近似形式代入到式（6）中，可得到临界阻尼因子的估计值为

将式（3）～（5）代入到式（2）中，可得

于是式（9）可简化为

因此，若已知 ut-Δt和ut，则根据式（11）即可求得ut+Δt。那么求解时仅需要按时间逐步推进，直至达到容许误差进而得到稳态解。式（11）是一个时间推进的显式迭代方程，时间步的选择是条件稳定的。根据Courant-Friedrichs-Lewy条件[20]，可得

式中：β为控制计算稳定性和收敛速度的因子（β≤1.0）。对于常应变单元，β值接近 1.0，但对于高阶单元，β值需大大减小。β值越小，数值计算的稳定性越好，但达到精度所需的迭代次数越多，计算时间越长，收敛越慢；l为单元上相邻节点之间的最小距离；Vc是介质中的弹性压缩波速，将Vc的表达式代入式（12），可得

其中：ρ为密度；v为泊松比；E为弹性模量。在式（13）的基础上，引入一种换算方法，实质上是为了对网格均匀化[7,21]。即在单元的基础上由式（13）计算虚拟质量密度，以使压缩波穿越单元（任意尺寸和刚度）所需的时间相同。在虚拟动力分析中，Δt可任意选定，若假定固定时间步长Δt=1，由式（13）即可求出每个单元的虚拟质量密度为

利用这种方法处理虚拟动力响应分析，变形波和不平衡力在迭代过程中能够均匀地传递到有限元网格模型的整个区域，而不需考虑有限元网格离散时的不均匀性以及介质非线性引起的单元刚度的不一致性。在动态松弛有限元法中惟一需调整的参数是β，其他参数在平衡迭代过程中自适应调整。在对砂土及砂土地基的有限元计算中[14－15,22－26]，β设为0.6可兼顾显式动态松弛计算的稳定性和运算速度。

此外，动态松弛法在追踪应力-应变空间中整个平衡路径的方法有3种：①荷载控制法；②直接位移控制法；③弧长法。在本研究中，采用直接位移控制法对整个平衡路径进行追踪。位移控制法中，在每个载荷步某些位移分量是固定的，其他的位移分量是需要求解的。固定位移分量可通过下式表示：

其他的位移分量根据式（11）计算。位移控制的动态松弛法的流程图如图1所示，其中Δus为各载荷步之间的增量位移。

图1 位移控制的动态松弛法的流程图Fig.1 Flow chart for the direct displacement control solution

2.2 收敛准则

计算时，收敛状态由以下3个标准来判定。

（1）整体残余力相对标准：

（2）整体内能增量（残余力在对应位移增量上所作的功）标准：

（3）两次相邻迭代（n和n+1）间的整体残余力绝对标准：

2.3 回归映射法

在有限元分析中，在给定增量位移后，高斯点上的应力需要进行更新，相应的算法称为应力更新算法或本构积分算法（即材料本构模型的数值算法）。利用动态松弛法迭代时，对应变和应力增量进行更新。在这里，材料模型取砂土的3要素弹黏塑性本构模型[12－15]为例进行说明。

应变增量 dε(n)由动态松弛法第n步迭代平衡位置的位移计算得出

式中：dε(n)为第n步迭代时的应变增量；du(n)为第n步迭代时的位移增量。则无黏性应力增量 dσf(n)可通过下式计算：

式中：dσf(n)为第n步迭代时的无黏性应力增量；Dep为弹塑性矩阵。若利用式（20）进行应力更新，则每次均需计算弹塑性矩阵，使计算变得异常复杂，不利于在有限元程序中实现。

在动态松弛有限元法中，对于无黏性弹塑性方程的积分，采用回归映射法进行应力更新，它是一种近似的一阶向后欧拉积分算法[27]。在利用回归映射法进行应力更新时，由总应变的增量驱动弹性预测状态，而由无黏性参数的增量驱动无黏性修正状态。因此，在弹性预测阶段，不可恢复应变和内变量保持固定；而在无黏性修正阶段，总应变是不变的。假定是通过弹性模型求得的应力（弹性预测）。当超出更新无黏性屈服面，应该回到更新无黏性屈服面对应的应力。与回归路径相关的弹性应变dεe为

式中：De为弹性矩阵。

在无黏性修正阶段，总应变保持不变，即总应变增量dε=0。根据dε=dεe+dεir，弹性与不可恢复应变增量存在如下关系

假定流动法则为不相关联的，无黏性势函数为G，则由正交法则可得

相应地不可恢复应变为

同时，硬化参量κ也发生相应的变化，即

假定，硬化参量的增量dκ与无黏性因子λ存在如下关系

于是，式（26）可改写为

将式（24）和（28）代入到无黏性屈服函数F(σf, κ)=0中，可得

对式（29）进行Taylor展开，并忽略高阶项，则λ可表示为

在回归过程中，通过式（30）可以迭代计算无黏性修正值，回归路径如图2所示。每个迭代步，无黏性应力和硬化参量根据式（24）和（28）同步更新，直至达到规定的容差范围内，即 F(, κ) ≈ 0。B

在砂土三要素弹黏塑性本构模型中，无黏性势函数G为Drucker-Prager形式，硬化参量为修正的不可恢复应变能Wir*

图2 回归映射法示意图Fig.2 Schematic diagram illustrating return mapping algorithm

根据式（23）和（31），可得硬化参量的增量

对于Drucker-Prager形式的无黏性势函数，式（32）可改写为

即

于是，式（30）可以改写为

在无黏性应力回映迭代的同时，根据3要素模型中的速率相关的非线性黏性模型计算黏性应力，进而求得总应力，使之满足屈服条件。根据非线性3要素模型进行应力更新的具体流程如图3所示，图中变量的意义请参见文献[12－15]。

2.4 缩减积分

在有限元分析中，需要进行单元体积内或面积内的积分。对于等参元来说，通常情况下是无法进行显式积分的，必须借助于数值积分。高斯求积是最常用的一种数值积分方法。在很多情况下，缩减积分[28]常常可以得到较精确积分更好的结果。但缩减积分会使单元产生沙漏（hourglass）或零能模式，这会导致系统刚度矩阵奇异，从而影响了求解结果。因此，在实际分析中，必须防止沙漏或零能模式的出现。在动态松弛有限元法中，采用等参元及缩减积分，并通过弹性刚度法[29]控制沙漏，即在任何土体单元开始形成沙漏时，将一很小的弹性刚度（真实材料刚度的 0.05%）附加到非线性系统中作为抵抗沙漏的节点力，这在砂土地基承载力有限元计算中已得到验证[15, 25－26]。

图3 根据非线性3要素模型进行应力更新的流程图Fig.3 Flow chart for stress updating according to nonlinear three-component model

3 有限元模拟

本文所提的动态松弛有限元法已成功应用于无加筋砂土、土工格栅加筋砂土及加筋砂土地基的变形和强度特征的数值计算[14－15,22－26]，但砂土多是采用弹塑性本构模型，未能考虑速率相关的黏性特性。此处通过对密实砂土平面应变压缩试验的典型应力-应变关系进行弹黏塑性有限元模拟，验证动态松弛有限元法在求解材料非线性问题方面的优越性，尤其是峰值应力状态之后的应变软化特性。

试验砂为饱和的Toyoura砂（D50=0.21 mm，Uc=1.2，Gs=2.65，emax=0.98，emin=0.62），试样尺寸（宽×高×长）为 96 mm×120 mm×62 mm（σ2方向）。初始孔隙比e0为0.665，相对密度Dr为87.39 %，初始围压σ3为30 kPa。随后以轴向应变速率为0.04%/min经行平面应变压缩试验。对此试验分别采用单个单元和多个单元进行有限元模拟，网格划分如图4所示，尺寸同试样σ2平面上的尺寸。材料模型采用3要素弹黏塑性本构模型[12－15]。有限元计算时，先施加初始围压30 kPa，再在网格上下两端节点上施加竖直方向的位移速度。图5为有限元计算结果，其中应力比R=σ1/σ3（σ1和σ3分别为竖直和水平方向的有效应力）和竖直应变εv均取整个有限元网格的平均值。可以看出，不管是单个单元还是多个单元，均能对密实砂土典型应力-应变关系进行很好的模拟，尤其是峰值后应变软化阶段，体现了动态松弛有限元法的特点。

图4 有限元网格Fig.4 FEM meshes

图5 密实砂土平面应变压缩试验典型应力-应变关系的有限元模拟Fig.5 FEM simulation of stress-strain relation of dense sand from PSC test

4 结 论

（1）针对砂土材料应变软化引起的强非线性问题的数值解法进行探讨，阐述了动态松弛有限元法的相关关键问题。通过中心差分法导出了动态松弛有限元基本控制方程及相关参数的确定方法。

（2）对于无黏性弹塑性本构方程的积分使用了一阶向后欧拉积分的回归映射法，并在无黏性应力回映迭代的同时，根据3要素模型中的速率相关的非线性黏性模型计算黏性应力，进而对总应力进行更新。

（3）有限元中采用等参元及缩减积分，提高了高度非线性材料的伪平衡中的求解界限。为了避免任何可能的沙漏模式，在动态松弛有限元法中使用了反沙漏方案。

（4）通过单个单元和多个单元的对密实砂土典型应力-应变关系进行有限元数值计算，表明动态松弛有限元法可以很好地模拟砂土材料应变软化引起的强非线性问题，尤其是峰值后的软化破坏特性。

[1] ZIENKIEWICZ O C, TAYLOR R L. The finite element method[M]. New York: McGraw-Hill, 1991.

[2] OTTER J R H. Computations for prestressed concrete reactor pressure vessels using dynamic relaxation[J].Nuclear Structure Engineering, 1965, 1(1): 61－75.

[3] DAY A S. An introduction to dynamic relaxation[J]. The Engineer, 1965, 219: 218－221.

[4] WELSH A K. Discussion on dynamic relaxation[C]//Proceedings of the Institution of Civil Engineers. [S. l.]:[s. n.], 1967.

[5] CASSEL A C, KINSEY P J, SEFTON D J. Cylindrical shell analysis by dynamic relaxation[C]//Proceedings of the Institution of Civil Engineers. [S. l.]: [s. n.], 1968.

[6] RUSHTON K R. Dynamic-relaxation solutions of elastic-plate problems[J]. The Journal of Strain Analysis for Engineering Design, 1968, 3(1): 23－32.

[7] UNDERWOOD P G. Dynamic relaxation - A review[C]//Proceedings of Computational Methods for Transient Analysis. North Holland: Amsterdam, 1983.

[8] PAPADRAKAKIS M. A method for the automatic evaluation of the dynamic relaxation parameters[J].Computer Methods in Applied Mechanics and Engineering, 1981, 25(1): 35－48.

[9] TANAKA T, KAWAMOTO O. Three-dimensional finite element collapse analysis for foundations and slopes using dynamic relaxation[C]//Proceedings of the 6th International Conference on Numerical Methods in Geomechanics. Rotterdam: A. A. Balkema, 1988: 1213－1218.

[10] TANAKA T, KAWAMOTO O. Numerical modelling for softening with localization[M]. California: Springer-Verlag, 1990.

[11] SIDDIQUEE M S A. Finite element analysis of settlement and bearing capacity of footing on sand[D].Tokyo: The University of Tokyo, 1994.

[12] TATSUOKA F, ISHIHARA M, DI BENEDETTO H, et al.Time-dependent shear deformation characteristics of geomaterials and their simulation[J]. Soils and Foundations, 2002, 42(2): 103－129.

[13] 彭芳乐, 李福林, 李建中, 等. 加载速率变化条件下砂土的黏塑特性及本构模型[J]. 岩石力学与工程学报,2008, 27(8): 1576－1585.PENG Fang-le, LI Fu-lin, LI Jian-zhong, et al.Viscoplastic behaviors and constitutive modeling of sands under change of loading rates[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(8): 1576－1585.

[14] 彭芳乐, 李福林, 白晓宇, 等. 考虑应力路径和加载速率效应砂土的弹黏塑性本构模型[J]. 岩石力学与工程学报, 2009, 28(5): 929－938.PENG Fang-le, LI Fu-lin, BAI Xiao-yu, et al. Elastoviscoplastic constitutive model of sandy soil considering stress path and loading rate[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 28(5): 929－938.

[15] PENG Fang-le, SIDDIQUE M S A, TATSUOKA F, et al.Strain energy-based elasto-viscoplastic constitutive modelling of sand for numerical simulation[J]. Soils and Foundations, 2009, 49(4): 611－629.

[16] CUNDALL P A. Explicit finite difference methods in geomechanics[C]//Proceedings of the EF Conference on Numerical Methods in Geomechanics, Vol.1. New York:ASCE, 1976: 132－150.

[17] BICĆANIĆ N, HINTON E. Spurious modes in twodimensional isoparametric elements[J]. International Journal of Numerical Methods in Engineering, 1979,14(10): 1545－1557.

[18] PICA A, HINTON E. Transient and pseudo-transient analysis of Mindlin plates[J]. International Journal of Numerical Methods in Engineering, 1980, 15(2): 189－208.

[19] BUNCE J W. A note on the estimation of critical damping in dynamic relaxation[J]. International Journal of Numerical Methods in Engineering, 1972, 4(2): 301－303.

[20] COURANT R, FRIEDRICHS K, LEWY H. On the partial difference equations of mathematical physics[J].Mathematische Annalen, 1928, 100(1): 32－74.

[21] KEY S W, STONE C M, KRIEG R D. Dynamic relaxation applied to the quasi-static, large deformation,inelastic response of axisymmetric solids[C]//Proceedings of the Europe-US Workshop on Nonlinear Finite Element Analysis in Structural Mechanics. Berlin: Springer, 1981:585－620.

[22] PENG Fang-le, KOTAKE N, TATSUOKA F, et al. Plane strain compression behaviour of geogrid- reinforced sand and its numerical analysis[J]. Soils and Foundations,2000, 40(3): 55－74.

[23] 彭芳乐, 小竹望, 龙冈文夫. 土工格栅加筋砂土的变形与破坏机理解析[J]. 岩土力学, 2004, 25(6): 843－849.PENG Fang-le, KOTAKE N, TATSUOKA F. Numerical analysis of deformation and failure mechanism of geogrid-reinforced sand[J]. Rock and Soil Mechanics,2004, 25(6): 843－849.

[24] 彭芳乐, 龙冈文夫, 廖少明, 等. 土工格栅加筋砂土的二维等价有限元解析方法[J]. 同济大学学报(自然科学版), 2005, 33(2): 143－148.PENG Fang-le, TATSUOKA F, LIAO Shao-ming, et al.Equivalent two dimensional finite element method of plane strain compression on geogrid reinforced sand[J].Journal of Tongji University (Natural Science), 2005,33(2): 143－148.

[25] 彭芳乐, 龙冈文夫. 加筋砂土地基中加筋宽板效果的数值解析研究(英文)[J]. 岩石力学与工程学报, 2005,24(2): 268－277.PENG Fang-le, TATSUOKA F. Numerical study for wide-slab effect on reinforced sandy ground[J]. Chinese Journal of Rock Mechanics and Engineering, 2005,24(2): 268－277.

[26] 彭芳乐, 希迪克, 龙冈文夫. 砂土地基中加筋深度效果研究(英文)[J]. 岩石力学与工程学报, 2005, 24(9): 1512－1521.PENG Fang-le, SIDDIQUE M S A, TATSUOKA F. Deepfooting reinforcing mechanism on reinforced sandy ground[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(9): 1512－1521.

[27] ORTIZ M, SIMO J C. An analysis of a new class of integration algorithms for elastoplastic constitutive relations[J]. International Journal for Numerical Methods in Engineering, 1986, 23(3): 353－366.

[28] ZIENKIEWICZ O C, TAYLOR R L, TOO J M. Reduced integration technique in general analysis of plates and shells[J]. International Journal for Numerical Methods in Engineering, 1971, 3(2): 275－290.

[29] FLANAGAN D P, BELYTSCHKO T. A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[J]. International Journal for Numerical Methods in Engineering, 1981, 17(5): 679－706.