基于量纲分析法的地运动折合位移势规律研究

何 增, 浦锡锋, 王海兵, 王智环,2, 田 宙

(1.西北核技术研究所，陕西 西安 710024；2.清华大学 工程物理系，北京 100084)

0 引言

折合位移势在应力波理论、地下爆炸现象学和地震学等领域具有广泛的应用背景[1]。李孝兰[2]较为系统地介绍了空腔解耦爆炸的基础理论知识，指出折合位移势是一个与距离无关的函数，因而在地震波运动的分析中极为方便。朱号锋等[3]利用一维球对称有限差分数值计算源程序模拟了硬岩中地下强爆炸的震源函数，得到了地下爆炸的折合位移势及其源频谱。肖卫国等[4-5]开展了不同介质和不同方式地下爆炸地震耦合效应的实验测量和数值模拟工作。卢强等[6]基于标准线性固体模型，给出了球面应力波的折合位移势在Laplace域的理论解，并指出折合位移势的稳态值的依赖关系。

可以看出，当前对折合位移势的研究在理论、实验和数值均有涉及。由于折合位移势对载荷特征、本构关系和材料参数的依赖极为复杂，而量纲分析是对复杂问题进行初步分析的有效方法，在爆炸力学等领域得到了广泛应用[7]。如李丽萍等[8]采用量纲分析建立了爆炸冲击波效应靶模型，并开展了实验研究；赵传荣等[9]对影响冲击波压力峰值和脉冲宽度的因素进行量纲分析，建立了经验模型并开展了实例分析；钟巍等[10]基于量纲分析推导了爆炸冲击波作用后钢化玻璃碎片质量与飞散距离的函数关系式，并通过实测数据验证了公式的合理性。

为了快速估算和探索规律，对折合位移势的稳态值进行了量纲分析。通过观察地运动模拟结果的分布特点，给出了一定条件下折合位移势的5个经验估算公式。采用最小二乘法确定了拟合参数，计算结果表明公式的误差约为10%。最后，利用所得公式讨论了空腔半径、压力峰值、作用时间、介质密度、屈服强度和弹性模量对折合位移势的影响规律。

1 数值模拟的基本情况

采用Lagrange结构动力学计算程序[11-12]，计算了无限大弹塑性介质中球形空腔表面作用三角衰减的载荷，得到了介质中应力波衰减的情况，积累了大量地运动的模拟结果。球形空腔用半径R表征；压力载荷用压力峰值p0和作用时间T表征；采用弹性-线性硬化塑性模型描述介质的本构关系，模型参数主要包括：密度ρ、弹性模量E、泊松比μ、屈服强度Y、塑性硬化模量Et。

本文主要关心折合位移势ψ，表1列出了计算中空腔半径R、压力峰值p0、作用时间t和屈服强度Y的取值情况。介质密度ρ、弹性模量E、泊松比μ、塑性硬化模量Et均固定。

表1 数值计算中的参数取值

图1 3种屈服强度下的折合位移势曲线

图1给出了空腔半径112.74 cm、压力峰值0.3 GPa、作用时间1 ms情形下的折合位移势曲线。可以观察到，折合位移势在10.0～17.5 ms之间基本保持不变，故数据处理中统一取15 ms时刻的折合位移势为稳态值，记作ψ∞。进一步将折合位移势的稳态区域放大可以发现：对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，因此只有当ψ较大时才能认为在反射波到达前保持稳定。在分析中，将忽略ψ∞<100 cm3(本质上是ψ∞/R3<6.979×10-5，见下文)的模拟结果，并在不引起混淆时简称ψ∞为折合位移势。

2 数值模拟的基本情况

2.1 无量纲量的确定

基于线性代数的程序化无量纲量求解方法进行量纲分析[7,10]。决定折合位移势ψ∞的主要物理量包括：空腔半径R；表征压力载荷的压力峰值p0和作用时间t；表征介质特性的密度ρ、弹性模量E、泊松比μ、屈服强度Y、塑性硬化模量Et。该问题所涉及物理量的单位和量纲总结在表2中。

表2 相关物理量的单位和量纲

由表2可知，该问题共有9个变量，涉及3个基本量纲(L，M和T)，故可形成6个无量纲量。根据问题的特点，选择球形空腔半径R、载荷作用时间t和介质弹性模量E组成基本量对其他物理量进行无量纲化，即参考量纲矩阵A与原始量纲矩阵B分别为

(1)

(2)

根据转换矩阵的计算公式M=BA-1，可得

(3)

(4)

2.2 无量纲表达式的确定

(5)

在数值模拟中，由于弹性模量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的计算结果。可以明显观察到，y与x1呈正相关关系，其物理含义是折合位移势随压力峰值的增大而增大。进一步，当p0= 0时，折合位移势应取零，即要求f1(x1=0)=0。

图2 y与x1之间的依赖关系f1(固定Y= 0.105 GPa，R = 112.74 cm)

f12和f13的参数也有类型现象，这表明y与x2呈正相关关系，但增速低于正比。

表3 3种形式f1的参数拟合结果

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的计算结果(p0≤0.1 GPa的折合位移势过小，故舍掉)，其中空腔半径R和作用时间t取表1中所有值。可以明显观察到，y与x2呈正相关关系，但增速低于正比，与上文的结论一致。进一步，当t= 0时，折合位移势应取零，即要求

f2(x2=0)=0

(7)

图3 y与x2之间的依赖关系f2(x2)(固定Y= 0.105 GPa)

表4 3种形式f2的参数拟合结果

