本發(fā)明涉及一種離散元虛擬壓縮試驗(yàn)數(shù)值仿真方法，具體地說，涉及一種基于離散-連續(xù)耦合的離散元虛擬壓縮試驗(yàn)數(shù)值仿真方法。本發(fā)明屬于巖土工程材料虛擬數(shù)值試驗(yàn)仿真領(lǐng)域，實(shí)現(xiàn)對材料介質(zhì)室內(nèi)壓縮試驗(yàn)的計(jì)算機(jī)真實(shí)仿真模擬。

背景技術(shù)：

近些年來隨著計(jì)算機(jī)和數(shù)值計(jì)算方法的發(fā)展，虛擬數(shù)值試驗(yàn)技術(shù)在巖土力學(xué)研究領(lǐng)域得到了廣泛的應(yīng)用。虛擬數(shù)值試驗(yàn)作為室內(nèi)試驗(yàn)的一種替代或者輔助工具，其在巖土介質(zhì)材料復(fù)雜力學(xué)行為研究方面具有得天獨(dú)厚的優(yōu)勢。顆粒離散元虛擬數(shù)值試驗(yàn)作為一種常用的虛擬數(shù)值試驗(yàn)，特別適用于巖土材料這類非連續(xù)介質(zhì)的模擬，尤其是對材料介質(zhì)細(xì)觀力學(xué)機(jī)理的研究。

目前，顆粒離散元虛擬壓縮試驗(yàn)主要是通過無摩擦剛性墻體對顆粒試樣施加側(cè)向圍壓，這種圍壓加載方式屬于剛性加載，與室內(nèi)試驗(yàn)的柔性加載方式不同(即，通過試樣外層橡膠皮套來對試樣施加側(cè)向圍壓，試驗(yàn)中允許試樣側(cè)向自由變形)，故這種類型的虛擬壓縮試驗(yàn)無法實(shí)現(xiàn)對室內(nèi)壓縮試驗(yàn)的真實(shí)模擬。同時(shí)，由于虛擬數(shù)值試驗(yàn)是采用剛性加載方式來施加圍壓，試驗(yàn)過程中限制了試樣側(cè)向的自由變形，對試樣的力學(xué)響應(yīng)也會產(chǎn)生了一定的影響。

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

為解決常規(guī)離散元虛擬壓縮試驗(yàn)中無法施加側(cè)向柔性應(yīng)力邊界條件的難題，本發(fā)明的目的是提供一種基于離散-連續(xù)耦合的離散元虛擬壓縮試驗(yàn)數(shù)值仿真方法，該方法能夠在離散元虛擬壓縮試驗(yàn)中實(shí)現(xiàn)側(cè)向圍壓的柔性加載，達(dá)到對室內(nèi)壓縮試驗(yàn)的真實(shí)模擬。

為實(shí)現(xiàn)上述目的，本發(fā)明采用以下技術(shù)方案：一種基于離散-連續(xù)耦合的離散元虛擬壓縮試驗(yàn)數(shù)值仿真方法，包括如下步驟：

(1)根據(jù)壓縮試驗(yàn)試樣尺寸，建立常規(guī)離散元虛擬壓縮試驗(yàn)?zāi)Ｐ停?/p>

(2)利用應(yīng)力伺服控制原理，將試樣模型伺服到初始固結(jié)應(yīng)力狀態(tài)；

(3)利用有限元單元替代側(cè)向剛性墻體，并建立離散-連續(xù)耦合虛擬數(shù)值試驗(yàn)?zāi)Ｐ停?/p>

(4)基于離散-連續(xù)耦合計(jì)算原理，開展離散-連續(xù)耦合虛擬壓縮試驗(yàn)數(shù)值仿真模擬。

其中，步驟(1)的具體步驟包括：

(1.1)根據(jù)試樣尺寸，建立無摩擦剛性墻體模型；

(1.2)根據(jù)試樣顆粒級配要求，建立試樣顆粒模型。

步驟(2)的具體步驟包括：

(2.1)確定當(dāng)前時(shí)步側(cè)向墻體和豎向墻體的平均應(yīng)力值，計(jì)算公式如下：

