高度非線性彈塑性模型的改進顯式算法數(shù)值實現(xiàn)

耿大將， 郭培軍,2， 周順華

(1. 同濟大學(xué) 道路與交通工程教育部重點實驗室，上海 201804； 2. 麥克馬斯特大學(xué) 土木工程系,漢密爾頓 L8S4L7)

天然形成的結(jié)構(gòu)性軟土與重塑軟土有不同的物理力學(xué)性質(zhì)，天然軟土的結(jié)構(gòu)性會隨土體變形的發(fā)展逐漸散失，土體強度會逐漸降低，甚至?xí)霈F(xiàn)應(yīng)變軟化效應(yīng)，而且天然結(jié)構(gòu)性軟土一般是各向異性的[1-3].因此，在對處于天然結(jié)構(gòu)性軟土中的工程行為進行數(shù)值模擬時用重塑軟土的本構(gòu)模型是不盡合理的[4-6]，應(yīng)該用能夠反映軟土的結(jié)構(gòu)性的本構(gòu)模型.鑒于此，國內(nèi)外眾多學(xué)者和研究人員提出了各種不同的描述原狀結(jié)構(gòu)性軟黏土力學(xué)性質(zhì)的本構(gòu)模型[7-11].為了更好地反映原狀結(jié)構(gòu)性軟土的力學(xué)特性，一般情況下，這些模型有很強的非線性性.模型建立完成后，必須先對其進行數(shù)值算法實現(xiàn)，才能將其用于工程數(shù)值計算中.因此，數(shù)值實現(xiàn)是從模型理論到工程應(yīng)用的關(guān)鍵環(huán)節(jié).對于高度非線性的彈塑性本構(gòu)模型，數(shù)值實現(xiàn)的難度比較大.本文以考慮軟土結(jié)構(gòu)性的SANICLAY模型[7]為例，對高度非線性彈塑性本構(gòu)模型的顯式算法實現(xiàn)展開探討.

常微分方程組的求解方法一般可分為隱式算法與顯式算法兩大類.與之對應(yīng)的彈塑性模型的數(shù)值實現(xiàn)方法也可以分為隱式算法與顯式算法兩大類.當(dāng)采用隱式算法時，在塑性修正步一般需要采用Newton迭代法求解回退映射非線性方程組.對于高度非線性的彈塑性本構(gòu)模型，Newton迭代法所對應(yīng)的Jacobian矩陣很難求出，即便求出形式也較復(fù)雜.而且Newton迭代法需要保證Jacobian矩陣是可逆的，否則無法進行求解，對于高度非線性的本構(gòu)模型，這點也是很難保證的.例如，對于考慮軟土結(jié)構(gòu)性的SANICLAY模型，在應(yīng)變硬化過渡到應(yīng)變軟化的階段，Jacobian矩陣不可逆.因此，隱式算法用于高度非線性彈塑性本構(gòu)模型數(shù)值實現(xiàn)的難度比較大.與隱式算法相比較，顯式算法一般不需要迭代求解復(fù)雜的非線性方程組，因此單個增量步的計算量小，但其計算精度低，并且存在誤差的積累.

對于高度非線性彈塑性本構(gòu)模型，如考慮土體結(jié)構(gòu)性的SANICLAY模型，當(dāng)采用傳統(tǒng)的顯式算法對模型進行數(shù)值實現(xiàn)時，計算精度難以控制.在求解本構(gòu)方程所對應(yīng)的常微分方程組時，為了提高顯式算法的計算精度，Sloan等[12]、Cabot等[13]、Gonzalez等[14]將向前Euler法與多步法結(jié)合使用.多步法本質(zhì)上是將一個增量步分為多個增量步進行求解，實際上是減小了增量步步長，這樣必然會大幅降低計算效率，而且多步法并不能完全消除傳統(tǒng)的全顯式算法所具有的誤差積累和精度較低的缺點.在實際應(yīng)用中，為了提高計算精度，往往會將增量步步長取得很小，這無疑會降低計算效率.

