本發(fā)明涉及巖土工程計(jì)算機(jī)輔助設(shè)計(jì)技術(shù)領(lǐng)域，特別涉及一種土體正交異性本構(gòu)模型的數(shù)值分析方法。

背景技術(shù)：

土作為一種非連續(xù)摩擦型散粒體工程材料，除表現(xiàn)為非線性非彈性、壓硬性、剪脹性、應(yīng)力-應(yīng)變與應(yīng)力歷史和應(yīng)力路徑相關(guān)性等諸多特性外，在工程實(shí)踐中，還特別表現(xiàn)出土體的各向異性。土之所以表現(xiàn)出各向異性是由于兩種原因：一種是固有各向異性，主要是天然土顆粒在沉積過(guò)程中或人工土在填筑過(guò)程中水平方向以及垂直方向的排列不同而引起的力學(xué)性狀和參數(shù)的不同；另外一種是應(yīng)力誘發(fā)各向異性，主要是由于土體在各個(gè)方向上所受的應(yīng)力不同，導(dǎo)致土顆粒在空間排列上不同，因而形成各向異性。由于土的成層特性，可以把土體視為正交異性，即平行于沉積面的各個(gè)方向的力學(xué)特性大致相同，而在垂直于沉積面的方向則表現(xiàn)出不同的特性。

已有研究表明，土體的各向異性對(duì)其破壞特性有顯著影響，例如，在研究模型條形基礎(chǔ)的承載力時(shí)，在加載方向和沉積面垂直時(shí)得到的最大承載力和在加載方向與沉積面平行時(shí)得到的最小承載力差別可以達(dá)到34％。對(duì)于黏土，強(qiáng)度的各向異性同樣很明顯。因此，對(duì)土體正交異性的合理描述具有重要的工程意義?，F(xiàn)有的本構(gòu)模型多把土體看作連續(xù)介質(zhì)，以經(jīng)典彈性、彈塑性理論為理論基礎(chǔ)，一些大型通用有限元軟件，如ANASYS和ABAQUS等，一般只有各向同性彈性、各向同性彈塑性各向異性彈性和各向異性彈塑性等簡(jiǎn)單的本構(gòu)關(guān)系模型，無(wú)法準(zhǔn)確反映考慮土體正交異性本構(gòu)條件下的應(yīng)力-應(yīng)變特性。

技術(shù)實(shí)現(xiàn)要素：

發(fā)明目的：本發(fā)明針對(duì)現(xiàn)有研究的不足，提供了一種土體正交異性本構(gòu)模型的數(shù)值分析方法。

技術(shù)方案：本發(fā)明提供了一種土體正交異性本構(gòu)模型的數(shù)值分析方法，其特征在于：(1)土體是正交異性材料，其本構(gòu)模型是基于摩爾-庫(kù)倫屈服準(zhǔn)則建立的；(2)土體正交異性本構(gòu)模型包括彈性正交異性和塑性正交異性兩個(gè)方面，并通過(guò)雙曲線擬合對(duì)摩爾-庫(kù)倫屈服方程進(jìn)行了修正；(3)土體的正交異性本構(gòu)模型是通過(guò)編譯用戶材料子程序UMAT，并嵌入到ABAQUS軟件中，最終實(shí)現(xiàn)土體正交異性本構(gòu)模擬的數(shù)值計(jì)算。

所述特征(1)具體為：土體為正交異性材料，即平行于沉積面的各個(gè)方向的力學(xué)特性大致相同，而在垂直于沉積面的方向則表現(xiàn)出不同的特性。本發(fā)明的土體正交異性本構(gòu)模型是基于摩爾-庫(kù)倫屈服準(zhǔn)則而建立的；

所述特征(2)具體為：土體正交異性模型同時(shí)考慮了彈性正交異性和塑性正交異性。彈性正交異性模型包含5個(gè)未知量，分別為沉積方向的彈性模量E_v和泊松比υ_vh、垂直于沉積方向的彈性模量E_h和泊松比υ_h、剪切模量G_vh，其他參數(shù)可由這5個(gè)未知量之間的相互關(guān)系得到。塑性本構(gòu)模型中假定土體內(nèi)摩擦角為常數(shù)，考慮土體粘聚力c存在各向異性效應(yīng)，且沉積面和垂直于沉積面上的粘聚力c_v、c_h存在關(guān)聯(lián)性，即：c_θ＝c_h+(c_v-c_h)sin²α。由于摩爾-庫(kù)倫屈服準(zhǔn)則在π平面存在六個(gè)棱角(奇異點(diǎn))，使得計(jì)算變繁和收斂緩慢，本發(fā)明采用雙曲線擬合對(duì)摩爾-庫(kù)倫屈服面進(jìn)行修正，在棱角處得到光滑連續(xù)的處理；