式中，s為顆粒作用于側(cè)向或豎向兩墻體的平均應(yīng)力，f為顆粒作用于墻體上的作用力，l為當(dāng)前時(shí)步側(cè)向或豎向墻體間的相對距離(即，試樣的寬度或高度)，下標(biāo)x和y分別表示模型側(cè)向和豎向方向，上標(biāo)(i)中i表示墻體的編號，1代表下部墻體，2右側(cè)墻體，3表示上部墻體，4表示左側(cè)墻體；

(2.2)確定下一時(shí)步墻體的移動速度，計(jì)算公式如下：

式中，v為墻體移動速度，其下標(biāo)x和y與上標(biāo)(i)代表含義與步驟(2.1)相同，和為側(cè)向和豎向墻體初始固結(jié)應(yīng)力，g為增益系數(shù)，按如下公式計(jì)算：

式中，kⁿ為顆粒與墻體間的法向接觸剛度，n為顆粒與側(cè)向或豎向墻體的接觸數(shù)目，其下標(biāo)x和y代表含義均與步驟(2.1)相同，δt為當(dāng)前計(jì)算時(shí)步，α為松弛因子，取值范圍為(0，1]，通?？扇?.5；

(2.3)離散元循環(huán)迭代一個(gè)時(shí)步，更新墻體和顆粒的位置；

(2.4)重復(fù)執(zhí)行上述步驟(2.1)至(2.3)，直到當(dāng)前墻體應(yīng)力值與初始固結(jié)應(yīng)力值滿足以下條件：

式中，ε為容差，取值范圍為[0.001，0.005]。

步驟(3)的具體步驟包括：

(3.1)刪除左右側(cè)墻體，并在原左右側(cè)墻體位置處分別建立一列有限元單元；

(3.2)根據(jù)建立的有限元單元，沿左右列單元分別建立一組耦合控制墻體；

(3.3)確定和存儲每個(gè)控制墻體頂點(diǎn)與有限元單元節(jié)點(diǎn)的對應(yīng)關(guān)系。

步驟(4)的具體步驟包括：

(4.1)在有限元單元外側(cè)邊界施加試驗(yàn)側(cè)向圍壓；

(4.2)設(shè)置上下墻體的軸向試驗(yàn)加載速度；

(4.3)在動力求解模式下，將有限元和離散元設(shè)置為相同的分析時(shí)步；

(4.4)開始離散-連續(xù)耦合虛擬壓縮試驗(yàn)仿真模擬，每隔一定計(jì)算時(shí)步記錄試驗(yàn)的軸向偏應(yīng)力與軸向應(yīng)變；

(4.5)當(dāng)軸向應(yīng)變達(dá)到指定目標(biāo)值，停止試驗(yàn)加載。

更進(jìn)一步，步驟(4.4)中離散-連續(xù)耦合虛擬數(shù)值壓縮試驗(yàn)是基于有限元(fem)和離散元(dem)相互耦合計(jì)算實(shí)現(xiàn)的，其中fem是作為服務(wù)器(server)，dem作為客戶端(client)，兩者耦合同步進(jìn)行，在每一步耦合計(jì)算過程中，fem先計(jì)算循環(huán)一步，更新單元節(jié)點(diǎn)位置，并將耦合節(jié)點(diǎn)位置發(fā)送給dem，dem更新控制墻體的位置并循環(huán)計(jì)算一步，更新顆粒的位置，接著將控制墻體上的作用力發(fā)送給fem，通過作用力映射原理獲得耦合節(jié)點(diǎn)上的作用力；接著fem再循環(huán)計(jì)算一步，更新單元節(jié)點(diǎn)位置，將耦合節(jié)點(diǎn)速度再次發(fā)送給dem，dem再次更新控制墻體的位置并循環(huán)計(jì)算一步，更新離散域內(nèi)顆粒的位置，接著將控制墻體的作用力再次發(fā)送給fem；依次類推，耦合逐步計(jì)算下去；

更進(jìn)一步的，步驟(4.4)中作用力映射原理核心是將顆粒作用于控制墻體上力和力矩轉(zhuǎn)為對應(yīng)的耦合節(jié)點(diǎn)上等效節(jié)點(diǎn)力；

