王海兵， 張海波， 田宙， 歐卓成， 周剛

(1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 100081; 2.西北核技術(shù)研究所, 陜西 西安 710024)

?

巖石動(dòng)力學(xué)計(jì)算中的網(wǎng)格效應(yīng)及機(jī)理研究

王海兵1,2， 張海波2， 田宙2， 歐卓成1， 周剛2

(1.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 100081; 2.西北核技術(shù)研究所, 陜西 西安 710024)

巖石動(dòng)力學(xué)計(jì)算中，網(wǎng)格尺寸對(duì)數(shù)值計(jì)算結(jié)果的可靠性有重要影響。采用數(shù)值實(shí)驗(yàn)的方法，對(duì)巖石爆炸應(yīng)力波傳播數(shù)值計(jì)算中的網(wǎng)格尺寸效應(yīng)及其敏感性機(jī)理進(jìn)行了研究。研究結(jié)果表明：合適的網(wǎng)格尺寸要根據(jù)載荷特征和波傳播介質(zhì)的屬性來(lái)決定；當(dāng)一個(gè)載荷波長(zhǎng)內(nèi)的網(wǎng)格數(shù)達(dá)到16個(gè)以上時(shí)，計(jì)算得到的各物理量的波形和峰值基本趨于穩(wěn)定；計(jì)算還給出了各物理量與網(wǎng)格密度的關(guān)系；隨著爆心距的增加，物理量對(duì)網(wǎng)格尺寸的敏感性降低，其機(jī)理是載荷中的高頻成分逐漸衰減、載荷的波長(zhǎng)變大，模型所需的網(wǎng)格尺寸變大；時(shí)間步長(zhǎng)系數(shù)對(duì)計(jì)算結(jié)果的影響也非常明顯，當(dāng)時(shí)間步長(zhǎng)系數(shù)取0.05時(shí)，位移穩(wěn)態(tài)值趨于收斂值。

兵器科學(xué)與技術(shù)； 數(shù)值計(jì)算； 網(wǎng)格尺寸； 時(shí)間步長(zhǎng)； 敏感性機(jī)理

0　引言

巖石中的爆炸動(dòng)力學(xué)問題廣泛應(yīng)用于國(guó)家經(jīng)濟(jì)建設(shè)、國(guó)防軍事以及科研領(lǐng)域等方面。巖石中的爆炸動(dòng)力學(xué)過程非常復(fù)雜，通常很難進(jìn)行精確地解析求解，數(shù)值計(jì)算是常用的有效手段。該方面的計(jì)算主要涉及三方面的物理過程：一是炸藥的起爆及爆炸產(chǎn)物與沖擊波在爆室內(nèi)的傳播；二是爆炸產(chǎn)物及沖擊波與圍巖的相互作用；三是圍巖的動(dòng)態(tài)響應(yīng)。與這三方面物理過程相對(duì)應(yīng)的三類計(jì)算問題為：爆炸載荷的求解問題、爆炸載荷與圍巖的能量耦合問題以及應(yīng)力波在圍巖中的傳播問題。計(jì)算這些問題常用的數(shù)值方法是有限差分方法和有限元方法。與有限差分相比，有限元方法具有靈活的網(wǎng)格劃分、邊界條件和邊界形狀的處理等優(yōu)點(diǎn)，現(xiàn)在已成為爆炸動(dòng)力學(xué)中的重要數(shù)值方法。通常來(lái)講，數(shù)值計(jì)算的精度依賴于描述材料行為所采用的數(shù)學(xué)模型，如本構(gòu)關(guān)系和狀態(tài)方程等。由于目前描述爆炸作用下材料響應(yīng)的數(shù)學(xué)模型還不完善，且都是近似的，因此數(shù)值分析的精度一般不高于近似方程的精度[1]。如果不考慮數(shù)學(xué)模型的近似程度，單從數(shù)值計(jì)算的角度講，一種算法的穩(wěn)定性與精度決定于[2]：1)積分格式與其所用的參數(shù)；2)時(shí)間步長(zhǎng)；3)計(jì)算模型所劃分的網(wǎng)格尺寸。

由于爆炸波持續(xù)時(shí)間非常短并且能量在不同的網(wǎng)格之間傳輸，數(shù)值模擬結(jié)果對(duì)有限單元的尺寸非常敏感，數(shù)值結(jié)果的精度強(qiáng)烈地依賴于所使用的網(wǎng)格尺寸。通常需要開展網(wǎng)格尺寸收斂性測(cè)試以獲得合適的網(wǎng)格尺寸。爆炸沖擊波數(shù)值模擬中網(wǎng)格尺寸效應(yīng)研究的核心問題就在于如何針對(duì)具體的爆炸問題確定合適的網(wǎng)格尺寸，以在確保數(shù)值模擬精度的同時(shí)盡可能地減少網(wǎng)格數(shù)量和提高計(jì)算效率。然而，由于各類介質(zhì)的物理屬性差異，使得針對(duì)某一種介質(zhì)開展網(wǎng)格尺寸效應(yīng)分析確定的網(wǎng)格尺寸在介質(zhì)改變時(shí)，其適用性具有較大的局限性。采用同一網(wǎng)格尺寸模擬不同介質(zhì)中的沖擊波傳播問題時(shí)，計(jì)算精度可能存在較大區(qū)別。因此，合適的網(wǎng)格尺寸不僅與所模擬的爆炸情況相關(guān)，還與爆炸沖擊波的傳播介質(zhì)相關(guān)[3]。當(dāng)爆炸載荷作用在爆室壁圍巖時(shí)，就涉及流體與固體(簡(jiǎn)稱流固)耦合問題，研究表明，流固耦合算法的選取及關(guān)鍵耦合參數(shù)的設(shè)置以及爆炸流體和圍巖固體網(wǎng)格尺寸的相對(duì)大小，對(duì)計(jì)算結(jié)果都有很大的影響。限于篇幅限制，本文關(guān)于網(wǎng)格尺寸對(duì)流固耦合效果的影響暫不展開討論，主要研究爆炸沖擊波在空氣中及巖石內(nèi)傳播的網(wǎng)格尺寸效應(yīng)。

