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

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

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

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

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

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

1 SANICLAY模型簡介

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

1.1 彈性關系

模型的彈性關系依賴于球應力p和孔隙比e，具體為

(1)

C=λ′δδ+2GI

(2)

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

1.2 塑性勢函數及流動法則

塑性勢函數

(3)

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

其中

χα為模型參數.將塑性勢函數ψ對應力σ求導可得塑性流動方向r=?ψ/?σ.

1.3 屈服函數

屈服函數

(4)

圖1 屈服函數與塑性勢函數

1.4 硬化規(guī)則

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

(5)

(6)

(7)

(8)

(9)

1.5 一致性條件

與屈服函數相應的一致性條件

(10)

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

(11)

其中，塑性模量

(12)

2 顯式算法數值實現

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

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

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

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

(13)

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

(14)

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

2.2 改進顯式算法

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

(15)

其中，

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

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

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

(16)

其中，上標(k-1)和(k)表示修正次數.

圖2 切平面算法原理

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

圖3 改進顯式算法框圖

3 不同算法對比

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

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

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

表1 計算結果對比

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

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

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

圖5 改進顯式算法1的計算結果

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

圖6 改進顯式算法2的計算結果

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

4 多單元分析

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

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

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

圖8 Mises應力云圖

5 結論

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

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

(3) 改進顯式算法可以有效對高度非線性彈塑性本構模型進行數值實現.

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

參考文獻:

[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] 孫德安, 陳波. 結構性軟土力學特性的試驗研究[J]. 土木工程學報, 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] 祝恩陽, 姚仰平.結構性土UH模型[J]. 巖土力學, 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.