基于Hoek-Brown準則的應(yīng)變軟化模型有限元數(shù)值實現(xiàn)研究

金俊超，佘成學(xué)，尚朋陽

(武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，湖北，武漢 430072)

自1980年Hoek-Brown(H-B)準則被首次提出以來，由于其具備綜合考慮了巖塊和結(jié)構(gòu)面的強度以及巖體結(jié)構(gòu)等多種因素的影響，能夠反映巖體的非線性破壞特征，較好地解決了Mohr-Coulomb(M-C)、Drucker-Prager(D-P)強度準則受拉破壞處理上的困難，并已建立了一套比較完善的參數(shù)確定方法等特點，在巖體工程中得到了廣泛使用[1?5]。

然而，當前僅有少量有限元軟件中集成了H-B準則本構(gòu)模型，這限制了H-B準則在數(shù)值計算中的應(yīng)用。雖然也有學(xué)者[6?12]通過二次開發(fā)等方法，將H-B準則理想彈塑性模型嵌入ABAQUS、GeoFBA3D等有限元軟件中，但上述研究并未實現(xiàn)H-B準則應(yīng)變軟化問題的求解。盡管文獻[13]采用顯示差分算法，實現(xiàn)了H-B準則應(yīng)變軟化求解，但該方法并不適用于有限元計算。而雖然文獻[14]給出了考慮應(yīng)變軟化的H-B準則本構(gòu)積分，但文中對于具體的有限元軟化求解過程未給出明確說明，考慮到經(jīng)典塑性力學(xué)對于軟化速率是有限制的[15]，該方法可能無法求解峰后軟化速率較大的情況。因此，有必要進一步研究基于H-B準則的應(yīng)變軟化模型有限元數(shù)值實現(xiàn)方法。

關(guān)于應(yīng)變軟化問題有限元求解，文獻[16-17]指出可將應(yīng)變軟化過程轉(zhuǎn)化為一系列的應(yīng)力脆性跌落與塑性流動過程，軟化求解問題就轉(zhuǎn)化為一系列脆塑性求解問題，有效避免了經(jīng)典塑性力學(xué)對于軟化速率的限制[15]，并提出了基于最小主應(yīng)力不變跌落的M-C準則應(yīng)變軟化模型有限元實現(xiàn)方法。 但是上述研究主要存在下述兩點不足有待解決：1) 采用的脆塑性計算方法的合理性有待分析，最小主應(yīng)力不變跌落能否正確描述巖石的變形破壞，是否應(yīng)引入其他脆塑性計算方法來進行軟化求解，還需要針對不同類型破壞現(xiàn)象，進行分析論證；2) 未解決H-B準則應(yīng)變軟化模型的有限元求解問題。

鑒于上述存在的不足，本文首先分析當前不同脆塑性計算方法的合理性，確定合理的脆塑性計算方法。在此基礎(chǔ)上，推導(dǎo)給出H-B準則脆塑性及理想彈塑性的隱式本構(gòu)積分算法，并采用一系列脆性跌落-塑性流動，將H-B準則應(yīng)變軟化模型嵌入有限元軟件ABAQUS中。然后，通過與應(yīng)變軟化圓隧圍巖的力學(xué)響應(yīng)規(guī)律的解析解對比，檢驗所建模型的正確性。最后，將模型應(yīng)用于某薄上覆蓋巖層高內(nèi)水壓輸水隧洞工程，分析運行期圍巖結(jié)構(gòu)的穩(wěn)定性。

1 不同脆塑性計算方法合理性分析

針對脆塑性求解問題，國內(nèi)外諸多學(xué)者進行了研究，根據(jù)假設(shè)不同，有以下3種脆塑性計算方法：偏應(yīng)力等比例跌落[18-19]、最小主應(yīng)力不變跌落[16-17]及塑性位勢跌落[15,20]。

巖石的變形破壞可分為單軸拉伸破壞、拉剪破壞、純剪破壞及壓剪破壞，見圖1。基于M-C準則及H-B準則，分析上述3種方法對單軸拉伸破壞、單軸壓縮破壞及二向純剪破壞計算是否正確，若單軸壓縮破壞、單軸拉伸破壞及二向純剪破壞可正確計算，則拉剪、壓剪等復(fù)合破壞類型也可正確計算。

