巖石爆破理論及LS-DYNA數(shù)值模擬方法

巖石在爆破荷載作用下表現(xiàn)的性質(zhì)十分復(fù)雜，在實(shí)際的鉆爆開挖中容易引發(fā)安全事故。埋入巖體中的炸藥發(fā)生化學(xué)爆炸后，巖石將在爆炸近區(qū)發(fā)生壓剪破壞形成粉碎區(qū)，在爆炸遠(yuǎn)區(qū)發(fā)生張拉破壞形成裂隙區(qū)。非線性動(dòng)力分析有限元軟件LS-DYNA以其高效、準(zhǔn)確等優(yōu)點(diǎn)被廣泛應(yīng)用于爆破、沖擊等領(lǐng)域。對(duì)于巖石爆破這類大變形問(wèn)題，其中的ALE算法能有效避免因網(wǎng)格畸變而導(dǎo)致的計(jì)算難題。它結(jié)合了Lagrange算法與Euler算法的優(yōu)點(diǎn)，可用于解決流體-固體耦合問(wèn)題。

本人針對(duì)各種工況下的爆破方法錄制了大量視頻課程，如對(duì)爆破數(shù)值模擬感興趣的朋友可以購(gòu)買學(xué)習(xí)。

1.1 巖石爆破理論

1.1.1 無(wú)限域巖石的爆破

現(xiàn)有的爆破理論認(rèn)為，埋入無(wú)限巖石中的炸藥發(fā)生化學(xué)爆炸后，將在巖石中形成以裝藥為中心由近及遠(yuǎn)的不同程度的破壞區(qū)域，依次稱為粉碎區(qū)、裂隙區(qū)和彈性震動(dòng)區(qū)，如圖1所示。

?

圖1 無(wú)限域巖石單孔爆破

Fig 1 Single hole blasting of infinite field rock

?

巖石在粉碎區(qū)內(nèi)受到的爆炸載荷加載率最高，該處應(yīng)力波的強(qiáng)度遠(yuǎn)大于巖石的動(dòng)態(tài)抗壓強(qiáng)度，因此巖石發(fā)生粉碎性破壞并形成空腔。另外，粉碎區(qū)距離炸藥較近，還會(huì)受到爆炸氣體的高溫高壓作用。隨著爆炸應(yīng)力波向遠(yuǎn)處傳播，應(yīng)力波強(qiáng)度逐漸減小，當(dāng)其強(qiáng)度小于巖石動(dòng)態(tài)抗壓強(qiáng)度后，巖石內(nèi)將不再形成粉碎區(qū)。然后，巖石徑向受壓時(shí)會(huì)伴隨有環(huán)向拉應(yīng)力，由于巖石的動(dòng)態(tài)抗拉強(qiáng)度遠(yuǎn)小于其動(dòng)態(tài)抗壓強(qiáng)度，因而巖石會(huì)在環(huán)向拉應(yīng)力的作用下發(fā)生張拉破壞，并形成徑向拉伸裂紋。在裂隙區(qū)外，隨著應(yīng)力波繼續(xù)衰減，此時(shí)巖石中不再出現(xiàn)明顯的破壞，可認(rèn)為該區(qū)域的巖石只產(chǎn)生了彈性震動(dòng)，因而稱為彈性震動(dòng)區(qū)。

具體實(shí)現(xiàn)過(guò)程可參考本人課程《（SCI）LS-DYNA的巖石單孔及雙孔爆破裂紋擴(kuò)展模擬》和《LS-DYNA的三維巖石爆破裂紋擴(kuò)展模擬》。

1.1.2 臨自由面巖石的爆破

如前所述，埋入無(wú)限巖石中的炸藥爆炸后，將在巖石中依次形成粉碎區(qū)、裂隙區(qū)和彈性震動(dòng)區(qū)。當(dāng)埋入巖石中的炸藥臨近自由面時(shí)，爆炸應(yīng)力波會(huì)在自由面附近發(fā)生反射作用。此時(shí)，徑向壓縮波會(huì)在自由面處反射為徑向拉伸波，使得自由面附近的巖石內(nèi)萌生更多的環(huán)向拉伸裂紋。而根據(jù)炸藥到巖石自由面距離的不同，還會(huì)在自由表面附近引起巖石的破裂、鼓包和拋擲，進(jìn)一步形成漏斗狀的炸坑，稱為爆破漏斗。

