共軛梯度法在彈塑性模型數(shù)值實(shí)現(xiàn)中的應(yīng)用

代 寧，耿大將，郭培軍，周順華，狄宏規(guī)

（1.同濟(jì)大學(xué)道路與交通工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 201804；2.同濟(jì)大學(xué)上海市軌道交通結(jié)構(gòu)耐久與系統(tǒng)安全重點(diǎn)實(shí)驗(yàn)室，上海 201804；3.新城控股集團(tuán)股份有限公司上海第二分公司，上海 201800；4.麥克馬斯特大學(xué)土木工程系，漢密爾頓L8S4L7）

彈塑性本構(gòu)模型的數(shù)值算法主要包括顯式算法與隱式算法兩大類。通常，顯式算法［1-3］計(jì)算精度較低，實(shí)際應(yīng)用少，而隱式算法計(jì)算精度較高，能實(shí)現(xiàn)自校正，不會(huì)發(fā)生誤差傳遞，實(shí)際應(yīng)用普遍。隱式算法一般采用彈性預(yù)測(cè)-塑性修正算法，在塑性修正步一般會(huì)用牛頓-最近點(diǎn)投影法（Newton-CPPM）求解回退映射非線性方程組，而Newton-CPPM 需確保Jacobian 矩陣是可逆矩陣，否則無(wú)法進(jìn)行求解。同時(shí)，Newton-CPPM 是一種局部收斂性算法，迭代初值選取的合理性直接影響算法收斂性。因此，Jacobian 矩陣奇異與Jacobian 矩陣不收斂是用傳統(tǒng)隱式算法進(jìn)行高度非線性彈塑性本構(gòu)模型數(shù)值實(shí)現(xiàn)所遇到的主要問(wèn)題。

針對(duì)該問(wèn)題，國(guó)內(nèi)外學(xué)者開展了許多研究并取得了一定進(jìn)展。通常，在采用具有局部收斂性的Newton-CPPM 時(shí)，只要所選取的初值十分接近真實(shí)解，該算法最后總可以收斂到問(wèn)題的真解。因此，為改進(jìn)算法的迭代初值，Bicanic 等［4］和Stupkiewicz等［5］采用輔助投影面法初步解決了迭代初值的合理選取問(wèn)題。然而，針對(duì)如何有效地建立輔助投影面這個(gè)問(wèn)題，其并未做出很好的解答。Valentini等［6］、Jia［7］、Majid 等［8］、Penasa 等［9］、Wang 等［10］、Lee等［11］采用多步法，即先將一個(gè)增量步分為多個(gè)增量步，再對(duì)方程進(jìn)行求解，最終實(shí)現(xiàn)了迭代初值的改進(jìn)。然而，該方法顯然會(huì)大幅增加計(jì)算量，從而降低計(jì)算效率。Hernandez 等［12］采用顯式方法實(shí)現(xiàn)了對(duì)部分狀態(tài)變量迭代初值的改進(jìn)，然而該方法仍未解決對(duì)于改進(jìn)哪些狀態(tài)變量的迭代初值就可以有效改進(jìn)算法收斂性的問(wèn)題。因此，為改善收斂性，耿大將等［13］引入同倫算法以改善收斂性并避免Jacobian 矩陣奇異。為了避免Jacobian 矩陣求解時(shí)巨大的計(jì)算量，Homel等［14-15］和Sharifian等［16］采用分階段迭代算法來(lái)求解高度非線性回退映射方程組，首先求得一致性參數(shù)，再將一致性參數(shù)視為已知量，最終求解其余狀態(tài)變量。值得注意的是，單獨(dú)應(yīng)用分階段迭代算法求解方程組仍存在結(jié)果不收斂問(wèn)題，對(duì)于高度非線性的本構(gòu)模型，不收斂問(wèn)題更加顯著。此外，在求解非線性回退映射方程組方面，Bilotta 等［17］、Placidi［18］、Armoro［19］也做了大量工作，主要采用優(yōu)化方法或搜索技術(shù)，而這類方法的計(jì)算特點(diǎn)是計(jì)算過(guò)程復(fù)雜、計(jì)算量較大，故該方法在彈塑性本構(gòu)模型數(shù)值實(shí)現(xiàn)中的應(yīng)用十分有限。