圖1 基于M-C準則和H-B準則的巖石破壞分類Fig.1 Rock failure types based on the M-C criterion and the H-B criterion

1.1 偏應(yīng)力等比例跌落不合理性論證

偏應(yīng)力等比例跌落假設(shè)應(yīng)力由峰值屈服面徑向跌落至殘余屈服面上，有且僅有各向應(yīng)力偏量按同一比例β衰減[18-19]。若設(shè)峰值強度面應(yīng)力張量為σp，則殘余強度面應(yīng)力張量為σr為：

1) 單軸拉伸破壞不可行論證

當巖石單軸拉伸破壞時，峰值強度面應(yīng)力張量σp為(應(yīng)力以拉為正)：

聯(lián)立式(1)和式(2)，可得到殘余強度面應(yīng)力張量σr為：

將σp和σr代入屈服函數(shù)，得到：

式中：上標M-C表示M-C準則屈服函數(shù)；上標H-B表示H-B準則屈服函數(shù)；下標p表示峰值強度面；下標r表示殘余強度面。

求解式(4a)和式(4b)可以發(fā)現(xiàn)，不管是采用M-C準則還是H-B準則，殘余強度系數(shù)β均不為1，即說明當巖石發(fā)生單軸拉伸破壞后，殘余階段不滿足單軸拉伸的應(yīng)力狀態(tài)，與試驗結(jié)果不符，證明偏應(yīng)力等比例跌落不能正確計算巖石單軸拉伸破壞。

2) 單軸壓縮破壞不可行論證

當巖石單軸壓縮破壞時，峰值強度面應(yīng)力張量σp為：

聯(lián)立式(1)和式(5)，得到殘余強度面應(yīng)力張量σr：

將σp和σr代入屈服函數(shù)，得到：

求解式(7a)和式(7b)可以發(fā)現(xiàn)，不管是采用M-C準則還是H-B準則，殘余強度系數(shù)β均不為1，即說明當巖石發(fā)生單軸壓縮破壞后，其殘余階段不滿足單軸壓縮的應(yīng)力狀態(tài)，與試驗結(jié)果不符，證明偏應(yīng)力等比例跌落不能正確計算巖石單軸拉伸破壞。

3) 二向純剪破壞可行論證

當巖石二向純剪破壞時，峰值強度面應(yīng)力張量σp為：

聯(lián)立式(1)和式(8)，得到殘余強度面應(yīng)力張量σr：

將σp和σr代入屈服函數(shù)，將和代入屈服函數(shù)，即可求解殘余強度系數(shù)β，限于篇幅，不再詳細展開。不管是采用M-C準則還是H-B準則，殘余階段均滿足二向純剪的應(yīng)力狀態(tài)，與試驗結(jié)果相符，證明偏應(yīng)力等比例跌落可正確計算巖石二向純剪破壞過程。

綜合上述分析可知，偏應(yīng)力等比例跌落存在不足，不能正確計算巖石單軸拉伸破壞和單軸壓縮破壞，進一步拉剪、壓剪等復(fù)合破壞類型也可能計算錯誤。

1.2 最小主應(yīng)力不變跌落不合理性論證

最小主應(yīng)力不變跌落針對以壓應(yīng)力為正的情況，假設(shè)應(yīng)力跌落過程中最小主應(yīng)力不變，應(yīng)力從峰值強度跌落至殘余強度時Lode參數(shù)不變[16―17]。對于以拉應(yīng)力為正的情況，有：

1) 單軸拉伸破壞不可行論證

對于單軸拉伸破壞，其峰值強度面應(yīng)力張量σp見式(2)，通過聯(lián)立式(2)和式(10)，得到殘余強度面應(yīng)力張量σr：

將σp和σr代入屈服函數(shù)，得到：

求解式(12a)和式(12b)發(fā)現(xiàn)無解，證明最小主應(yīng)力不變跌落不可正確計算巖石單軸拉伸破壞。

2) 單軸壓縮破壞可行論證

對于單軸壓縮破壞，其峰值強度面應(yīng)力張量σp見式(5)，通過聯(lián)立式(5)和式(10)，得到殘余強度面應(yīng)力張量σr：

將σp和σr代入屈服函數(shù)，即可求得量值，限于篇幅，不再詳細展開。不管是采用M-C準則還是H-B準則，殘余階段均滿足單軸壓縮的應(yīng)力狀態(tài)，與試驗結(jié)果相符，證明最小主應(yīng)力不變跌落可正確計算巖石單軸壓縮破壞。

