彭芳樂(lè)，李福林，白曉宇

（1.同濟(jì)大學(xué) 地下建筑與工程系，上海 200092；2.同濟(jì)大學(xué) 巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200092；3. 中國(guó)礦業(yè)大學(xué) 力學(xué)與建筑工程學(xué)院，江蘇 徐州 221116）

1 引 言

砂土材料，尤其是密實(shí)砂土在峰值后具有顯著的應(yīng)變軟化特性。軟化引起的材料強(qiáng)非線性問(wèn)題在現(xiàn)有的有限元分析中仍然非常難以處理。目前，大部分有限元計(jì)算程序仍使用增量加載的步驟和修正Newton-Raphson迭代來(lái)消除不平衡荷載[1]。該法與Newton-Raphson法的區(qū)別在于剛度矩陣不再需要在每個(gè)增量步進(jìn)行重組或重新分解，但修正Newton-Raphson帶來(lái)的高效卻因?yàn)槭諗啃缘慕档投窒?，特別是在后期迭代時(shí)。Newton型求解方法雖然取得了很多重要的成果（如共扼牛頓法，割線牛頓法，準(zhǔn)牛頓法等），但由于這種方法需要大量的存儲(chǔ)空間以及計(jì)算復(fù)雜，不是非常適合材料的強(qiáng)非線性問(wèn)題的求解。而動(dòng)態(tài)松弛法由于每次求解過(guò)程中不需要形成整體剛度矩陣，且要求內(nèi)存空間少，迭代計(jì)算簡(jiǎn)單，從而為材料的強(qiáng)非線性問(wèn)題的求解提供了一條有效的數(shù)值計(jì)算途徑。動(dòng)態(tài)松弛首先出現(xiàn)在20世紀(jì)60年代中期Otter[2]和Day[3]的論文中。Welsh[4]、Cassel等[5]引入虛擬質(zhì)量對(duì)動(dòng)態(tài)松弛法做了改進(jìn)。Rushton[6]將這種方法用于處理非線性問(wèn)題。Underwood[7]提出了可適應(yīng)的動(dòng)態(tài)松弛方法，Papadrakakis[8]提出如何自動(dòng)選取動(dòng)態(tài)松弛迭代參數(shù)的方法。Tanaka和 Kawamoto[9－10]將動(dòng)態(tài)松弛法應(yīng)用到三維彈塑性有限元中，成功地預(yù)測(cè)了土體結(jié)構(gòu)的破壞荷載。Siddiquee[11]首先提出動(dòng)態(tài)松弛法的相同路徑追蹤算法，用來(lái)求解材料的非線性問(wèn)題。

為解決砂土應(yīng)變軟化引起的強(qiáng)非線性問(wèn)題，本文提出一種動(dòng)態(tài)松弛有限元法。文中對(duì)動(dòng)態(tài)松弛有限元技術(shù)進(jìn)行詳細(xì)闡述，給出了動(dòng)態(tài)松弛有限元法的基本控制方程、直接位移控制的平衡路徑追蹤以及應(yīng)力更新的回歸映射算法等。最后將動(dòng)態(tài)松弛有限元法應(yīng)用于砂土平面應(yīng)變壓縮試驗(yàn)的數(shù)值計(jì)算，驗(yàn)證了動(dòng)態(tài)松弛有限元法在求解砂土應(yīng)變軟化特性引起的材料強(qiáng)非線性問(wèn)題方面的優(yōu)勢(shì)。在應(yīng)力更新算法和數(shù)值分析中，材料模型采用砂土的3要素彈黏塑性本構(gòu)模型[12－15]。

2 動(dòng)態(tài)松弛有限元求解技術(shù)

2.1 動(dòng)態(tài)松弛法

非線性有限元的平衡控制方程為

式中：P為內(nèi)力矢量；Pinit為初始應(yīng)力引起的節(jié)點(diǎn)力矢量；F為外力矢量；B為幾何矩陣；N為有限元中離散的單元個(gè)數(shù)；σ為單元高斯點(diǎn)上的應(yīng)力矢量；Ve為每個(gè)單元的體積。

動(dòng)態(tài)松弛法的核心思想就是將上述靜力學(xué)問(wèn)題看作動(dòng)力學(xué)問(wèn)題來(lái)處理，即通過(guò)求解動(dòng)力學(xué)運(yùn)動(dòng)方程的穩(wěn)態(tài)響應(yīng)來(lái)獲得上述平衡控制方程的解。動(dòng)態(tài)松弛法在式（1）中引入假想力，即慣性力和阻尼力

