二维颗粒结构中非线性波的传播和脉冲衰减①

王超科 周振江 程相乐 D.Michael McFarland 卢奂采

(浙江工业大学特种装备制造与先进加工技术教育部重点实验室 杭州 310014)

0 引言

有序颗粒介质对冲击能量捕获和吸收[1-3]和冲击振动抑制[4-5]等,具有潜在的和重要的研究价值。Nesterenko[6]在1983 年发现了脉冲在一维均匀颗粒链中以一种特殊的孤立波形式传播,即Nesterenko孤立波,可以在较长距离内稳定地传播,并且集中了脉冲激励的大部分能量[7-12]。Nesterenko[13]及此后的研究人员[14-15]还发现了如声真空和呼吸子等强非线性动力学现象。

通过改变颗粒的大小、质量和材质等,以形成复合颗粒链[16-20]是进一步研究颗粒链中的波传播特性最常用的方法。仿真和实验研究发现,由轻、重2种颗粒构成的二段式一维复合颗粒链,当孤立波从轻颗粒子链传播到重颗粒子链,通过轻-重颗粒交界面时,大部分能量会在界面处被反射回来,反射波在轻颗粒子链形成一组幅值较小的子波,而透射波在重颗粒子链以小幅值孤立波的形式传播。反之,孤立波从重颗粒子链传播到轻颗粒子链,通过重-轻颗粒交界面时,大部分能量会在界面处被透射。Hong[21]研究了脉冲在由大小相同,质量和材料不同的3 种球形颗粒子链构成无预压缩的一维非均匀颗粒链中的传播,发现部分脉冲会被短时约束在轻珠子链中,称之为颗粒容器,并通过时间滞后的方式释放,使得高幅值孤立波被分散为许多低幅值孤立波,从而降低了冲击力。

众多学者对二维耦合颗粒链进行了研究,发现脉冲波在二维耦合均匀颗粒链中传播时,会出现脉冲等分[22-24]、定向能量转移[25]等现象。另外,Zhang等人[24]在2015 年对二维耦合颗粒结构作了系统的研究,研究了由2 列相同的、交错耦合的一维均匀颗粒链构成的二维耦合颗粒链模型,实验验证了新型混合孤立波-剪切波的存在,并阐明了颗粒链间的能量相互交换是由混合波引起的、并最终会达到能量均分的物理现象。而在他们的研究中,并未对二维耦合颗粒链在冲击能量的吸收和降低传递力方面作具体的考察。

通过二维非均匀颗粒链对脉冲能量的局部约束和延时传递来降低传递力,是颗粒介质动力学研究的重要题目。本文提出了由材质相同、半径不同的3 种不同质量的球形颗粒,即重珠、轻珠和耦合珠,有序排列为2 个不同的二维耦合颗粒链系统、二维均匀颗粒链系统和二维非均匀颗粒链系统,运用分子动力学方法仿真,比较研究了这2 种颗粒链系统中的脉冲波传播特性、能量传递现象以及上述机理对颗粒链末端传递力的影响。

1 数学模型

在图1 所示的二维均匀颗粒链系统中,颗粒间在纵向和横向均无预紧力作用,其中,上、下2 列颗粒链均由N颗直径为R1的重珠紧密排列组成。由初始速度为V0冲击激励的颗粒链为激励链,另一条链则为吸收链。激励链与吸收链之间通过嵌入(N-1) 颗直径为R2的耦合珠提供耦合性,耦合珠与相邻重珠满足相切的几何关系。系统的上侧、下侧和右侧为刚性壁边界条件。

图1 二维均匀颗粒链系统

颗粒之间的接触类型和颗粒与边界之间的接触类型,主要为点与点接触和点与面接触,相邻颗粒间的接触力可以用赫兹接触定律[26-27]来描述,表达式为

其中,k是赫兹接触系数,δ是2 个颗粒之间的变形量。不同几何形状的颗粒间的接触定律和参数详情见文献[26,27]。

当2 个球形颗粒接触时,式(1)中的接触系数为