3) 二向純剪破壞不可行論證

對于二向純剪破壞，其峰值強度面應(yīng)力張量σp見式(8)，通過聯(lián)立式(8)和式(10)，得到殘余強度面應(yīng)力張量σr：

將σp和σr代入屈服函數(shù)，得到：

求解式(15a)和式(15b)可以發(fā)現(xiàn)，不管是采用M-C準則還是H-B準則，均有說明當巖石發(fā)生二向純剪破壞后，殘余階段不滿足二向純剪的應(yīng)力狀態(tài)，明顯與試驗結(jié)果不符，證明最小主應(yīng)力不變跌落計算不能正確計算巖石二向純剪破壞。

綜合上述分析可知，最小主應(yīng)力不變跌落存在不足，不能正確計算巖石單軸拉伸破壞及二向純剪破壞，進一步拉剪、壓剪等復(fù)合破壞類型也可能計算錯誤。

1.3 塑性位勢跌落合理性論證

塑性位勢跌落假設(shè)脆性材料仍滿足Il’yushin公設(shè)，跌落過程產(chǎn)生的塑性應(yīng)變增量pΔε的方向仍然滿足塑性位勢理論[15,20]，為：

式中，Δλ為塑性乘子，塑性勢函數(shù)g取M-C準則。

1) 單軸拉伸破壞可行論證

對于單軸拉伸破壞，其峰值強度面應(yīng)力張量σp見式(2)，由于應(yīng)力位于π平面張拉棱線上，造成塑性位勢求導(dǎo)的數(shù)值奇異，為此本文將單軸拉伸破壞問題退化為一維問題進行分析，僅有拉伸方向存在峰值應(yīng)力此時塑性勢函數(shù)退化為：

聯(lián)立式(16)和式(17)可知，只有拉伸方向產(chǎn)生塑性應(yīng)變對應(yīng)的殘余應(yīng)力

式中，E為彈性模量。

2) 單軸壓縮破壞可行論證

對于單軸壓縮破壞，其峰值強度面應(yīng)力張量σp見式(5)，由于應(yīng)力位于π平面壓棱線上，造成塑性位勢求導(dǎo)的數(shù)值奇異，為此本文將單軸壓縮破壞問題退化為一維問題進行分析，僅有壓縮方向存在峰值應(yīng)力此時塑性勢函數(shù)退化為：

聯(lián)立式(16)和式(19)可知，只有壓縮方向產(chǎn)生塑性應(yīng)變對應(yīng)的殘余應(yīng)力

3) 二向純剪破壞可行論證

對于二向純剪破壞，其峰值強度面應(yīng)力張量σp見式(8)。假設(shè)巖石純剪破壞過程不產(chǎn)生塑性體應(yīng)變，即剪脹角為0°[9-10]，聯(lián)立式(8)和式(16)可知，只有1方向和3方向產(chǎn)生塑性應(yīng)變，且：

對應(yīng)的殘余強度面應(yīng)力張量σr：

式中，G為剪切模量。

將σp和σr代入屈服函數(shù)，即可求得量值，限于篇幅，不再詳細展開。不管是采用M-C準則還是H-B準則，殘余階段均滿足二向純剪的應(yīng)力狀態(tài)，與事實相符，證明塑性位勢跌落可正確計算二向純剪破壞。

綜合上述分析結(jié)果可知，塑性位勢跌落可正確計算巖石單軸拉伸破壞、單軸壓縮破壞及二向純剪破壞，至于拉剪、壓剪等復(fù)合破壞類型也可正確計算，而偏應(yīng)力等比例跌落和最小主應(yīng)力不變跌落均存在不足。

2 H-B準則應(yīng)變軟化模型有限元數(shù)值實現(xiàn)

以拉應(yīng)力為正，壓應(yīng)力為負，以應(yīng)力不變量表示的H-B屈服準則為：

式中：I1為第一應(yīng)力不變量；J2為第二偏應(yīng)力不變量；J3為第三偏應(yīng)力不變量；

通過建立強度參數(shù)隨塑性參數(shù)的演化方程，即可實現(xiàn)應(yīng)變軟化過程的模擬，具體的演化模型研究可參見文獻[21-23]。