嘗試用共軛梯度法［20-22］對(duì)傳統(tǒng)隱式算法進(jìn)行改進(jìn)，以避免在高度非線性彈塑性本構(gòu)模型的隱式算法實(shí)現(xiàn)中出現(xiàn)Jacobian矩陣奇異和不收斂問(wèn)題。

1 傳統(tǒng)隱式算法

彈塑性本構(gòu)模型的數(shù)值實(shí)現(xiàn)主要包括狀態(tài)變量的更新與一致切線剛度的計(jì)算兩部分內(nèi)容。

1.1 狀態(tài)變量的更新

在隱式算法對(duì)彈塑性本構(gòu)模型進(jìn)行狀態(tài)變量更新的過(guò)程中，第（n+1）個(gè)增量步所對(duì)應(yīng)的基本方程如下所示：

式中：σ和Δε分別為應(yīng)力和應(yīng)變?cè)隽?；C為彈性剛度；Δλ為一致性參數(shù)；r為塑性流動(dòng)方向；ξ和分別為硬化內(nèi)變量和相應(yīng)的硬化方向；φ為屈服函數(shù)。以應(yīng)力σ為自變量，對(duì)塑性勢(shì)函數(shù)求偏導(dǎo)即可得出塑性流動(dòng)方向r。硬化內(nèi)變量可以有多個(gè)，具體數(shù)量取決于具體的本構(gòu)模型。

狀態(tài)變量更新的過(guò)程實(shí)際上就是在滿足Δλ≥0的前提下對(duì)以σn+1、ξn+1和Δλ為自變量的回退映射非線性方程組式（1）進(jìn)行求解的過(guò)程。目前，多采用彈性預(yù)測(cè)-塑性修正算法進(jìn)行求解，求解過(guò)程主要包括以下三個(gè)步驟。

（1）彈性預(yù)測(cè)。

帶下標(biāo)tr 的量表示的是經(jīng)彈性預(yù)測(cè)后得出的狀態(tài)變量。

（2）塑性檢驗(yàn)。如果φtr，n+1≤0，這說(shuō)明材料處于彈性狀態(tài)，否則應(yīng)該進(jìn)行塑性修正。

（3）塑性修正。一般采用Newton-CPPM 求解以σn+1、ξn+1和Δλ為未知量的非線性方程組式（1），迭代初值一般取為彈性預(yù)測(cè)后所得狀態(tài)變量的值，如下所示：

在率無(wú)關(guān)彈塑性力學(xué)框架內(nèi)建立或改造本構(gòu)模型主要可以通過(guò)以下兩種途徑實(shí)現(xiàn)：提高彈性關(guān)系、屈服函數(shù)、塑性勢(shì)函數(shù)或硬化規(guī)則的復(fù)雜程度，適當(dāng)增加內(nèi)變量。如對(duì)于結(jié)構(gòu)性軟土的彈塑性本構(gòu)模型，為反映土體結(jié)構(gòu)性的影響，祝恩陽(yáng)等［23］、Suebsuk等［24］和Dafalias 等［25-26］均在經(jīng)典修正劍橋模型［27］（MCC）的基礎(chǔ)上進(jìn)行了改進(jìn)。