1　網(wǎng)格尺寸對(duì)爆炸場(chǎng)載荷計(jì)算的影響

爆炸載荷的確定依賴于爆炸場(chǎng)的求解，也直接影響著外圍結(jié)構(gòu)響應(yīng)的計(jì)算和安全性評(píng)估的可靠性，因此準(zhǔn)確的爆炸載荷計(jì)算非常重要。一般情況下，爆室或容器內(nèi)的爆炸屬于近區(qū)爆炸，在爆炸近區(qū)，爆轟產(chǎn)物的密度遠(yuǎn)高于空氣密度，產(chǎn)物對(duì)爆炸載荷有著顯著的影響。對(duì)此，國(guó)內(nèi)外學(xué)者已開展了諸多研究，研究結(jié)果表明，在爆炸近區(qū)網(wǎng)格尺寸對(duì)計(jì)算結(jié)果的精度影響非常敏感。胡八一等[4]指出，要想模擬出來(lái)的反射沖擊波壓力峰值及波形與實(shí)測(cè)波形一致，最為關(guān)鍵的一點(diǎn)是爆室內(nèi)或容器內(nèi)部網(wǎng)格尺寸至少應(yīng)小于或等于1 mm. Chapman[5]在用類似的商業(yè)軟件計(jì)算沖擊波超壓時(shí)也發(fā)現(xiàn)當(dāng)網(wǎng)格由10 mm細(xì)化到1 mm時(shí)，超壓峰值提高了24.6%. 文獻(xiàn)[6-7]的計(jì)算結(jié)果顯示，采用LS-DYNA等商業(yè)軟件計(jì)算得到的計(jì)算結(jié)果小于經(jīng)驗(yàn)公式值，受軟件和計(jì)算機(jī)硬件容量的限制，只能提供定性的計(jì)算結(jié)果，很難提供準(zhǔn)確的定量計(jì)算。雷鳴等[8]采用Autodyn研究了網(wǎng)格大小對(duì)空中化爆自由場(chǎng)參數(shù)計(jì)算結(jié)果的影響，結(jié)果表明網(wǎng)格大小對(duì)質(zhì)量守恒沒有影響，對(duì)能量守恒有一定的影響，爆炸產(chǎn)物界面第一次擴(kuò)張是在大約0.9 m/kg1/3爆心距位置附近，在比例距離1 m/kg1/3內(nèi)每個(gè)單元取0.5 mm，在更遠(yuǎn)的區(qū)域每個(gè)單元5 mm計(jì)算得到的值與標(biāo)準(zhǔn)實(shí)驗(yàn)值較符合較好。姚成寶等[9]利用LS-DYNA計(jì)算空中爆炸時(shí)發(fā)現(xiàn)，網(wǎng)格尺寸在0.1～5 mm時(shí)計(jì)算結(jié)果差別不大，在網(wǎng)格尺寸小于1 mm時(shí)計(jì)算結(jié)果趨于收斂，當(dāng)網(wǎng)格的密度尺寸小于5 mm/kg1/3時(shí)，計(jì)算結(jié)果基本趨于收斂，此時(shí)得到的結(jié)果基本滿足精度要求。王海兵等[10]采用二維方法研究了容器內(nèi)的爆炸載荷，結(jié)果指出：在比距離小于1.0 m/kg1/3時(shí)，2 mm的網(wǎng)格尺寸計(jì)算結(jié)果趨于收斂；在比距離大于1.0 m/kg1/3、網(wǎng)格尺寸取5 mm時(shí)，計(jì)算得到的峰值超壓基本滿足精度。Shi等[11]詳細(xì)比較了網(wǎng)格尺寸對(duì)空中爆炸自由場(chǎng)諸參量的影響，圖1～圖3給出了其主要結(jié)果。圖1給出了不同網(wǎng)格尺寸在比距離R為1 m/kg1/3處的正超壓時(shí)程曲線比較，從圖中可看出，隨著網(wǎng)格尺寸的增加，壓力p的上升率減緩、波形變得平緩、波幅拉寬，峰值降低；但波前到達(dá)時(shí)間和正壓作用時(shí)間與網(wǎng)格尺寸無(wú)關(guān)。圖2為峰值超壓與網(wǎng)格尺寸和比距離的關(guān)系，從圖中看出，當(dāng)比距離小于0.8 m/kg1/3時(shí)，即使5 mm的網(wǎng)格尺寸也不能準(zhǔn)確計(jì)算峰壓；而當(dāng)比距離超過2 m/kg1/3時(shí)，100 mm的網(wǎng)格尺寸也能準(zhǔn)確地計(jì)算峰壓。圖3給出了正壓沖量I與網(wǎng)格尺寸和比距離的關(guān)系，從圖中看出，正壓沖量與網(wǎng)格尺寸基本無(wú)關(guān)，100 mm的網(wǎng)格尺寸已能達(dá)到計(jì)算收斂?？偨Y(jié)以上學(xué)者的研究結(jié)果，可以得到以下結(jié)論：當(dāng)爆炸近區(qū)的網(wǎng)格尺寸取毫米量級(jí)時(shí)，可使計(jì)算獲得的所有力學(xué)參量符合實(shí)測(cè)；隨著網(wǎng)格尺寸的增加，壓力的上升率減緩、波形變得平緩、波幅拉寬，峰值降低；超壓峰值對(duì)網(wǎng)格的敏感性很高，隨著比距離的增加，超壓峰值對(duì)網(wǎng)格的敏感降低；超壓峰值和壓力上升率是網(wǎng)格密度的正相關(guān)量，而波前到達(dá)時(shí)間、正壓作用時(shí)間、正壓沖量等是網(wǎng)格尺寸的無(wú)關(guān)量。