具體實(shí)現(xiàn)過(guò)程可參考本人課程《LS-DYNA巖石爆破漏斗模擬》和《LS-DYNA的三維爆破SPH-FEM耦合法》。

1.2 LS-DYNA軟件介紹

1976年，美國(guó)Lawrence Livermore實(shí)驗(yàn)室的John Hallquist博士發(fā)布了最早的DYNA軟件。20世紀(jì)80年代，該軟件被法國(guó)ESI公司商品化。1989年，John Hallquist博士開始經(jīng)營(yíng)自己的公司Livermore Software Technology Corporation(LSTC)，并推廣了如今的DYNA軟件版本：LS-DYNA。

在此之后，LS-DYNA陸續(xù)推出了930版、936版、940版、950版、960版和970版，并且ANSYS公司已收購(gòu)了LS-DYNA軟件的使用權(quán)。如今，LS-DYNA軟件已是一款功能強(qiáng)大的幾何非線性（大位移、大轉(zhuǎn)動(dòng)和大應(yīng)變）、材料非線性（140多種材料動(dòng)態(tài)模型）以及摩擦和接觸分離等界面狀態(tài)的非線性程序。它以Lagrange算法為主，兼有ALE和Euler算法；以顯式求解為主，兼有隱式求解功能；以結(jié)構(gòu)分析為主，兼有熱分析和流固耦合功能；以非線性動(dòng)力分析為主，兼有靜力分析功能。綜合來(lái)看，它是一款軍用和民用相結(jié)合的通用結(jié)構(gòu)分析非線性有限元程序。

LS-DYNA顯示動(dòng)力分析采用了中心差分法，結(jié)構(gòu)系統(tǒng)各節(jié)點(diǎn)在第n個(gè)時(shí)間步結(jié)束時(shí)刻tn的加速度向量為：

?????????????????????????????????????????????? ???????????????????????????????????????? (2)

式中：P為第n個(gè)時(shí)間步結(jié)束時(shí)刻結(jié)構(gòu)上所施加的節(jié)點(diǎn)外力向量；Fint為tn時(shí)刻的內(nèi)力矢量，表達(dá)式為：

???????????????????????????????????????????? ?????????????????????????????????????? (3)

式中：右邊三項(xiàng)依次為單元應(yīng)力場(chǎng)等效節(jié)點(diǎn)力、沙漏阻力和接觸力矢量。

根據(jù)中心差分法的基本思路，加速度可由速度的一階中心差分求出，而速度可由位移的一階中心差分計(jì)算得到。因而有如下表達(dá)式：

????????????????????????????????? ??????????????????????????? (4)

???????????????????????????????????????????? ?????????????????????????????????????? (5)

因此，節(jié)點(diǎn)速度向量可由程序計(jì)算出的加速度結(jié)合差分公式表示，而節(jié)點(diǎn)位移向量可由節(jié)點(diǎn)速度向量結(jié)合差分公式表示，即：

???????????????????????????????????? ?????????????????????????????? (6)

??????????????????????????????????????????????? ????????????????????????????????????????? (7)

新的幾何構(gòu)型由初始構(gòu)型x0和位移增量u相加得到，即：

???????????????????????????????????????????????????????? ????????????????????????????????????????????????? (8)

為了保證數(shù)值計(jì)算過(guò)程中的穩(wěn)定性，應(yīng)當(dāng)滿足自由度解在若干時(shí)間步后仍是有限值。通過(guò)理論分析可知，保證數(shù)值計(jì)算穩(wěn)定的臨界時(shí)間步長(zhǎng)應(yīng)滿足：

??????????????????????????????????????????????????????? ????????????????????????????????????????????????? (9)

式中：ωn為系統(tǒng)中各單元的最高階固有振動(dòng)頻率，一般由系統(tǒng)中的最小單元決定。系統(tǒng)中最小單元的振動(dòng)特征值方程為：

???????????????????????????????????????????????????????? ???????????????????????????????????????????????? (10)

由式(10)可計(jì)算出的最大特征值即為臨界時(shí)間步長(zhǎng)中的各單元最高階固有振動(dòng)頻率。為了保證數(shù)值計(jì)算時(shí)的收斂，LS-DYNA程序采用了變步長(zhǎng)積分法，每時(shí)刻的積分步長(zhǎng)由當(dāng)前構(gòu)形網(wǎng)格中的最小單元來(lái)確定。而在具體的計(jì)算中，還可在LS-DYNA程序的默認(rèn)時(shí)間步長(zhǎng)前乘一個(gè)小于1的時(shí)間步調(diào)整系數(shù)。