動(dòng)態(tài)松弛法是一種顯式算法。將中心差分法應(yīng)用于式（2），可得t時(shí)刻的速度矢量和加速度矢量，分別為

式中：Δt為時(shí)間增量；ut-Δt、ut和ut+Δt分別為t-Δt、t和t+Δt時(shí)刻的位移矢量。

顯式算法的每一步求解均為矩陣乘法，系數(shù)矩陣均可以轉(zhuǎn)換成對(duì)角陣。為了保證動(dòng)態(tài)松弛法的可靠性，質(zhì)量矩陣m選用對(duì)角陣，則由如下關(guān)系可得阻尼矩陣為

其中，α為阻尼因子，它是確定動(dòng)態(tài)松弛法求解時(shí)能否快速達(dá)到穩(wěn)態(tài)解的最關(guān)鍵因素。如果任意給定一個(gè)阻尼因子，就可能需要大量的計(jì)算時(shí)間來(lái)獲得臨界阻尼值。因此，許多研究者建議了多種確定臨界阻尼因子的方法[6,8,16－18]?？紤]到動(dòng)態(tài)松弛法中阻尼力是假想力的特點(diǎn)，只要能夠使系統(tǒng)振蕩趨于收斂，任意估計(jì)的臨界阻尼因子均可接受。在本研究中，采用Rayleigh商確定臨界阻尼因子[19]

式中：x為特征向量；K為整體剛度矩陣；M為整體質(zhì)量矩陣。這里，x可用t時(shí)刻的速度矢量t近似。M用對(duì)角陣m來(lái)近似可達(dá)到理想效果，如集中質(zhì)量。K可用對(duì)角切線剛度矩陣Klt來(lái)近似取值[7]，具體形式如下

式中：Pt和Pt-Δt分別為t和t-Δt時(shí)刻的內(nèi)力矢量。將上述這些近似形式代入到式（6）中，可得到臨界阻尼因子的估計(jì)值為

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

于是式（9）可簡(jiǎn)化為

因此，若已知 ut-Δt和ut，則根據(jù)式（11）即可求得ut+Δt。那么求解時(shí)僅需要按時(shí)間逐步推進(jìn)，直至達(dá)到容許誤差進(jìn)而得到穩(wěn)態(tài)解。式（11）是一個(gè)時(shí)間推進(jìn)的顯式迭代方程，時(shí)間步的選擇是條件穩(wěn)定的。根據(jù)Courant-Friedrichs-Lewy條件[20]，可得

式中：β為控制計(jì)算穩(wěn)定性和收斂速度的因子（β≤1.0）。對(duì)于常應(yīng)變單元，β值接近 1.0，但對(duì)于高階單元，β值需大大減小。β值越小，數(shù)值計(jì)算的穩(wěn)定性越好，但達(dá)到精度所需的迭代次數(shù)越多，計(jì)算時(shí)間越長(zhǎng)，收斂越慢；l為單元上相鄰節(jié)點(diǎn)之間的最小距離；Vc是介質(zhì)中的彈性壓縮波速，將Vc的表達(dá)式代入式（12），可得

其中：ρ為密度；v為泊松比；E為彈性模量。在式（13）的基礎(chǔ)上，引入一種換算方法，實(shí)質(zhì)上是為了對(duì)網(wǎng)格均勻化[7,21]。即在單元的基礎(chǔ)上由式（13）計(jì)算虛擬質(zhì)量密度，以使壓縮波穿越單元（任意尺寸和剛度）所需的時(shí)間相同。在虛擬動(dòng)力分析中，Δt可任意選定，若假定固定時(shí)間步長(zhǎng)Δt=1，由式（13）即可求出每個(gè)單元的虛擬質(zhì)量密度為

利用這種方法處理虛擬動(dòng)力響應(yīng)分析，變形波和不平衡力在迭代過(guò)程中能夠均勻地傳遞到有限元網(wǎng)格模型的整個(gè)區(qū)域，而不需考慮有限元網(wǎng)格離散時(shí)的不均勻性以及介質(zhì)非線性引起的單元?jiǎng)偠鹊牟灰恢滦?。在?dòng)態(tài)松弛有限元法中惟一需調(diào)整的參數(shù)是β，其他參數(shù)在平衡迭代過(guò)程中自適應(yīng)調(diào)整。在對(duì)砂土及砂土地基的有限元計(jì)算中[14－15,22－26]，β設(shè)為0.6可兼顧顯式動(dòng)態(tài)松弛計(jì)算的穩(wěn)定性和運(yùn)算速度。