所述特征(3)具體為：取塑性勢(shì)函數(shù)與屈服函數(shù)的表達(dá)式一致，采用隱式后退歐拉算法，對(duì)摩爾-庫(kù)倫屈服函數(shù)一次求導(dǎo)、勢(shì)函數(shù)一次求導(dǎo)、勢(shì)函數(shù)二次求導(dǎo)，通過(guò)UMAT接口，將編寫有土體正交異性本構(gòu)模型的UMAT子程序嵌入到ABAQUS軟件中實(shí)現(xiàn)土體正交異性本構(gòu)模擬的數(shù)值計(jì)算。

有益效果：本發(fā)明克服了ABAQUS自帶的摩爾-庫(kù)倫屈服準(zhǔn)則無(wú)法同時(shí)考慮材料彈性正交異性和塑性各向異性的缺陷，通過(guò)編譯UMAT子程序?qū)⑼馏w正交異性的本構(gòu)模型嵌入到ABAQUS軟件中，實(shí)現(xiàn)了土體正交異性本構(gòu)模擬的數(shù)值計(jì)算。本發(fā)明所提供的數(shù)值分析方法能夠全面反映土體沉積所引起的正交異性的特點(diǎn)，為合理描述土體正交異性的特性提供了參考，也對(duì)指導(dǎo)工程實(shí)際問(wèn)題具有重要的科學(xué)價(jià)值。

附圖說(shuō)明:

圖1為本發(fā)明的數(shù)值方法流程圖；

圖2為土體粘聚力隨方向的各向異性；

圖3為用雙曲線近似擬合摩爾庫(kù)倫屈服函數(shù)；

圖4為p-q平面摩爾庫(kù)倫準(zhǔn)則的雙曲線逼近。

圖中，c_v、c_h分別為沉積面和垂直于沉積面上土體的粘聚力，c_θ表示土體在一個(gè)與沉積面成θ角的平面上的黏聚力；p為等效壓應(yīng)力；q為Mises等效應(yīng)力。

具體實(shí)施方式

以下結(jié)合附圖詳細(xì)敘述本發(fā)明的具體實(shí)施方式。本發(fā)明的保護(hù)范圍并不僅僅局限于本實(shí)施方式的描述。

一種土體正交異性本構(gòu)模型的數(shù)值分析方法，如圖1所示流程，具體步驟如下：

(1)根據(jù)問(wèn)題的實(shí)際情況，提取結(jié)構(gòu)幾何參數(shù)，建立ABAQUS有限元模型；

(2)在輸入文件中，使用關(guān)鍵詞“USER MATERIAL”，表示定義用戶材料屬性，并根據(jù)地質(zhì)勘測(cè)和室內(nèi)試驗(yàn)結(jié)果，將材料參數(shù)輸入到數(shù)值模型中；

(3)在ABAQUS的LOAD模塊中，對(duì)模型的邊界條件進(jìn)行設(shè)置，并施加初始荷載；

(4)選取合適的單元類型劃分網(wǎng)格，網(wǎng)格劃分前應(yīng)以劃分均勻合格的有限元網(wǎng)格為準(zhǔn)則對(duì)幾何模型進(jìn)行整理，并進(jìn)行網(wǎng)格尺寸敏感性分析，以確定同時(shí)滿足計(jì)算精度和計(jì)算效率的最優(yōu)網(wǎng)格劃分方式；

(5)編譯UMAT用戶子程序，并通過(guò)ABAQUS主程序接口嵌入到ABAQUS有限元模型中，提交模型并使其運(yùn)算完成。

所述的一種土體正交異性本構(gòu)模型的數(shù)值分析方法，其特征在于：(1)土體是正交異性材料，其本構(gòu)模型是基于摩爾-庫(kù)倫屈服準(zhǔn)則建立的；(2)土體正交異性本構(gòu)模型包括彈性正交異性和塑性正交異性兩個(gè)方面，并通過(guò)雙曲線擬合對(duì)摩爾-庫(kù)倫屈服方程進(jìn)行了修正；(3)土體的正交異性本構(gòu)模型是通過(guò)編譯用戶材料子程序UMAT，并嵌入到ABAQUS軟件中，最終實(shí)現(xiàn)土體正交異性本構(gòu)模擬的數(shù)值計(jì)算。

所述特征(1)具體為：土體為正交異性材料，即平行于沉積面的各個(gè)方向的力學(xué)特性大致相同，而在垂直于沉積面的方向則表現(xiàn)出不同的特性。摩爾-庫(kù)倫屈服準(zhǔn)則是目前巖土工程領(lǐng)域應(yīng)用最廣的強(qiáng)度準(zhǔn)則之一，并具有兩個(gè)其他屈服準(zhǔn)則無(wú)法比擬的優(yōu)點(diǎn)：第一，它的所有參數(shù)都具有直接的物理意義；第二，所有參數(shù)均能夠由常規(guī)的實(shí)驗(yàn)方法測(cè)得，因此摩爾庫(kù)倫屈服準(zhǔn)則的使用特別簡(jiǎn)單，本發(fā)明的土體正交異性本構(gòu)模型是基于摩爾-庫(kù)倫屈服準(zhǔn)則而建立的；

