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

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

0 引言

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

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

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

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

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

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

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

圖1 3種屈服強度下的折合位移勢曲線

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

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

2.1 無量綱量的確定

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

表2 相關物理量的單位和量綱

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

(1)

(2)

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

(3)

(4)

2.2 無量綱表達式的確定

(5)

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

(6)

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

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

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

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

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

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

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

f2(x2=0)=0

(7)

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

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

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

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

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

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

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

(8)

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

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

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

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

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

假設在所關注的范圍內(nèi)，各自變量對因變量的影響是彼此獨立的，即x1、x2和x3并沒有耦合在一起。此外，根據(jù)簡潔性原則，可考慮正常形式(x,y)和對數(shù)形式(lnx, lny)2種類型的表達式。綜上，考慮各自變量的函數(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)的候選表達式有2個，分別為

(11)

(12)

(13)

(14)

(15)

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

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

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

di(a)=f(Xi)-Yi

(16)

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

(17)

為此，要求D(a)對a的偏導數(shù)為零，即

(18)

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

(19)

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

(20)

3 結(jié)果與討論

3.1 y≥yc的結(jié)果

選取如下計算條件下的模擬結(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。取無量綱位移勢y≥yc=6.979×10-5的模擬結(jié)果進行擬合，共計203組。

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

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

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

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

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

圖6 無量綱折合位移勢的模擬值與計算值的比較

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

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

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

(21)

(22)

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

(23)

3.3.1 半徑的影響

對lnψ∞關于空腔半徑R求偏導數(shù)，得

(24)

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

3.3.2 載荷的影響

(1)壓力峰值的影響。對lnψ∞關于壓力峰值p0求偏導數(shù)，得

(25)

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

(2)作用時間的影響。對lnψ∞關于作用時間t求偏導數(shù)，得

(26)

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

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

3.3.3 介質(zhì)的影響

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

(27)

這表明折合位移勢隨介質(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條件下折合位移勢隨介質(zhì)密度的變化曲線及適用范圍。

(2)屈服強度的影響。對lnψ∞關于屈服強度Y求偏導數(shù)，得

(28)

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

(3)彈性模量的影響。對lnψ∞關于彈性模量E求偏導數(shù)，得

(29)

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

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

4 結(jié)論

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