假設(shè)某個(gè)控制墻體上的作用力和力矩分別為(fx，fy，)和m，與該控制墻體頂點(diǎn)相對應(yīng)的兩個(gè)耦合節(jié)點(diǎn)為節(jié)點(diǎn)a和節(jié)點(diǎn)b，其坐標(biāo)分別為(xa，ya)和(xb，yb)，則對應(yīng)耦合節(jié)點(diǎn)上的等效節(jié)點(diǎn)力可由如下公式計(jì)算：

式中，(fxa，fya)和(fxb，fyb)分別為施加在節(jié)點(diǎn)a和節(jié)點(diǎn)b上的等效節(jié)點(diǎn)力，γ為一個(gè)系數(shù)，可由下式計(jì)算：

有益效果：與現(xiàn)有技術(shù)相比，本發(fā)明的一種基于離散-連續(xù)耦合的離散元虛擬壓縮試驗(yàn)數(shù)值仿真方法可以實(shí)現(xiàn)側(cè)向圍壓的柔性施加，能夠真實(shí)地模擬室內(nèi)壓縮試驗(yàn)。本發(fā)明方法原理簡單、數(shù)值計(jì)算效率高且穩(wěn)定，可以適用于分析巖土材料(如，土體、土石混合體、堆石料)宏細(xì)觀力學(xué)特性。

附圖說明

圖1是本發(fā)明基于離散-連續(xù)耦合的離散元虛擬壓縮試驗(yàn)數(shù)值仿真方法總的流程圖；

圖2是本發(fā)明實(shí)施例的虛擬試樣模型示意圖；

圖3是根據(jù)圖2建立的離散元虛擬壓縮試驗(yàn)?zāi)Ｐ蛨D；

圖4是圖3應(yīng)力伺服后的虛擬壓縮試驗(yàn)?zāi)Ｐ蛨D；

圖5是根據(jù)圖4建立的有限元單元模型圖；

圖6是沿圖5中有限元單元邊建立的控制墻體模型圖；

圖7是有限元單元耦合節(jié)點(diǎn)與墻體頂點(diǎn)對應(yīng)關(guān)系示意圖；

圖8是虛擬壓縮試驗(yàn)加載條件示意圖；

圖9是虛擬壓縮試驗(yàn)有限元與離散元耦合計(jì)算示意圖；

圖10是虛擬壓縮獲得的軸向應(yīng)變與軸向偏應(yīng)力曲線圖；

圖11是壓縮試驗(yàn)后的虛擬試驗(yàn)?zāi)Ｐ蛨D；

圖12a是壓縮試驗(yàn)后單元節(jié)點(diǎn)位移矢量圖；

圖12b是壓縮試驗(yàn)后試樣顆粒位移矢量圖。

具體實(shí)施方式

下面結(jié)合具體實(shí)施例，進(jìn)一步闡明本發(fā)明，應(yīng)理解這些實(shí)施例僅用于說明本發(fā)明而不用于限制本發(fā)明的范圍，在閱讀了本發(fā)明之后，本領(lǐng)域技術(shù)人員對本發(fā)明的各種等價(jià)形式的修改均落于本申請所附權(quán)利要求所限定的范圍。

如圖1所示，本發(fā)明公開了一種基于離散-連續(xù)耦合的離散元虛擬壓縮試驗(yàn)數(shù)值仿真方法，包括以下步驟：

s1：根據(jù)壓縮試驗(yàn)試樣尺寸，建立常規(guī)離散元虛擬壓縮試驗(yàn)?zāi)Ｐ汀?/p>

具體的步驟為：

s1.1：根據(jù)試樣尺寸，建立無摩擦剛性墻體模型；

圖2所示為虛擬試樣，其高為h，寬度為w。根據(jù)試樣模型圖建立壓縮試驗(yàn)墻體模型，如圖3所示，壓縮試驗(yàn)墻體模型是由四個(gè)無摩擦的剛性墻體組成，根據(jù)空間位置分別標(biāo)記為上下墻體和左右墻體，每個(gè)墻體均為一條線段，線段端點(diǎn)即為墻體的頂點(diǎn)，且上下墻體(豎向墻體)和左右墻體(側(cè)向墻體)間的相對距離即為試樣的寬度h和高度w。

s1.2：根據(jù)試樣顆粒級配要求，建立試樣顆粒模型；