無(wú)論采用哪種改進(jìn)模式均會(huì)導(dǎo)致回退映射方程組式（1）非線性程度的提高。式（1）的第一式中包含彈性剛度Cn+1和塑性流動(dòng)方向rn+1，當(dāng)彈性剛度Cn+1的表達(dá)式中包含應(yīng)力σn+1時(shí)，顯然第一式的非線性程度會(huì)提高；當(dāng)塑性勢(shì)函數(shù)的復(fù)雜程度提高或內(nèi)變量增加時(shí)，相應(yīng)的塑性流動(dòng)方向rn+1的復(fù)雜程度也會(huì)提高，所含內(nèi)變量的數(shù)量也會(huì)增加，同樣會(huì)導(dǎo)致第一式非線性程度的提高。第二式主要反映的是硬化規(guī)則的影響，包含硬化方向，當(dāng)硬化規(guī)則的復(fù)雜程度提高后，顯然第二式的非線性程度會(huì)提高，當(dāng)內(nèi)變量的數(shù)目增加時(shí)，第二式中所含等式的數(shù)量會(huì)做相應(yīng)的增加，顯然也會(huì)造成第二式非線性程度的提高。第三式主要體現(xiàn)的是屈服函數(shù)，屈服函數(shù)變復(fù)雜時(shí)，第三式的非線性程度顯然會(huì)提高。

當(dāng)要進(jìn)行數(shù)值實(shí)現(xiàn)的彈塑性本構(gòu)模型非線性程度較高時(shí)，相應(yīng)的回退映射方程組式（1）的非線性程度也會(huì)較高。對(duì)于高度非線性的彈塑性本構(gòu)模型，采用隱式的彈性預(yù)測(cè)-塑性修正算法進(jìn)行狀態(tài)變量更新容易出現(xiàn)Jacobian 矩陣奇異和不收斂問(wèn)題。這是因?yàn)閺椥灶A(yù)測(cè)-塑性修正算法主要是基于Newton-CPPM，而Newton-CPPM 是一種局部收斂方法，只有當(dāng)初值接近于問(wèn)題的真解時(shí)收斂性才較好，同時(shí)Newton-CPPM 必須保證Jacobian 矩陣是可逆的。對(duì)于高度非線性的回退映射方程組式（1），這兩點(diǎn)都很難保證。結(jié)合方程組式（1）的構(gòu)成和Newton-CPPM的特點(diǎn)可以看出，可能導(dǎo)致Jacobian矩陣奇異或不收斂的原因有以下幾個(gè)方面：

（1）彈性預(yù)測(cè)點(diǎn)離真解較遠(yuǎn)。例如，對(duì)于同時(shí)可以考慮硬化和軟化的模型，當(dāng)真解在軟化區(qū)，而根據(jù)彈性預(yù)測(cè)得出的迭代初值在硬化區(qū)，這種情況下可能導(dǎo)致Jacobian矩陣奇異或不收斂。

（2）組成回退映射方程組式（1）的方程數(shù)量級(jí)相差較大。例如，考慮軟土結(jié)構(gòu)性的Saniclay 模型［26］，包含表示旋轉(zhuǎn)硬化的內(nèi)變量，數(shù)量級(jí)為10-1左右，包含表示結(jié)構(gòu)性的內(nèi)變量，數(shù)量級(jí)為100～101，包含表示各向同性硬化的內(nèi)變量，數(shù)量級(jí)可以達(dá)到102以上，將數(shù)量級(jí)相差較大的方程組成方程組然后求解，可能導(dǎo)致Jacobian矩陣奇異或不收斂。

（3）彈性剛度較為復(fù)雜。例如，修正劍橋模型的彈性剛度依賴于靜水壓力，屬于壓力依賴型模型，可能導(dǎo)致不收斂。

（4）硬化規(guī)則較為復(fù)雜導(dǎo)致局部迭代過(guò)程中硬化內(nèi)變量的改變過(guò)大。例如，修正劍橋模型中代表屈服面大小的硬化內(nèi)變量在迭代過(guò)程中可能由正值轉(zhuǎn)化為負(fù)值，導(dǎo)致不收斂。

