何 增, 浦錫鋒, 王海兵, 王智環(huán),2, 田 宙

(1.西北核技術(shù)研究所，陜西 西安 710024；2.清華大學(xué) 工程物理系，北京 100084)

0 引言

折合位移勢(shì)在應(yīng)力波理論、地下爆炸現(xiàn)象學(xué)和地震學(xué)等領(lǐng)域具有廣泛的應(yīng)用背景[1]。李孝蘭[2]較為系統(tǒng)地介紹了空腔解耦爆炸的基礎(chǔ)理論知識(shí)，指出折合位移勢(shì)是一個(gè)與距離無(wú)關(guān)的函數(shù)，因而在地震波運(yùn)動(dòng)的分析中極為方便。朱號(hào)鋒等[3]利用一維球?qū)ΨQ有限差分?jǐn)?shù)值計(jì)算源程序模擬了硬巖中地下強(qiáng)爆炸的震源函數(shù)，得到了地下爆炸的折合位移勢(shì)及其源頻譜。肖衛(wèi)國(guó)等[4-5]開(kāi)展了不同介質(zhì)和不同方式地下爆炸地震耦合效應(yīng)的實(shí)驗(yàn)測(cè)量和數(shù)值模擬工作。盧強(qiáng)等[6]基于標(biāo)準(zhǔn)線性固體模型，給出了球面應(yīng)力波的折合位移勢(shì)在Laplace域的理論解，并指出折合位移勢(shì)的穩(wěn)態(tài)值的依賴關(guān)系。

可以看出，當(dāng)前對(duì)折合位移勢(shì)的研究在理論、實(shí)驗(yàn)和數(shù)值均有涉及。由于折合位移勢(shì)對(duì)載荷特征、本構(gòu)關(guān)系和材料參數(shù)的依賴極為復(fù)雜，而量綱分析是對(duì)復(fù)雜問(wèn)題進(jìn)行初步分析的有效方法，在爆炸力學(xué)等領(lǐng)域得到了廣泛應(yīng)用[7]。如李麗萍等[8]采用量綱分析建立了爆炸沖擊波效應(yīng)靶模型，并開(kāi)展了實(shí)驗(yàn)研究；趙傳榮等[9]對(duì)影響沖擊波壓力峰值和脈沖寬度的因素進(jìn)行量綱分析，建立了經(jīng)驗(yàn)?zāi)Ｐ筒㈤_(kāi)展了實(shí)例分析；鐘巍等[10]基于量綱分析推導(dǎo)了爆炸沖擊波作用后鋼化玻璃碎片質(zhì)量與飛散距離的函數(shù)關(guān)系式，并通過(guò)實(shí)測(cè)數(shù)據(jù)驗(yàn)證了公式的合理性。

為了快速估算和探索規(guī)律，對(duì)折合位移勢(shì)的穩(wěn)態(tài)值進(jìn)行了量綱分析。通過(guò)觀察地運(yùn)動(dòng)模擬結(jié)果的分布特點(diǎn)，給出了一定條件下折合位移勢(shì)的5個(gè)經(jīng)驗(yàn)估算公式。采用最小二乘法確定了擬合參數(shù)，計(jì)算結(jié)果表明公式的誤差約為10%。最后，利用所得公式討論了空腔半徑、壓力峰值、作用時(shí)間、介質(zhì)密度、屈服強(qiáng)度和彈性模量對(duì)折合位移勢(shì)的影響規(guī)律。

1 數(shù)值模擬的基本情況

采用Lagrange結(jié)構(gòu)動(dòng)力學(xué)計(jì)算程序[11-12]，計(jì)算了無(wú)限大彈塑性介質(zhì)中球形空腔表面作用三角衰減的載荷，得到了介質(zhì)中應(yīng)力波衰減的情況，積累了大量地運(yùn)動(dòng)的模擬結(jié)果。球形空腔用半徑R表征；壓力載荷用壓力峰值p0和作用時(shí)間T表征；采用彈性-線性硬化塑性模型描述介質(zhì)的本構(gòu)關(guān)系，模型參數(shù)主要包括：密度ρ、彈性模量E、泊松比μ、屈服強(qiáng)度Y、塑性硬化模量Et。