采用全顯式算法進行高度非線性彈塑性本構(gòu)模型的數(shù)值實現(xiàn)時，為提高計算效率和計算精度，本文采用四階Dormand and Prince Runge-Kutta法[15](DPRK法)代替?zhèn)鹘y(tǒng)的全顯式算法中的向前Euler法，并且與切平面法[16-17]相結(jié)合，形成了改進顯式算法.最后，對傳統(tǒng)的全顯式算法、改進顯式算法和隱式算法在計算收斂性、計算效率和計算精度方面進行了對比，得出了一些有益的結(jié)論.

1 SANICLAY模型簡介

在原始SANICLAY彈塑性本構(gòu)模型[18]的基礎(chǔ)上，Taiebat等[7]引入了表征土體強度和旋轉(zhuǎn)硬化的內(nèi)變量，以考慮隨著變形的發(fā)展，結(jié)構(gòu)性的逐漸散失，強度的逐漸降低和各向異性的演化.該模型的屈服函數(shù)和塑性勢函數(shù)是不同的，具體如下.

1.1 彈性關(guān)系

模型的彈性關(guān)系依賴于球應(yīng)力p和孔隙比e，具體為

(1)

C=λ′δδ+2GI

(2)

λ′表示拉梅常數(shù);δ與I分別表示二階和四階等同張量;K、G分別表示體變模量和剪切模量;μ表示泊松比;k表示e-lnp空間內(nèi)回彈線的斜率.

1.2 塑性勢函數(shù)及流動法則

塑性勢函數(shù)

(3)

式中:s表示偏應(yīng)力；α為表征塑性勢函數(shù)旋轉(zhuǎn)的張量；pα確保塑性勢面過當(dāng)前應(yīng)力點；M*=SfM;Sf表征摩擦結(jié)構(gòu)性;M通過應(yīng)力羅德角θα在拉壓臨界狀態(tài)線的斜率Me和Mc間插值確定

其中

χα為模型參數(shù).將塑性勢函數(shù)ψ對應(yīng)力σ求導(dǎo)可得塑性流動方向r=?ψ/?σ.

1.3 屈服函數(shù)

屈服函數(shù)

(4)

圖1 屈服函數(shù)與塑性勢函數(shù)

1.4 硬化規(guī)則

模型有Si、Sf、p0、α和β共5個內(nèi)變量，對應(yīng)的演化規(guī)則為

(5)

(6)

(7)

(8)

(9)

1.5 一致性條件

與屈服函數(shù)相應(yīng)的一致性條件

(10)

將硬化規(guī)則式(5)～(9)代入可以得到一致性參數(shù)

(11)

其中，塑性模量

(12)

2 顯式算法數(shù)值實現(xiàn)

考慮土體結(jié)構(gòu)性的SANICLAY模型考慮的因素較全面，使得模型比較復(fù)雜，導(dǎo)致模型具有很強的非線性，數(shù)值實現(xiàn)的難度較大.結(jié)合本文第一部分可以看出，在采用顯式算法對考慮土體結(jié)構(gòu)性的SANICLAY模型進行數(shù)值實現(xiàn)的過程中，可能導(dǎo)致計算精度差的因素包括：①彈性關(guān)系依賴于球應(yīng)力和孔隙比；②屈服函數(shù)和塑性勢函數(shù)的曲率較大；③模型考慮了各向異性——旋轉(zhuǎn)硬化；④硬化內(nèi)變量較多，且硬化規(guī)則具有較高的非線性性.總之，一般情況下，會導(dǎo)致本構(gòu)模型的非線性程度增加的因素都會導(dǎo)致顯式算法的精度降低.

2.1 傳統(tǒng)的顯式算法