考慮到通常采用剪脹角來衡量巖體的剪脹變形，以應(yīng)力不變量表示的M-C準則塑性勢函數(shù)為：

采用一系列脆性跌落-塑性流動，進行H-B應(yīng)變軟化模型的有限元求解，其中數(shù)值實現(xiàn)的關(guān)鍵在于H-B脆塑性及理想彈塑性本構(gòu)積分的實現(xiàn)。

盡管基于D-P、M-C等準則的塑性位勢跌落本構(gòu)積分算法已有相關(guān)研究[15,20]，但并未有學(xué)者給出基于H-B準則的塑性位勢跌落本構(gòu)積分算法；另外，雖然關(guān)于H-B準則理想彈塑性本構(gòu)積分算法也已有大量研究[6-12]，但其塑性函數(shù)通常采用H-B準則，少有采用M-C準則的。為此，2.1節(jié)推導(dǎo)給出基于M-C準則塑性勢函數(shù)的，H-B準則脆塑性及理想彈塑性的隱式積分算法。

2.1 H-B準則脆塑性隱式本構(gòu)積分算法

假定第n+1增量步開始時應(yīng)力為σn，應(yīng)變增量為Δεn，塑性應(yīng)變?yōu)橛嬎銖椥栽囂綉?yīng)力σtr：

其中，D為彈性剛度矩陣。

若應(yīng)力在峰值強度屈服面fp內(nèi)，即則σn+1=σtr；若應(yīng)力位于或超出峰值強度屈服面fp，即如圖2所示，則發(fā)生脆塑性應(yīng)力跌落，到殘余強度屈服面，初次回拉應(yīng)力σC為：

式中：下標B表示對σB的計算結(jié)果；由于發(fā)生脆塑性應(yīng)力跌落后，巖石進入理想彈塑性變形階段，取A=0。

圖2 隱式本構(gòu)積分算法示意Fig.2 The diagram of implicit integration algorithm

通常，初次回拉應(yīng)力σC不會在殘余強度屈服面上，為此采用隱式向后歐拉算法，通過多次迭代計算使σC準確位于殘余強度屈服面fr上。以下省略公式中下標C表示為當前構(gòu)型的應(yīng)力狀態(tài)，對式(26)進行微分，并化簡得到：

式中：I為單位矩陣；為相較初始Δεn的改變量；為相較初始σC的改變量；為相較Δλ的改變量。

根據(jù)屈服函數(shù)f的一致性條件，為了使當前應(yīng)力σ落在殘余強度屈服面上，應(yīng)使=0，即：

將式(29)代入式(28)，并化簡得到：

將式(30)代入式(28)中，即可得到一致性切線模量Dct：

通過多次回拉計算，最終滿足下述方程組，使σC準確位于殘余強度屈服面上：

在隱式應(yīng)力更新過程中，涉及到屈服函數(shù)f及塑性勢函數(shù)g對應(yīng)力σ的求導(dǎo)，具體計算公式可參見文末附錄。

2.2 H-B準則理想彈塑性隱式本構(gòu)積分算法

H-B準則理想彈塑性的隱式本構(gòu)積分算法，與脆塑性隱式本構(gòu)積分算法相近，也可按照式(25)～式(32)進行理想彈塑性應(yīng)力更新，區(qū)別僅在于式中屈服函數(shù)取峰值強度面屈服函數(shù)，不再進行贅述。

采用FORTRAN語言，將上述H-B準則脆塑性及理想彈塑性隱式本構(gòu)積分算法編入有限元軟件ABAQUS中，并按照一系列應(yīng)力跌落-塑性流動，實現(xiàn)H-B準則應(yīng)變軟化模型求解。

3 數(shù)值算例驗證

H-B準則常見的應(yīng)用之一是判別山嶺隧道開挖后圍巖的穩(wěn)定性，因此其數(shù)值實現(xiàn)的驗證常基于如圖3(a)所示的一類經(jīng)典彈塑性力學(xué)問題，即計算處于平面應(yīng)變狀態(tài)的應(yīng)變軟化圓隧圍巖的力學(xué)響應(yīng)規(guī)律，并與解析解對比，有限元模型見圖3(b)。

圖3 圍巖力學(xué)響應(yīng)理論模型及有限元模型Fig.3 Example of a circular tunnel excavation and finite element model