圖1　　　　不同網(wǎng)格尺寸在比距離1 m/kg1/3處的正超壓時(shí)程曲線比較[11]Fig.1　　　　Comparison of pressure-time histories obtained with different mesh sizes at scaled distance of 1 m/kg1/3[11]

圖2　峰壓與網(wǎng)格尺寸和比距離的關(guān)系[11]Fig.2　　　　Variation of the positive incident peak pressure with mesh size and scaled distance[11]

圖3　正壓沖量與網(wǎng)格尺寸和比距離的關(guān)系[11]Fig.3　　　　Variation of normalized positive incident impulse with mesh size and scaled distance[11]

在實(shí)際應(yīng)用中，要使計(jì)算達(dá)到工程使用價(jià)值，計(jì)算模型的尺度至少要達(dá)到1～10 m的范圍，在如此大的范圍內(nèi)，若將所有的力學(xué)參量都計(jì)算準(zhǔn)確，需要將單元尺寸細(xì)化到毫米量級(jí)。對(duì)于三維計(jì)算，即使采取1/8對(duì)稱模型，網(wǎng)格數(shù)量也將達(dá)到109～1012個(gè)量級(jí)左右。這樣的計(jì)算規(guī)模即使采用目前最為先進(jìn)的巨型工作站進(jìn)行計(jì)算，也絕非易事。但通常關(guān)注的效應(yīng)是爆炸作用下的結(jié)構(gòu)響應(yīng)，而決定結(jié)構(gòu)響應(yīng)的參量是正壓沖量。前文的網(wǎng)格收斂性計(jì)算說明，正壓沖量是網(wǎng)格尺寸的不敏感量，100 mm的網(wǎng)格尺寸已經(jīng)能基本滿足計(jì)算精度要求。所以，當(dāng)計(jì)算規(guī)模較大時(shí)，雖然由于網(wǎng)格尺寸限制，不能對(duì)超壓峰值進(jìn)行準(zhǔn)確計(jì)算，但仍然能夠準(zhǔn)確計(jì)算正壓沖量，進(jìn)而能夠獲得準(zhǔn)確的爆炸效應(yīng)。文獻(xiàn)[11]還提出一種網(wǎng)格尺寸效應(yīng)的數(shù)值修正方法，其原理是利用網(wǎng)格尺寸的不敏感量來(lái)修正網(wǎng)格尺寸敏感量，從而提高大網(wǎng)格尺寸情況下有限元模擬結(jié)果的精度。具體方法是通過沖量等效的原理來(lái)修正超壓峰值，修正后，即使采用200 mm的網(wǎng)格尺寸也能得到較為精確的超壓峰值。

2　網(wǎng)格尺寸對(duì)巖石中應(yīng)力波傳播計(jì)算的影響

關(guān)于網(wǎng)格尺寸對(duì)巖石動(dòng)力學(xué)計(jì)算的影響，已有學(xué)者開展了相關(guān)的研究。梁正召等[12]討論了非均勻性巖石中網(wǎng)格尺寸對(duì)裂紋擴(kuò)展即巖石破裂的影響，研究認(rèn)為，隨著網(wǎng)格尺寸的減小，巖石強(qiáng)度逐漸降低，并趨于穩(wěn)定；在均質(zhì)巖石中，網(wǎng)格尺寸主要和結(jié)構(gòu)特征相關(guān)；在非均質(zhì)巖石中，網(wǎng)格尺寸還與材料的非均勻性和細(xì)觀特征尺度密切相關(guān)，引入單元統(tǒng)計(jì)性的非均勻性后，單元尺寸必須接近細(xì)觀特征尺度才能保證計(jì)算結(jié)果的穩(wěn)定性和可靠性。崔煥平等[13]提出了把極限拉應(yīng)變作為網(wǎng)格尺寸的函數(shù)，以保證混凝土開裂單位面積吸收的能量唯一，從而有效地消除混凝土非線性有限元分析中的網(wǎng)格尺寸效應(yīng)。門建兵等[14]研究了網(wǎng)格對(duì)混凝土侵徹?cái)?shù)值模擬的影響，研究認(rèn)為彈丸半徑方向要?jiǎng)澐?個(gè)網(wǎng)格，靶板在彈丸半徑尺寸方向上劃分6個(gè)網(wǎng)格得到的計(jì)算結(jié)果比較理想。Shayanfar等[15]討論了網(wǎng)格尺寸對(duì)載荷- 位移、載荷- 應(yīng)變、裂紋模式和最終載荷的影響，并開發(fā)了一個(gè)簡(jiǎn)單的程序來(lái)計(jì)算混凝土的最終拉伸應(yīng)變，以此來(lái)消除計(jì)算結(jié)果對(duì)網(wǎng)格尺寸的依賴，并將其嵌入到非線性有限元程序中。以上的研究大多都是通過網(wǎng)格的收斂性測(cè)試來(lái)獲得合適的網(wǎng)格尺寸，然而對(duì)于大規(guī)模計(jì)算來(lái)講，很難對(duì)模型進(jìn)行整體的網(wǎng)格收斂性測(cè)試；另外，這些研究并未考慮載荷的特征，而爆炸沖擊型載荷包含大量的高頻分量，其輸入荷載的波形與模型網(wǎng)格劃分密切相關(guān)，如果不同時(shí)考慮荷載的頻譜特性，在確定時(shí)間步長(zhǎng)和網(wǎng)格劃分時(shí)往往帶有較大的盲目性，得到的結(jié)果也不具備一般通用性。因此，本文以輸入載荷的特性為出發(fā)點(diǎn)，來(lái)研究網(wǎng)格尺寸對(duì)巖石中應(yīng)力波傳播計(jì)算的影響。