1.3 爆破數(shù)值模擬方法

為避免孔壁周圍的巖石單元發(fā)生大變形而出現(xiàn)負(fù)體積的錯(cuò)誤，采用流固耦合方法。另外，巖體的預(yù)應(yīng)力和重復(fù)爆破過(guò)程分別由Dynain文件和完全重啟動(dòng)技術(shù)實(shí)現(xiàn)。

1.3.1 流固耦合方法

在LS-DYNA軟件中主要有四種數(shù)值方法模擬巖石爆破過(guò)程，包括共節(jié)點(diǎn)方法、接觸方法、施加爆炸應(yīng)力波荷載方法和流固耦合方法。具體實(shí)現(xiàn)過(guò)程如下：

（1）共節(jié)點(diǎn)方法

巖石和炸藥在建模時(shí)采用共節(jié)點(diǎn)，如圖2所示。該方法的計(jì)算速度很快，但缺點(diǎn)也非常明顯，即炸藥單元的畸變會(huì)引起巖石單元也產(chǎn)生較大的變形，同時(shí)由于采用了共節(jié)點(diǎn)，炸藥單元的滑移變形會(huì)受到限制，進(jìn)而引起附加的虛假滑移剛度，會(huì)對(duì)計(jì)算精度產(chǎn)生影響。

?

圖2 共節(jié)點(diǎn)方法

Fig 2 Common node method

?

（2）接觸方法

巖石與炸藥分開建模，并在兩個(gè)部件之間設(shè)置侵蝕接觸或滑移接觸，如圖3所示。該方法存在的問(wèn)題是炸藥單元會(huì)產(chǎn)生較大變形，容易出現(xiàn)求解過(guò)程突然終止且不報(bào)錯(cuò)的現(xiàn)象，并且網(wǎng)格畸變也會(huì)對(duì)計(jì)算精度產(chǎn)生很大的影響。另外，接觸的設(shè)置有可能導(dǎo)致計(jì)算不收斂。該方法往往需要借助小型重啟動(dòng)技術(shù)，在炸藥爆炸過(guò)程基本完成后通過(guò)*DELETE_PART關(guān)鍵字刪除炸藥PART，并使用小型重啟動(dòng)繼續(xù)計(jì)算。

?

圖3 接觸方法

Fig 3 Contact method

?

（3）施加爆炸應(yīng)力波荷載方法

該方法不需要建立炸藥模型，只需建立巖石模型即可。典型的爆炸應(yīng)力波荷載曲線如圖4所示，可以根據(jù)荷載峰值和加載持續(xù)時(shí)間描述不同藥量和爆速的炸藥。該方法計(jì)算效率較大，但最明顯的缺點(diǎn)是無(wú)法描述高溫高壓的爆炸氣體作用于孔壁的過(guò)程和爆生氣體促進(jìn)裂紋擴(kuò)展的準(zhǔn)靜態(tài)作用。

具體實(shí)現(xiàn)過(guò)程可參考本人發(fā)布的課程《LS-DYNA施加爆炸應(yīng)力波曲線法模擬巖石爆破裂紋擴(kuò)展》。

圖4 爆炸應(yīng)力波荷載曲線

Fig 4 Load curve of explosion stress wave

?

圖5 流固耦合方法

Fig 5 Fluid-structure coupling method

?

（4）流固耦合方法

巖石采用Lagrange算法，空氣和炸藥采用ALE或Euler算法，并通過(guò)*CONSTRAINED_LAGRANGE_IN_SOLID關(guān)鍵字實(shí)現(xiàn)流固耦合，如圖5所示。其中，Lagrange算法的單元作為從面，ALE或Euler算法的單元作為主面。在前處理建模時(shí)，巖石應(yīng)和空氣、炸藥采用分離式建模，而空氣和炸藥之間采用共節(jié)點(diǎn)建模，并將空氣和炸藥放在同一個(gè)ALE組里。需注意的是，應(yīng)在炮孔周圍設(shè)置更大范圍的空氣PART，這部分空氣單元應(yīng)與巖石單元相互重合，作為耦合域。

