细观非连续介质的应力波传播研究

袁良柱; 陈美多; 谢雨珊; 陆建华; 王鹏飞; 徐松林

doi:10.11883/bzycj-2023-0365

细观非连续介质的应力波传播研究

doi: 10.11883/bzycj-2023-0365

袁良柱^1,,
陈美多¹,
谢雨珊¹,
陆建华¹,
王鹏飞¹,
徐松林^{1, 2, ,}

1.
中国科学技术大学中国科学院材料力学行为和设计重点实验室，安徽合肥 230027
2.
中国地震局地震预测研究所高压物理与地震科技联合实验室，北京 100036

基金项目: 国家自然科学基金 (11672286, 11872361)；高压物理与地震科技联合实验室室开放基金 (2019HPPES01)；中石油与中科院重大战略合作项目(2015A-4812)；中央高校基本科研业务费专项资金(WK2480000008)

详细信息

作者简介:
袁良柱（1998－　），男，博士研究生，ylzustcedu@mail.ustc.edu.cn

通讯作者:
徐松林（1971－　），男，博士，研究员，博士生导师，slxu99@ustc.edu.cn

中图分类号: O347.3
计量
- 文章访问数: 351
- HTML全文浏览量: 110
- PDF下载量: 150
- 被引次数: 0
出版历程
- 收稿日期: 2023-10-09
- 修回日期: 2023-12-14
- 网络出版日期: 2024-02-04
- 刊出日期: 2024-09-20

Investigation on stress wave propagation in mesoscopic discontinuous medium

1.
CAS Key Laboratory for Mechanical Behavior and Design of Materials, University of Science and Technology of China, Hefei 230027, Anhui, China
2.
United Laboratory of High-Pressure Physics and Earthquake Science, Institute of Earthquake Forecasting, China Earthquake Administration, Beijing 100036, China

摘要

摘要: 固体介质，如岩石、混凝土、贝壳和多孔材料等均具有细观非连续、宏观连续的特性，揭示这种细观非连续性对材料动力学响应的影响规律，对于材料设计、安全防护等具有重要意义。从广义Taylor公式出发，推导了分数阶定义下的非连续介质的一维波动方程，引入等效分数阶简化了控制方程。利用有限差分法得到了控制方程的数值解，结果表明：控制方程中的等效分数阶阶数越小，计算得到的波形衰减的程度越大。为了验证方程的可靠性，并进一步研究非连续介质的波传播规律，在考虑多孔材料、岩石等介质的结构特征的基础上，基于ABAQUS软件建立了随机多孔介质模型。分析发现：多孔介质的波传播受到介质细观非连续程度、材料属性和输入波脉宽的影响，但对应的等效分数阶阶数只与介质细观非连续程度相关，因此，其可以作为评价非连续介质动态响应的一个依据。等效分数阶阶数随着孔隙率的增加而减小，孔洞相对数量分布大致相同的情况下，其统计关系近似呈线性关系。研究结果可为研究多孔材料、贝壳等细观非连续介质的波动传播提供新思路。
- 分数阶导数 /
- 细观非连续 /
- 波传播 /
- 有限差分 /
- 等效分数阶阶数
Abstract: Solid mediums, like rocks, concretes, shells and porous materials, etc., has the characteristics of microscopic discontinuity and macroscopic continuity. It is of great significance for material design, safety protection and other fields to reveal the influence of the meso-discontinuity on the dynamic response of the material. In this paper, based on the generalized Taylor’s formula under fractional definition, the governing equation of 1-D wave propagation in discontinuous medium is derived. Equivalent fractional order is introduced and the simplified form of the governing equation is presented for easily calculating. By using the finite difference method, the numerical solution of the governing equation is obtained. The influence of equivalent fractional order on wave propagation are analyzed. By the time domain analysis, the smaller the equivalent fractional order, the greater the degree of attenuation of the calculated waveform. By the frequency domain analysis, both high frequency wave and low frequency wave exhibit attenuation, and the attenuation of high frequency wave is higher than that of low frequency wave, which makes the pulse duration of the wave being larger. It is obvious that the equivalent fractional order has a certain relationship with the spatial structure of discontinuous medium. Based on the structural characteristics of some meso-discontinuous medium, e.g., porous materials and rocks, a randomly distributed pores model is established by using ABAQUS to verify the reliability of the governing equation and study the wave propagation of meso-discontinuous medium. The effects of porosity, material properties and input waves on wave propagation are analyzed. The degree of wave attenuation is positively related to the porosity of the medium, and negatively related to the wave velocity and the pulse duration of input wave. However, the equivalent fractional order is only related to the porosity and pore distribution of the discontinuous medium. When the spatial structure of the discontinuous medium remains unchanged, the corresponding equivalent fractional order does not change with the material property and the pulse duration of the input wave. By the randomly distributed pores model with various porosities, it is found that the equivalent fractional order decreases with the increase of porosity. Under the same porosity, the heterogeneity of pore distribution will result in different waveforms, while with the increase of porosity, this difference becomes more obvious, but the corresponding equivalent fractional order only has little difference. The statistical relation between equivalent fractional order and porosity is approximately linear when the pore distribution is almost the same. Compared with the randomly distributed pores medium, the statistical relation between equivalent fractional order and porosity of discontinuous medium with uniform distribution of different porosity shifts upward, indicating that the attenuation effect of random structure on wave is higher than that of uniform structure. This paper provides a new approach to investigate wave propagation in meso-discontinuous medium such as porous materials, rocks, shells, etc. It can be used as a basis to evaluate the dynamic response of discontinuous medium.
- fractional derivative /
- meso-discontinuous /
- wave propagation /
- finite difference method /
- equivalent fractional order

HTML全文

非连续介质广泛存在于日常生活中：泡沫铝、纸蜂窝等多孔材料^[1-4]因具有轻质高强、缓冲减震等优点，被大量应用于航空航天、交通运输等领域；岩石和土壤内部多含裂纹、孔洞等^[5-6]；而贝壳、碳纳米管纤维、肌肉组织等，虽然空间结构在宏观上是连续的，但其细微观结构却是非连续的^[7-8]。