2.1典型的爆炸載荷特征

以1 t TNT爆炸為例，讀取爆心距3.06 m位置處的壓力波形如圖4所示。從圖4中的載荷等效曲線確定外載荷的周期約為300 μs.

圖4　1 t TNT爆炸在3.06 m位置處的壓力波形Fig.4　　　　Pressure waveform at 3.06 m from the distance of 1 t TNT explosive source

2.2巖石中的波速及波長(zhǎng)

彈性波傳播的速度不僅僅取決于波的類型，而且還與介質(zhì)周圍的條件有關(guān)。彈性縱波在單軸應(yīng)變(相當(dāng)于無(wú)限介質(zhì))條件下傳播的速度為

(1)

式中：ρ、Λ和μ分別是介質(zhì)的密度、拉梅常數(shù)和剪切模量；cim是彈性波在無(wú)限介質(zhì)中的傳播速度。

對(duì)于一維桿的波傳播問題，根據(jù)波動(dòng)方程，有

(2)

式中：E是彈性模量；codb是彈性波在一維桿中的傳播速度?！盁o(wú)限介質(zhì)”和“桿”等模型都是一些極限情況的近似，所謂“無(wú)限介質(zhì)”是指在物理上介質(zhì)的尺度遠(yuǎn)遠(yuǎn)大于所傳播波的波長(zhǎng)，“桿”是指介質(zhì)的橫向尺度遠(yuǎn)遠(yuǎn)小于波的波長(zhǎng)。通常研究的介質(zhì)其尺度各不相同，研究的波長(zhǎng)也從地震波的幾百千米至實(shí)驗(yàn)室內(nèi)超聲波波長(zhǎng)的幾個(gè)毫米乃至幾微米不等。因此，在討論波速及波傳播問題時(shí)，需要研究介質(zhì)的尺度與波長(zhǎng)λ大小的相互關(guān)系。取花崗巖參數(shù)：密度為2.62 g/cm3，泊松比υ為0.24，彈性模量為58.7 GPa. 再根據(jù)彈性力學(xué)關(guān)系式：

(3)

(4)

得到彈性縱波在花崗巖無(wú)限介質(zhì)中的傳播速度cim=5 138 m/s，在花崗巖中一維桿的波速codb=4 733 m/s.

因此，2.1節(jié)中的爆炸載荷在花崗巖無(wú)限介質(zhì)中傳播的波長(zhǎng)為

λim=cimT=1.541 4 m，

(5)

式中：T為載荷加載周期。

在花崗巖介質(zhì)一維桿中傳播時(shí)的波長(zhǎng)為

λodb=codbT=1.419 9 m.

(6)

由此看到，應(yīng)力波的波長(zhǎng)由外載荷和傳播介質(zhì)的屬性共同決定。

2.3巖石中波傳播計(jì)算所需的網(wǎng)格尺寸

下面研究當(dāng)網(wǎng)格尺寸取多少時(shí)，巖石中的波傳播計(jì)算結(jié)果趨于穩(wěn)定。

建立一維桿模型，計(jì)算桿中的彈性波傳播問題。在桿的端部施加半周期正弦載荷F=Asinωt，其中，A為振幅，ω為角頻率，頻率f=1/T=ω/2π. 取T=300 μs，A=10 MPa，低于巖石強(qiáng)度，彈性加載。輸入的加載波形如圖5所示。

圖5　輸入的壓力加載波形Fig.5　Loading waveform of input pressure

區(qū)別于大多數(shù)文獻(xiàn)中通過逐漸減小網(wǎng)格尺寸加密網(wǎng)格數(shù)來(lái)獲得網(wǎng)格收斂性結(jié)果的方法，本文將網(wǎng)格尺寸和載荷波長(zhǎng)相關(guān)聯(lián)，研究一個(gè)波長(zhǎng)長(zhǎng)度內(nèi)不同網(wǎng)格數(shù)對(duì)計(jì)算各物理量的影響。取網(wǎng)格尺寸Δl=λ/n，其中n為網(wǎng)格密度數(shù)，表示一個(gè)波長(zhǎng)長(zhǎng)度內(nèi)的網(wǎng)格數(shù)，λ取花崗巖一維桿中的波長(zhǎng)1.419 9 m.n分別取2、4、8、16、20、40和60.

2.3.1網(wǎng)格尺寸對(duì)壓力波形的影響

圖6給出了一個(gè)波長(zhǎng)范圍內(nèi)取不同網(wǎng)格數(shù)計(jì)算得到的壓力波形比較。從圖6中可以看出：壓力的波形和峰值對(duì)網(wǎng)格密度非常敏感，當(dāng)一個(gè)波長(zhǎng)內(nèi)的網(wǎng)格達(dá)到16個(gè)時(shí)，壓力波形和峰值基本趨于穩(wěn)定；隨著網(wǎng)格密度的增大，壓力波形的首峰上升率增大、振蕩次數(shù)減少、振蕩頻率減??；壓力波形以最終穩(wěn)態(tài)值為中心進(jìn)行振蕩，不同網(wǎng)格密度計(jì)算得到的最終壓力穩(wěn)態(tài)值相同。

圖6　　　　一個(gè)波長(zhǎng)范圍內(nèi)取不同網(wǎng)格數(shù)計(jì)算得到的壓力波形比較Fig.6　Comparison of pressure waveforms obtained from different grid number within a wavelength range

2.3.2網(wǎng)格尺寸對(duì)速度波形的影響