彈塑性本構(gòu)模型的顯式算法數(shù)值實現(xiàn)主要是進行狀態(tài)變量的更新.考慮土體結(jié)構(gòu)性的SANICLAY模型的狀態(tài)變量主要包括應(yīng)力σ、應(yīng)變ε、彈性應(yīng)變εe、塑性應(yīng)變εp、孔隙比e、表征各向同性結(jié)構(gòu)性的內(nèi)變量Si、表征摩擦結(jié)構(gòu)性的內(nèi)變量Sf、各向同性硬化內(nèi)變量p0、塑性勢函數(shù)和屈服函數(shù)中表征旋轉(zhuǎn)硬化的內(nèi)變量α與β.

彈塑性本構(gòu)模型狀態(tài)變量的更新一般采用應(yīng)變驅(qū)動模式，所對應(yīng)的基本問題為：給定滿足本構(gòu)模型的狀態(tài)變量的初值σ(t0)、e(t0)、Si(t0)、Sf(t0)、p0(t0)、α(t0)、β(t0)以及應(yīng)變歷史ε(t)，t∈[t0,T]，求σ(t)、e(t)、Si(t)、Sf(t)、p0(t)、α(t)和β(t)，使其滿足如下本構(gòu)方程：

(13)

當(dāng)采用向前Euler法對方程(13)進行離散可以得到，在第n+1個增量步對應(yīng)狀態(tài)變量的更新方程為

(14)

式(14)即為傳統(tǒng)的全顯式算法所對應(yīng)狀態(tài)變量的更新方程.

2.2 改進顯式算法

考慮到傳統(tǒng)的全顯式算法所對應(yīng)的向前Euler法只有一階精度，為提高顯式算法的計算效率，本文用四階DPRK法代替向前Euler法，對方程(13)進行離散，可以得到第n+1個增量步對應(yīng)狀態(tài)變量的更新方程變?yōu)?/p>

(15)

其中，

式中：Δσi、Δξi分別為應(yīng)力和硬化內(nèi)變量的第i次修正量，并非表示應(yīng)力增量Δσ和內(nèi)變量增量Δξ的分量；σni、ξni分別為計算過程中所采用的修正后的應(yīng)力和內(nèi)變量，并非表示應(yīng)力σn和內(nèi)變量ξn的分量.

理論上說，顯式算法經(jīng)過以上改進后計算效率可以得到提高，但是其所具有的誤差積累和精度較低的缺點仍無法消除，這在后文中的分析可以得到說明.造成顯式算法精度低和誤差積累的原因可以從顯式算法的原理得到解釋.顯式算法是基于已知的狀態(tài)變量的值直接外推下一增量步的狀態(tài)變量的值，在外推過程中并沒有對下一增量步的狀態(tài)變量進行有效的限制.例如，外推得出的下一增量步的狀態(tài)變量不一定能滿足屈服函數(shù)的一致性條件.因此，通過顯示算法更新后的狀態(tài)變量會存在誤差.同時，在本增量步的狀態(tài)變量存在誤差的情況下，顯式算法會直接用本增量步的有誤差的狀態(tài)變量外推下一增量步的狀態(tài)變量，顯然這樣會造成誤差的積累.因此，提高顯式算法精度的關(guān)鍵就在于提高每個增量步的計算精度，以減小誤差的積累和傳播.

在第n+1個增量步，采用式(15)所對應(yīng)的改進顯式算法進行狀態(tài)變量的更新，更新后的狀態(tài)變量會存在誤差.如圖2所示，根據(jù)更新后的狀態(tài)變量得出的屈服函數(shù)值Φn+1>0，這與傳統(tǒng)的彈塑性理論是相違背的，因此需要對更新后的狀態(tài)變量進行修正.本文嘗試采用切平面算法對更新后的狀態(tài)變量進行修正，以提高顯式算法的計算精度.具體采用的修正方程為

(16)

其中，上標(biāo)(k-1)和(k)表示修正次數(shù).

圖2 切平面算法原理