（5）應(yīng)力空間中，屈服面或塑性勢(shì)面的曲率較大或存在尖角。例如，考慮軟土結(jié)構(gòu)性的Saniclay模型的屈服面和塑性勢(shì)面曲率比較大，Mohr-Coulomb模型的屈服面存在尖角，可能導(dǎo)致Jacobian 矩陣奇異或不收斂。

1.2 一致切線剛度的計(jì)算

當(dāng)狀態(tài)變量在彈性或彈塑性狀態(tài)時(shí)，相應(yīng)的一致切線剛度D顯然是不同的。

（1）彈性狀態(tài)時(shí)，對(duì)于σn+1和εn+1，兩者均滿足方程σn+1-σn-Cn+1∶Δεn+1=0，將σn+1和Δεn+1同時(shí)視為變量，對(duì)該式微分，即可確定一致切線剛度。

（2）彈塑性狀態(tài)時(shí)，σn+1和εn+1滿足回退映射方程組式（1），將σn+1、ξn+1、Δλ、Δεn+1視作變量，對(duì)該方程組微分，得到微分方程組，然后消去dξn+1與dΔλ即可得到dσn+1和dΔεn+1間的關(guān)系，從而確定一致切線剛度。

2 改進(jìn)隱式算法

彈塑性本構(gòu)模型數(shù)值實(shí)現(xiàn)的關(guān)鍵在于狀態(tài)變量的更新，而狀態(tài)變量的更新關(guān)鍵在于塑性修正步。傳統(tǒng)的隱式算法在塑性修正步中需要采用Newton-CPPM 求解回退映射非線性方程組，但不能保證迭代初值位于收斂域內(nèi)，更不能保證Jacobian 矩陣可逆，對(duì)于高度非線性的彈塑性本構(gòu)模型更是如此。

為了避免由于組成回退映射方程組的方程數(shù)量級(jí)相差較大而造成迭代過(guò)程的不收斂，將式（1）修正為

式中：ρξ和ρφ表示權(quán)重系數(shù)。當(dāng)?shù)踔等閺椥灶A(yù)測(cè)結(jié)果時(shí)，權(quán)重系數(shù)的取值要能保證組成方程組（4）的各個(gè)方程對(duì)應(yīng)的殘差位于同一數(shù)量級(jí)。顯然，方程組（1）和（4）的解是完全相同的。

即便對(duì)高度非線性的方程組式（1）做如式（4）所示的修正，仍無(wú)法完全避免迭代求解過(guò)程的不收斂（這是由Newton-CPPM 的局部收斂性所決定的），更無(wú)法避免Jacobian 矩陣奇異問(wèn)題。為此，將方程組式（4）簡(jiǎn)記為

式中：x=(σn+1，ξn+1，Δλ)T。在解存在且唯一的情況下，方程組式（5）的解等價(jià)于如下最小化問(wèn)題的解：

共軛梯度法對(duì)應(yīng)的迭代公式為

上標(biāo)（k+1）和k分別表示第（k+1）次和第k次迭代。步長(zhǎng)h(k)由下式得出：

搜索方向d(k)由下式得出：

式中：?F為函數(shù)F(x)的梯度；βk-1為計(jì)算因子；gk為函數(shù)F(x)在x=x(k)處的梯度。

采用改進(jìn)隱式算法求解回退映射方程組式（1）的過(guò)程如圖1 所示。求解過(guò)程中，首先考慮數(shù)量級(jí)問(wèn)題，將方程組式（1）修正為式（4），然后將彈性預(yù)測(cè)結(jié)果作為塑性修正步的迭代初值，用Newton-CPPM 求解方程組式（4），求解過(guò)程中出現(xiàn)Jacobian矩陣奇異或不收斂問(wèn)題時(shí)，用共軛梯度法求解方程組式（4）對(duì)應(yīng)的最小化問(wèn)題式（6），將求得結(jié)果作為改進(jìn)后的迭代初值，再次采用Newton-CPPM 求解方程組式（4）。由于共軛梯度法是大范圍收斂算法，因此可以有效避免Newton-CPPM所存在的不收斂問(wèn)題。從共軛梯度法的求解過(guò)程可以看出，整個(gè)過(guò)程不需要進(jìn)行矩陣求逆運(yùn)算，因此也不會(huì)存在Jacobian矩陣奇異問(wèn)題。綜上，改進(jìn)隱式算法可以同時(shí)解決不收斂或Jacobian矩陣奇異問(wèn)題。