基于H-B準則，根據(jù)表1中計算參數(shù)(上標p表示峰值強度參數(shù)，上標r表示殘余強度參數(shù))，遵循強度參數(shù)和剪脹角隨塑性剪應(yīng)變η分段線性演化規(guī)律，Lee等[24]提出將塑性區(qū)按相等的應(yīng)力增量劃分，圍壓逐漸遞減的差分方法；Wang等[25]及Zhang等[26]提出了將軟化區(qū)劃分成若干個彈脆性模型，從彈塑性交界面逐步計算至洞壁，給出了應(yīng)變軟化圓隧圍巖力學(xué)響應(yīng)規(guī)律解析解，見圖4。為了檢驗?zāi)Ｐ偷暮侠硇裕瑢⒈疚耐ㄟ^ABAQUS有限元計算的圓隧收斂位移及應(yīng)力分布也繪于圖4。

表1 文獻[24-26]中應(yīng)變軟化圓隧圍巖計算參數(shù)Table 1 Calculation parameters of tunnel rock mass material exhibiting strain softening behavior in reference [24-26]

觀察圖4(a)可以發(fā)現(xiàn)，本文計算的圍巖特征曲線與Lee等[24]、Wang等[25]及Zhang等[26]解析解十分接近；圖4(b)中，本文計算的pi=0時的圍巖徑向及切向應(yīng)力分布與Lee等[24]及Wang等[25]的解析解也具有較好的可比性，說明本文所述的數(shù)值實現(xiàn)工作在整體模型層次是正確有效的，擴展了ABAQUS功能。

圖4(b)中切向應(yīng)力分布解析解與本文有限元解存在一定的差異，一方面，可能是因為解析解在計算過程中，每一個微元的切向應(yīng)力取跌落前和跌落后的平均值，而本文采用有限元計算結(jié)果；另一方面，也可能與采用的脆塑性計算步數(shù)有關(guān)，出于計算效率考慮，本文在計算過程中，設(shè)置了100個脆性跌落-塑性流動計算步，脆塑性計算步越多，峰后模擬曲線與試驗結(jié)果越接近。

圖4 不同方法計算結(jié)果對比Fig.4 Comparisons of calculation results by different methods

4 工程應(yīng)用

4.1 工程概況

某盾構(gòu)輸水隧洞工程采用雙線形式布置，隧洞設(shè)計斷面為圓形，盾構(gòu)外徑6 m，兩洞室平行布置，中心間距12 m。混凝土管片襯砌內(nèi)徑為5.4 m，厚0.3 m，寬1.5 m。

選取埋深42.5 m的某一洞段作為研究對象，該洞段上覆弱風(fēng)化巖層厚3 m，上部為土層。根據(jù)規(guī)范[27]，進行運行期承載能力極限狀態(tài)分析，內(nèi)水壓力水頭為196 m，外水壓力水頭為18 m，地應(yīng)力按自重應(yīng)力場考慮。

4.2 計算參數(shù)及有限元模型

不考慮襯砌透水，土體及管片襯砌采用線彈性模型，圍巖采用所建H-B準則應(yīng)變軟化模型及理想彈塑性模型，管片與管片之間、管片與圍巖之間接縫采用庫倫摩擦模型，接縫面之間可以傳遞壓應(yīng)力和小于臨界切向力的剪應(yīng)力，不能傳遞拉應(yīng)力。具體力學(xué)參數(shù)見表2，假定圍巖強度參數(shù)及剪脹角隨塑性剪應(yīng)變分段線性變化，臨界塑性剪應(yīng)變?yōu)?.2×10?3。

表2 圍巖、管片襯砌及土層力學(xué)參數(shù)Table 2 Calculation parameters of surrounding rock, segment lining and soil

有限元模型沿洞軸向取4環(huán)管片長度，四周邊界均取3倍洞徑，頂部超出部分土層作為均布壓應(yīng)力施加在上表面。與洞軸向垂直、平行的邊界面均沿其面法向施加水平鏈桿約束，底面施加垂直鏈桿約束，有限元模型網(wǎng)格見圖5。

4.3 計算結(jié)果及分析

計算步驟為：① 建立初始地應(yīng)力平衡；② 左、右隧洞同時開挖；③ 施作管片襯砌，施加內(nèi)水壓力及外水壓力。