該方法能有效避免因網(wǎng)格畸變而導(dǎo)致的計(jì)算難題，大大提高了計(jì)算精度，被廣泛應(yīng)用于巖石爆破和沖擊領(lǐng)域。本人發(fā)布的大多數(shù)課程均采用流固耦合方法。

1.3.2 預(yù)應(yīng)力巖體爆破模擬方法

由于深部巖體處在預(yù)應(yīng)力的狀態(tài)下，其圍壓的施加屬于靜力分析范疇，而爆炸是超動(dòng)態(tài)過(guò)程，屬于動(dòng)力分析范疇。因此，這兩個(gè)步驟需要分開進(jìn)行求解。目前，巖石地應(yīng)力的施加主要有三種方法：動(dòng)力松弛法、Dynain文件法和顯隱式序列法。其中，Dynain文件法兼具準(zhǔn)確、高效和操作方便等優(yōu)點(diǎn)。該方法的主要分析過(guò)程如下：

（1）建立爆破數(shù)值模型，并輸出原始K文件；

（2）基于原始K文件，將空氣、炸藥和堵塞PART刪除，定義圍壓載荷曲并施加于巖石模型外邊界，并通過(guò)*INTERFACE_SPRINGBACK_LSDYNA關(guān)鍵字輸出包含巖石應(yīng)力狀態(tài)的Dynain文件，進(jìn)而將該K文件遞交求解，完成靜力分析過(guò)程；

（3）基于原始K文件，將巖石PART刪除，并通過(guò)*INCLUDE關(guān)鍵字將預(yù)應(yīng)力巖石模型導(dǎo)入，最后將該K文件遞交求解，完成動(dòng)力分析過(guò)程。

具體實(shí)現(xiàn)過(guò)程可參考本人發(fā)布的課程《LS-DYNA地應(yīng)力下的隧道爆破三維模擬》和《（SCI復(fù)現(xiàn)）LS-DYNA地應(yīng)力掏槽爆破》。

1.3.3 重啟動(dòng)模擬方法

對(duì)于巖石循環(huán)或重復(fù)起爆問(wèn)題，需要借助重啟動(dòng)技術(shù)進(jìn)行分析。在LS-DYNA軟件中，重啟動(dòng)分析主要分為以下三種類型：

（1）簡(jiǎn)單重啟動(dòng)

當(dāng)求解過(guò)程被意外終止時(shí)，用戶可直接根據(jù)生成的d3dump01文件進(jìn)行簡(jiǎn)單重啟動(dòng)分析，不需要對(duì)K文件進(jìn)行任何修改。此時(shí)，求解將從意外終止的時(shí)刻開始，繼續(xù)求解至最初預(yù)定的求解完成時(shí)刻。

（2）小型重啟動(dòng)

當(dāng)需要對(duì)分析的求解設(shè)置和模型進(jìn)行修改，并且仍要繼續(xù)求解運(yùn)算時(shí)，可采用小型重啟動(dòng)技術(shù)。比如，求解達(dá)到預(yù)定計(jì)算時(shí)間結(jié)束后，發(fā)現(xiàn)爆炸仍在繼續(xù)進(jìn)行，則用戶可以修改計(jì)算時(shí)長(zhǎng)，并將新的K文件和上次求解生成的d3dump01文件遞交給求解器進(jìn)行小型重啟動(dòng)分析。

（3）完全重啟動(dòng)

當(dāng)需要對(duì)數(shù)值模型進(jìn)行較大修改，如需要增減PART、增減接觸信息、修改荷載條件和修改控制參數(shù)等，則需要用到完全重啟動(dòng)分析。該技術(shù)可以在繼承前一次分析結(jié)果的基礎(chǔ)上進(jìn)行全新的分析。

循環(huán)爆破、循環(huán)掘進(jìn)或重復(fù)侵徹模擬的具體實(shí)現(xiàn)過(guò)程可參考本人發(fā)布的課程《LS-DYNA巖石循環(huán)爆破/重復(fù)起爆-完全重啟動(dòng)技術(shù)》、《LS-DYNA彈體多次重復(fù)侵徹靶體-完全重啟動(dòng)》和《LS-DYNA流固耦合法模擬隧道循環(huán)掘進(jìn)爆破（建炸藥、完全重啟動(dòng)）-循環(huán)進(jìn)尺3次》。