本文主要關(guān)心折合位移勢(shì)ψ，表1列出了計(jì)算中空腔半徑R、壓力峰值p0、作用時(shí)間t和屈服強(qiáng)度Y的取值情況。介質(zhì)密度ρ、彈性模量E、泊松比μ、塑性硬化模量Et均固定。

表1 數(shù)值計(jì)算中的參數(shù)取值

圖1 3種屈服強(qiáng)度下的折合位移勢(shì)曲線

圖1給出了空腔半徑112.74 cm、壓力峰值0.3 GPa、作用時(shí)間1 ms情形下的折合位移勢(shì)曲線?？梢杂^察到，折合位移勢(shì)在10.0～17.5 ms之間基本保持不變，故數(shù)據(jù)處理中統(tǒng)一取15 ms時(shí)刻的折合位移勢(shì)為穩(wěn)態(tài)值，記作ψ∞。進(jìn)一步將折合位移勢(shì)的穩(wěn)態(tài)區(qū)域放大可以發(fā)現(xiàn)：對(duì)Y=0.03 GPa，位移勢(shì)的范圍在606.3 cm3<ψ< 617.3 cm3，波動(dòng)幅度為0.90%；對(duì)Y= 0.3 GPa，位移勢(shì)的范圍在22.98 cm3<ψ<31.87 cm3，波動(dòng)幅度為16.2%；對(duì)Y= 0.5 GPa，位移勢(shì)的范圍在-3.210 cm3<ψ< 6.291 cm3，波動(dòng)幅度為308%?？梢钥闯觯灰苿?shì)的變化幅度Δψ≈10 cm3，因此只有當(dāng)ψ較大時(shí)才能認(rèn)為在反射波到達(dá)前保持穩(wěn)定。在分析中，將忽略ψ∞<100 cm3(本質(zhì)上是ψ∞/R3<6.979×10-5，見(jiàn)下文)的模擬結(jié)果，并在不引起混淆時(shí)簡(jiǎn)稱ψ∞為折合位移勢(shì)。

2 數(shù)值模擬的基本情況

2.1 無(wú)量綱量的確定

基于線性代數(shù)的程序化無(wú)量綱量求解方法進(jìn)行量綱分析[7,10]。決定折合位移勢(shì)ψ∞的主要物理量包括：空腔半徑R；表征壓力載荷的壓力峰值p0和作用時(shí)間t；表征介質(zhì)特性的密度ρ、彈性模量E、泊松比μ、屈服強(qiáng)度Y、塑性硬化模量Et。該問(wèn)題所涉及物理量的單位和量綱總結(jié)在表2中。

表2 相關(guān)物理量的單位和量綱

由表2可知，該問(wèn)題共有9個(gè)變量，涉及3個(gè)基本量綱(L，M和T)，故可形成6個(gè)無(wú)量綱量。根據(jù)問(wèn)題的特點(diǎn)，選擇球形空腔半徑R、載荷作用時(shí)間t和介質(zhì)彈性模量E組成基本量對(duì)其他物理量進(jìn)行無(wú)量綱化，即參考量綱矩陣A與原始量綱矩陣B分別為

(1)

(2)

根據(jù)轉(zhuǎn)換矩陣的計(jì)算公式M=BA-1，可得

(3)

(4)

2.2 無(wú)量綱表達(dá)式的確定

(5)

在數(shù)值模擬中，由于彈性模量E、泊松比μ、塑性硬化模量Et保持不變，因此式(5)可簡(jiǎn)化為

(6)

2.2.1y與x1之間的依賴關(guān)系y=f1(x1)

為了研究y與x1之間的依賴關(guān)系f1(x1)，需固定x2與x3。圖2給出了Y=0.105 GPa，R=112.74 cm，t分別取0.5 ms、1.0 ms和2.0 ms的計(jì)算結(jié)果?？梢悦黠@觀察到，y與x1呈正相關(guān)關(guān)系，其物理含義是折合位移勢(shì)隨壓力峰值的增大而增大。進(jìn)一步，當(dāng)p0= 0時(shí)，折合位移勢(shì)應(yīng)取零，即要求f1(x1=0)=0。