应力波在具有非连续结构的材料中传播时会发生衰减^[9-11]，对具有空间结构特征的材料力学响应的研究^[5,11-17]很多。多孔材料在受到冲击时，其孔隙会被逐渐压实，造成其冲击过程中的局部密实化现象^[18-19]。Reid等^[12]最早提出了一维刚性-塑性-锁定(rigid, perfectly plastic, locking, R-PP-L)冲击波模型来描述这种材料在应力波作用下的动态响应。在此之后，根据工程环境中多种复杂的工况，发展了更多的本构模型，如：理想弹性-理想塑性-刚性(elastic, perfectly-plastic, rigid, E-PP-R)模型^[13]、刚性-非线性塑性硬化(rigid, non-linear plastic hardening, R-NLPH)模型^[14]、率无关-刚性-塑性硬化(rigid, plastic hardening, R-PH)和率相关-刚性-塑性硬化(dynamic, rigid, plastic hardening, D-R-PH)模型^[15]等。大量的数值模拟工作验证了这些模型在解释多孔材料受到冲击时动态响应方面的可靠性^[19-23]。与上述多孔材料不同，岩石类材料属于脆性材料，其内部多存在裂纹、孔洞等不连续结构^[5-6]。刘永贵等^[5]将岩体裂隙作为内边界处理，研究了岩体内部裂纹、孔洞等不连续结构对弹性波的散射。在此基础上，徐松林等^[16]、谭子翰等^[17]讨论了岩体中的谐波入射角、裂纹排列方向、孔隙率等对弹性波散射的影响。而贝壳、碳纳米管、肌肉组织等材料在宏观上可以看作连续结构，但在微观上具有十分规律的空间结构特征^[7-8,11,24]。例如贝壳珍珠层和碳纳米管纤维，其结构均类似砖泥结构，常常用剪滞(shear lag)模型研究其力学响应^[24-25]。

分数阶本构因具有较强的时间记忆效应，近年来被广泛应用于讨论岩土、沥青等多孔材料的蠕变、松弛特性^[26-30]等。汪成贵等^[26]建立了考虑钙质砂颗粒破碎和状态相关特性的分数阶塑性边界面本构模型，能够模拟不同初始密实度和围压条件下钙质砂三轴排水试验结果和颗粒破碎影响下的状态相关行为。颜可珍等^[29]利用改进的五参数分数阶导数模型较好地表征了沥青混合料的动态黏弹函数，如存储模量、损失模量、相位角、动态模量、损耗因子等。袁良柱等^[30]利用分数阶本构推导了复杂介质的波传播控制方程，以此为基础分析了牡蛎壳微观结构对波传播的影响。这些理论在处理材料时往往将介质视为连续介质，再用较复杂的本构模型对其力学响应进行解释，目前尚缺乏从非连续角度出发的理论。连续介质力学在处理问题时，常常应用Taylor展开公式进行增量形式的表达。Zaid等^[31]根据Caputo形式的分数阶导数提出了广义Taylor公式：

$f\left( x \right) = f\left( a \right) + \frac{{{{\left( {x - a} \right)}^\alpha }}}{{\varGamma \left( {\alpha + 1} \right)}}\;{}_a^{\mathrm{C}}D_x^\alpha f\left( a \right) + \frac{{{{\left( {x - a} \right)}^{2\alpha }}}}{{\varGamma \left( {2\alpha + 1} \right)}}\;{}_a^{\mathrm{C}}D_x^{2\alpha }f\left( a \right) + \cdots$

(1)

