剪脹角的案例
【FLAC系列】使用fish實現狀態相關摩爾本構庫倫
而狀態相關摩爾庫倫是眾多修正摩爾庫倫模型的一種，能夠模擬出摩擦角和剪脹角參數隨著狀態參數（孔隙比與臨界孔隙比差值）的變化而變化的現象。
??本文將基于“Gao L, Guo N, Yang Z X, et al. MPM modeling of pile installation in sand: Contact improvement and quantitative analysis[J]. Computers and Geotechnics, 2022, 151: 104943.”文章中使用的狀態相關摩爾庫倫本構模型，僅用Fish函數實現其二次開發。
??本文包括以下內容：1、介紹狀態相關摩爾庫倫（MC）；2、狀態相關摩爾庫倫的關鍵方程組；3、楓丹白露砂的狀態相關摩爾庫倫的標定；4、基于fish嵌入FLAC的上述狀態相關MC開發。
1、狀態相關摩爾庫倫簡介
??狀態相關MC采納了臨界狀態的概念，認為砂土受剪切達到臨界狀態時，處于一種“流動狀態”，即剪脹角為0。而我們知道，常規的MC模型擁有一個固定的剪脹角。另外一個方面在于，密砂實際上擁有峰值強度（對應峰值摩擦角）和殘余強度（對應于殘余摩擦角），而常規MC僅有一個摩擦角，無法模擬出從峰值強度到殘余強度的軟化過程，而狀態相關MC是可以的。
??因此，適用性來說，摩爾庫倫模型僅適用于小變形下的強度分析，而狀態相關MC可以模擬砂土在大變形下的力學行為。
2、狀態相關摩爾庫倫的關鍵方程組
??首先是臨界狀態線參數：
??其次是定義狀態變量：
??然后假設剪脹D和狀態變量線性相關：
??根據物理定義求剪脹角：
??同樣的步驟給出摩擦角的定義：
??以上具體內容可以看上面提到的文獻，有給出很詳細的參考文獻資料。
展開 基于Runge-Kutta算法的硬化土模型二次開發
土體處于塑性階段時，采用各向同性硬化法則與非關聯的流動法則反映土體的剪脹性，同時引入帽蓋屈服面反映土體的壓縮硬化，形成了雙硬化本構模型[5]。
圖1 硬化土模型彈性階段應力-應變關系
1.1 剪切屈服面
HS模型的剪切屈服面可用如下公式表示：
式中：Fs為剪切屈服函數；qa為極限偏應力；E50為對應50%強度時的割線模量；q為剪應力；Eur為卸載再加載模量；γ p為塑性剪切應變。
式中：E50ref為對應參考圍壓σref時的E50模量；σ3為第三主應力；c為黏聚力；φ為摩擦角；m為與土體性質有關的冪指數。
式中：Eurref為對應參考圍壓σref時的Eur模量。
式中：σ1為第一主應力；f為處于破壞時的狀態；Rf為破壞比。
1.2 塑性勢函數
HS模型剪切屈服面采用的是不相適應的流動規則，其剪切塑性勢函數如下：
式中：Qs為剪切塑性勢函數；ψm為機動剪脹角，由于不允許負剪脹角的存在，當ψm<0時，取0。
式中：φm為機動摩擦角；φcv為臨界摩擦角。
式中：φ為土體固有剪脹角。
1.3 壓縮屈服面
HS模型中為體現土體壓縮硬化特性加入了帽蓋屈服面，其屈服函數如下：
式中：F為壓縮屈服函數；δ為土體計算參數；σ2為第二主應力；M為摩擦常數；Pc為前期固結應力。
壓縮屈服面采用的是相適應的流動規則，其硬化定律如下：
式中：dp為硬化參數增量；Eoedref為切線模量；dεvp為塑性應變增量。
展開 Abaqus通過USDFLD子程序進行泥巖的應變軟化模擬
對于三軸壓縮試驗，η用塑性剪切應變來表示
塑性屈服準則采用Mohr-Coulomb準則，則粘聚力和內摩擦角的參數演化可以用下式表示
在巖石的塑性變形過程中會產生比較明顯的剪脹現象，而用來描述這一現象的較常用的力學參數就是剪脹角 Ψ，
上述模型可以通過USDFLD子程序進行實現，流程圖如下
有限元模型如下圖所示
計算得到的應力云圖及不同圍壓下的載荷位移響應如下所示
參考文獻：張力偉,賈善坡,鄒江濤,舒婧曦.泥巖的峰后軟化力學模型.中國科技論文,2016,11(21):2456-2461
有關于abaqus子程序開發的相關問題可以通過公眾號聯系我們.
公眾號: 320科技工作室
展開 Abaqus通過USDFLD子程序進行泥巖的應變軟化模擬
對于三軸壓縮試驗，η用塑性剪切應變來表示
塑性屈服準則采用Mohr-Coulomb準則，則粘聚力和內摩擦角的參數演化可以用下式表示
在巖石的塑性變形過程中會產生比較明顯的剪脹現象，而用來描述這一現象的較常用的力學參數就是剪脹角 Ψ，
上述模型可以通過USDFLD子程序進行實現，流程圖如下
有限元模型如下圖所示
計算得到的應力云圖及不同圍壓下的載荷位移響應如下所示
參考文獻：張力偉,賈善坡,鄒江濤,舒婧曦.泥巖的峰后軟化力學模型.中國科技論文,2016,11(21):2456-2461
有關于abaqus子程序開發的相關問題可以聯系扣扣1653004885或者關注cae320公眾號
展開 
懸臂式基坑開挖模擬
土體的彈性模量為66000kPa，泊比=0.2，土體重度=20kN/m3，黏聚力=0，摩擦角=30°，剪脹角=0，水平壓力系數K=2。墻體彈性模量E=28GPa，泊松比v=0.15。墻與土之間的摩擦角取為30°(摩擦系數為0577)。
圖1 模型示意圖
2. 模型的建立與求解
建立100mX100m的二維變形體部件，分隔開挖區域的幾何形狀。
圖2 建立土體表面
定義名為geo的Geostatic分析步，接受默認選項；隨后建立名為remove的Static，general通用分析步，起始時間0.1，最大步長0.2，采用非對稱算法。
載荷、邊界條件應在initial或geo分析步就激活生效。
在geo分析步對土體加體力-20，模擬重力荷載。
圖3 設置初始應力
網格單元類型為CPE4，選擇合適的網格密度。
圖4 有限元網格
處理器：i5-8700
內存：16G
計算時間：30s
軟件版本：ABAQUS 2020
3.
展開 有限元模擬三軸固結排水試驗，abaqus 初始應力場設置（一）
有限元模擬三軸固結排水試驗
模型概況
土體試樣尺寸：高 8 cm，直徑 4 cm；
土體力學參數：彈性模量 10MPa，泊松比 0.3，粘聚力 10 kPa，內摩擦角 30°;
試驗荷載：圍壓 100kPa；
試驗類型：等應變式三軸試驗，豎向應變為 10%；
模擬的目標
1、等壓固結完成時的應力狀態
2、獲得三軸試驗剪切破壞時的豎向應力
模型注意事項
1、簡化為軸對稱問題
2、彈性階段采用線彈性本構模型，塑性階段采用莫爾-庫倫本構模型
3、將固結完成后的應力狀態作為初始狀態
4、不考慮等壓固結的變形
5、采用 abaqus 的 Geostatic 分析步模擬等壓固結完成后的應力狀態
6、采用軸對稱應力單元 CAX4 ，只劃分一個單元
7、剪脹角采用 abaqus 默認的最小值 0.1°
有限元模型
注：斜體樣式只劃分一個單元，單元類型 ：4節點線性軸對稱應力單元
豎向應力與豎向應變關系
得到土體試樣剪切破壞時的豎向應力為 334.6kPa，與理論計算結果一致。
土體試樣的初始應力場設置
初始應力的設置需要滿足平衡條件：等效節點荷載要和外部荷載、邊界條件平衡。如果達不到平衡，將不能得到一個位移為零的初始狀態。此時所產生的應力場也不是所施加的初始應力場。
在本例中，等壓固結完成后的應力場為：三個方向的主應力都為 100kPa。在初始步設置初始應力如下：
在 Geostatic 分析步定義邊界條件為：對稱軸處 X 方向位移為零，底部 Y 方向位移為零。在頂面和右側施加圍壓 100kPa。得到的初始應力場如下：
對應的土體試樣位移云圖如下，可以判斷 Geostatic 分析步未產生位移：
展開 專業論文 | 土壓力作用下的隧道受力ABAQUS建模分析
創建材料，命名為soil，設置密度為1750，彈性模量隨溫度變化為
Temp=1 E=3E7 V=0.32
Temp=1 E=1.5E7 V=0.32
摩爾庫倫布模型：摩擦角30剪脹角0.1
3.裝配部件（Assembly）
點擊創建整體，選取上面創建的兩種部件tunnel和soil，由于坐標完備，可直接裝備完全。
4.定義分析步
編輯命名分析步：Geostat,選擇geostatic，在增量中改初始和最大增量為0.1，0.1。其余保留默認，確定。
編輯命名分析步：Reducation stiffness,選擇，static，general在增量中改初始和最大增量為0.1，0.1，other欄目中選擇matrix storage選擇Unsymmetric。其余保留默認，確定。
5.定義接觸
設置接觸命名為Add shotcrete，分析步在shotcrete，綁定隧道和土體。
設置接觸命名為remove shotcrete，分析步在geostat，綁定隧道和土體。
6.定義荷載，邊界條件
邊界條件：設置為初始分析步，對土體左右兩側U1方向進行固定，對土體底部進行U2固定。
荷載條件：選擇geo分析步點擊gravity，選取整體，設置重力大小為-9.81。
7.
展開 螺栓連接的彈塑性變形分析 附線性隨動強化彈塑性理論基礎下載
非關聯流動:
– 對摩擦材料，通常需要非關聯流動法則 (在 Drucker-Prager 模型中， 剪脹角與內摩擦角不同)。
強化準則:
? 強化準則描述屈服面如何隨塑性變形的結果而變化 (大小、中心、 形狀)。
? 強化準則決定如果繼續加載或卸載, 材料將何時再次屈服。
– 這與呈現無硬化– 即屈服面保持固定的彈性-理想塑性材料完全不同。
? 等向強化 指屈服面在塑性流動期間均勻擴張。 ‘等向’ 一詞指屈服面的均勻擴張，和 ‘各向同性’ 屈服準則(即材料取向)不同。
等向強化適用于大應變、比例加載情況。不適與循環加載。工程數據模塊提供了雙線性和多線性等向強化彈塑性模型。
對線性隨動強化, 屈服面在塑性流動過程中進行剛體平移。
屈服后最初的各向同性塑性行為不再各向同性 (隨動強化是各向異性強化的一種形式)
彈性區等于 2 倍的初始屈服應力，這稱為包辛格效應。
展開 土壓力作用下的隧道受力ABAQUS建模分析
創建材料，命名為soil，設置密度為1750，彈性模量隨溫度變化為
Temp=1 E=3E7 V=0.32
Temp=1 E=1.5E7 V=0.32
摩爾庫倫布模型：摩擦角30 剪脹角0.1
保存
(3) 裝配部件（Assembly）
點擊創建整體，選取上面創建的兩種部件tunnel和soil，由于坐標完備，可直接裝備完全。
(4) 定義分析步
編輯命名分析步：Geostat,選擇geostatic，在增量中改初始和最大增量為0.1，0.1。其余保留默認，確定。
編輯命名分析步：Reducation stiffness,選擇，static，Genral在增量中改初始和最大增量為0.1，0.1，other欄目中選擇matrix storage 選擇Unsymmetric。其余保留默認，確定。
(5) 定義接觸
設置接觸命名為Add shotcrete，分析步在shotcrete，綁定隧道和土體。
設置接觸命名為remove shotcrete，分析步在geostat，綁定隧道和土體。
確認
(6) 定義荷載，邊界條件
邊界條件；設置為初始分析步，對土木左右兩側U1方向進行固定，對土體底部進行U2固定。
荷載條件：，選擇geo分析步點擊gravity，選取整體，設置重力大小為-9.81。
(7) 設置初始應力
如下所示
圖片 1.png
選擇整個隧道，如下設置
(8) 劃分網格
對隧道進行網格劃分
對土體進行網格劃分，
對靠近管道位置的土體進行細分
網格單元類型
(9) 提交任務
4.
展開 隧道開挖三維模擬-1
土體采用莫爾庫倫模型，凝聚力5kPa，摩擦角25度，剪脹角12.5度，地下水位于地表，有效重度10kN/m3，水平土壓力系數1.0，彈性模量隨深度向下線性增加，即E=4000（1+0.2*depth）。
襯砌混凝土采用彈性模型模擬，參數如圖5所示。
圖4
圖5
4、建模過程
（1）通過拉伸建立三維實體，partition出隧道位置（開挖土體位置），為劃分網格方便，可在隧道周圍分割出一個矩形，便于網格過渡（圖6）。為便于后期操作，將每步開挖的土體創建集合r1、r2等。
圖6
（2）在Property模塊中創建土體材料，注意選擇參數與溫度相關（圖7）。實現彈性模量的分布也可以通過參數隨場變量變化的功能實現，但是場變量的定義不能在cae中進行。方便起見，這里利用溫度分布（可以在cae中設置）來實現。圖7中的溫度0對應模型底面，60對應模型頂面。建立相應section并分配給土體。
圖7
（3）創建襯砌混凝土材料，并創建板shell截面屬性（圖8），輸入厚度0.168。
圖8
（4）創建蒙皮。在模型中利用createdisplay group（按鈕或tools菜單下）將開挖土體移除，執行special-skin-create將第一段開挖土體后余下土體的表面創建蒙皮（skin）（圖9）。
圖9
（5）執行Assignsection，將板截面屬性分配給蒙皮。注意將選擇對象改為skin，在屏幕底端順便建立集合s1（圖10）。依次建立所有的襯砌集合。將所有襯砌建為集合sall。
圖10
（6）拼裝實體后，進入step模塊，建立geostatic分析步，用于初始應力平衡（時間步長可選自動）。
展開 開源XFEM程序：PhiPsi介紹
由于擴展有限元法一般采用水平集函數
和
描述裂縫面，此外
和
還在塑性力學中用來表示內摩擦角和剪脹角，故取名PhiPsi。