此外，動(dòng)態(tài)松弛法在追蹤應(yīng)力-應(yīng)變空間中整個(gè)平衡路徑的方法有3種：①荷載控制法；②直接位移控制法；③弧長(zhǎng)法。在本研究中，采用直接位移控制法對(duì)整個(gè)平衡路徑進(jìn)行追蹤。位移控制法中，在每個(gè)載荷步某些位移分量是固定的，其他的位移分量是需要求解的。固定位移分量可通過(guò)下式表示：

其他的位移分量根據(jù)式（11）計(jì)算。位移控制的動(dòng)態(tài)松弛法的流程圖如圖1所示，其中Δus為各載荷步之間的增量位移。

圖1 位移控制的動(dòng)態(tài)松弛法的流程圖Fig.1 Flow chart for the direct displacement control solution

2.2 收斂準(zhǔn)則

計(jì)算時(shí)，收斂狀態(tài)由以下3個(gè)標(biāo)準(zhǔn)來(lái)判定。

（1）整體殘余力相對(duì)標(biāo)準(zhǔn)：

（2）整體內(nèi)能增量（殘余力在對(duì)應(yīng)位移增量上所作的功）標(biāo)準(zhǔn)：

（3）兩次相鄰迭代（n和n+1）間的整體殘余力絕對(duì)標(biāo)準(zhǔn)：

2.3 回歸映射法

在有限元分析中，在給定增量位移后，高斯點(diǎn)上的應(yīng)力需要進(jìn)行更新，相應(yīng)的算法稱為應(yīng)力更新算法或本構(gòu)積分算法（即材料本構(gòu)模型的數(shù)值算法）。利用動(dòng)態(tài)松弛法迭代時(shí)，對(duì)應(yīng)變和應(yīng)力增量進(jìn)行更新。在這里，材料模型取砂土的3要素彈黏塑性本構(gòu)模型[12－15]為例進(jìn)行說(shuō)明。

應(yīng)變?cè)隽?dε(n)由動(dòng)態(tài)松弛法第n步迭代平衡位置的位移計(jì)算得出

式中：dε(n)為第n步迭代時(shí)的應(yīng)變?cè)隽?；du(n)為第n步迭代時(shí)的位移增量。則無(wú)黏性應(yīng)力增量 dσf(n)可通過(guò)下式計(jì)算：

式中：dσf(n)為第n步迭代時(shí)的無(wú)黏性應(yīng)力增量；Dep為彈塑性矩陣。若利用式（20）進(jìn)行應(yīng)力更新，則每次均需計(jì)算彈塑性矩陣，使計(jì)算變得異常復(fù)雜，不利于在有限元程序中實(shí)現(xiàn)。

在動(dòng)態(tài)松弛有限元法中，對(duì)于無(wú)黏性彈塑性方程的積分，采用回歸映射法進(jìn)行應(yīng)力更新，它是一種近似的一階向后歐拉積分算法[27]。在利用回歸映射法進(jìn)行應(yīng)力更新時(shí)，由總應(yīng)變的增量驅(qū)動(dòng)彈性預(yù)測(cè)狀態(tài)，而由無(wú)黏性參數(shù)的增量驅(qū)動(dòng)無(wú)黏性修正狀態(tài)。因此，在彈性預(yù)測(cè)階段，不可恢復(fù)應(yīng)變和內(nèi)變量保持固定；而在無(wú)黏性修正階段，總應(yīng)變是不變的。假定是通過(guò)彈性模型求得的應(yīng)力（彈性預(yù)測(cè)）。當(dāng)超出更新無(wú)黏性屈服面，應(yīng)該回到更新無(wú)黏性屈服面對(duì)應(yīng)的應(yīng)力。與回歸路徑相關(guān)的彈性應(yīng)變dεe為

式中：De為彈性矩陣。

在無(wú)黏性修正階段，總應(yīng)變保持不變，即總應(yīng)變?cè)隽縟ε=0。根據(jù)dε=dεe+dεir，彈性與不可恢復(fù)應(yīng)變?cè)隽看嬖谌缦玛P(guān)系