圖7給出了一個(gè)波長(zhǎng)范圍內(nèi)取不同網(wǎng)格數(shù)計(jì)算得到的速度v波形比較。從圖7中可以看出：隨著網(wǎng)格密度的增大，速度峰值逐漸增大，速度波形的首峰上升率增大；當(dāng)網(wǎng)格密度達(dá)到16時(shí)，速度峰值和波形基本趨于穩(wěn)定值；隨著網(wǎng)格密度的增大，速度波形的振蕩(次數(shù)和周期)逐漸減少。

圖7　一個(gè)波長(zhǎng)范圍不同網(wǎng)格數(shù)計(jì)算得到的速度波形比較Fig.7　Comparison of velocity waveforms obtained from different grid numbers within a wavelength range

2.3.3網(wǎng)格尺寸對(duì)位移波形的影響

圖8給出了一個(gè)波長(zhǎng)范圍內(nèi)取不同網(wǎng)格數(shù)時(shí)計(jì)算得到的位移D波形比較。從圖8中可以看出：隨著網(wǎng)格密度的增大，位移峰值逐漸減小、首峰上升率增大、位移波形的振蕩周期逐漸減小、振蕩次數(shù)逐漸減少；當(dāng)每個(gè)波長(zhǎng)內(nèi)的網(wǎng)格數(shù)量達(dá)到16時(shí)，位移峰值和波形基本趨于穩(wěn)定；比較不同網(wǎng)格密度的位移波形可以看到，雖然波形不同、振蕩周期不同，但位移的穩(wěn)態(tài)值基本相同(振蕩波形取穩(wěn)態(tài)值為波峰或波谷間的平均值)。由此可認(rèn)為，位移的波形和峰值受網(wǎng)格尺寸影響明顯，但位移穩(wěn)態(tài)值與網(wǎng)格尺寸無(wú)關(guān)，穩(wěn)態(tài)位移值是網(wǎng)格無(wú)關(guān)量。為證明該結(jié)論的正確性，下面進(jìn)行進(jìn)一步計(jì)算驗(yàn)證。

圖8　　　　一個(gè)波長(zhǎng)范圍不同網(wǎng)格數(shù)計(jì)算得到的位移波形比較Fig.8　Comparison of displacement waveforms obtained from different grid numbers within a wavelength range

首先進(jìn)一步減小網(wǎng)格密度(增大網(wǎng)格尺寸)，比較位移波形和穩(wěn)態(tài)值。計(jì)算中增加桿長(zhǎng)、增加計(jì)算時(shí)間，將網(wǎng)格密度進(jìn)一步減小到一個(gè)波長(zhǎng)一個(gè)網(wǎng)格和兩個(gè)波長(zhǎng)一個(gè)網(wǎng)格，比較不同網(wǎng)格密度計(jì)算得到的位移波形，如圖9所示。圖9中標(biāo)出了不同網(wǎng)格密度計(jì)算得到的位移穩(wěn)態(tài)值，從圖中看到，不同網(wǎng)格密度計(jì)算得到的位移穩(wěn)態(tài)值完全相同。

圖9　一個(gè)波長(zhǎng)范圍內(nèi)不同網(wǎng)格數(shù)計(jì)算得到的位移波形比較Fig.9　Comparison of displacement waveforms obtained from different grid numbers within a wavelength range

另外，為了驗(yàn)證加載波類型對(duì)結(jié)論的影響，將加載波由半周期正弦波改為單脈沖方波，對(duì)不同網(wǎng)格密度進(jìn)行計(jì)算驗(yàn)證，計(jì)算得到了同樣的結(jié)論：當(dāng)一個(gè)波長(zhǎng)內(nèi)的網(wǎng)格數(shù)從1個(gè)增加到32個(gè)時(shí)，計(jì)算得到的穩(wěn)態(tài)位移值完全相同；但網(wǎng)格密度不同，計(jì)算到達(dá)穩(wěn)態(tài)位移所需的時(shí)間不同，網(wǎng)格密度越大，計(jì)算達(dá)到穩(wěn)態(tài)位移的時(shí)間越短。該驗(yàn)證說明，獲得的上述結(jié)論和加載波的類型無(wú)關(guān)。

從能量的角度可以解釋穩(wěn)態(tài)位移計(jì)算結(jié)果的網(wǎng)格無(wú)關(guān)性。對(duì)各種不同的網(wǎng)格密度模型輸入的外力載荷相同，因此對(duì)模型做的總功相同。而功等于力與位移的乘積，如果計(jì)算中能量守恒，輸入的外載荷又相同，那么得到的穩(wěn)態(tài)位移值一定相等。

2.3.4巖石中爆炸波傳播計(jì)算所需的網(wǎng)格尺寸

根據(jù)2.3.3節(jié)得到的結(jié)論，當(dāng)一個(gè)波長(zhǎng)內(nèi)的網(wǎng)格數(shù)達(dá)到16個(gè)時(shí)，計(jì)算波形基本趨于穩(wěn)定，計(jì)算結(jié)果基本收斂，可以推算花崗巖無(wú)限介質(zhì)中爆炸波傳播計(jì)算所需的網(wǎng)格尺寸為

(7)

花崗巖介質(zhì)中一維桿波傳播計(jì)算所需的網(wǎng)格尺寸為

(8)

根據(jù)計(jì)算得到的網(wǎng)格尺寸，測(cè)試圍巖自由場(chǎng)計(jì)算精度，結(jié)果如圖10所示，圖中，vmax為速度峰值，d為爆心距。從圖10中看到，10 cm的網(wǎng)格尺寸，基本能滿足計(jì)算精度要求。

圖10　　　　1 t TNT填實(shí)爆炸時(shí)的不同網(wǎng)格尺寸自由場(chǎng)速度峰值比較Fig.10　　　　Comparison of free field peak velocities with different mesh size of 1 t TNT tamped explosion

3　網(wǎng)格尺寸敏感性機(jī)理分析