綜合本節(jié)內(nèi)容可以看出，改進顯式算法對傳統(tǒng)的全顯式算法做了計算效率和計算精度兩方面的改進.如圖3所示，為提高計算效率而引入了DPRK法，為改善計算精度對更新后的狀態(tài)變量采用切平面算法進行了修正.

圖3 改進顯式算法框圖

3 不同算法對比

為明確文中的改進顯式算法的計算精度和計算效率，同時也為后續(xù)研究工作中采用哪種算法提出建議，筆者將顯式算法、改進顯式算法和隱式算法進行對比.

本文進行單個單元的計算分析以模擬原狀土的三軸不排水壓縮試驗.具體在數(shù)值計算中采用的是8節(jié)點3維實體單元.在考慮軟土結(jié)構(gòu)性的SANICLAY模型中，考慮土體結(jié)構(gòu)性的內(nèi)變量主要包括各向同性結(jié)構(gòu)性內(nèi)變量Si和摩擦結(jié)構(gòu)性內(nèi)變量Sf，具體的模型參數(shù)取值參考文獻[7].軸向初始應(yīng)力σa=74.667 kPa，周向初始應(yīng)力σr=37.667 kPa.采用應(yīng)變控制式比例加載εa/εr=2/(-1)，加載至軸向應(yīng)變εa=0.150.在具體計算中，每種算法都分別考慮了增量步步長為10-1、10-2、10-3、10-4的情況.

計算結(jié)果及對比如表1所示.表1中，傳統(tǒng)顯式算法指的是與傳統(tǒng)的向前Euler法相對應(yīng)的顯式算法，改進顯式算法1指的是用四階的DPRK法代替向前Euler法所得出的顯式算法；改進顯式算法2指的是用四階的DPRK法代替向前Euler法，然后再采用切平面算法進行修正后所得出的顯式算法.改進顯式算法1和改進顯式算法2的主要區(qū)別在于是否采用切平面算法進行精度改善.表中的收斂性可以直接反映出相應(yīng)算法的收斂性，計算時間可以間接反映計算效率的高低，通過對比q-εq曲線可以反映算法的計算精度.當(dāng)增量步步長為10-1時，傳統(tǒng)顯式算法、改進顯式算法1、改進顯式算法2會得出錯誤的結(jié)果，而隱式算法會不收斂.

表1 計算結(jié)果對比

結(jié)合圖4和表1可以看出，對傳統(tǒng)顯式算法，當(dāng)增量步步長較大(如10-1和10-2)時，會得出錯誤的結(jié)果，只有在增量步步長較小時算法才會穩(wěn)定.與隱式算法相比，即使在增量步步長取值很小(如10-4)，計算時間很長的情況下，傳統(tǒng)顯式算法仍然存在誤差積累，計算精度仍然較差.為了盡可能提高傳統(tǒng)顯式算法的計算精度，其增量步步長必須取得相當(dāng)小，比如小于等于10-4，對應(yīng)的計算時間為397.30 s，計算效率會很低.

圖4 傳統(tǒng)顯式算法的計算結(jié)果

結(jié)合圖5和表1可以看出，對只做計算效率改進的顯式算法1，當(dāng)增量步步長特別大(如10-1)時會得出錯誤的結(jié)果，隨著增量步步長的減小，計算時間大幅增加，計算趨于收斂.當(dāng)增量步步長較大(如10-2)和步長較小(如10-4)時，計算精度差別不大.與傳統(tǒng)顯式算法相比，當(dāng)增量步長相同時，改進顯式算法1的計算效率要稍低，這是因為與改進顯式算法1相對應(yīng)的是DPRK法，而與傳統(tǒng)顯式算法相對應(yīng)的是向前Euler法，向前Euler法顯然比DPRK法計算量要小.當(dāng)增量步長可以發(fā)生變化時，改進顯式算法1在增量步長較大(如10-2)時即可達到與傳統(tǒng)顯式算法在增量步長較小(如10-4)時相同的精度，相應(yīng)的計算時間分別為3.40 s和397.30 s.所以，整體上在變步長情況下，改進顯式算法1比傳統(tǒng)顯式算法計算效率要高很多.與隱式算法相比較，改進顯式算法1與傳統(tǒng)顯式算法的計算精度一樣，都比較差.