3 隱式算法比較

理論上，改進(jìn)隱式算法可以同時(shí)解決傳統(tǒng)隱式算法在高度非線性彈塑性本構(gòu)模型數(shù)值實(shí)現(xiàn)中所存在的不收斂和Jacobian 矩陣奇異問(wèn)題，但實(shí)際效果如何還有待檢驗(yàn)?？紤]軟土結(jié)構(gòu)性的Saniclay 模型具有較高的非線性，屬于壓力依賴性模型，該模型的塑性勢(shì)面和屈服面的曲率較大，含有多達(dá)五個(gè)的硬化內(nèi)變量，而且硬化內(nèi)變量的數(shù)量級(jí)相差較大，硬化規(guī)則的非線性程度較高，該模型也可以同時(shí)考慮硬化和軟化效應(yīng)。因此，以該模型為例對(duì)傳統(tǒng)隱式算法和改進(jìn)隱式算法進(jìn)行計(jì)算精度、收斂性和計(jì)算效率的對(duì)比。

圖1 改進(jìn)隱式算法Fig.1 Improved implicit algorithm

對(duì)比分析建立在單元體分析基礎(chǔ)上，計(jì)算中所采用的初始應(yīng)力狀態(tài)和加載應(yīng)變路徑分別如表1和表2 所示。表2 中，εa和εr分別表示軸向應(yīng)變和徑向應(yīng)變。計(jì)算中建立單個(gè)三維八節(jié)點(diǎn)實(shí)體單元，本構(gòu)模型參數(shù)取值參考文獻(xiàn)［13，26］。

如表3和表4所示，每種初始應(yīng)力狀態(tài)對(duì)應(yīng)的每種變應(yīng)變路徑均采用兩種算法進(jìn)行數(shù)值計(jì)算，并且每種算法均考慮了10-1、10-2、10-3和10-4等四種增量步長(zhǎng)In。表3、4 中，×表示結(jié)果不收斂，√表示結(jié)果收斂。對(duì)每種情況下的收斂性Cov、計(jì)算時(shí)間t、廣義剪應(yīng)力q的終值qt進(jìn)行統(tǒng)計(jì)對(duì)比。qt-εq曲線如圖2～5 所示。圖2～5 中，表示廣義剪應(yīng)變，其中J2ε表示偏應(yīng)變張量第二分量。

表1 初始應(yīng)力狀態(tài)Tab.1 Initial stress state

通過(guò)對(duì)應(yīng)力-應(yīng)變曲線的對(duì)比發(fā)現(xiàn)，對(duì)于特定的初始應(yīng)力狀態(tài)、特定的應(yīng)變路徑和特定的步長(zhǎng)，在收斂情況下，傳統(tǒng)隱式算法（K0）和改進(jìn)隱式算法（K1）總能得到近乎完全相同的結(jié)果，如圖2～5 所示。這是因?yàn)樵诔跏紤?yīng)力狀態(tài)、應(yīng)變路徑及步長(zhǎng)相同的情況下，兩種算法所對(duì)應(yīng)的基本方程都是式（1），在該方程的解存在且唯一的情況下，精度控制相同時(shí)，不論采用哪種算法，計(jì)算結(jié)果理應(yīng)一致。

表2 應(yīng)變路徑Tab.2 Strain path

表3 初始應(yīng)力狀態(tài)K0計(jì)算結(jié)果對(duì)比Tab.3 Comparison of calculation results under initial stress state K0

表4 初始應(yīng)力狀態(tài)K1計(jì)算結(jié)果對(duì)比Tab.4 Comparison of calculation results under initial stress state K1