假定流動(dòng)法則為不相關(guān)聯(lián)的，無(wú)黏性勢(shì)函數(shù)為G，則由正交法則可得

相應(yīng)地不可恢復(fù)應(yīng)變?yōu)?/p>

同時(shí)，硬化參量κ也發(fā)生相應(yīng)的變化，即

假定，硬化參量的增量dκ與無(wú)黏性因子λ存在如下關(guān)系

于是，式（26）可改寫為

將式（24）和（28）代入到無(wú)黏性屈服函數(shù)F(σf, κ)=0中，可得

對(duì)式（29）進(jìn)行Taylor展開，并忽略高階項(xiàng)，則λ可表示為

在回歸過(guò)程中，通過(guò)式（30）可以迭代計(jì)算無(wú)黏性修正值，回歸路徑如圖2所示。每個(gè)迭代步，無(wú)黏性應(yīng)力和硬化參量根據(jù)式（24）和（28）同步更新，直至達(dá)到規(guī)定的容差范圍內(nèi)，即 F(, κ) ≈ 0。B

在砂土三要素彈黏塑性本構(gòu)模型中，無(wú)黏性勢(shì)函數(shù)G為Drucker-Prager形式，硬化參量為修正的不可恢復(fù)應(yīng)變能Wir*

圖2 回歸映射法示意圖Fig.2 Schematic diagram illustrating return mapping algorithm

根據(jù)式（23）和（31），可得硬化參量的增量

對(duì)于Drucker-Prager形式的無(wú)黏性勢(shì)函數(shù)，式（32）可改寫為

即

于是，式（30）可以改寫為

在無(wú)黏性應(yīng)力回映迭代的同時(shí)，根據(jù)3要素模型中的速率相關(guān)的非線性黏性模型計(jì)算黏性應(yīng)力，進(jìn)而求得總應(yīng)力，使之滿足屈服條件。根據(jù)非線性3要素模型進(jìn)行應(yīng)力更新的具體流程如圖3所示，圖中變量的意義請(qǐng)參見文獻(xiàn)[12－15]。

2.4 縮減積分

在有限元分析中，需要進(jìn)行單元體積內(nèi)或面積內(nèi)的積分。對(duì)于等參元來(lái)說(shuō)，通常情況下是無(wú)法進(jìn)行顯式積分的，必須借助于數(shù)值積分。高斯求積是最常用的一種數(shù)值積分方法。在很多情況下，縮減積分[28]常?？梢缘玫捷^精確積分更好的結(jié)果。但縮減積分會(huì)使單元產(chǎn)生沙漏（hourglass）或零能模式，這會(huì)導(dǎo)致系統(tǒng)剛度矩陣奇異，從而影響了求解結(jié)果。因此，在實(shí)際分析中，必須防止沙漏或零能模式的出現(xiàn)。在動(dòng)態(tài)松弛有限元法中，采用等參元及縮減積分，并通過(guò)彈性剛度法[29]控制沙漏，即在任何土體單元開始形成沙漏時(shí)，將一很小的彈性剛度（真實(shí)材料剛度的 0.05%）附加到非線性系統(tǒng)中作為抵抗沙漏的節(jié)點(diǎn)力，這在砂土地基承載力有限元計(jì)算中已得到驗(yàn)證[15, 25－26]。

圖3 根據(jù)非線性3要素模型進(jìn)行應(yīng)力更新的流程圖Fig.3 Flow chart for stress updating according to nonlinear three-component model

3 有限元模擬

本文所提的動(dòng)態(tài)松弛有限元法已成功應(yīng)用于無(wú)加筋砂土、土工格柵加筋砂土及加筋砂土地基的變形和強(qiáng)度特征的數(shù)值計(jì)算[14－15,22－26]，但砂土多是采用彈塑性本構(gòu)模型，未能考慮速率相關(guān)的黏性特性。此處通過(guò)對(duì)密實(shí)砂土平面應(yīng)變壓縮試驗(yàn)的典型應(yīng)力-應(yīng)變關(guān)系進(jìn)行彈黏塑性有限元模擬，驗(yàn)證動(dòng)態(tài)松弛有限元法在求解材料非線性問(wèn)題方面的優(yōu)越性，尤其是峰值應(yīng)力狀態(tài)之后的應(yīng)變軟化特性。