圖5 改進顯式算法1的計算結(jié)果

結(jié)合圖6和表1可以看出，對做計算效率和計算精度改進的顯式算法2，當(dāng)增量步步長較大(如等于10-2)時，即可以得到收斂而且穩(wěn)定的結(jié)果.與傳統(tǒng)顯式算法和改進的顯式算法1相比較，當(dāng)增量步長相同時，改進顯式算法2的計算效率要低，這是因為與改進顯式算法2相對應(yīng)的是DPRK法，并且進行了精度修正.當(dāng)增量步長可以發(fā)生變化時，改進顯式算法2在增量步長較大(如10-2)時即可達到比傳統(tǒng)顯式算法在增量步長較小(如10-4)時高很多的計算精度，相應(yīng)的計算時間分別為3.70 s和397.30 s.所以，整體上在變步長情況下，改進顯式算法2也比傳統(tǒng)顯式算法計算效率要高很多.與傳統(tǒng)顯式算法相比較，改進顯式算法2的計算精度有了極大的改善，基本不存在誤差的積累，但是計算精度比隱式算法仍然略差.

圖6 改進顯式算法2的計算結(jié)果

綜合上述分析可以看出，與傳統(tǒng)顯式算法相比較，改進顯示算法在計算效率和計算精度方面均有較大的改善.計算效率的提高是因為與傳統(tǒng)顯式算法對應(yīng)的向前Euler法具有一階精度，而與改進顯式算法對應(yīng)的DPRK法具有四階精度.為了達到相同的計算精度，對于高階的算法，增量步步長可以取得大一些，這樣總的計算時間就會降低，計算效率就會提高.計算精度的提高是因為在每個增量步都采用切平面法對更新后的狀態(tài)變量進行了修正，使其較嚴(yán)格地滿足了屈服函數(shù)，這樣即可避免誤差的傳遞和積累.

4 多單元分析

本文的算法比較是建立在簡單應(yīng)變路徑下的單個單元分析基礎(chǔ)之上的.然而，實際工程計算往往需要的是多單元分析的結(jié)果.多單元分析中有眾多材料點，基本上每個材料點的應(yīng)變路徑或應(yīng)力路徑都是較為復(fù)雜的.因此，有必要對改進顯式算法在多單元分析中的適用性進行驗證分析.

以某隧道開挖工程為例進行多單元計算分析，所建立的分層土的幾何模型和模型參數(shù)取值如圖7所示，所采用的SANICLAY模型參數(shù)取值參考文獻[7] .隧道開挖采用斷面收縮法模擬，本次模擬中斷面收縮率取為8‰，最后計算結(jié)果如圖8所示.這說明改進顯式算法不但可用于單個單元的簡單應(yīng)變路徑計算，而且用于工程實際中的多單元復(fù)雜應(yīng)變路徑或應(yīng)力路徑分析也是可行的.

圖7 幾何模型及參數(shù)取值(尺寸單位：m)

圖8 Mises應(yīng)力云圖

5 結(jié)論

(1) 與隱式算法相比，傳統(tǒng)顯式算法雖然較為簡單，但其存在誤差的積累，計算精度較差，建議盡量不使用.

(2) 與傳統(tǒng)顯式算法相比，改進顯式算法的計算效率和計算精度都得到大幅提高.但是，與隱式算法相比較，其計算精度仍然略差.

(3) 改進顯式算法可以有效對高度非線性彈塑性本構(gòu)模型進行數(shù)值實現(xiàn).