其中,Rs1、Es1、vs1和Rs2、Es2、vs2分别为2 个颗粒的半径、弹性模量、泊松比。当颗粒与刚性壁接触时,假设颗粒的半径为Rs1,刚性壁的曲率半径为Rs2且Rs2→∞,式(1)中的接触系数为

式(4)中,α为重珠和重珠之间的接触系数,η为重珠和耦合珠之间的接触系数,均由式(2)可得。由于边界条件不同,激励链首尾颗粒的运动控制方程,不同于式(4)。激励链首尾颗粒、吸收链和耦合链中颗粒的运动控制方程,不在此赘述,详见附录A。

通过四阶龙格-库塔法,求解颗粒的控制方程,得到每个颗粒的位移和速度时间序列。求解过程中引入以下无量纲计算:

其中,ζ是时间t和初始速度激励V0的无量纲化比例因子

2 仿真计算和分析

2.1 二维均匀颗粒链系统

在如图1 所示的二维均匀颗粒链系统中,设定激励链和吸收链均由N=20 颗重珠组成,所有颗粒的材料参数相同,即弹性模量E=193 MPa,密度ρ=7930 kg/m3,泊松比v=0.3。重珠的半径为R1=12.5 mm,耦合珠的半径为5.18 mm。激励的无量纲初始速度为=1 ×10-4。

图2中给出了激励链和吸收链中颗粒的水平方向和垂直方向速度分量的时间序列响应,不同颜色的曲线代表激励链和吸收链中不同颗粒的速度曲线。从图中可以看出,有2 种不同类型波的传播,分别沿水平方向和垂直方向传播。

图2 激励链与吸收链中颗粒的速度时间序列

如图2(a)所示,激励链和吸收链中颗粒的水平速度分量均形成了孤立脉冲波(solitary pulse waves,简称P 波)。这2 个孤立波都类似于在一维均匀颗粒链中产生的Nesterenko 孤立波,是由颗粒链水平方向的动能传递形成的。由图2(b)所示,由于耦合珠将激励链和吸收链耦合起来,2 珠链之间存在着剧烈的能量交换。两个孤立波在前12 颗重珠中传播时的幅值发生着变化,激励链中的孤立波幅值逐渐减小,吸收链中的孤立波幅值逐渐增大,到达约为第13 颗重珠(图中线宽较粗的曲线)时发生脉冲等分[22-24],即激励链和吸收链中传播的孤立波的幅值趋于相同。

垂直速度分量波形如图2(b)所示,为剪切波(shear waves,简称S 波)。颗粒的横向高频率振荡(与x轴传播的P波的近零频率形成对比)是由其与固定刚性边界之间的赫兹接触相互作用产生的,导致颗粒的运动形成强烈的非线性剪切型波。图3给出了激励链和吸收链中颗粒的水平和垂直速度分量在空间和时间上的变化,从中能更好地看到非线性剪切波的形成。如图3 所示,这些非线性剪切波是近场运动,即它们在空间上靠近激励施加点的地方,并在最初的瞬变过程中形成,同时P 波在2 条链中传播[29-30]。因此,非线性剪切波是由颗粒的水平和垂直运动的非线性耦合在运动的初始高能阶段(即非线性效应更明显时)产生的,并且不能传递到系统的远场。

图3 激励链与吸收链中颗粒的速度时空图

图4 给出了激励链和吸收链中颗粒的动能随时间和空间的变化,从中能更清楚地观察到能量交换和脉冲等分现象。从图中可看出,在激励施加之后,由于耦合珠的存在,图4(a)激励链中颗粒的动能逐渐减小,图4(b)吸收链中颗粒的动能逐渐增大,最终2 链中颗粒的能量大小近似相同,即能量从激励链转移到吸收链,最终导致2 条链之间的能量均衡分配。

图4 激励链与吸收链中颗粒的动能时空图

在二维均匀颗粒链系统中,脉冲能量在初始瞬变后发生等分,将在整个孤立波传播过程中保持不变,使得每条颗粒链独立产生等幅值的孤立波传播现象。在施加冲击之后存在一个初始阶段,即脉冲能量从初始瞬变到等分这一阶段,在此期间产生混合型波的传播,从而在2 条链之间产生剧烈的能量交换,随后脉冲在2 条链中各自沿x方向等幅值传播,并且未衰减地传播到远场。