隧洞開挖完成后未出現(xiàn)塑性變形，由圍巖承受全部開挖應(yīng)力釋放荷載。施加管片襯砌，內(nèi)、外水壓力分別作用在管片內(nèi)外側(cè)。鑒于有限元網(wǎng)格及荷載對稱，僅以右洞為例，分析圍巖的變形破壞。

1) 圍巖拉應(yīng)力對比分析

觀察圖6所示圍巖拉應(yīng)力對比可以發(fā)現(xiàn)，采用應(yīng)變軟化模型計算的隧洞頂部及底部圍巖拉應(yīng)力，明顯低于采用理想彈塑性模型的計算結(jié)果，這主要是因為在考慮應(yīng)變軟化情況下，隨著圍巖塑性屈服，其承載能力逐漸下降，應(yīng)力發(fā)生轉(zhuǎn)移。

圖5 有限元模型網(wǎng)格Fig.5 The finite element model

圖6 圍巖拉應(yīng)力對比Fig.6 Tensile stress distribution of surrounding rocks

2) 圍巖塑性區(qū)對比分析

圖7中采用應(yīng)變軟化模型和采用理想彈塑性模型計算的圍巖塑性開裂區(qū)范圍差異明顯，采用理想彈塑性模型計算的頂部圍巖塑性區(qū)深度為1.6 m，未貫穿上覆蓋巖層，認為輸水隧洞是安全穩(wěn)定的；而采用應(yīng)變軟化模型計算的頂部圍巖塑性區(qū)貫穿上覆蓋巖層，整體結(jié)構(gòu)接近失穩(wěn)破壞。顯然，以理想塑性計算的結(jié)果是偏危險的，而應(yīng)考慮巖石峰后塑性軟化效應(yīng)。

圖7 圍巖塑性區(qū)對比Fig.7 Plastic zone distribution of surrounding rocks

從考慮峰后軟化效應(yīng)的計算結(jié)果來看，在給定內(nèi)水壓力作用下，即使圍巖與管片襯砌共同工作，由于上覆巖層過薄，難于滿足安全要求。處理的方法可考慮在管片襯砌內(nèi)側(cè)增加鋼筋混凝土襯砌，或增加鋼板襯砌。

5 結(jié)論

(1) 對當前不同脆塑性計算方法的合理性進行分析，發(fā)現(xiàn)塑性位勢跌落可正確計算巖石不同類型破壞，而偏應(yīng)力等比例跌落和最小主應(yīng)力不變跌落均存在不足。

(2) 推導(dǎo)給出基于塑性位勢跌落的H-B準則脆塑性隱式積分算法，及H-B準則理想彈塑性隱式積分算法，并采用一系列脆性跌落與塑性流動，將H-B準則應(yīng)變軟化模型嵌入有限元軟件ABAQUS中。

(3) 比較應(yīng)變軟化圓隧圍巖收斂位移及應(yīng)力分布的解析解與本文有限元解，發(fā)現(xiàn)二者吻合良好，說明本文所述的數(shù)值實現(xiàn)工作在整體模型層次是正確有效的，擴展了ABAQUS功能。

(4) 對某薄上覆蓋層高內(nèi)水壓輸水隧洞的計算結(jié)果表明，相較理想彈塑性模型，所建應(yīng)變軟化模型可更合理地描述圍巖塑性破壞及應(yīng)力調(diào)整過程，反映隧洞頂部圍巖塑性區(qū)貫通引起的整體結(jié)構(gòu)失穩(wěn)破壞現(xiàn)象，為工程選擇襯砌方案提供依據(jù)。

附錄：

選取σ1＞σ2＞σ3區(qū)間進行分析，當應(yīng)力位于π平面棱線或頂點時，可進行單獨處理，具體可參考文獻[6, 14]。由于脆塑性隱式本構(gòu)積分算法與理想彈塑性本構(gòu)積分算法，屈服函數(shù)f及塑性勢函數(shù)g的求導(dǎo)相同，故統(tǒng)一表示如下：

1) 屈服函數(shù)f一次求導(dǎo)的表達式如下：

其中：

2) 塑性勢函數(shù)g一次求導(dǎo)的表達式如下：

3) 塑性勢函數(shù)g二次求導(dǎo)的表達式如下：

其中：