利用離散元模型構(gòu)建方法(動力方法或幾何方法)在墻體模型圍成的試樣區(qū)域內(nèi)生成試樣顆粒模型。如圖3所示，建立的顆粒模型為一系列不同大小的圓形顆粒，且顆粒間相互接觸。

s2：利用應(yīng)力伺服控制原理，將試樣模型伺服到初始固結(jié)應(yīng)力狀態(tài)。具體的步驟為：

s2.1：根據(jù)公式(1)確定當(dāng)前時(shí)步側(cè)向墻體和豎向墻體的平均應(yīng)力值；

式中，s為顆粒作用于側(cè)向或豎向兩墻體的平均應(yīng)力，f為顆粒作用于墻體上的作用力，l為當(dāng)前時(shí)步側(cè)向或豎向墻體間的相對距離(即，試樣的寬度或高度)，下標(biāo)x和y分別表示模型側(cè)向和豎向方向，上標(biāo)(i)中i表示墻體的編號，1代表下部墻體，2右側(cè)墻體，3表示上部墻體，4表示左側(cè)墻體。

s2.2：根據(jù)公式(2)確定下一時(shí)步墻體的移動速度；

式中，v為墻體移動速度，其下標(biāo)x和y與上標(biāo)(i)代表含義與步驟(2.1)相同，和為側(cè)向和豎向墻體初始固結(jié)應(yīng)力，g為增益系數(shù)，按如下公式計(jì)算：

式中，kⁿ為顆粒與墻體間的法向接觸剛度，n為顆粒與側(cè)向或豎向墻體的接觸數(shù)目，其下標(biāo)x和y代表含義均與步驟(2.1)相同，δt為當(dāng)前計(jì)算時(shí)步，α為松弛因子，取值范圍為(0，1]，通?？扇?.5；

s2.3：離散元循環(huán)迭代一個(gè)時(shí)步，更新墻體和顆粒的位置；

s2.4：重復(fù)執(zhí)行上述步驟s2.1至步驟s2.3，直到當(dāng)前墻體應(yīng)力值與初始固結(jié)應(yīng)力值小于指定容差。

通常，當(dāng)前墻體應(yīng)力值與初始固結(jié)應(yīng)力值滿足以下條件：

式中，ε為指定容差，取值范圍為[0.001，0.005]；

在離散元中，應(yīng)力伺服的原理是通過反復(fù)調(diào)整墻體的位置使試樣內(nèi)部應(yīng)力逐步逼近目標(biāo)應(yīng)力值。如圖4所示，通過步驟s2將圖2中的顆粒試樣模型伺服到指定的固結(jié)應(yīng)力狀態(tài)，此過程相當(dāng)于室內(nèi)試驗(yàn)的固結(jié)過程；在應(yīng)力伺服后，試樣尺寸發(fā)生了變化，其中高度變?yōu)閔′，寬度變?yōu)閣′。

s3：利用有限元單元替代側(cè)向剛性墻體，建立離散-連續(xù)耦合虛擬試驗(yàn)數(shù)值模型。具體的步驟為：

s3.1：刪除左右側(cè)墻體，并在原左右側(cè)墻體位置處分別建立一列有限元單元；

如圖5所示，刪除圖4中左右側(cè)墻體，取而代之在左右側(cè)墻體位置建立兩列相同數(shù)目和大小有限元單元，每個(gè)有限元單元均為一個(gè)規(guī)則的四邊形單元。

s3.2：根據(jù)建立的有限元單元，沿左右列單元分別建立一組耦合控制墻體；

如圖6所示，沿圖5中左右列有限元單元邊分別建立一組耦合控制墻體；其中，左側(cè)控制墻體與左列有限元單元的右側(cè)邊重合，右側(cè)控制墻體與右列有限元單元的左側(cè)邊重合，故控制墻體頂點(diǎn)與相應(yīng)有限元單元頂點(diǎn)重合。

s3.3：確定和存儲每個(gè)控制墻體頂點(diǎn)與有限元單元節(jié)點(diǎn)的對應(yīng)關(guān)系。

如圖7所示，根據(jù)控制墻體與有限元單元邊的建立關(guān)系，確定和存儲每個(gè)控制墻體頂點(diǎn)與有限元單元節(jié)點(diǎn)的對于關(guān)系，用于后續(xù)步驟s4中耦合數(shù)據(jù)的傳遞與交換。