對(duì)傳統(tǒng)隱式算法，當(dāng)采用平面應(yīng)變加載模式時(shí)，即使增量步減小到10-4時(shí)，仍然不收斂；當(dāng)采用不排水加載模式時(shí)，增量步減小至10-4時(shí)，才會(huì)達(dá)到收斂。這說(shuō)明，對(duì)于高度非線性彈塑性本構(gòu)模型，傳統(tǒng)隱式算法對(duì)步長(zhǎng)有很強(qiáng)的依賴性且收斂性很差。圖6 為初始應(yīng)力狀態(tài)K0 和平面應(yīng)變路徑下增量步為10-2時(shí)某一增量步迭代誤差Δ隨迭代次數(shù)n的變化。從圖6 可以看出，改進(jìn)隱式算法收斂性明顯優(yōu)于傳統(tǒng)隱式算法。這是因?yàn)閭鹘y(tǒng)隱式算法在求解彈塑性本構(gòu)模型中的高度非線性方程組時(shí)，對(duì)于局部收斂的Newton-CPPM，較大的增量步使彈性預(yù)測(cè)所得的迭代初值很難位于收斂域內(nèi)，因此極易出現(xiàn)不收斂現(xiàn)象。

圖2 初始應(yīng)力狀態(tài)K0下平面應(yīng)變加載所得qt-εq曲線Fig.2 qt-εq curve under initial stress state K0 and plane strain loading

圖3 初始應(yīng)力狀態(tài)K0下不排水加載所得qt-εq曲線Fig.3 qt-εq curve under initial stress state K0 and undrained loading

對(duì)于改進(jìn)隱式算法，當(dāng)采用平面應(yīng)變加載模式時(shí)，增量步為10-2時(shí)即可收斂；當(dāng)采用不排水加載模式時(shí)，增量步為10-1時(shí)即可收斂。顯然，改進(jìn)隱式算法基本不依賴于步長(zhǎng)且收斂性好。原因在于，改進(jìn)隱式算法中高度非線性方程組的求解問(wèn)題轉(zhuǎn)化為最小化問(wèn)題，通過(guò)采用具有全局收斂性的共軛梯度法求解最小化問(wèn)題，使得結(jié)果收斂性較好。

圖4 初始應(yīng)力狀態(tài)K1下平面應(yīng)變加載所得qt-εq曲線Fig.4 qt-εq curve under initial stress state K1 and plane strain loading

圖5 初始應(yīng)力狀態(tài)K1下不排水加載所得qt-εq曲線Fig.5 qt-εq curve under initial stress state K1 and undrained loading

圖6 兩種算法收斂性對(duì)比Fig.6 Comparison of convergence between two algorithms

在同樣的初始應(yīng)力狀態(tài)、加載路徑及增量步長(zhǎng)下，當(dāng)兩種算法均收斂時(shí)，如初始應(yīng)力狀態(tài)K0 不排水加載條件下，當(dāng)增量步為10-4時(shí)，改進(jìn)隱式算法采用傳統(tǒng)的Newton-CPPM可以成功求解高度非線性方程組，而不需要采用共軛梯度法求解相應(yīng)的最小化問(wèn)題，此時(shí)改進(jìn)隱式算法則退化為傳統(tǒng)隱式算法，因此兩種算法計(jì)算效率理應(yīng)相同。為了節(jié)約計(jì)算時(shí)間，實(shí)際工程計(jì)算中常采用變步長(zhǎng)方法，因此有必要對(duì)變步長(zhǎng)情況下兩種算法的計(jì)算效率進(jìn)行比較。如初始應(yīng)力狀態(tài)K0不排水加載時(shí)，傳統(tǒng)隱式算法在增量步為10-4時(shí)可以收斂，對(duì)應(yīng)計(jì)算時(shí)間為310.70 s；改進(jìn)隱式算法在增量步為10-1時(shí)可以收斂，對(duì)應(yīng)計(jì)算時(shí)間為0.20 s。顯然，整體上改進(jìn)隱式算法比傳統(tǒng)隱式算法計(jì)算效率要高得多。

