針對砂土應變軟化強非線性問題的動態松弛有限元法研究

彭芳樂，李福林，白曉宇

（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.