s4：基于離散-連續(xù)耦合計(jì)算原理，開展離散-連續(xù)耦合虛擬壓縮試驗(yàn)數(shù)值仿真模擬。具體步驟為：

s4.1：在有限元單元外側(cè)邊界施加試驗(yàn)側(cè)向圍壓；

如圖7所示，在左右列有限元單元外側(cè)施加側(cè)向圍壓，模擬室內(nèi)試驗(yàn)的側(cè)向圍壓。

s4.2：設(shè)置上下墻體的軸向試驗(yàn)加載速度；

如圖8所示，虛擬壓縮試驗(yàn)是通過加載速度方式進(jìn)行軸向加載的，即對上部墻體施加一個(gè)恒定的向下的移動速度，對下部墻體施加一個(gè)相同大小恒定的向上的移動速度；同時(shí)將有限元單元底部和頂部8個(gè)節(jié)點(diǎn)的速度設(shè)置為與相應(yīng)墻體相同的移動速度，即上部4個(gè)節(jié)點(diǎn)的速度與上部墻體的移動速度一致，下部4個(gè)節(jié)點(diǎn)的速度與下部墻體移動速度一致。

s4.3：在動力求解模式下，將有限元和離散元設(shè)置為相同的分析時(shí)步；

為保證有限元與離散元同步進(jìn)行，將有限元和離散元分析時(shí)步設(shè)置為同一時(shí)步大小。

s4.4：開始離散-連續(xù)耦合虛擬壓縮試驗(yàn)仿真模擬，每隔一定計(jì)算時(shí)步記錄試驗(yàn)的軸向偏應(yīng)力與軸向應(yīng)變；

本發(fā)明中，有限元與離散元耦合是同步進(jìn)行的，即有限元計(jì)算一個(gè)時(shí)步后，緊接著離散元計(jì)算一個(gè)時(shí)步，每一個(gè)計(jì)算時(shí)步中進(jìn)行數(shù)據(jù)間的相互傳遞和交換，整個(gè)耦合計(jì)算過程如圖9所示。本發(fā)明離散-連續(xù)耦合虛擬數(shù)值壓縮試驗(yàn)是基于有限元(fem)和離散元(dem)相互耦合計(jì)算實(shí)現(xiàn)的，其中fem是作為服務(wù)器(server)，dem作為客戶端(client)，兩者耦合同步進(jìn)行，在每一步耦合計(jì)算過程中，fem先計(jì)算循環(huán)一步，更新單元節(jié)點(diǎn)位置，并將耦合節(jié)點(diǎn)位置發(fā)送給dem，dem更新控制墻體的位置并循環(huán)計(jì)算一步，更新顆粒的位置，接著將控制墻體上的作用力發(fā)送給fem，通過作用力映射原理獲得耦合節(jié)點(diǎn)上的作用力；接著fem再循環(huán)計(jì)算一步，更新單元節(jié)點(diǎn)位置，將耦合節(jié)點(diǎn)速度再次發(fā)送給dem，dem再次更新控制墻體的位置并循環(huán)計(jì)算一步，更新離散域內(nèi)顆粒的位置，接著將控制墻體的作用力再次發(fā)送給fem；依次類推，耦合逐步計(jì)算下去。

在耦合逐步計(jì)算過程中，作用力映射原理核心是將顆粒作用于控制墻體上的力和力矩轉(zhuǎn)為對應(yīng)的耦合節(jié)點(diǎn)上等效節(jié)點(diǎn)力；

假設(shè)某個(gè)控制墻體上的作用力和力矩分別為(fx，fy)和m，與該控制墻體頂點(diǎn)相對應(yīng)的兩個(gè)耦合節(jié)點(diǎn)為節(jié)點(diǎn)a和節(jié)點(diǎn)b，其坐標(biāo)分別為(xa，ya)和(xb，yb)，則對應(yīng)耦合節(jié)點(diǎn)上的等效節(jié)點(diǎn)力可由如下公式計(jì)算：

式中，(fxa，fya)和(fxb，fyb)分別為施加在節(jié)點(diǎn)a和節(jié)點(diǎn)b上的等效節(jié)點(diǎn)力，γ為一個(gè)系數(shù)，可由下式計(jì)算：

