1 引言

由于巖土材料問題的物理不穩定性，非線性材料的路徑依賴性以及應力-應變響應的非線性特征，選擇合適的本構模型會變得非常困難。《本構模型(Constitutive Models)選擇》《IMASS---FLAC3D和3DEC新的本構模型(1)》《FLAC3D 7.0 新特性簡介(P3)---新的本構模型》《FLAC2D---過去，現在和將來》。迄今為止，FLAC3D/3DEC已經內置了超過35種本構模型，當模擬一個問題時，我們不可能試驗所有的本構模型。大多數情況下，總是從最簡單的各向同性的彈性模型(Isotropic Elastic Model)入手。彈性模型(block zone cmodel assign elastic)只需要兩個材料參數：體積模量和剪切模量(block zone prop  dens=0.0026 bulk=4000 shear=3000)，運行速度最快。在解決全尺寸的邊值問題之前使用彈性模型測試能夠對模擬的問題呈現出一個big picture，粗略判斷出應力集中的位置，了解模型的預期響應，有助于重新定義網格的密度以及改變材料的本構模型。

對于一般的彈塑性問題，可以直接使用Mohr-Coulomb模型進行測試。Mohr-Coulomb模型是塑性模型組(Plastic Model Group)里最簡單的模型。輸入參數除了體積模量和剪切模量外，只需提供材料的密度，粘結力，內摩擦角和抗拉強度即可。

block zone property density 2.5E3 bulk 1.19E10 shear 1.1E10 friction 44 cohesion 2.72E5 tension 2E5

下面的例子演示了一個由Mohr-Coulomb材料組成的壓縮試驗。這個問題的物理意義是模擬一個礦柱的屈服行為。

2 模擬步驟

(1) 幾何形狀。幾何形狀必須在一定程度上代表著真實的物理問題，主要考慮以下幾個方面：[1] 模型應加入多少細節(節理，斷層，層面等)來表示地質結構，或者是否使用離散斷裂網絡DFN表示隨機的巖體結構；[2] 模型邊界的位置將如何影響模型結果？邊界設置的太小，計算結果不準確，而邊界設置的太大，計算時間將大大增加；[3] 如果使用可變形塊體，在感興趣的區域內，需要多大的網格密度才能得到精確的解答？這三個方面決定了合適模型的尺寸。block create 可以直接創建6種基元，它們分別是：polyhedron, brick, drum, prism, tunnel, tetrahedron, 在本例中使用drum產生一個圓柱體。

block create drum center-1 (0,0,0) center-2 (0,0,2) ...                     radius-1 1 radius-2 1 edges 16

(2) 劃分網格。網格劃分是一門藝術，沒有萬能的準則，在許多情況下可能需要不斷調節尺寸以獲得最優的結果。一些有限元軟件，例如Rocscience RS2和MIDAS GTS NX(巖土工程有限元分析軟件MIDAS GTS NX 2021 V1.1)中，確實內置了一些網格劃分準則幫助用戶進行網格劃分，然而在FLAC3D/3DEC中沒有這樣顯式的準則，不過提供了非常多的選項供用戶自行選擇，在本問題中，使用最大邊長max-edge作為單元劃分準則。

block zone generate-new max-edge 0.25

(3) 本構模型。如引言中所述，本構模型使用Mohr-Coulomb模型，同時設置該模型相應的材料參數。

block zone cmodel assign mohr-coulombblock zone property dens 2500 bulk 1.19e10 shear 1.1e10block zone property cohesion 2.72e5 friction 44 tension 2e5

(4) 邊界條件。為了模擬一個礦柱的受力狀態，在模型的頂部和底部施加數值相等的速度邊界，類似于單軸抗壓強度試驗。

block gridpoint apply velocity (0, 0, 1e-3) range position-z 0 block gridpoint apply velocity (0, 0, -1e-3) range position-z 2

(5) 監測點布置。布置3個監測點。第一個點監測模型底部中心點的垂直位移，第二個點檢測模型內部中點的垂直應力，第三個點檢測模型側表面中點的垂直應力。

block history displacement-z position (0,0,0)block history stress-zz position (0,0,1)block history stress-zz position (1,0,1)

3 模擬結果

下圖所示的是在運行一定時步后(model step 20000)的位移矢量圖和位移等值線圖。

下圖所示的是兩個垂直應力監測點隨礦柱頂板位移增加的變化。可以看出，礦柱內部的應力隨著位移的增加而增加，而礦柱表面的應力增加不是非常顯著。

4 塑性結果

四面體單元計算出的塑性結果不如六面體單元或高階單元計算的塑性結果準確。如果需要非常精確的塑性計算，應該使用六面體單元。由于3DEC是四面體單元，因此為了得到精確的塑性計算結果，需要啟用節點混合離散化(Nodal Mixed Discretization, NMD)算法。NMD算法在大多數情況下能夠產生更精確的塑性計算結果，因此為四面體網格提供了一個更精確的解決方案。默認狀態下，這個命令是關閉的(off), 需要使用下述命令打開:

block zone nodal-mixed-discretization on

比較結果顯示，NMD on對最終位移的影響不大，但對應力有顯著影響。NMD on使得在相同位移情況下塑性應力的值增大，因而這個NMD on導致計算結果和設計更保守一些。