從2.3.3節(jié)和2.3.4節(jié)的分析看到，對(duì)于網(wǎng)格敏感量(以壓力p為例)，爆心距越小計(jì)算所需的網(wǎng)格尺寸越小，說明在不同的爆心距，物理量對(duì)網(wǎng)格的依賴程度不一樣。文獻(xiàn)[16]在研究空氣和水下爆炸時(shí)指出，采用同一網(wǎng)格尺寸，在炸藥量較大的情況下可以獲得較為精確的數(shù)值結(jié)果，而炸藥量較小時(shí)與經(jīng)驗(yàn)公式計(jì)算值偏離較大，即炸藥當(dāng)量均對(duì)數(shù)值計(jì)算的精度有較大影響——炸藥當(dāng)量越大，可使用的網(wǎng)格尺寸越大。由此產(chǎn)生疑問：當(dāng)量對(duì)計(jì)算網(wǎng)格產(chǎn)生影響的原因是什么？進(jìn)而還產(chǎn)生疑問，為什么計(jì)算相同的物理量在不同爆心距對(duì)網(wǎng)格的尺寸要求不同？為此進(jìn)行了網(wǎng)格尺寸的敏感性機(jī)理分析。

首先分析不同炸藥當(dāng)量產(chǎn)生的差別。從爆炸縮比規(guī)律來(lái)講，壓力波的作用時(shí)間和當(dāng)量的三分之一次方呈正比。所以爆炸當(dāng)量越大，產(chǎn)生的爆炸載荷波長(zhǎng)越長(zhǎng)。而之前的研究結(jié)果表明，當(dāng)一個(gè)波長(zhǎng)內(nèi)的網(wǎng)格數(shù)目大于等于16個(gè)時(shí)，基本能夠準(zhǔn)確模擬各種物理量。波長(zhǎng)越長(zhǎng)，相應(yīng)的模擬所需的網(wǎng)格尺寸也就越大，因此，對(duì)于同一網(wǎng)格尺寸，爆炸當(dāng)量越大，模擬的結(jié)果越精確。

為研究物理量在不同爆心距對(duì)網(wǎng)格尺寸的依賴程度不同，本文采用1 mm大小的網(wǎng)格尺寸計(jì)算得到TNT空中爆炸時(shí)不同比例爆心距R的速度波形(如圖11所示)，并對(duì)其進(jìn)行了頻譜分析(如圖12所示)。從圖12看到，隨著爆心距的增大，速度波形中的高頻分量逐漸減小，而低頻分量相對(duì)逐漸增多。通常來(lái)講，爆炸載荷是由多種頻率的波疊加而成，高頻波波長(zhǎng)短，模擬所需的網(wǎng)格尺寸較小，而高頻波衰減較快，隨著爆心距的增加，高頻波成分逐漸減小，低頻波相對(duì)成分逐漸增加，因此隨著爆心距的增大，加載波的波長(zhǎng)逐漸增加，模擬所需的網(wǎng)格尺寸逐漸增大，因此，隨著爆心距的增大，物理量對(duì)網(wǎng)格尺寸的依賴性逐漸降低。

圖11　TNT爆炸時(shí)不同比例爆心距的速度波形Fig.11　　　　Velocity waveforms at different scaled distances from explosive source during TNT explosion

圖12　TNT爆炸時(shí)不同比例爆心距的速度波形頻譜分析Fig.12　Spectral analysis of velocity waveforms at different scaled distances from explosive source during TNT explosion

圖13為對(duì)某次化爆實(shí)驗(yàn)實(shí)測(cè)的不同爆心距位置的地震速度波形進(jìn)行頻譜分析得到的速度譜圖像。從圖中看到，隨著爆心距的增加，振幅峰值前移，說明隨著爆心距的增加高頻分量耗散的較多，低頻份額相對(duì)增加。由此導(dǎo)致整體的爆炸波載荷頻率降低，波長(zhǎng)增大，從而導(dǎo)致對(duì)計(jì)算網(wǎng)格尺寸的要求也越來(lái)越低，所以會(huì)出現(xiàn)隨著爆心距的增加，計(jì)算網(wǎng)格尺寸的依賴性降低的現(xiàn)象。從理論上講，根據(jù)不同爆心距位置沖擊波的優(yōu)勢(shì)頻率和波長(zhǎng)確定合適的網(wǎng)格尺寸更接近問題的本質(zhì)。然而沖擊波的高頻成分特別豐富，且高頻分量衰減較快，如何確定不同爆心距處的優(yōu)勢(shì)波長(zhǎng)尚存在一定困難。雖然目前可以通過網(wǎng)格漸變的劃分方法來(lái)定性描述物理量隨著爆心距的增大對(duì)網(wǎng)格的依賴性降低的現(xiàn)象，但如何準(zhǔn)確地決定網(wǎng)格漸變率仍然有待深入研究。

圖13　巖石中爆炸時(shí)不同爆心距的速度譜分析Fig.13　Spectral analysis of velocity waveforms at different scaled distances from explosive source during explosion in rock

4　時(shí)間步長(zhǎng)系數(shù)對(duì)計(jì)算結(jié)果的影響

在動(dòng)力學(xué)問題中，時(shí)間步長(zhǎng)大致等于以聲速通過單元最小特征長(zhǎng)度所需的時(shí)間，計(jì)算中一般采用顯式方法進(jìn)行求解。在顯式積分方法中，為保證算法的穩(wěn)定性，積分步長(zhǎng)是由最小單元的尺寸控制的。中心差分法是最常用的顯式方法。中心差分法的臨界時(shí)間步長(zhǎng)取決于單元的特征長(zhǎng)度。在小變形問題中，單元的特征長(zhǎng)度不變，因此在整個(gè)求解過程中可以采用相同的時(shí)間積分步長(zhǎng)。但在大變形問題中，單元的特征長(zhǎng)度不斷減小，因此中心差分法的臨界時(shí)間步長(zhǎng)也不斷減小，需要采用變步長(zhǎng)的中心差分法。LS-DYNA3D程序采用的就是變步長(zhǎng)增量解法。對(duì)于三維實(shí)體單元，其極限時(shí)間步長(zhǎng)采用如(9)式算法[17]：

