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

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

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

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

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

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

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

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

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

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

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

（1）彈性預測。

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

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

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

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

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

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

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

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

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

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

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

1.2 一致切線剛度的計算

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

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

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

2 改進隱式算法

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

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

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

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

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

共軛梯度法對應的迭代公式為

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

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

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

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

3 隱式算法比較

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

圖1 改進隱式算法Fig.1 Improved implicit algorithm

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

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

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

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

表2 應變路徑Tab.2 Strain path

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

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

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

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

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

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

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

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

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

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

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

4 算例

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

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

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

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

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

5 結(jié)論

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

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

作者貢獻聲明:

樂 云：分析當前研究現(xiàn)狀，指導研究方向，審查研究內(nèi)容和結(jié)果的合理性。

雪克來提·亥依熱特：負責具體的研究工作，收集并分析案例，組織驗證實驗，整理論文思路并撰寫論文。

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

虞 濤：提供案例，組織實驗。