式中： ${}_a^{\mathrm{C}}D_x^\alpha f\left( x \right)$ 为函数f(x)的Caputo形式的分数阶导数，即 ${}_a^{\mathrm{C}}D_x^\alpha f\left( x \right) = {1 / {\varGamma \left( {1 - \alpha } \right)}}\int_a^x {f'\left( \tau \right){{\left( {x - \tau } \right)}^{ - \alpha }}{\text{d}}\tau }$ ，其中a为积分下限，C表示分数阶导数为Caputo形式，α为分数阶阶数（0<α<1）， $\varGamma \left( {\textit{z}} \right) = \int_0^\infty {{{\text{e}}^{ - t}}{t^{ {\textit{z}} - 1}}{\text{d}}t}$ 为伽马函数。由此可以相应地对非连续介质进行分数阶增量形式的描述。

为此，本文中首先从分数阶定义的广义Taylor公式出发，推导一维非连续介质的弹性波传播方程，利用有限差分法得到其数值解。然后，基于多孔泡沫和岩石等材料的结构特征，利用ABAQUS建立随机多孔数值模型，分析介质细观非连续程度、材料属性和输入波脉宽对介质波传播的影响。最后，配合有限差分法定量分析等效分数阶阶数与介质细观非连续程度、材料属性和输入波脉宽的具体关系。

1. 控制方程

1.1 广义Taylor公式

对于连续介质来说，在考察其微元体x+dx处的物理量ψ(x+dx)（ψ为任一物理量，如应力σ、位移u和速度v等）时，通常考虑用Taylor公式的前两项来得到，即：

$\psi \left( {x + {\text{d}}x} \right) = \psi \left( x \right) + \psi '\left( x \right){\text{d}}x$

(2)

对于非连续介质来说，由于其空间结构的非连续性，经典Taylor展开公式不再满足。Zaid等^[31]根据Caputo形式的分数阶导数提出了广义Taylor公式，取式(1)的前两项，则非连续介质微元体x+dx处的物理量ψ(x+dx)可表示为：

$\psi \left( {x + {\text{d}}x} \right) = \psi \left( x \right) + \frac{{{{\left( {{\text{d}}x} \right)}^\alpha }}}{{\varGamma \left( {\alpha + 1} \right)}}\;{}_x^{\mathrm{C}}D_{x + {\text{d}}x}^\alpha \psi \left( x \right) = \psi \left( x \right) + \frac{{{\text{d}}{x^*}}}{{\varGamma \left( {\alpha + 1} \right)}}\;{}_x^{\mathrm{C}}D_{x + {\text{d}}x}^\alpha \psi \left( x \right)$

(3)

1.2 一维细观非连续介质的动力学方程

如图1所示的半无限长非连续介质，其左端受外力σ₀。在不考虑其横向应力、变形等，将其作为一维问题考虑的情况下，取x处的微元体，则t时刻微元体左、右两侧的应力分别为σ(x，t)和σ(x+dx，t)。根据牛顿第二定律，应有：

图 1 非连续介质微元体受力和变形的示意图

Figure 1. Schematic diagram of force and deformation of discontinuous medium

下载: 全尺寸图片幻灯片

$\left[{\sigma \left( {x + {\text{d}}x,t} \right) - \sigma \left( {x,t} \right)} \right]S = {\text{d}}M\left( x \right)\frac{{{\partial ^2}u\left( {x,t} \right)}}{{\partial {t^2}}}$

(4)

式中：S为介质的横截面积，u为微元体的位移。dM(x)为微元体的质量，其表达式可写为：

${\text{d}}M\left( x \right) = \rho _0^*S{\text{d}}{x^*}$

(5)

式中： $\rho _0^*$ 定义为介质的分数阶密度，其量纲为[ML^−2−α]。

根据式(3)、(5)，式(4)可改写为：

${}_x^{\mathrm{C}}D_{x + {\text{d}}x}^\alpha \sigma \left( {x,t} \right) = \varGamma \left( {\alpha + 1} \right)\rho _0^*\frac{{{\partial ^2}u\left( {x,t} \right)}}{{\partial {t^2}}}$

(6)

同样地，微元体的变形可表示为：

${\varepsilon ^*}\left( {x,t} \right) = \frac{{u\left( {x + {\text{d}}x,t} \right) - u\left( {x,t} \right)}}{{{{\left( {{\text{d}}x} \right)}^\alpha }}} = \frac{{{}_x^{\mathrm{C}}D_{x + {\text{d}}x}^\alpha u\left( {x,t} \right)}}{{\varGamma \left( {\alpha + 1} \right)}}$

(7)

式中：ε^*为分数阶表达的细观非连续介质应变，其量纲为[L^1−α]。

介质的本构关系采用拟线性表达，即：

$\sigma = {E^*}{\varepsilon ^*}$

(8)

式中：E^*定义为介质的分数阶弹性模量，其量纲为[MT⁻²L^α−2]。

由式(6)～(8)可以得到细观非连续介质的动力学控制方程为：

${}_x^{\mathrm{C}}D_{x + {\text{d}}x}^{2\alpha }\sigma \left( {x,t} \right) = {\left[ {\varGamma \left( {\alpha + 1} \right)} \right]^2}\frac{{\rho _0^*}}{{{E^*}}}\frac{{{\partial ^2}\sigma \left( {x,t} \right)}}{{\partial {t^2}}}$

(9)

1.3 控制方程的简化

Caputo分数阶导数有以下性质：

${}_a^{\mathrm{C}}D_x^\alpha \left[ {\lambda f\left( x \right) + \mu g\left( x \right)} \right] = \lambda \;{}_a^{\mathrm{C}}D_x^\alpha f\left( x \right) + \mu\; {}_a^{\mathrm{C}}D_x^\alpha g\left( x \right)$

(10a)

$\underbrace {{}_a^{\mathrm{C}}D_x^\alpha\; {}_a^{\mathrm{C}}D_x^\alpha \ldots {}_a^{\mathrm{C}}D_x^\alpha }_nf\left( x \right) = {}_a^{\mathrm{C}}D_x^{n\alpha }f\left( x \right)$

(10b)

式中：λ和μ为任一常数，n为正整数。

考虑如图2所示的情况，介质在受力变形后，各处的空间结构发生改变，或者非连续介质本身的空间结构具有非均质性，其各处对应的分数阶阶数不同，根据式(6)～(8)，对于第i个微元体来说，其控制方程为：

图 2 非连续介质的非均质性

Figure 2. Heterogeneity of discontinuous medium

下载: 全尺寸图片幻灯片

${}_{x - {\text{d}}x} {}^{\mathrm{C}}D_x^{2{\alpha _i}}\sigma \left( {x,t} \right) = \varGamma \left( {{\alpha _i} + 1} \right)\rho _0^*\frac{{{\partial ^2}{}_{x - {\text{d}}x} {}^{\mathrm{C}}D_x^{{\alpha _i}}u\left( {x,t} \right)}}{{\partial {t^2}}}$

(11a)

${}_{x - {\text{d}}x} {}^{\mathrm{C}}D_x^{{\alpha _i}}\sigma \left( {x,t} \right) = {E^*}\frac{{{}_{x - {\text{d}}x} {}^{\mathrm{C}}D_x^{2{\alpha _i}}u\left( {x,t} \right)}}{{\varGamma \left( {{\alpha _i} + 1} \right)}}$

(11b)

式中：α_i为第i个微元体的分数阶阶数。

由于每一个微元体的控制方程对应的分数阶阶数均有差异，不利于计算，考虑用一等效分数阶阶数 $\bar \alpha$ 来代表介质整体的分数阶阶数。将第1个微元体至第i个微元体的控制方程相加，得到：

${}_0^{\mathrm{C}}D_x^{\bar \alpha }\sigma \left( {x,t} \right) = \varGamma \left( {\bar \alpha + 1} \right)\rho _0^*\frac{{{\partial ^2}u\left( {x,t} \right)}}{{\partial {t^2}}}$

(12a)

$\sigma \left( {x,t} \right) = {E^*}\frac{{{}_0^{\mathrm{C}}D_x^{\bar \alpha }\;u\left( {x,t} \right)}}{{\varGamma \left( {\bar \alpha + 1} \right)}}$

(12b)

联立式(12a)～(12b)，可以得到方程(9)的简化形式：

${}_0^{\mathrm{C}}D_x^{2\bar \alpha }\sigma \left( {x,t} \right) = {\left[ {\varGamma \left( {\bar \alpha + 1} \right)} \right]^2}\frac{{\rho _0^*}}{E^{\text{*}}}\frac{{{\partial ^2}\sigma \left( {x,t} \right)}}{{\partial {t^2}}}$

(13)

2. 控制方程的数值解

2.1 分数阶方程的有限差分格式

方程(13)中的分数阶导数为Caputo定义形式下的，可以使用L1公式进行差分，对于阶数满足0<α<1.0的情况，有：

${\left. {{}_0^{\mathrm{C}}D_x^\alpha f\left( x \right)} \right|_{x = {x_j}}} = \frac{{\Delta {x^{ - \alpha }}}}{{\varGamma \left( {2 - \alpha } \right)}}\left[ {a_0^{\left( \alpha \right)}f\left( {{x_j}} \right) - \sum\limits_{k = 1}^{j - 1} {\left( {a_{j - k - 1}^{\left( \alpha \right)} - a_{j - k}^{\left( \alpha \right)}} \right)f\left( {{x_k}} \right) - a_{j - 1}^{\left( \alpha \right)}f\left( {{x_0}} \right)} } \right]$

(14)

式中： $a_l^{\left( \alpha \right)} = {\left( {l + 1} \right)^{1 - \alpha }} - {l^{1 - \alpha }}$ ，Δx为空间步长。Gao等^[32]证明了式(14)具有(2−α)阶精度。

对于阶数满足1.0<γ<2.0的情况，令g(x)= $f' (x)$ 和α=γ−1，则：

${}_0^{\mathrm{C}}D_x^\gamma f\left( x \right) = \frac{1}{{\varGamma \left[ {1 - \left( {\gamma - 1} \right)} \right]}}\int_0^x {\frac{{g'\left( {x'} \right)}}{{{{\left( {x - x'} \right)}^{\gamma - 1}}}}{\text{d}}x'} = {}_0^{\mathrm{C}}D_x^\alpha g\left( x \right)$

(15)

根据式(15)，方程(13)的差分格式为：

$\begin{split} \frac{{\Delta {x^{1 - 2\bar \alpha }}}}{{\varGamma \left( {3 - 2\bar \alpha } \right)}}& \Bigg\{ {\frac{{\sigma _{j + 1}^n - \sigma _j^n}}{{\Delta x}} - \sum\limits_{k = 1}^{j - 1} {\left[ {2{{\left( {j - k} \right)}^{2 - 2\bar \alpha }} - {{\left( {j - k - 1} \right)}^{2 - 2\bar \alpha }} - {{\left( {j - k + 1} \right)}^{2 - 2\bar \alpha }}} \right]\frac{{\sigma _{k + 1}^n - \sigma _k^n}}{{\Delta x}}} -} \\ & { \left[ {{j^{2 - 2\bar \alpha }} - {{\left( {j - 1} \right)}^{2 - 2\bar \alpha }}} \right]\frac{{\sigma _1^n - \sigma _0^n}}{{\Delta x}}} \Bigg\} = {\left[ {\varGamma \left( {\bar \alpha + 1} \right)} \right]^2}\frac{{\rho _0^*}}{E^{\text{*}}}\frac{{\sigma _j^{n + 1} - 2\sigma _j^n + \sigma _j^{n - 1}}}{{\Delta {t^2}}} \end{split}$

(16)

式中：Δt为时间步长。

2.2 一维细观非连续介质中的波动传播

方程(13)表明了非连续介质的波传播与等效分数阶阶数 $\bar \alpha$ 有关，图3为 $\bar \alpha$ 分别等于0.999、0.900和0.800时，波传播至0.05、0.10和0.15 m处的波形，其中，输入的边界条件为幅值为100 MPa、脉宽为100 μs的单脉宽正弦波， $\rho _0^*$ 和 $E^{\text{*}}$ 的数值分别取为2700 kg/m^2+α和70 kg/(m^2−α·s²)。

图 3 不同分数阶阶数对应的波传播

Figure 3. Wave propagation with different fractional order

下载: 全尺寸图片幻灯片

从图3(a)可以看出，当 $\bar \alpha$ <1时，应力波的幅值会随着传播距离的增大发生明显的衰减，应力波的脉宽随着传播距离的增大而变大，这符合波在非连续介质中的传播规律^[9]。当 $\bar \alpha$ 逐渐变小时，波的衰减程度变得更加剧烈，通过图3(b)的频率分析可知，高频波和低频波均有损耗，即波在非连续介质中传播时，各频率的波均会反射一部分，使得非连续介质中波的幅值会发生衰减。进一步分析发现，各个 $\bar \alpha$ 计算得到的高频波衰减的程度明显高于低频波，使得非连续介质中波的脉宽会变大。

在图3(a)中，当 $\bar \alpha$ =0.999时，波在传播过程中几乎不会发生变化，注意到当 $\bar \alpha$ =1.0时，方程(13)变为下式中的弹性波传播方程：

$\frac{{{\partial ^2}\sigma \left( {x,t} \right)}}{{\partial {x^2}}} = \frac{{\rho _0^*}}{E^{\text{*}}}\frac{{{\partial ^2}\sigma \left( {x,t} \right)}}{{\partial {t^2}}}$

(17)

即 $\bar \alpha$ =1.0对应了连续介质的波传播。

3. 细观多孔介质的波传播规律

3.1 有限元模型

为了分析等效分数阶阶数与介质非连续结构之间的联系，利用ABAQUS/ EXPLICIT建立了如图4(a)所示的二维有限元模型(finite element model，FEM)。模型由输入杆、输出杆和试样三个部分组成，其中输入杆、输出杆为长1500 mm、高50 mm的细长杆，试样为由多种大小孔洞（直径分别为10、4和1 mm三种尺寸的孔洞）组成的多孔结构（长度可变、高50 mm），通过改变孔洞的数量可获得不同细观非连续程度的介质。模型中，试样由大约40000个以上的S4R壳单元（单元边长为0.2 mm）组成，材料采用简单的线弹性本构，即密度为2700 kg/m³、模量为70 GPa。输入杆、输出杆和试样的材料属性设置为同一种材料，并将输入杆、输出杆和试样绑定，以避免由于杆和试样材料属性不同带来的波的反射。

图 4 细观多孔介质有限元模型及相应的波形

Figure 4. FEM of porous medium and corresponding waveform

下载: 全尺寸图片幻灯片

设定输入波为幅值100 MPa、脉宽100 μs的半脉宽正弦波。本文中，输入波形取离试样500 mm处的单元进行监测，由于试样内部孔洞的随机分布，波在传到输出杆时，沿着杆的横截面的应力是非均匀分布的，因此，输出波形统一地取离试样50 mm后界面上所有单元的波形再平均得到。图4(b)显示了孔隙率为11.03%，长度分别为0.05、0.10和0.15 m的试样的波形，分数阶密度和模量分别取为2700 kg/m^2+α和70 kg/(m^2−α·s²)，通过波形对比可以确定其对应的等效分数阶阶数为0.95，数值模拟和理论计算的波形重合度良好（数值模拟中波形下降沿之后出现的波动来自于杆的横向惯性效应），表明空间结构相同的多孔介质对应的等效分数阶阶数是相同的。

3.2 孔隙率对波传播的影响

孔隙率P是块状材料中孔隙体积与材料在自然状态下总体积的百分比，它是评价介质密实程度的一个重要指标。改变孔洞数量，得到孔隙率分别为11.03%、26.06%、43.46%和60.83%的试样，如图5(a)所示，同样设定输入波为幅值100 MPa、脉宽100 μs的半脉宽正弦波，相应的波形如图5(b)所示。随着孔隙率的增加，介质的细观非连续程度增大，对应的波形幅值减小、脉宽增大。这4种孔隙率对应的等效分数阶阶数分别为0.95、0.85、0.70、0.60，数值模拟和理论计算的结果对比如图5 (b)所示。值得注意的是，图5(a)中介质的内部孔洞实际上是非均匀分布的，其分数阶阶数沿其长度方向也是非均匀的。当介质的孔隙率较小（孔隙率为11.03%、26.06%）时，介质内部孔洞对波传播的影响较小，理论计算的结果与模拟结果一致性较好；当介质的孔隙率较大（孔隙率为43.46%、60.83%）时，介质内部孔洞对波传播的影响显著，其非均匀分布的作用不可忽略，理论计算结果也与模拟结果在下降沿段出现差异。观察图5(b)中孔隙率与等效分数阶阶数之间的关系可以发现，分数阶阶数是随着介质孔隙率的增大而减小的，这表明等效分数阶阶数可以在一定程度上反映非连续介质的空间结构特性。

图 5 不同孔隙率的细观多孔介质的波传播及对应的分数阶阶数

Figure 5. Wave propagation of porous medium with different porosities and corresponding fractional order

下载: 全尺寸图片幻灯片

在孔隙率相同的情况下，各个尺寸的孔洞在介质中的分布差异也会带来波传播情况的不同。当孔隙率较小时，波的衰减程度很小，因此，孔洞分布对波衰减的作用还不明显，当孔隙率较大时，这种由于孔洞分布引起的结构差异开始对波形产生影响。图6显示了不同孔隙率下、不同孔洞分布的介质波传播的波形，不同孔隙率下的波形均因孔洞分布的不一致产生差异，并且随着孔隙率的升高，这种差异越来越明显，但与之相对应的等效分数阶阶数的差异不大。

图 6 各孔隙率下不同孔洞分布的波形

Figure 6. Waveforms of different hole distribution under different porosity

下载: 全尺寸图片幻灯片

为进一步分析孔洞分布不均对波传播的影响，对5种不同孔隙率（19.58%、30.83%、39.25%、43.27%和58.84%）下均匀分布孔洞的细观模型进行了数值模拟，波传播的波形及其对应的分数阶阶数如图7所示。对比图5和图7可以发现，孔隙分布不均会对输出波形结果产生影响，在同样的孔隙率下，均匀分布孔隙输出的波形比随机分布孔隙输出的波形衰减更小。理论模型和数值模拟结果下降沿的偏差，随机分布孔洞出现在孔隙率为43.46%时，均匀分布孔洞出现在孔隙率为58.84%时，因此，高孔隙率下理论模型和数值模拟结果的偏差一部分来自于细观数值模型非均匀分布的影响。另一方面，由于试样内部孔洞的不均匀分布，其分数阶阶数沿长度方向是不均匀的，即使是均匀分布孔隙的数值模型，其分数阶阶数也并不是处处相等，因此，高孔隙率下理论模型和数值模拟结果偏差的另一部分来自于理论模型分数阶阶数均匀化的影响。

图 7 均匀分布的细观多孔介质的波传播及对应的分数阶阶数

Figure 7. Wave propagation of porous medium with uniformly distributed pores and corresponding fractional order

下载: 全尺寸图片幻灯片

为了探讨孔隙率与等效分数阶阶数的具体关系，统计了多种孔隙率下对应的等效分数阶阶数，如表1所示，相应的统计关系如图8所示，细观多孔介质的等效分数阶阶数近似与介质的孔隙率成线性关系，通过拟合可以得到其具体关系为：

表 1 不同孔隙率的细观多孔介质对应的等效分数阶阶数

Table 1. Equivalent fractional order of porous medium with different porosity

编号	孔隙率/%	等效分数阶阶数	编号	孔隙率/%	等效分数阶阶数	编号	孔隙率/%	等效分数阶阶数
1	11.02	0.95	8	34.77	0.81	15	67.69	0.54
2	17.37	0.95	9	43.46	0.70	16	73.80	0.50
3	17.37	0.95	10	43.46	0.71	17	19.58	0.95
4	26.06	0.85	11	52.13	0.68	18	30.83	0.90
5	26.07	0.85	12	52.13	0.66	19	39.25	0.83
6	34.76	0.78	13	60.83	0.60	20	43.27	0.83
7	34.77	0.78	14	60.83	0.58	21	58.84	0.71

下载: 导出CSV

| 显示表格

图 8 等效分数阶阶数与细观多孔介质孔隙率之间的统计关系

Figure 8. Statistical relationship between the equivalent fractional order and the porosity of the porous medium

下载: 全尺寸图片幻灯片

$\bar \alpha = - 0.65P + 1$

(18)

式中：P为孔隙率， $\bar \alpha$ 为孔隙率对应的等效分数阶阶数。相较于随机分布的细观多孔介质，孔隙率不同的均匀分布的细观多孔介质的分数阶阶数与孔隙率之间的统计关系整体向上偏移，表明非均匀结构对波的衰减作用要高于均匀分布的结构。当孔隙率较大时，由于孔洞分布不均匀，因此产生的波形差异较大，对应的等效分数阶阶数差异也较大，因此式(18)不适用于孔隙率较大的情况。此外，当孔隙率达到73.80%时，对应的等效分数阶阶数约为0.50，此时，数值模拟和理论计算的波形也出现了较大的偏差。

3.3 材料属性和输入波脉宽对波传播的影响

材料属性是影响介质力学响应的重要物理参数。为研究材料属性对于非连续介质波传播的影响，选取了4种材料（钢、铝、铜和环氧树脂）进行数值模拟，其材料参数见表2。采用相同的试样(孔隙率为43.46%)进行模拟以避免介质细观非连续程度带来的影响，同样设定输入波为幅值100 MPa、脉宽100 μs的半脉宽正弦波，得到波传播至0.05 m处的波形，如图9所示。由图9可知，材料属性对波在非连续介质中的传播会产生影响，呈现波速越小，波衰减程度越大的趋势。

表 2 不同材料的物理参数

Table 2. Physical parameters of different materials

材料	密度/（kg·m⁻³）	弹性模量/GPa	波速/（m·s⁻¹）
钢	7800	210	5189
铝	2700	70	5092
铜	8960	120	3660
环氧树脂	1200	3	1581

下载: 导出CSV

| 显示表格

图 9 材料属性对波传播的影响

Figure 9. Effect of material properties on wave propagation

下载: 全尺寸图片幻灯片

输入波的频率或者说脉宽是影响非连续介质波传播的一个重要因素^[5,16-17]。选取50、100和200 μs等3种脉宽的输入波进行模拟，同样地，采用相同的试样(孔隙率为43.46%)进行模拟以避免介质细观非连续程度带来的影响，同样设定输入波为幅值100 MPa、脉宽100 μs的半脉宽正弦波，得到了波传播至0.05 m处的波形，如图10所示。由图10可知，输入波的脉宽对波的传播会有影响，脉宽越小，波衰减的程度越大，但对于3种脉宽的输入波，其对应的等效分数阶阶数相同。

图 10 输入波脉宽对波传播的影响

Figure 10. Effect of pulse duration of input wave on wave propagation

下载: 全尺寸图片幻灯片

图9～10说明了材料属性和输入波脉宽均会对波在非连续介质中的传播产生影响，但介质对应的等效分数阶阶数与材料属性和输入波关系较小，在相同的等效分数阶阶数下，计算得到的波形与模拟得到的波形相匹配，说明细观多孔介质的等效分数阶阶数只与介质的细观非连续程度相关。

4. 结　论

从分数阶定义的广义Taylor公式出发，推导了一维细观非连续介质中的波传播控制方程，利用有限差分法得到了控制方程的数值解，基于ABAQUS软件建立的二维随机多孔数值模型验证了所建立方程分析问题的可靠性，以此为基础研究了孔隙率、材料属性和输入波脉宽等因素对细观多孔介质波传播的影响，建立了控制方程中的等效分数阶阶数与孔隙率之间的联系，得到的主要结论如下。

(1)得到了一维细观非连续介质中的波传播控制方程。控制方程中的等效分数阶阶数等于1.0时，对应着连续介质的波传播，小于1.0时，对应着空间非连续介质的波传播。

(2)非连续介质的波传播与介质细观非连续程度、材料属性和输入波脉宽相关，具体表现为介质细观非连续程度越高、材料波速越小、输入波脉宽越大，应力波幅值衰减越厉害。

(3)控制方程中的等效分数阶阶数与介质细观非连续程度有关，与材料属性和输入波脉宽无关。以孔隙率作为评价介质细观非连续程度的标准，孔隙率越大，应力波幅值衰减得越厉害，相应的等效分数阶阶数也越小。孔隙率和等效分数阶阶数的统计关系近似为线性关系。

图 1 非连续介质微元体受力和变形的示意图

Figure 1. Schematic diagram of force and deformation of discontinuous medium

下载: 全尺寸图片幻灯片

图 2 非连续介质的非均质性

Figure 2. Heterogeneity of discontinuous medium

下载: 全尺寸图片幻灯片

图 3 不同分数阶阶数对应的波传播

Figure 3. Wave propagation with different fractional order

下载: 全尺寸图片幻灯片

图 4 细观多孔介质有限元模型及相应的波形

Figure 4. FEM of porous medium and corresponding waveform

下载: 全尺寸图片幻灯片

图 5 不同孔隙率的细观多孔介质的波传播及对应的分数阶阶数

Figure 5. Wave propagation of porous medium with different porosities and corresponding fractional order

下载: 全尺寸图片幻灯片

图 6 各孔隙率下不同孔洞分布的波形

Figure 6. Waveforms of different hole distribution under different porosity

下载: 全尺寸图片幻灯片

图 7 均匀分布的细观多孔介质的波传播及对应的分数阶阶数

Figure 7. Wave propagation of porous medium with uniformly distributed pores and corresponding fractional order

下载: 全尺寸图片幻灯片

图 8 等效分数阶阶数与细观多孔介质孔隙率之间的统计关系

Figure 8. Statistical relationship between the equivalent fractional order and the porosity of the porous medium

下载: 全尺寸图片幻灯片

图 9 材料属性对波传播的影响

Figure 9. Effect of material properties on wave propagation

下载: 全尺寸图片幻灯片

图 10 输入波脉宽对波传播的影响

Figure 10. Effect of pulse duration of input wave on wave propagation

下载: 全尺寸图片幻灯片

表 1 不同孔隙率的细观多孔介质对应的等效分数阶阶数

Table 1. Equivalent fractional order of porous medium with different porosity

编号	孔隙率/%	等效分数阶阶数	编号	孔隙率/%	等效分数阶阶数	编号	孔隙率/%	等效分数阶阶数
1	11.02	0.95	8	34.77	0.81	15	67.69	0.54
2	17.37	0.95	9	43.46	0.70	16	73.80	0.50
3	17.37	0.95	10	43.46	0.71	17	19.58	0.95
4	26.06	0.85	11	52.13	0.68	18	30.83	0.90
5	26.07	0.85	12	52.13	0.66	19	39.25	0.83
6	34.76	0.78	13	60.83	0.60	20	43.27	0.83
7	34.77	0.78	14	60.83	0.58	21	58.84	0.71

下载: 导出CSV

表 2 不同材料的物理参数

Table 2. Physical parameters of different materials

材料	密度/（kg·m⁻³）	弹性模量/GPa	波速/（m·s⁻¹）
钢	7800	210	5189
铝	2700	70	5092
铜	8960	120	3660
环氧树脂	1200	3	1581

下载: 导出CSV

参考文献(32)

[1]	王鹏飞, 徐松林, 胡时胜. 不同温度下泡沫铝压缩行为与变形机制探讨 [J]. 振动与冲击, 2013, 32(5): 16–19. DOI: 10.13465/j.cnki.jvs.2013.05.009. WANG P F, XU S L, HU S S. Compressive behavior and deformation mechanism of aluminum foam under different temperature [J]. Journal of Vibration and Shock, 2013, 32(5): 16–19. DOI: 10.13465/j.cnki.jvs.2013.05.009.
[2]	王鹏飞, 徐松林, 李志斌, 等. 高温下轻质泡沫铝动态力学性能实验 [J]. 爆炸与冲击, 2014, 34(4): 433–438. DOI: 10.11883/1001-1455(2014)04-0433-06. WANG P F, XU S L, LI Z B, et al. An experimental study on dynamic mechanical property of ultra-light aluminum foam under high temperatures [J]. Explosion and Shock Waves, 2014, 34(4): 433–438. DOI: 10.11883/1001-1455(2014)04-0433-06.
[3]	WANG P F, XU S L, HU S S. Experimental and numerical study of the effect of micro-structure on the rate-sensitivity of cellular foam [J]. Mechanics of Advanced Materials and Structures, 2016, 23(8): 888–895. DOI: 10.1080/15376494.2015.1047477.
[4]	SHAN J F, XU S L, ZHOU L J, et al. Dynamic fracture of aramid paper honeycomb subjected to impact loading [J]. Composite Structures, 2019, 223: 110962. DOI: 10.1016/j.compstruct.2019.110962.
[5]	刘永贵, 徐松林, 席道瑛, 等. 节理玄武岩体弹性波频散效应研究 [J]. 岩石力学与工程学报, 2010, 29(1): 3314–3320. DOI: CNKI:SUN:YSLX.0.2010-S1-106. LIU Y G, XU S L, XI D Y, et al. Dispersion effect of elastic wave in jointed basalts [J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(1): 3314–3320. DOI: CNKI:SUN:YSLX.0.2010-S1-106.
[6]	席道瑛, 徐松林, 宛新林. 高孔隙岩石局部压缩屈服与帽盖模型 [J]. 地球物理进展, 2015, 30(4): 1926–1934. DOI: 10.6038/pg20150454. XI D Y, XU S L, WAN X L. The cap model and local compressive yield process of high porous rock [J]. Progress in Geophysics, 2015, 30(4): 1926–1934. DOI: 10.6038/pg20150454.
[7]	WEI X D, MOHAMMAD N, HORACIO D E. Optimal length scales emerging from shear load transfer in natural materials: application to carbon-based nanocomposite design [J]. ACS Nano, 2012, 6(3): 2333–2344. DOI: 10.1021/nn204506d.
[8]	MENG X S, ZHOU L C, LEI L, et al. Deformable hard tissue with high fatigue resistance in the hinge of bivalve Cristaria plicata [J]. Science, 2023, 380(6651): 1252–1257. DOI: 10.1126/science.ade2038.
[9]	胡亚峰, 刘建青, 顾文斌, 等. PVDF应力测试技术及其在多孔材料爆炸冲击实验中的应用 [J]. 爆炸与冲击, 2016, 36(5): 655–662. DOI: 10.11883/1001-1455(2016)05-0655-08. HU Y F, LIU J Q, GU W B, et al. Stress-testing method by PVDF gauge and its application in explosive test of porous material [J]. Explosion and Shock Waves, 2016, 36(5): 655–662. DOI: 10.11883/1001-1455(2016)05-0655-08.
[10]	孙晓旺, 李永池, 叶中豹, 等. 新型空壳颗粒材料在人防工程中应用的实验研究 [J]. 爆炸与冲击, 2017, 37(4): 643–648. DOI: 10.11883/1001-1455(2017)04-0643-06. SUN X W, LI Y C, YE Z B, et al. Experimental study of a novel shelly cellular material used in civil defense engineering [J]. Explosion and Shock Waves, 2017, 37(4): 643–648. DOI: 10.11883/1001-1455(2017)04-0643-06.
[11]	YUAN L Z, MIAO C H, XU S L, et al. Stress-wave propagation in multilayered and density-graded viscoelastic medium [J]. International Journal of Impact Engineering, 2023, 173: 104415. DOI: 10.1016/j.ijimpeng.2022.104415.
[12]	REID S R, PENG C. Dynamic uniaxial crushing of wood [J]. International Journal of Impact Engineering, 1997, 19(5/6): 531–570. DOI: 10.1016/S0734-743X(97)00016-X.
[13]	LOPATNIKOV S L, GAMA B A, JAHIRUL H M, et al. Dynamics of metal foam deformation during Taylor cylinder–Hopkinson bar impact experiment [J]. Composite Structures, 2003, 61(1/2): 61–71. DOI: 10.1016/S0263-8223(03)00039-4.
[14]	HANSSEN A G, HOPPERSTAD O S, LANGSETH M, et al. Validation of constitutive models applicable to aluminium foams [J]. International Journal of Mechanical Sciences, 2002, 44(2): 359–406. DOI: 10.1016/S0020-7403(01)00091-1.
[15]	ZHENG Z, WANG C, Yu J, et al. Dynamic stress–strain states for metal foams using a 3 D cellular model [J]. Journal of the Mechanics and Physics of Solids, 2014, 72: 93–114. DOI: 10.1016/j.jmps.2014.07.013.
[16]	徐松林, 刘永贵, 席道瑛, 等. 弹性波在含双裂纹岩体中的传播分析 [J]. 地球物理学报, 2012, 55(3): 944–952. DOI: 10.6038/j.issn.0001-5733.2012.03.024. XU S L, LIU Y G, XI D Y, et al. Analysis of propagation of elastic wave in rocks with double-crack model [J]. Chinese Journal of Geophysics, 2012, 55(3): 944–952. DOI: 10.6038/j.issn.0001-5733.2012.03.024.
[17]	谭子翰, 徐松林, 刘永贵, 等. 含多种尺寸缺陷岩体中的弹性波散射 [J]. 应用数学和力学, 2013, 34(1): 38–48. DOI: 10.3879/j.issn.1000-0887.2013.01.005. TAN Z H, XU S L, LIU Y G, et al. Scattering of elastic waves by multi-size defects in rock mass [J]. Applied Mathematics and Mechanics, 2013, 34(1): 38–48. DOI: 10.3879/j.issn.1000-0887.2013.01.005.
[18]	章超, 徐松林, 王鹏飞, 等. 不同冲击速度下泡沫铝变形和应力的不均匀性 [J]. 爆炸与冲击, 2015, 35(4): 567–575. DOI: 10.11883/1001-1455(2015)04-0567-09. ZHANG C, XU S L, WANG P F, et al. Deformation and stress nonuniformity of aluminum foam under different impact speeds [J]. Explosion and Shock Waves, 2015, 35(4): 567–575. DOI: 10.11883/1001-1455(2015)04-0567-09.
[19]	刘冕, 王根伟, 宋辉, 等. 负梯度泡沫金属中的局部密实化现象 [J]. 高压物理学报, 2020, 34(4): 044204. DOI: 10.11858/gywlxb.20190866. LIU M, WANG G W, SONG H, et al. Phenomenon of local densification in negative graded metal foam [J]. Chinese Journal of High Pressure Physics, 2020, 34(4): 044204. DOI: 10.11858/gywlxb.20190866.
[20]	范东宇, 苏彬豪, 彭辉, 等. 多孔泡沫牺牲层的动态压溃及缓冲吸能机理研究 [J]. 力学学报, 2022, 54(6): 1630–1640. DOI: 10.6052/0459-1879-22-047. FAN D Y, SU B H, PENG H, et al. Research on dynamic crushing and mechanism of mitigation and energy absorption of cellular sacrificial layers [J]. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(6): 1630–1640. DOI: 10.6052/0459-1879-22-047.
[21]	常白雪, 郑志军, 赵凯, 等. 具有恒定冲击荷载的梯度泡沫金属材料设计 [J]. 爆炸与冲击, 2019, 39(4): 041101. DOI: 10.11883/bzycj-2018-0431. CHANG B X, ZHENG Z J, ZHAO K, et al. Design of gradient foam metal materials with a constant impact load [J]. Explosion and Shock Waves, 2019, 39(4): 041101. DOI: 10.11883/bzycj-2018-0431.
[22]	CHANG B X, ZHENG Z J, ZHANG Y L, et al. Crashworthiness design of graded cellular materials: an asymptotic solution considering loading rate sensitivity [J]. International Journal of Impact Engineering, 2020, 143: 103611. DOI: 10.1016/j.ijimpeng.2020.103611.
[23]	蔡正宇, 丁圆圆, 王士龙, 等. 梯度多胞牺牲层的抗爆炸分析 [J]. 爆炸与冲击, 2017, 37(3): 396–404. DOI: 10.11883/1001-1455(2017)03-0396-09. CAI Z Y, DING Y Y, WANG S L, et al. Anti-blast analysis of graded cellular sacrificial cladding [J]. Explosion and Shock Waves, 2017, 37(3): 396–404. DOI: 10.11883/1001-1455(2017)03-0396-09.
[24]	WANG D Y, WANG P F, WU Y F, et al. Temperature and rate-dependent plastic deformation mechanism of carbon nanotube fiber: experiments and modeling [J]. Journal of the Mechanics and Physics of Solids, 2023, 173: 105241. DOI: 10.1016/j.jmps.2023.105241.
[25]	LIU J J, ZHU W Q, YU Z L, et al. Dynamic shear-lag model for understanding the role of matrix in energy dissipation in fiber-reinforced composites [J]. Acta Biomaterialia, 2018, 74: 270–279. DOI: 10.1016/j.actbio.2018.04.031.
[26]	汪成贵, 束善治, 肖杨, 等. 考虑钙质砂颗粒破碎的分数阶边界面本构模型 [J]. 岩土工程学报, 2023, 45(6): 1162–1170. DOI: 10.11779/CJGE20220229. WANG C G, SHU S Z, XIAO Y, et al. Fractional-order bounding surface model considering breakage of calcareous sand [J]. Chinese Journal of Geotechnical Engineering, 2023, 45(6): 1162–1170. DOI: 10.11779/CJGE20220229.
[27]	刘泉声, 罗慈友, 彭星新, 等. 软岩现场流变试验及非线性分数阶蠕变模型 [J]. 煤炭学报, 2020, 45(4): 1348–1356. DOI: 10.13225/j.cnki.jccs.2019.0479. LIU Q S, LUO C Y, PENG X X, et al. Research on field rheological test and nolinear fractional derivative creep model of weak rock mass [J]. Journal of China Coal Society, 2020, 45(4): 1348–1356. DOI: 10.13225/j.cnki.jccs.2019.0479.
[28]	孙逸飞, 沈扬, 刘汉龙. 粗粒土的分数阶应变率及其与分形维度的关系 [J]. 岩土力学, 2018, 39(1): 297–317. DOI: 10.16285/j.rsm.2017.1320. SUN Y F, SHEN Y, LIU H L. Fractional strain rate and its relation with fractal dimension of granular soils [J]. Rock and Soil Mechanics, 2018, 39(1): 297–317. DOI: 10.16285/j.rsm.2017.1320.
[29]	颜可珍, 杨胜丰, 黎国凯, 等. 沥青混合料动态黏弹性分数阶导数模型 [J]. 中国公路学报, 2022, 35(5): 12–22. DOI: 10.19721/j.cnki.1001-7372.2022.05.002. YAN K Z, YANG S F, LI G K, et al. Fractional derivative model for dynamic viscoelasticity of asphalt mixtures [J]. China Journal of Highway and Transport, 2022, 35(5): 12–22. DOI: 10.19721/j.cnki.1001-7372.2022.05.002.
[30]	袁良柱, 陆建华, 苗春贺, 等. 基于分数阶模型的牡蛎壳动力学特性研究 [J]. 爆炸与冲击, 2023, 43(1): 011101. DOI: 10.11883/bzycj-2022-0318. YUAN L Z, LU J H, MIAO C H, et al. Dynamic properties of oyster shells based on a fractional-order model [J]. Explosion and Shock Waves, 2023, 43(1): 011101. DOI: 10.11883/bzycj-2022-0318.
[31]	ZAID M O, NABIL T S. Generalized Taylor’s formula [J]. Applied Mathematics and Computation, 2007, 186: 286–293. DOI: 10.1016/j.amc.2006.07.102.
[32]	GAO G H, SUN Z H, ZHANG H W. A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications [J]. Journal of Computational Physics, 2014, 259: 33–50. DOI: 10.1016/j.jcp.2013.11.017.