試驗(yàn)砂為飽和的Toyoura砂（D50=0.21 mm，Uc=1.2，Gs=2.65，emax=0.98，emin=0.62），試樣尺寸（寬×高×長(zhǎng)）為 96 mm×120 mm×62 mm（σ2方向）。初始孔隙比e0為0.665，相對(duì)密度Dr為87.39 %，初始圍壓σ3為30 kPa。隨后以軸向應(yīng)變速率為0.04%/min經(jīng)行平面應(yīng)變壓縮試驗(yàn)。對(duì)此試驗(yàn)分別采用單個(gè)單元和多個(gè)單元進(jìn)行有限元模擬，網(wǎng)格劃分如圖4所示，尺寸同試樣σ2平面上的尺寸。材料模型采用3要素彈黏塑性本構(gòu)模型[12－15]。有限元計(jì)算時(shí)，先施加初始圍壓30 kPa，再在網(wǎng)格上下兩端節(jié)點(diǎn)上施加豎直方向的位移速度。圖5為有限元計(jì)算結(jié)果，其中應(yīng)力比R=σ1/σ3（σ1和σ3分別為豎直和水平方向的有效應(yīng)力）和豎直應(yīng)變?chǔ)舦均取整個(gè)有限元網(wǎng)格的平均值?？梢钥闯?，不管是單個(gè)單元還是多個(gè)單元，均能對(duì)密實(shí)砂土典型應(yīng)力-應(yīng)變關(guān)系進(jìn)行很好的模擬，尤其是峰值后應(yīng)變軟化階段，體現(xiàn)了動(dòng)態(tài)松弛有限元法的特點(diǎn)。

圖4 有限元網(wǎng)格Fig.4 FEM meshes

圖5 密實(shí)砂土平面應(yīng)變壓縮試驗(yàn)典型應(yīng)力-應(yīng)變關(guān)系的有限元模擬Fig.5 FEM simulation of stress-strain relation of dense sand from PSC test

4 結(jié) 論

（1）針對(duì)砂土材料應(yīng)變軟化引起的強(qiáng)非線性問(wèn)題的數(shù)值解法進(jìn)行探討，闡述了動(dòng)態(tài)松弛有限元法的相關(guān)關(guān)鍵問(wèn)題。通過(guò)中心差分法導(dǎo)出了動(dòng)態(tài)松弛有限元基本控制方程及相關(guān)參數(shù)的確定方法。

（2）對(duì)于無(wú)黏性彈塑性本構(gòu)方程的積分使用了一階向后歐拉積分的回歸映射法，并在無(wú)黏性應(yīng)力回映迭代的同時(shí)，根據(jù)3要素模型中的速率相關(guān)的非線性黏性模型計(jì)算黏性應(yīng)力，進(jìn)而對(duì)總應(yīng)力進(jìn)行更新。

（3）有限元中采用等參元及縮減積分，提高了高度非線性材料的偽平衡中的求解界限。為了避免任何可能的沙漏模式，在動(dòng)態(tài)松弛有限元法中使用了反沙漏方案。

（4）通過(guò)單個(gè)單元和多個(gè)單元的對(duì)密實(shí)砂土典型應(yīng)力-應(yīng)變關(guān)系進(jìn)行有限元數(shù)值計(jì)算，表明動(dòng)態(tài)松弛有限元法可以很好地模擬砂土材料應(yīng)變軟化引起的強(qiáng)非線性問(wèn)題，尤其是峰值后的軟化破壞特性。

[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] 彭芳樂(lè), 李福林, 李建中, 等. 加載速率變化條件下砂土的黏塑特性及本構(gòu)模型[J]. 巖石力學(xué)與工程學(xué)報(bào),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] 彭芳樂(lè), 李福林, 白曉宇, 等. 考慮應(yīng)力路徑和加載速率效應(yīng)砂土的彈黏塑性本構(gòu)模型[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 彭芳樂(lè), 小竹望, 龍岡文夫. 土工格柵加筋砂土的變形與破壞機(jī)理解析[J]. 巖土力學(xué), 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] 彭芳樂(lè), 龍岡文夫, 廖少明, 等. 土工格柵加筋砂土的二維等價(jià)有限元解析方法[J]. 同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版), 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] 彭芳樂(lè), 龍岡文夫. 加筋砂土地基中加筋寬板效果的數(shù)值解析研究(英文)[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 彭芳樂(lè), 希迪克, 龍岡文夫. 砂土地基中加筋深度效果研究(英文)[J]. 巖石力學(xué)與工程學(xué)報(bào), 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.