(9)

式中：Q是體黏性系數(shù)C0和C1的函數(shù)，其表達(dá)式為

(10)

圖14　不同時(shí)間步長(zhǎng)系數(shù)計(jì)算得到的空腔半徑比較Fig.14　　　　Comparison of calculated cavity radii with different time step coefficients

5　結(jié)論

本文通過對(duì)巖石中爆炸應(yīng)力波傳播數(shù)值計(jì)算中的網(wǎng)格尺寸效應(yīng)及物理量對(duì)網(wǎng)格尺寸的敏感性機(jī)理的研究，得到了以下結(jié)論：

1) 合適的網(wǎng)格尺寸應(yīng)根據(jù)載荷特征和波傳播介質(zhì)的屬性來(lái)決定。

2) 網(wǎng)格大小對(duì)壓力、速度和位移波形的影響較大；當(dāng)一個(gè)波長(zhǎng)內(nèi)的網(wǎng)格個(gè)數(shù)達(dá)到16個(gè)以上時(shí)，計(jì)算得到的壓力峰值、速度峰值和位移峰值等參量基本趨于穩(wěn)定。

3) 壓力峰值、壓力上升率、速度峰值等參量是網(wǎng)格密度的正相關(guān)量；壓力脈寬、波形的振蕩次數(shù)、振蕩頻率、波形穩(wěn)定時(shí)間和位移峰值等參量是網(wǎng)格密度的負(fù)相關(guān)量；而波前到達(dá)時(shí)間、正壓作用時(shí)間、正壓沖量、位移穩(wěn)態(tài)值和爆室內(nèi)的穩(wěn)定壓力等參量是網(wǎng)格無(wú)關(guān)量。

4) 隨著爆心距的增加，物理量對(duì)網(wǎng)格尺寸的依賴性降低，其機(jī)理是隨著爆心距的增加，載荷中的高頻成分逐漸衰減、載荷的波長(zhǎng)變大，模擬所需的網(wǎng)格尺寸變大。

5) 步長(zhǎng)系數(shù)對(duì)計(jì)算結(jié)果的影響也非常明顯，當(dāng)時(shí)間步長(zhǎng)系數(shù)取0.05時(shí)，位移穩(wěn)態(tài)值趨于收斂值。

本文通過對(duì)網(wǎng)格尺寸這一關(guān)鍵計(jì)算因素對(duì)計(jì)算結(jié)果影響的研究，了解了哪些物理量是計(jì)算的確定量，哪些物理量是計(jì)算的不確定量；網(wǎng)格尺寸對(duì)各種物理量影響規(guī)律的獲得，對(duì)計(jì)算模型的修正及對(duì)計(jì)算結(jié)果的可靠性把握起了很大的提升作用，使大規(guī)模的定量計(jì)算應(yīng)用于實(shí)際工程問題成為可能。

References)

[1]張雄, 王天舒. 計(jì)算動(dòng)力學(xué)[M]. 北京: 清華大學(xué)出版社, 2007.

ZHANG Xiong, WANG Tian-shu. Computational dynamics[M]. Beijing: Tsinghua University Press, 2007. (in Chinese)

[2]陸金甫, 關(guān)治. 偏微分方程數(shù)值解法[M]. 北京:清華大學(xué)出版社, 2004.LU Jin-fu, GUAN Zhi. Numerical solution of partial differential equation[M]. Beijing: Tsinghua University Press, 2004. (in Chinese)

[3]李順波, 東兆星, 齊燕軍, 等. 爆炸沖擊波在不同介質(zhì)中傳播衰減規(guī)律的數(shù)值模擬[J]. 振動(dòng)與沖擊, 2009, 28(7):115-117.

LI Shun-bo, DONG Zhao-xing, QI Yan-jun, et al. Numerical simulation for spread decay of blasting shock wave in different media[J]. Journal of Vibration and Shock, 2009, 28(7):115-117.(in Chinese)

[4]胡八一, 柏勁松, 劉大敏, 等. 爆炸容器的工程設(shè)計(jì)方法及其應(yīng)用[J]. 壓力容器, 2000, 17(2):39-41.HU Ba-yi, BAI Jin-song, LIU Da-min, et al. The engineering design method of explosion-containment vessel and its application[J]. Pressure Vessel Technology, 2000, 17(2): 39-41. (in Chinese)

[5]Chapman T C. Blast wave simulation using AUTODYN 2D: a parametric study[J]. International Journal of Impact Engineering, 1995, 16(5): 777-787.

[6]胡八一, 李平, 張振宇, 等.爆炸塔內(nèi)壁特征點(diǎn)的反射壓力數(shù)值模擬[J]. 計(jì)算力學(xué)學(xué)報(bào), 2009, 26(4): 573-578.

HU Ba-yi, LI Ping, ZHANG Zhen-yu, et al. Numerical simulation of characteristic points reflective pressure on the inner surface of the explosion chamber[J]. Chinese Journal of Computational Mechanics, 2009, 26(4):573-578.(in Chinese)

[7]楊鑫, 石少卿, 程鵬飛, 等. 爆炸沖擊波在空氣中傳播規(guī)律的經(jīng)驗(yàn)公式對(duì)比及數(shù)值模擬[J]. 四川建筑, 2007, 27(5): 71-73.

YANG Xin, SHI Shao-qing, CHENG Peng-fei, et al. Numerical simulation and comparison with empirical formula on the law of shock wave generated by TNT explosions in air[J]. Sichuan Architecture, 2007, 27(5): 71-73.(in Chinese)