本文對高度非線性彈塑性本構(gòu)模型的顯式算法數(shù)值實現(xiàn)做了一些研究工作.結(jié)合本文研究可以看出，當(dāng)對計算精度要求特別高時，對高度非線性彈塑性本構(gòu)模型的隱式算法實現(xiàn)，尤其是在不需要求Jacobian矩陣的情況下的數(shù)值實現(xiàn)，是極具價值的研究方向.

參考文獻:

[1] CALLISTO L, RAMPELLO S. Shear strength and small-strain stiffness of a natural clay under general stress conditions[J]. Geotechnique, 2002, 52(8):547.

[2] 孫德安, 陳波. 結(jié)構(gòu)性軟土力學(xué)特性的試驗研究[J]. 土木工程學(xué)報, 2011(增刊2):65.

SUN Dean, CHEN Bo. Experimental study on the mechanical behavior of structural soft clay[J]. China Civil Engineering Journal, 2011(Suppl.2):65.

[3] JUNG Y H, FINNO R J, CHO W. Stress-strain responses of reconstituted and natural compressible Chicago glacial clay[J]. Engineering Geology, 2012, 129:9.

[4] LIU W Z, SHI M L, MIAO L C,etal. Constitutive modeling of the destructuration and anisotropy of natural soft clay[J]. Computers and Geotechnics, 2013, 51:24.

[5] PANAYIDES S, ROUAINIA M, WOOD D M. Influence of degradation of structure on the behaviour of a full-scale embankment[J]. Canadian Geotechnical Journal, 2012, 49(3):344.

[6] ROCCHI G, VACIAGO G, FONTANA M,etal. Understanding sampling disturbance and behaviour of structured clays through constitutive modelling[J]. Soils and Foundations, 2013, 53(2):315.

[7] TAIEBAT M, DAFALIAS Y F, PEEK R. A destructuration theory and its application to SANICLAY model[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2010, 34(10):1009.

[8] SUEBSUK J, HORPIBULSUK S, LIU M D. Modified structured cam clay: a generalised critical state model for destructured, naturally structured and artificially structured clays[J]. Computers and Geotechnics, 2010, 37:956.

[9] 祝恩陽, 姚仰平.結(jié)構(gòu)性土UH模型[J]. 巖土力學(xué), 2015(11):3101.

ZHU Enyang, YAO Yangping. A UH constitutive model for structured soils[J]. Rock and Soil Mechanics, 2015(11):3101.

[10] CALLISTO L, GAJO A, WOOD D M. Simulation of triaxial and true triaxial tests on natural and reconstituted Pisa clay[J]. Geotechnique, 2002, 52(9):649.

[11] LIU M D, CARTER J P, AIREY D W. Sydney soil model. I: theoretical formulation[J]. International Journal of Geomechanics, 2011, 11(3):211.

[12] SLOAN S W, ABBO A J, SHENG D C. Refined explicit integration of elastoplastic models with automatic error control[J]. Engineering Computations, 2001, 18(1/2):121.

[13] CABOT M L, SLOAN S W, SHENG D C,etal. Error behaviour in explicit integration algorithms with automatic sub-stepping[J]. International Journal for Numerical Methods in Engineering, 2016, 108:1030.

[14] GONZALEZ N A, GENS A. Evaluation of a constitutive model for unsaturated soils: stress variables and numerical implementation[C]∥ALONSO E, GENS A. In Fifth International Conference on Unsaturated Soils Barcelona. Abingdon: Taylor & Francis Group, 2011: 829-836.

[15] HAIRER E, NORSETT S P, WANNER G. Solving ordinary differential equations I[M]. 2nd ed. Berlin: Springer-Verlag, 1993.

[16] ORTIZ M, POPOV E P. Accuracy and stability of integration algorithms for elastoplastic constitutive relations[J]. International Journal for Numerical Methods in Engineering, 1985, 21(9):1561.

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

[18] DAFALIAS Y F, MANZARI M T, AKAISHI M. A simple anisotropic clay plasticity model[J]. Mechanics Research Communications, 2002, 29(4):241.