综上所述,由于二维均匀颗粒链系统的几何结构,在激励施加点附近,接近零频率的纵向传播的P波和横向高频震荡的S 波相互耦合共同形成一种混合波传播模式[24]。P 波构成了颗粒链系统沿水平方向的动量与能量传递的主要机制(类似于一维均匀颗粒链中的Nesterenko 孤立波),但是在垂直方向上,S 波提供了脉冲能量的部分高频散射的机制。因此,这些剪切波在减少最终传递到颗粒链系统远场的脉冲能量中起着重要作用。

传递到二维均匀颗粒链系统中激励链末端的脉冲能量衰减的主要原因包括:(1)由于耦合珠的存在,激励链和吸收链之间会发生能量交换,能量从激励链转移到吸收链;(2)一部分脉冲能量会通过颗粒横向高频振荡被约束在激励施加点附近。因此,与一维均匀颗粒链系统[6]相比,二维均匀颗粒链系统中激励链末端的传递力峰值降低到约100 N,如图5 所示。

图5 一维均匀颗粒链系统和二维均匀颗粒链系统中激励链末端传递力

2.2 二维非均匀颗粒链系统

在引言中提到,通过改变颗粒的大小、质量和材质,以形成复合颗粒链[16-20]系统,改变一维颗粒链中的波传播特性,以减弱颗粒链末端传递力峰值[21]。因此,本文引入轻珠子链构成二维非均匀颗粒链系统,如图6所示。本节考虑在冲击激励作用下,混合P 波和S 波的传播过程,对末端传递力进行了数值仿真,并与2.1 节中二维均匀颗粒链系统的结果进行对比。

图6 二维非均匀颗粒链系统示意图

设定二维非均匀颗粒链系统中激励链和吸收链均由N1=7,N2=6 和N3=N4=5 的4 段颗粒链构成,激励链和吸收链中的颗粒总数目为N=23,非均匀颗粒链的总长度和2.1 节中均匀颗粒链保持一致。所有颗粒,即重珠、轻珠和耦合珠的材料参数也与2.1 节中相同。重珠的半径为R1=12.5 mm,耦合珠的半径为≈5.18 mm,轻珠的半径为R3=6.25 mm。激励的无量纲初始速度仍为

从图7、8 和9 中仍可以看出,在激励施加点附近和初始高能阶段,存在P 波和S 波相互耦合的混合波传播模式[24],激励链和吸收链之间有强烈的能量交换现象。但是,如图7(a)和8(a)所示,不同于二维均匀颗粒链系统中脉冲传播方式,P 波从重珠传播到轻珠时,轻珠链中的波速提高,其波形和波幅都发生改变并继续向前传播,在轻-重珠交界位置发生反射。而在N2和N3区域,激励链和吸收链之间没有嵌入耦合珠,横向振动的幅值远小于存在耦合珠的N1和N4区域的颗粒振动,如图7(b)和8(b)所示。从能量角度分析,图9 所示的部分脉冲能量被局部约束在轻珠链中,并且延时传递到颗粒链末端,从而产生类似于一维非均匀颗粒容器的冲击缓释机理。

图8 激励链与吸收链中颗粒的速度时空图

图9 激励链与吸收链中颗粒的动能时空图

综上所述,除了激励链和吸收链之间发生能量交换之外,引起二维非均匀颗粒链系统中激励链末端的脉冲能量衰减的主要原因还包括:沿纵向传递的P 波在轻-重颗粒交界面发生反射,分解为一系列低幅值脉冲波,按照不同的波速、逐次释放到颗粒链末端;部分能量在通过轻珠链的过程,在x方向上发生多次反射,局部约束在轻珠链中,导致该部分能量的延时传递。因此,非均匀颗粒链系统能够进一步地分散脉冲能量和降低颗粒链的传递力。