s4.5：當(dāng)軸向應(yīng)變達(dá)到指定目標(biāo)值，停止試驗(yàn)加載。

下面以一個(gè)虛擬壓縮試驗(yàn)為例，說明本發(fā)明虛擬壓縮試驗(yàn)的具體實(shí)施過程：

1、根據(jù)壓縮試驗(yàn)的試樣尺寸，建立常規(guī)離散元虛擬壓縮試驗(yàn)?zāi)Ｐ汀?/p>

假設(shè)虛擬試樣代表的實(shí)際尺寸為5cm×10cm(w×h)。圖2所示為根據(jù)試樣尺寸建立的虛擬壓縮試驗(yàn)?zāi)Ｐ?，其中試樣顆粒模型是采用采用動力方法生成的，模型中顆?？倲?shù)為3314，顆粒的最小半徑為0.5mm，顆粒的最大半徑為0.8mm，顆粒的平均半徑為0.62mm。

2、將試樣模型伺服到初始固結(jié)應(yīng)力狀態(tài)。

實(shí)施例中壓縮試驗(yàn)的圍壓選取為500.0kpa，利用應(yīng)力伺服原理將上下墻體和左右墻體的應(yīng)力伺服到圍壓目標(biāo)值。應(yīng)力伺服終止的容差ε＝0.001，最終伺服結(jié)束后，上下墻體和左右墻體的平均應(yīng)力值分別為500.4kpa和499.8kpa。如圖4所示，伺服后的試樣尺寸變?yōu)?.2cm×10.1cm(w′×h′)。

3、建立離散-連續(xù)耦合虛擬數(shù)值試驗(yàn)?zāi)Ｐ汀?/p>

如圖5所示，根據(jù)伺服后試驗(yàn)?zāi)Ｐ?，刪除側(cè)向墻體，并在側(cè)向墻體位置建立兩列有限元單元，每列單元的數(shù)目均為20個(gè)，且每個(gè)單元的長度為0.01cm，高度為h′/20＝10.1/20＝0.505cm；然后根據(jù)建立的有限元單元，沿單元邊建立了兩組耦合墻體，每個(gè)墻體的長度為單元的高度，且墻體節(jié)點(diǎn)與單元節(jié)點(diǎn)在位置上是重疊的，如圖5所示；最后按照圖7所示的對應(yīng)關(guān)系，存儲單元節(jié)點(diǎn)與控制墻體頂點(diǎn)的對應(yīng)關(guān)系。

4、開展離散-連續(xù)耦合虛擬壓縮試驗(yàn)數(shù)值仿真模擬。

首先，如圖8中所示的圍壓和速度的施加方式，在左右側(cè)有限元單元外側(cè)施加500kpa圍壓，上下墻體及底部和頂部單元節(jié)點(diǎn)的移動速度v設(shè)為0.0005m/；其次將有限元和離散元設(shè)置為相同大小的分析時(shí)步，時(shí)步取為1.0e-5s；最后根據(jù)圖9所的耦合計(jì)算原理進(jìn)行虛擬壓縮試驗(yàn)耦合模擬，當(dāng)軸向應(yīng)變達(dá)到10％時(shí)，停止加載，并在模擬過程中記錄軸向應(yīng)變與軸向偏應(yīng)力。圖10為虛擬壓縮試驗(yàn)過程記錄的軸向應(yīng)變與軸向偏應(yīng)力曲線圖。圖11為壓縮試驗(yàn)后的虛擬試驗(yàn)?zāi)Ｐ蛨D，試樣高度最終被壓縮為h″＝h′-h′×10％＝9.1cm，圖12a和圖12b分別為壓縮試驗(yàn)后試驗(yàn)有限元節(jié)點(diǎn)位移矢量圖和試樣顆粒位移矢量圖。由圖11和圖12可看出，虛擬試驗(yàn)壓縮過程中側(cè)向有限元單元變形是不一致，實(shí)現(xiàn)了室內(nèi)壓縮試驗(yàn)側(cè)向圍壓柔性加載的目的，達(dá)到了本發(fā)明的最終目標(biāo)，且試樣壓縮后發(fā)生鼓脹破壞，這與室內(nèi)壓試驗(yàn)試樣破壞形式也相一致，也進(jìn)一步驗(yàn)證了本文明的優(yōu)勢。