圖2 y與x1之間的依賴關(guān)系f1(固定Y= 0.105 GPa，R = 112.74 cm)

f12和f13的參數(shù)也有類型現(xiàn)象，這表明y與x2呈正相關(guān)關(guān)系，但增速低于正比。

表3 3種形式f1的參數(shù)擬合結(jié)果

2.2.2y與x2之間的依賴關(guān)系y=f2(x2)

為了研究y與x2之間的依賴關(guān)系f2(x2)，需固定x1與x3。圖3給出了Y=0.105 GPa，p0分別取0.2 GPa、0.5 GPa、1.0 GPa和2.0 GPa的計(jì)算結(jié)果(p0≤0.1 GPa的折合位移勢(shì)過(guò)小，故舍掉)，其中空腔半徑R和作用時(shí)間t取表1中所有值?？梢悦黠@觀察到，y與x2呈正相關(guān)關(guān)系，但增速低于正比，與上文的結(jié)論一致。進(jìn)一步，當(dāng)t= 0時(shí)，折合位移勢(shì)應(yīng)取零，即要求

f2(x2=0)=0

(7)

圖3 y與x2之間的依賴關(guān)系f2(x2)(固定Y= 0.105 GPa)

表4 3種形式f2的參數(shù)擬合結(jié)果

另一方面，對(duì)模擬結(jié)果y和x2取對(duì)數(shù)，見(jiàn)圖4。圖4中還給出了拋物線函數(shù)lnf21=b1′+b2′lnx2+b2′(lnx2)2、冪函數(shù)lnf22=b1′(lnx2)b2′與指數(shù)函數(shù)lnf23=b1′+b2′eb3′ln x2的擬合情況，相應(yīng)的參數(shù)見(jiàn)表5?？梢钥闯觯?個(gè)函數(shù)的擬合度都很好，且冪函數(shù)的指數(shù)b2′和指數(shù)函數(shù)的系數(shù)b3′的數(shù)值比較穩(wěn)定。(lnx2,y)和(x2, lny)的擬合結(jié)果與(lnx2, lny)類似，根據(jù)形式簡(jiǎn)單原則，不再列出。

圖4 ln y與ln x2之間的依賴關(guān)系ln f2(固定Y=0.105 GPa)

表5 3種形式ln f2的參數(shù)擬合結(jié)果

2.2.3y與x3之間的依賴關(guān)系y=f3(x3)

為了研究y與x3之間的依賴關(guān)系f3(x3)，需固定x1與x2。圖5給出了R=112.74 cm，p0=2.0 GPa，t分別取0.5 ms、1.0 ms和2.0 ms的計(jì)算結(jié)果?？梢郧逦乜闯觯琹ny與x3符合線性關(guān)系

(8)

圖5 ln y與x3之間的依賴關(guān)系ln f3(固定Y= 2.0 GPa，R=112.74 cm)

其物理含義是折合位移勢(shì)隨著屈服強(qiáng)度的增大而指數(shù)衰減。表6列出了相應(yīng)的參數(shù)擬合結(jié)果。

表6 f3的參數(shù)擬合結(jié)果

2.2.4y=f(x1,x2,x3)的表達(dá)式

上文分別確定了f1、f2和f3的形式，總結(jié)一下：f1(x1)取拋物線函數(shù)與冪函數(shù)較佳；f2(x2)取指數(shù)函數(shù)(無(wú)b1′+b2′=0約束)，或者lnf2(lnx2)取拋物線函數(shù)、冪函數(shù)或指數(shù)函數(shù)較佳；f3(x3)取指數(shù)函數(shù)，即 lnf3(x3)取線性函數(shù)較佳。