由于非均匀颗粒结构导致的能量局部约束和延时传递,激励链末端的传递力的峰值减小,即传递力的波形由大幅值、小波长的脉冲波转变为一系列的小幅值、大波长的脉冲波。图10(a)、(b)分别给出了激励链中首、尾颗粒的能量分布和激励链末端的传递力,其能量衰减了约85%,激励链末端的传递力大约降低至80 N。与图11(a)、(b)中二维均匀颗粒链系统的结果相比,其能量大约减少了5%,传递力大约减小了20 N。因此,二维非均匀颗粒链系统能够将高幅值脉冲分散成一系列低幅值脉冲,从而降低激励链末端的传递力。

图10 二维非均匀颗粒链系统中激励链首尾颗粒的能量和末端传递力

图11 二维均匀颗粒链系统中激励链首尾颗粒的能量和末端传递力

3 结论

本文通过理论建模和数值仿真,研究了二维均匀颗粒链系统和非均匀颗粒链系统中的脉冲波传播特性、能量传递现象以及上述机理对减小颗粒链传递力的影响。结果表明,在冲击激励作用下,由于耦合珠的存在,二维均匀颗粒链系统中的脉冲波在激励施加点附近发生散射,激发出相互耦合的混合波,包括纵向传播的近零频率的孤立波(P 波)与横向高频剪切波(S 波)。在二维非均匀颗粒链系统中,这种强能量散射机制依旧存在,相比于前者,非均匀颗粒结构导致P 波从重珠传播到轻珠时,波形和波幅发生了改变,引起波速改变和时延的产生,在轻-重珠交界面位置发生反射。并且,非均匀颗粒结构能够将脉冲能量局部约束到轻珠链中,使得该部分能量延时传递到颗粒链末端,导致激励链末端传递力的峰值降低,力的波形由大幅值、小波长的脉冲波转变为一系列的小幅值、大波长的脉冲波。

与一维均匀颗粒链系统相比[6],二维均匀颗粒链系统能将激励链末端传递力的峰值降低到约100 N;而二维非均匀颗粒链系统对脉冲能量的削弱效果更好,激励链末端传递力的峰值进一步降低到约80 N。

本文通过引入轻珠链构成非均匀颗粒结构,将脉冲能量局部约束在轻珠链中,使得传递到颗粒链末端的脉冲能量逐级衰减,降低了末端的传递力的峰值。该研究结果对于二维颗粒结构在新型超结构和超材料的高效减振器应用方面具有重要的参考价值。

附录A 强非线性运动方程

在第1 节中已给出激励链中第i(2 ≤i≤N-1)颗颗粒的运动控制方程式(4)。在本节中将给出激励链首尾颗粒、吸收链和耦合链中颗粒的运动控制方程,列为式(A1)～(A6)。式(4)和(A1)～(A6)共同组成图1 所示的二维均匀颗粒链系统的控制方程。同第2 节中所述为颗粒质量为颗粒纵向位移分量为颗粒横向位移分量。其中,上标j=1、2、3 分别代表激励链、吸收链和耦合链,下标i代表颗粒的序号。当j=1、2 时,即对于激励链和吸收链;当j=3 时,即对于耦合链为颗粒的材料密度。

激励链中首颗粒的纵向位移和横向位移的运动方程,即i=1 时,

激励链中尾颗粒的纵向位移和横向位移的运动方程,即i=N时,有:

吸收链中首颗粒的纵向位移和横向位移的运动方程,即i=1 时,有:

吸收链中第i颗颗粒的纵向位移和横向位移的运动方程,当i=2、…、N-1 时,有:

吸收链中尾颗粒的纵向位移和横向位移的运动方程,即i=N时,有:

耦合链中第i颗颗粒的纵向位移和横向位移的运动方程,当i=1、…、N-1 时,有:

其中,α为重珠和重珠之间的接触系数,η为重珠和耦合珠之间的接触系数,二者的具体表达式可由式(2)获得。类似地,图6 所示的二维非均匀颗粒链系统中颗粒的运动控制方程也能根据赫兹接触定律和牛顿第二定律列出。