另一方面，对模拟结果y和x2取对数，见图4。图4中还给出了抛物线函数lnf21=b1′+b2′lnx2+b2′(lnx2)2、幂函数lnf22=b1′(lnx2)b2′与指数函数lnf23=b1′+b2′eb3′ln x2的拟合情况，相应的参数见表5。可以看出，3个函数的拟合度都很好，且幂函数的指数b2′和指数函数的系数b3′的数值比较稳定。(lnx2,y)和(x2, lny)的拟合结果与(lnx2, lny)类似，根据形式简单原则，不再列出。

图4 ln y与ln x2之间的依赖关系ln f2(固定Y=0.105 GPa)

表5 3种形式ln f2的参数拟合结果

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的计算结果。可以清晰地看出，lny与x3符合线性关系

(8)

图5 ln y与x3之间的依赖关系ln f3(固定Y= 2.0 GPa，R=112.74 cm)

其物理含义是折合位移势随着屈服强度的增大而指数衰减。表6列出了相应的参数拟合结果。

表6 f3的参数拟合结果

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

上文分别确定了f1、f2和f3的形式，总结一下：f1(x1)取抛物线函数与幂函数较佳；f2(x2)取指数函数(无b1′+b2′=0约束)，或者lnf2(lnx2)取抛物线函数、幂函数或指数函数较佳；f3(x3)取指数函数，即 lnf3(x3)取线性函数较佳。

假设在所关注的范围内，各自变量对因变量的影响是彼此独立的，即x1、x2和x3并没有耦合在一起。此外，根据简洁性原则，可考虑正常形式(x,y)和对数形式(lnx, lny)2种类型的表达式。综上，考虑各自变量的函数可写作

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 待定系数的求解

上述共计5个y=f(x)的表达式，每个表达式含5个待定系数a=(a1,…,a5)T。使用最小二乘法来确定a的具体数值。

为了同时描述2种表达式类型，用(X,Y)统一表示模拟情况：对类型一，(X,Y) = (x,y)；对类型二，(X,Y) = (lnx, lny)。使用下标i表示第i次的模拟值，则理论值和模拟值的差值记作

di(a)=f(Xi)-Yi

(16)

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

(17)

为此，要求D(a)对a的偏导数为零，即

(18)

式(18)是关于aj(j= 1,…,5)的非线性方程组，利用方程组的牛顿迭代法来解决该求根问题。记G(a)=(g1,…,g5)T，则牛顿迭代法的计算公式为

(19)

其中，G′(a)为向量函数G(a)的导数，即雅克比矩阵，其元素为

(20)

3 结果与讨论

3.1 y≥yc的结果

选取如下计算条件下的模拟结果确定参数：(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的模拟结果进行拟合，共计203组。

牛顿迭代法的关键是提供合适的初值，而上文已对此作了初步拟合。表7给出了待定参数的收敛值和平均相对误差。误差大于50%的个数均为15个。

表7 不同形式的f在y≥yc的拟合结果

以式(12)为例，图6(a)比较了无量纲折合位移势的模拟值与计算值。在y3×10-4时，拟合值与模拟值较为接近；在y10-4时，拟合值明显偏高。其他公式的结果是类似的。由于f1和f3的形式是固定的(不含式(11)的f1)，这也表明目前无法单纯根据模拟结果和量纲分析来确定f2的最佳形式。

3.2 y≥3yc的结果

选取y≥3yc的模拟结果(共计177组)重新进行最小二乘法拟合，相应结果如图6(b)和表8所示。误差大于50%的个数均为2个。与y≥yc的结果相比，y≥3yc最主要改变的是相对误差与误差超过50%的个数大幅减小，这也说明本文的表达式均不适合y值小的情形。

图6 无量纲折合位移势的模拟值与计算值的比较

表8 不同形式的f在y≥3yc的拟合结果

3.3 折合位移势的规律性

以式(12)为例讨论折合位移势的规律性。代入表8中的参数，可得折合位移势的经验公式为

(21)

(22)

为了便于讨论，对式(22)取对数，可得

(23)

3.3.1 半径的影响

对lnψ∞关于空腔半径R求偏导数，得

(24)

图7 折合位移势随空腔半径的变化曲线

3.3.2 载荷的影响

(1)压力峰值的影响。对lnψ∞关于压力峰值p0求偏导数，得

(25)

这表明折合位移势随压力峰值升高而增加。图8(a)展示了在R=150 cm、t=1.0 ms、ρ=2.46 g/cm3、E=210 GPa、Y=0.120 GPa条件下折合位移势随压力峰值的变化曲线及适用范围。

(2)作用时间的影响。对lnψ∞关于作用时间t求偏导数，得

(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 介质的影响

(1)密度的影响。对lnψ∞关于介质密度ρ求偏导数，得

(27)

这表明折合位移势随介质密度增加而减小。图9(a)展示了在R=150 cm、p0=0.5 GPa、t=1.0 ms、ρ=2.46 g/cm3、E= 210 GPa、Y=0.120 GPa条件下折合位移势随介质密度的变化曲线及适用范围。

(2)屈服强度的影响。对lnψ∞关于屈服强度Y求偏导数，得

(28)

这表明折合位移势随屈服强增强而减小。图9(b)展示了在R=150 cm、p0=0.5 GPa、t=1.0 ms、ρ=2.46 g/cm3、E= 210 GPa条件下折合位移势随屈服强度的变化曲线及适用范围。

(3)弹性模量的影响。对lnψ∞关于弹性模量E求偏导数，得

(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 折合位移势随介质参数的变化曲线

4 结论

本文只是使用量纲分析的手段对折合位移势的规律性进行初探，还可从以下2个方面进行拓展。首先，由于计算程序的限制，模型较为简单，个别参数的取值与典型值差距较大，后续可使用更复杂的模型、更准确的参数重新计算和分析。其次，对于地下爆炸过程涉及到的其他物理量，如折合速度势、准静态压力等，也可尝试利用量纲分析的方法进行处理。