假設(shè)在所關(guān)注的范圍內(nèi)，各自變量對(duì)因變量的影響是彼此獨(dú)立的，即x1、x2和x3并沒(méi)有耦合在一起。此外，根據(jù)簡(jiǎn)潔性原則，可考慮正常形式(x,y)和對(duì)數(shù)形式(lnx, lny)2種類型的表達(dá)式。綜上，考慮各自變量的函數(shù)可寫作

y=f(x1,x2,x3)=f1(x1)×f2(x2)×f3(x3)

(9)

或

lny=f(lnx1,lnx2,lnx3)=f1(lnx1)+f2(lnx2)+f3(lnx3)

(10)

(1)類型一：y=f(x1,x2,x3)。y=f(x1,x2,x3)的候選表達(dá)式有2個(gè)，分別為

(11)

(12)

(13)

(14)

(15)

2.2.5 待定系數(shù)的求解

上述共計(jì)5個(gè)y=f(x)的表達(dá)式，每個(gè)表達(dá)式含5個(gè)待定系數(shù)a=(a1,…,a5)T。使用最小二乘法來(lái)確定a的具體數(shù)值。

為了同時(shí)描述2種表達(dá)式類型，用(X,Y)統(tǒng)一表示模擬情況：對(duì)類型一，(X,Y) = (x,y)；對(duì)類型二，(X,Y) = (lnx, lny)。使用下標(biāo)i表示第i次的模擬值，則理論值和模擬值的差值記作

di(a)=f(Xi)-Yi

(16)

最小二乘法是使di的平方和D最小，其中

(17)

為此，要求D(a)對(duì)a的偏導(dǎo)數(shù)為零，即

(18)

式(18)是關(guān)于aj(j= 1,…,5)的非線性方程組，利用方程組的牛頓迭代法來(lái)解決該求根問(wèn)題。記G(a)=(g1,…,g5)T，則牛頓迭代法的計(jì)算公式為

(19)

其中，G′(a)為向量函數(shù)G(a)的導(dǎo)數(shù)，即雅克比矩陣，其元素為

(20)

3 結(jié)果與討論

3.1 y≥yc的結(jié)果

選取如下計(jì)算條件下的模擬結(jié)果確定參數(shù)：(1)R和T取表1中所有值，p0=0.1, 0.2, 0.5, 1.0, 2.0 GPa，Y= 0.105 GPa；(2)p0、T和Y取表1中所有值，R=112.74 cm。取無(wú)量綱位移勢(shì)y≥yc=6.979×10-5的模擬結(jié)果進(jìn)行擬合，共計(jì)203組。

牛頓迭代法的關(guān)鍵是提供合適的初值，而上文已對(duì)此作了初步擬合。表7給出了待定參數(shù)的收斂值和平均相對(duì)誤差。誤差大于50%的個(gè)數(shù)均為15個(gè)。

表7 不同形式的f在y≥yc的擬合結(jié)果

以式(12)為例，圖6(a)比較了無(wú)量綱折合位移勢(shì)的模擬值與計(jì)算值。在y3×10-4時(shí)，擬合值與模擬值較為接近；在y10-4時(shí)，擬合值明顯偏高。其他公式的結(jié)果是類似的。由于f1和f3的形式是固定的(不含式(11)的f1)，這也表明目前無(wú)法單純根據(jù)模擬結(jié)果和量綱分析來(lái)確定f2的最佳形式。

3.2 y≥3yc的結(jié)果

選取y≥3yc的模擬結(jié)果(共計(jì)177組)重新進(jìn)行最小二乘法擬合，相應(yīng)結(jié)果如圖6(b)和表8所示。誤差大于50%的個(gè)數(shù)均為2個(gè)。與y≥yc的結(jié)果相比，y≥3yc最主要改變的是相對(duì)誤差與誤差超過(guò)50%的個(gè)數(shù)大幅減小，這也說(shuō)明本文的表達(dá)式均不適合y值小的情形。

圖6 無(wú)量綱折合位移勢(shì)的模擬值與計(jì)算值的比較

表8 不同形式的f在y≥3yc的擬合結(jié)果

3.3 折合位移勢(shì)的規(guī)律性