綜合以上分析，與傳統(tǒng)隱式算法相比，改進(jìn)隱式算法收斂性好，整體計(jì)算效率高，計(jì)算精度基本相同。顯然，改進(jìn)改進(jìn)隱式算法能夠很好地解決傳統(tǒng)隱式算法所存在的Jacobian矩陣奇異和不收斂問(wèn)題。

4 算例

基坑開挖是土木工程領(lǐng)域常見的工程行為，以如圖7 所示的一小型無(wú)支護(hù)的基坑開挖為實(shí)例，基坑寬4.00 m，開挖深度4.00 m，考慮對(duì)稱性建立一半模型進(jìn)行計(jì)算?；油馏w分兩層，上部硬層采用Mohr-Coulomb 模型，下部軟土層采用考慮土體結(jié)構(gòu)性的Saniclay模型，具體模型參數(shù)的取值如表5所示，表5中各模型參數(shù)的物理意義可參考文獻(xiàn)［26］。計(jì)算過(guò)程主要包括兩個(gè)分析步，第一步為地應(yīng)力平衡，第二步進(jìn)行基坑開挖及鋼板樁施做，采用單元移除和激活功能實(shí)現(xiàn)開挖模擬和鋼板樁施做。

圖7 幾何模型（單位：m）Fig.7 Geometric model(unit：m)

表5 模型參數(shù)取值Tab.5 Parameter values of the model

當(dāng)采用改進(jìn)隱式算法進(jìn)行多單元數(shù)值計(jì)算時(shí)，在第二個(gè)分析步的最后一個(gè)增量步可以得到如圖8所示的位移結(jié)果，而采用傳統(tǒng)隱式算法時(shí)在第二個(gè)分析步的第一個(gè)增量步出現(xiàn)了不收斂現(xiàn)象。為明確不收斂的原因，對(duì)計(jì)算過(guò)程進(jìn)行實(shí)時(shí)監(jiān)控，發(fā)現(xiàn)第一個(gè)增量步中的傳統(tǒng)隱式算法無(wú)法對(duì)部分材料點(diǎn)的狀態(tài)變量進(jìn)行準(zhǔn)確且有效的更新，導(dǎo)致整體計(jì)算出現(xiàn)不收斂。由此可以表明，與傳統(tǒng)隱式算法相比，改進(jìn)隱式算法不僅適用于簡(jiǎn)單的單單元計(jì)算，還可以用于復(fù)雜的多單元分析，對(duì)于高度非線性的彈塑性本構(gòu)模型，改進(jìn)隱式算法顯然比傳統(tǒng)隱式算法更加可靠、有效。

圖8 地表沉降計(jì)算結(jié)果Fig.8 Calculation results of surface settlement

5 結(jié)論

（1）改進(jìn)隱式算法和傳統(tǒng)隱式算法的計(jì)算精度基本相同，但改進(jìn)隱式算法的收斂性和整體計(jì)算效率比傳統(tǒng)隱式算法要高很多。

（2）改進(jìn)隱式算法可以高效地解決傳統(tǒng)隱式算法所存在的不收斂和Jacobian矩陣奇異問(wèn)題。

作者貢獻(xiàn)聲明:

樂(lè) 云：分析當(dāng)前研究現(xiàn)狀，指導(dǎo)研究方向，審查研究?jī)?nèi)容和結(jié)果的合理性。

雪克來(lái)提·亥依熱特：負(fù)責(zé)具體的研究工作，收集并分析案例，組織驗(yàn)證實(shí)驗(yàn)，整理論文思路并撰寫論文。

李永奎：提供研究背景和主要研究方法。

虞 濤：提供案例，組織實(shí)驗(yàn)。