所述特征(2)具體為：土體正交異性模型同時(shí)考慮了彈性正交異性和塑性正交異性。在彈性階段，土體的應(yīng)力-應(yīng)變關(guān)系滿足正交異性本構(gòu)方程，如式(1)所示：

式中，E_v和E_h分別為沉積方向的垂直于沉積方向的彈性模量；υ_vh和υ_h分別為沉積方向的垂直于沉積方向的泊松比；G_h和G_vh分別為沉積方向的垂直于沉積方向的剪切模量。這些彈性參數(shù)之間存在一定的關(guān)系：

因此，彈性正交異性模型包含了5個(gè)未知量，其他的參數(shù)可由這5個(gè)未知量之間的相互關(guān)系得到。

在摩爾-庫(kù)倫屈服準(zhǔn)則中土的強(qiáng)度可以用c與兩個(gè)參數(shù)表示，研究表明，相對(duì)于內(nèi)摩擦角的各向異性，黏聚力c的各向異性效應(yīng)往往更為顯著，可以忽略內(nèi)摩擦角對(duì)于主應(yīng)力方向的依賴性，而認(rèn)為黏聚力依賴于主應(yīng)力方向。因此，這里假定內(nèi)摩擦角為常數(shù)，而考慮c的非均質(zhì)性與各向異性效應(yīng)。如圖2所示，正交異性土體在各個(gè)方向上的粘聚力表示為：

c_θ＝c_h+(c_v-c_h)sin²α (3)

式中，c_v、c_h分別為沉積面和垂直于沉積面上土體的粘聚力，c_θ表示土體在一個(gè)與沉積面成θ角的平面上的黏聚力。

若土體服從摩爾-庫(kù)倫屈服準(zhǔn)則，則其用應(yīng)力不變量表示的正交異性摩爾-庫(kù)倫屈服方程為：

在有限單元法分析中，由于摩爾庫(kù)倫屈服準(zhǔn)則在π平面上不是正六邊形，在拐角處的不連續(xù)往往會(huì)在數(shù)值計(jì)算中產(chǎn)生許多問(wèn)題。當(dāng)應(yīng)力狀態(tài)在這些奇異點(diǎn)處時(shí)，由于拐點(diǎn)的存在，屈服函數(shù)或勢(shì)函數(shù)存在不連續(xù)和梯度無(wú)法計(jì)算的情況，導(dǎo)致彈塑性應(yīng)力應(yīng)變矩陣無(wú)法求解。為了避免這種情況，本發(fā)明將以應(yīng)力不變量表示的摩爾庫(kù)倫準(zhǔn)則采用雙曲線近似模擬，如圖3所示。以圖中摩爾庫(kù)倫準(zhǔn)則線作為漸近線，畫出在p-q平面上相應(yīng)的曲線來(lái)修正摩爾-庫(kù)倫準(zhǔn)則，最終圓滑屈服面錐體頂點(diǎn)，消除奇異性與子午面拐點(diǎn)不可導(dǎo)的問(wèn)題，極大地提高了在有限元分析中彈塑性計(jì)算的穩(wěn)定性。圖中所對(duì)應(yīng)的參數(shù)l、d和h的物理意義分別為

修正后的正交異性摩爾庫(kù)倫屈服函數(shù)的表達(dá)式為：

由上式可看出，通過(guò)調(diào)整參數(shù)A可以更好的模擬摩爾庫(kù)倫屈服函數(shù)。在p-q面上，隨著參數(shù)A的改變，雙曲屈服函數(shù)與摩爾庫(kù)倫屈服函數(shù)的逼近關(guān)系如圖4所示。A的取值越小，雙曲屈服函數(shù)的近似效果越好，當(dāng)A＝1時(shí)，變回原摩爾庫(kù)倫屈服函數(shù)，不再具有消去奇異性的作用。對(duì)應(yīng)在子午面上，摩爾庫(kù)倫屈服準(zhǔn)則是六角錐體，錐體部分用雙曲屈服曲線逼近，使錐體變圓滑，讓程序在數(shù)值計(jì)算變得更加穩(wěn)定。

所述特征(3)具體為：取塑性勢(shì)函數(shù)與屈服函數(shù)的表達(dá)式一致，即

式中Ψ為土體剪脹角，K(θ)與屈服函數(shù)中的表達(dá)式相同，將摩擦角改寫為Ψ。

在有限元分析中，需要對(duì)勢(shì)函數(shù)和屈服函數(shù)求偏導(dǎo)，一是用于計(jì)算彈塑性應(yīng)力應(yīng)變關(guān)系，二是在應(yīng)力迭代算法中用于將應(yīng)力拉回屈服面。在UMAT子程序編寫中，采用隱式后退歐拉算法，對(duì)摩爾-庫(kù)倫屈服函數(shù)一次求導(dǎo)、勢(shì)函數(shù)一次求導(dǎo)、勢(shì)函數(shù)二次求導(dǎo)；通過(guò)UMAT接口，將編寫有土體正交異性本構(gòu)模型的UMAT子程序嵌入到ABAQUS軟件中實(shí)現(xiàn)土體正交異性本構(gòu)模擬的數(shù)值計(jì)算。