以式(12)為例討論折合位移勢(shì)的規(guī)律性。代入表8中的參數(shù)，可得折合位移勢(shì)的經(jīng)驗(yàn)公式為

(21)

(22)

為了便于討論，對(duì)式(22)取對(duì)數(shù)，可得

(23)

3.3.1 半徑的影響

對(duì)lnψ∞關(guān)于空腔半徑R求偏導(dǎo)數(shù)，得

(24)

圖7 折合位移勢(shì)隨空腔半徑的變化曲線

3.3.2 載荷的影響

(1)壓力峰值的影響。對(duì)lnψ∞關(guān)于壓力峰值p0求偏導(dǎo)數(shù)，得

(25)

這表明折合位移勢(shì)隨壓力峰值升高而增加。圖8(a)展示了在R=150 cm、t=1.0 ms、ρ=2.46 g/cm3、E=210 GPa、Y=0.120 GPa條件下折合位移勢(shì)隨壓力峰值的變化曲線及適用范圍。

(2)作用時(shí)間的影響。對(duì)lnψ∞關(guān)于作用時(shí)間t求偏導(dǎo)數(shù)，得

(26)

這表明折合位移勢(shì)隨作用時(shí)間增長(zhǎng)而增加。圖8(b)展示了在R=150 cm、p0= 0.5 GPa、ρ= 2.46 g/cm3、E=210 GPa、Y=0.120 GPa 條件下折合位移勢(shì)隨作用時(shí)間的變化曲線及適用范圍。

圖8 折合位移勢(shì)隨壓力載荷的變化曲線

3.3.3 介質(zhì)的影響

(1)密度的影響。對(duì)lnψ∞關(guān)于介質(zhì)密度ρ求偏導(dǎo)數(shù)，得

(27)

這表明折合位移勢(shì)隨介質(zhì)密度增加而減小。圖9(a)展示了在R=150 cm、p0=0.5 GPa、t=1.0 ms、ρ=2.46 g/cm3、E= 210 GPa、Y=0.120 GPa條件下折合位移勢(shì)隨介質(zhì)密度的變化曲線及適用范圍。

(2)屈服強(qiáng)度的影響。對(duì)lnψ∞關(guān)于屈服強(qiáng)度Y求偏導(dǎo)數(shù)，得

(28)

這表明折合位移勢(shì)隨屈服強(qiáng)增強(qiáng)而減小。圖9(b)展示了在R=150 cm、p0=0.5 GPa、t=1.0 ms、ρ=2.46 g/cm3、E= 210 GPa條件下折合位移勢(shì)隨屈服強(qiáng)度的變化曲線及適用范圍。

(3)彈性模量的影響。對(duì)lnψ∞關(guān)于彈性模量E求偏導(dǎo)數(shù)，得

(29)

在R= 150 cm、p0=0.5 GPa、t=1.0 ms、ρ=2.46 g/cm3、Y=0.120 GPa 條件下，當(dāng)EE*時(shí)，折合位移勢(shì)隨彈性模量增強(qiáng)而減小。圖9(c)展示了相同條件下折合位移勢(shì)隨彈性模量的變化曲線及適用范圍，并標(biāo)注了最大折合位移勢(shì)ψ∞(E*)=3 417 cm3的位置。

圖9 折合位移勢(shì)隨介質(zhì)參數(shù)的變化曲線

4 結(jié)論

本文只是使用量綱分析的手段對(duì)折合位移勢(shì)的規(guī)律性進(jìn)行初探，還可從以下2個(gè)方面進(jìn)行拓展。首先，由于計(jì)算程序的限制，模型較為簡(jiǎn)單，個(gè)別參數(shù)的取值與典型值差距較大，后續(xù)可使用更復(fù)雜的模型、更準(zhǔn)確的參數(shù)重新計(jì)算和分析。其次，對(duì)于地下爆炸過(guò)程涉及到的其他物理量，如折合速度勢(shì)、準(zhǔn)靜態(tài)壓力等，也可嘗試?yán)昧烤V分析的方法進(jìn)行處理。