[8]雷鳴, 田宙, 張海波, 等. 網(wǎng)格疏密對(duì)空中化爆自由場(chǎng)參數(shù)計(jì)算結(jié)果的影響[R]. 西安:西北核技術(shù)研究所, 2005.

LEI Ming, TIAN Zhou, ZHANG Hai-bo, et al. Study on the effects of mesh density to the results of freefield parameters about explosion in air[R]. Xi’an:Northwest Institute of Nuclear Technology, 2005.(in Chinese)

[9]姚成寶, 王宏亮, 張柏華, 等. TNT空中爆炸沖擊波傳播數(shù)值模擬及數(shù)值影響因素分析[J]. 現(xiàn)代應(yīng)用物理, 2014, 5(1): 39-44.

YAO Cheng-bao, WANG Hong-liang, ZHANG Bai-hua, et al. Numerical simulation of shock wave generated by TNT explosions in infinite air[J]. Modern Applied Physics, 2014, 5(1): 39-44.(in Chinese)

[10]王海兵, 曹淵, 張海波, 等. 5kg TNT當(dāng)量爆炸容器力學(xué)響應(yīng)數(shù)值模擬及參數(shù)敏感性分析[C]∥第二屆全國(guó)爆炸容器專題研討會(huì)論文集. 西安:西北核技術(shù)研究所, 2013.

WANG Hai-bing, CAO Yuan, ZHANG Hai-bo, et al. Numerical simulation on dynamic response of 5kg TNT equivalent explosion container and parameters sensitivity analysis[C]∥Proceedings of the Second National Symposium on Explosive Containers. Xi’an: Northwest Institute of Nuclear Technology, 2013. (in Chinese)

[11]Shi Y C, Li Z X, Hao H. Mesh size effect in numerical simulation of blast wave propagation and interaction with structures[J]. Transactions of Tianjin University, 2008, 14(6): 396-402.

[12]梁正召, 王述紅, 唐春安, 等. 非均勻性巖石破裂的網(wǎng)格效應(yīng)[J]. 巖石力學(xué)與工程學(xué)報(bào), 2005, 24(增刊1):5108-5112.

LIANG Zheng-zhao, WANG Shu-hong, TANG Chun-an, et al. Mesh effects on failure processes of heterogeneous rocks[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(S1):5108-5112.(in Chinese)

[13]崔煥平, 崔燕平, 王宗敏. 混凝土非線性有限元分析中的網(wǎng)格尺寸效應(yīng)[J]. 混凝土, 2007(6): 27-29.

CUI Huan-ping, CUI Yan-ping, WANG Zong-min. Mesh size effect in nonlinear finite element analysis of concrete[J]. Concrete, 2007(6): 27-29.(in Chinese)

[14]門建兵, 隋樹元, 蔣建偉, 等. 網(wǎng)格對(duì)混凝土侵徹?cái)?shù)值模擬的影響[J]. 北京理工大學(xué)學(xué)報(bào), 2005, 25(8):659-662.

MEN Jian-bing, SUI Shu-yuan, JIANG Jian-wei, et al. Mesh dependency for numerical simulation of concrete penetration[J]. Transactions of Beijing Institute of Technology, 2005, 25(8):659-662.(in Chinese)

[15]Shayanfar M A, Kheyroddin A, Mirza M S. Element size effects in nonlinear analysis of reinforced concrete members[J]. Computers and Structures, 1997, 62(2): 339-352.

[16]張社榮, 李宏璧, 王高輝, 等. 空中和水下爆炸沖擊波數(shù)值模擬的網(wǎng)格尺寸效應(yīng)對(duì)比分析[J]. 水利學(xué)報(bào), 2015, 46(3):298-306.

ZHANG She-rong, LI Hong-bi, WANG Gao-hui, et al. Comparative analysis of mesh size effects on numerical simulation of shock wave in air blast and underwater explosion[J]. Journal of Hydraulic Engineering, 2015, 46(3):298-306.(in Chinese)

[17]Hallquist J O. LS-DYNA theory manual[M]. Livermore, CA: Livermore Software Technology Corporation, 2006.

Mesh Size Effect and Its Mechanism Research in Numerical Calculation of Rock Dynamics

WANG Hai-bing1,2， ZHANG Hai-bo2， TIAN Zhou2， OU Zhuo-cheng1， ZHOU Gang2

(1.State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, China; 2.Northwest Institute of Nuclear Technology, Xi’an 710024, Shaanxi, China)

In rock dynamics calculation, the mesh size has an important influence on the reliability of numerically calculated results. A numerical experimental method is used to research the mesh size effect and its sensitivity mechanism in the numerical simulation of propagation of explosion stress wave in rock. The research results show that a proper mesh size should be specified both by the load characteristics and the property of wave propagation medium. When the number of mesh is up to 16 within one load wavelength, the waveforms and peak values of all calculated physical quantities trend to be stable.The relationships between the values of physical quantities and the different mesh densities are also presented. With the increase in distance from explosive source, the sensitivity of physical quantities to mesh size decreases. This is because high-frequency components attenuate gradually, and the wavelength of load becomes longer with the increase in distance from explosive source. The time step coefficient also has a great influence on computational results. When the time step coefficient equals to 0.05, the stable value of displacement tends to converge.

ordnance science and technology; numerical calculation; mesh size; time step; sensitivity mechanism

2016-03-07

國(guó)家自然科學(xué)基金項(xiàng)目(91330205)

王海兵(1980—)，男， 博士研究生。E-mail: wanghaibing@nint.ac.cn；

歐卓成(1961—)，男，教授，博士生導(dǎo)師。E-mail: zcou@bit.edu.cn

O383+.1

A

1000-1093(2016)10-1828-09

10.3969/j.issn.1000-1093.2016.10.009