近场动力学系统中波传播特性的探究1)

朱竞高 任晓丹

(同济大学土木工程学院,上海 200092)

引言

在载荷的作用下,固体的变形并不是瞬间完成的[1],而是波传播、演化的结果.从这种意义上来说,波的传播对固体的行为有着极其重要的影响,是固体力学中的基础问题[2].特别是对于断裂问题,波的传递更是影响着裂纹萌生、成核、分岔乃至破碎的全过程[3-5],是正确捕捉裂纹扩展的关键.

随着计算机技术的发展和固体力学研究的深入,数值计算方法纷纷涌现[6-8],为研究固体力学中的科学问题提供了更为便捷的途径.其中,近场动力学作为一类基于非局部思想的新固体力学方法[9],从邻域内粒子间的相互作用出发,考虑了微纳米尺度下的长程力,在模拟裂纹方面更具合理性.该方法采用积分形式的控制方程,利用微观键的失效来描述断裂,能够从根本上避免求解不连续问题的奇异性和复杂性,无需特殊处理即可实现对裂纹发展的捕捉[10-12],被广泛用于断裂问题的求解.经过发展,目前形成了键基[13-14]和态基[15-17]近场动力学两大分支.Zhou 等[18]和Zhu 等[19]还引入了剪切效应,以更准确地描述材料的力学性能.随着对材料非线性研究的深入,Gu 等[20]和Liu 等[21]建立了弹塑性模型,王涵等则发展了热黏塑性模型[22-23],在预测各类材料的断裂失效问题中取得重要的进展[24-27].

然而,近场动力学作为一种非局部系统,不可避免地会引入波色散问题.波色散会改变波的传播特性,不利于对固体行为的正确捕捉,进而影响对裂纹扩展的准确预测.Silling[13]在提出键基模型时,首先发现了近场动力学中的非线性色散行为.2016 年,随着Bazănt 等[28]从非局部性的角度对色散关系开展了进一步研究后,近场动力学中的色散问题逐渐引起了学者们的关注.Butt 等[29]对态基近场动力学中的色散关系进行了研究,提出调整邻域半径尺寸和影响函数来减弱波色散的方法.Wildman[30]分析了键基近场动力学中的色散关系,并实现了通过改变微弹模参数来匹配所需的色散方程.Gu 等[31]和Li 等[32]也开展了相关研究,探讨不同衰减函数对近场动力学中色散关系的影响.除了通过调整参数来减小波色散外,还有学者[33-34]从借鉴有限差分法中色散研究的角度出发,提出有限差分增强的近场动力学方法.但是,这些研究主要局限于低频成分,忽视了高频成分的色散行为,同时也缺乏对波间断性的考虑.

近年来,恐怖袭击和局部战争频发,极端载荷作用下工程结构的安全性问题引起广泛关注,成为国防安全领域的研究重点.极端载荷,如爆炸、冲击作用,往往伴随着激波的产生[35-37],激波的传播又影响着材料破碎和裂纹发展,是进行工程结构安全设计的关键.而激波因其高频、强间断的复杂性,给相关分析理论和计算方法带来了挑战.近场动学方法采用积分形式的控制方程,自然地适用于激波引起的材料破碎和裂纹发展模拟,但该方法对激波的准确捕捉却不可避免地会受到高频、强间断成分色散特性的制约.

因此,对近场动力学的波色散特性展开全面系统的研究,纳入对高频成分和波间断性的考虑,构成了文章的出发点.理论研究表明,相比于低频成分,高频成分还表现出振荡关系和零能模式,色散问题更为突出.高频域的色散关系会随波的传播方向发生变化,引起空间传播的方向性.数值计算则发现强间断性波的色散现象更为显著,呈现出Gibbs 不稳定性.说明色散不仅会对波的传播特性产生不利影响,还会引起数值不稳定性.为了解决上述问题,通过黏性引入耗散来抑制色散行为,并考虑对高频成分的选择性抑制,构造了相应的黏性力态.本文研究工作将为在近场动力学的框架下准确模拟波的传播过程,正确捕捉固体行为提供重要参考.

1 近场动力学方法

1.1 基本思想

近场动力学通过引入特征长度 δ,定义了邻域范围

如图1 所示,在邻域范围内,粒子间由键相连,键随着粒子运动而发生变形.对某一空间域 Ω,假设在某一时刻t,域内任一物质点x只与其邻域范围内的其他物质点x′之间存在相互作用力,则根据牛顿第二定律可得近场动力学的控制方程为[38]

图1 近场动力学模型示意图Fig.1 The schematic of peridynamics

其中,ρ为物质密度;u为物质点的位移;b为体力密度;f为力密度函数,包含了材料的本构关系.

根据本构关系不同,可将近场动力学分为键基近场动力学和态基近场动力学.键基近场动力学由于模型简单、计算方便,在断裂问题中的应用更为广泛.因此,本文的研究将基于键基模型进行展开.这里,对键基近场动力学作简要介绍.

1.2 键基近场动力学

相比于态基近场动力学,键基近场动力学的模型较为简化.该模型认为键与键之间相互独立,力密度函数只与单一键有关[13].若定义粒子间的相对位置和相对位移为

则未变形的键和变形后的键可以分别用 ξ和 η +ξ表示.由于采用了键间独立性的假定,粒子间的相互作用可视为中心弹簧,将微观势函数表示为

其中,s为键的相对伸长量,C (∥ξ∥)为微弹性模量.对势函数求导可得力密度函数的表达式

其中n为变形后键的方向

微弹性模量有多种函数形式,当取为常数c时,键基近场动力学就变成了微观弹脆性模型,将被用于后面的理论分析和数值计算.该模型根据应变能等效,可得不同维度下的微弹性模量

其中,E为弹性模量,A为一维杆的横截面积,B为二维板的厚度,二维情况只选取了平面应力问题作为代表.此时,力密度函数可写为

2 波的传播特性分析

这里,将从微观弹脆性模型入手,兼顾低频和高频成分,对近场动力学系统中波的传播特性展开全面的分析.采用小变形假设,微观弹脆性模型的控制方程可以写为[28]

2.1 一维谱分析

对于一维问题,式(9)可以写为

其中,A为一维模型的横截面积.为了研究波的传播特性,有必要在频域内展开分析.这里关注线弹性问题,因此,可对位移项进行傅里叶展开

其中,i 为虚数单位,满足 i2=−1 ;a为幅值;v为相速度;κ为 2 π长度内的波数.将式(11)代入式(10),可约去时间项得

进一步,考虑邻域的对称性,可约去奇数项得

若定义与波数有关的项为

则式(13)可视为相速度v的方程

求解可得相速度为

显然,式(16)中的虚数项为0,故近场动力学系统中没有数值耗散.当波色散产生不利影响时,如数值振荡、非物理变形等,系统本身无法对其进行抑制.随着不利因素的逐渐累积,数值计算的精度将会受到严重的影响,无法保证对固体行为的正确捕捉.这一方面说明了进行近场动力学系统中波色散研究的必要性,另一方面为抑制色散现象提供了方向.

进一步,考虑数值计算中的空间离散,对近场动力学系统的色散特性进行分析.根据邻域的对称性,经过数值离散,式(14)变为

其中,m为一维邻域内的粒子数,h为粒子间距,c为微弹性模量.采用式(7)中的微弹性模量,式(17)可写为

将式(18)代入式(16),可得一维近场动力学系统的色散关系

其中,v0为弹性理论波速

其中,E为弹性模量,ρ为密度.由于 κ δ和 κ ξi均与 κh有关,式(19)可视为Nh=κh/(2π)的函数,表示粒子间距内的波数.若取邻域半径为 δ=3.015h,式(19)中的色散关系将如图2 所示.

由图2 可知,在一维近场动力学系统中,相速度与理论波速的比值偏离局部理论的结果,表现出色散行为.局部理论各个频率成分的波的传播速度相同,在传播时会同时达到.而近场动力学中各个频率成分的波的传播速度不同,在传播时会分先后达到.对于传播速度小于理论波速的频率成分,在传播中会落后于其他波,表现为拖尾波.因此,定义临界点为与局部理论色散关系的交点,即的点.可以看出,大部分频率成分位于临界点以下,在波的传播中表现为拖尾波.

图2 一维近场动力学系统的色散关系Fig.2 The dispersion relation of 1D peridynamics

为区分色散行为和研究方便,将Nh≥1定义为高频域,0 ≤Nh<1定义为低频域.即高频域与低频域的划分与波长有关,高频域的波长范围为 λ ≤h,低频域的波长范围为 λ >h,其中h为粒子间距.显然,低频成分(虚线段)和高频成分(实线段)分别呈现出不同的色散特性.在高频域,波速不仅存在振荡的趋势,还在某些频率处等于0.根据文献[28-29],这些点的相速度为0,对应于零能模式点.另有学者将零能模式定义为群速度为0 的频散带,在这种定义下,首先有必要给出一维模型的频率表达式.基于式(19),有

显然,由于余弦函数具有周期性,频率将呈现在Nh取整区间内的周期性变化.为了更加直观地说明此周期性,将频率与Nh的关系绘制成图3.

图3 一维近场动力学系统的频率关系Fig.3 The frequency relation of 1D peridynamics

如图3 中蓝色线所示,频率随波数呈现周期性变化,并在波数取整处为0.与图2 不同的是,这种周期性在低频域(虚线)和高频域(实线)呈现相同的变化规律.由于频率呈现从0 到0 的周期性变化,对应的斜率应当呈现从正到负或从负到正的周期性变化.而频率的斜率正比于群速度,说明存在群速度为0 的点.为了进一步验证上述结论,基于频率的表达式(21),求导得到群速度与波数的关系,绘制图像进行表示.如图3 中红色线所示,群速度会在整波数区间内从正变为负,并在从正变为负的过程中出现零点.在这种情况下,零能模式点对应于图3 中红色线与y=0 水平线的交点.

综上所述,无论是采用何种形式的零能模式定义,一维近场动力学模型中均存在零能模式点.振荡趋势和零能模式均是非物理响应,会加重色散行为带来的不利影响.因此,相比于低频成分,高频成分的色散问题更为严重,应当给予关注.

2.2 二维谱分析

下面,考虑波的空间传播,对二维系统展开分析.此时,除方向x1外,还需引入方向x2,控制方程(9)变为

其中,B为二维模型的厚度,ξ为相对位置

同时,位移项的傅里叶展开变为

其中,a为与位移u共线的幅度向量;b为波传播方向的单位矢量;其他变量的定义同式(11).将式(24)代入式(22),可同时约去矢量项和时间项,有

同样地,根据邻域的对称性,可约去奇数项得

仍可将与波数有关的项用一个变量进行表示

此时,关于波速的方程可写为

解得相速度为

与一维的谱分析结果相同,相速度的虚数项为0,没有数值耗散,无法抑制色散带来的不利影响;实数项与波数有关,表现出色散行为.

进一步,考虑数值计算,对Nκ进行空间离散,有

其中,n为二维邻域内的粒子数,h为粒子间距,ξj·b为相对位置在波传播方向上的投影.仍采用式(7)中的微弹性模量,式(30)可写为

将式(31)代入式(29),可以得到二维近场动力学系统的色散关系

其中,v0为式(20)定义的弹性理论波速.显然,在二维问题中,波的色散行为不仅与频率有关,还与波的传播方向有关.

如图4 所示,对于对称邻域,对关于x2轴对称的传播方向b′,一定可以同样对称地找到关于x2轴对称的相对位置向量,使式(30)中的投影相等

图4 二维键基模型色散关系的对称性Fig.4 The symmetry of dispersion relation for 2D peridynamics

由于邻域内每个相对位置向量 ξj,都存在关于x2轴对称的相对位置向量满足投影相等关系,故最终求和得到的式(30)也相等.同样地,对于x1轴也满足此对称性

因此,可将波的传播方向缩小至 θ=0 ∼90◦范围内进行研究.显然,在 0◦∼90◦范围内,关于45◦方向存在对称性

因此,二维键基模型的色散行为将表现出沿45◦的对称性.

基于此对称性,对波传播方向与x1方向的夹角在 θ=0◦∼45◦之间的色散行为进行研究.这里,仍取Nh=κh/(2π)为横坐标,并取波的传播方向分别为θ=0◦,30◦,45◦,式(32)中的色散关系将如图5 所示.由图5 可知,当波沿x1方向传播时,二维近场动力学系统中色散关系的变化趋势与一维系统较为一致,即在低频成分和高频成分表现出不同的色散行为,并在高频域呈现出振荡趋势和零能模式.但随着传播方向改变,高频成分色散关系的振荡程度发生变化,零能模式点的位置也在移动.

图5 二维近场动力学系统的色散关系Fig.5 The dispersion relation of 2D peridynamics

为了更加直观地表征不同方向的频散行为,提取高频成分色散关系变化率的最大值作为振荡程度的衡量指标,将 0◦∼360◦范围内高频成分色散关系的振荡程度绘制雷达图进行对比分析.如图6 所示,对于对称邻域,高频成分色散关系的变化趋势呈现出关于 4 5◦方向的对称性,验证了前述分析的正确性.在 0◦∼45◦方向范围内,随着波的传播方向改变,振荡指标的大小发生变化,与图5 的分析结果相一致.

图6 不同方向频散行为的雷达图Fig.6 Radar chart of dispersion behavior in different directions

综合上述分析,相比于低频成分,高频成分的色散行为更为严重,表现出振荡趋势和零能模式.高频域的色散关系会随波的传播方向变化,呈现出沿45°的对称性.而近场动力学系统本身没有数值耗散,无法抑制色散带来的不利影响,会严重影响对固体行为的正确捕捉.

3 近场动力学方法的改进

事实上,作为一种非局部理论,近场动力学在高频段呈现频散行为是非局部效应的合理体现.对于非均质材料,频散行为是真实存在的,可以通过近场动力学来模拟非均质材料中的非线性色散关系,且近场作用半径与非均质材料的特征尺寸具有强相关性[29].但是,对于变异性不大的材料,通常可以认为是均质材料,近场动力学的频散行为就会带来不利影响.根据第2 章的研究,近场动力学的频散行为在高频域表现出振荡趋势.也就是说,理论上应该同时达到的不同频率成分的波,在近场动力学的计算中传播速度会出现差异,这种差异会在数值模拟中引起非物理变形和损伤.

对于极端载荷作用下的问题,一方面这种不利影响会加重,严重影响数值模拟的精度.另一方面,极端载荷作用往往伴随着裂纹的产生,而裂纹的扩展与波的传播紧密相关,异常的波传播过程必定会对裂纹扩展的预测产生不利影响.而极端载荷作用下工程结构的安全性问题是国防安全领域的研究重点,近场动力学方法的优势正是对于裂纹扩展捕捉的能力.因此,有必要抑制近场动力学中的频散行为,提高其对波传播的模拟精度,进而保证其对极端载荷作用下计算的可靠性,为国防安全邻域研究提供理论指导和技术支撑.本章将从引入数值耗散的角度出发,对近场动力学体系的色散行为进行抑制.并考虑高频域色散更为严重的特点,对高频成分进行选择性抑制.

3.1 黏性效应的研究

在固体力学中,通常认为黏性与数值耗散有关,为近场动力学中耗散机制的引入提供了方向.首先,有必要对黏性效应进行研究,证明黏性引入的有效性.这里,仍从微观弹脆性模型入手.考虑到黏性与应变率的相关性,可将黏性项取为力密度函数的率形式

其中,γ为无量纲系数.在小变形下,有

以一维问题为例,式(38)可化简为

其中,A为一维模型的横截面积.沿用第2 章的谱分析方法,采用式(11)中的傅里叶展开,式(39)可写为

其中,v′为黏性作用下的相速度.仍采用式(14)中对波数相关项的定义,式(40)可表示为相速度的方程

求解可得相速度

此时,虚数项不再为0

说明黏性项能够引入数值耗散,这与文献[40-43]的研究结论相一致,证明了黏性引入的有效性.

进一步,考虑高频域色散更为严重的特点,对黏性项在不同频率成分的表现进行研究.由于数值耗散的引入是为了抑制色散行为带来的不利影响,这里关注虚数项(43)与原始相速度(16)的比值.借助式(18),比值可以表示为

其中,v0为式(20)定义的弹性理论波速,γ为无量纲系数,δ为邻域半径,h为粒子间距.仍取Nh=κh/(2π)为横坐标,归一化后式(44)中的比值关系将如图7所示.由图7 可知,数值耗散对色散行为的抑制作用在低频域(虚线段)和高频域(实线段)呈现相同的周期变化,即黏性效应对低频和高频成分的表现相同.实际上,式(37)本质上是一种线性黏性,无法实现对不同频率成分的选择性抑制.

图7 数值耗散对不同频率成分的抑制作用Fig.7 The dissipation effect on different frequencies

3.2 黏性效应的引入

因此,可将黏性效应引入近场动力学,抑制波色散带来的不利影响.并有必要在黏性项的构造中,考虑对高频成分的选择性抑制.基于近场动力学的非局部性和键基模型的键间独立性假设,将在键层次以非局部作用的方式引入黏性效应.相应地,将控制方程定义为

其中,f为力密度函数,是传统力态;fvs为黏性力态.同时,考虑到在键基模型中,力密度函数的方向为变形后键的方向,故满足线动量守恒

和角动量守恒

为了保持上述守恒关系,取黏性力态的方向与力密度函数的方向相同.此时,黏性力态可以分解为标量项和方向项

其中,fvs为黏性力态的标量形式,n为变形后键的方向.显然,这种引入方式不仅保留了传统近场动力学模型的理论框架,而且只需要给出黏性力态标量形式的定义.

下面,将针对体积变形,这一常见的固体变形模式,构造黏性力态的标量表达式.参照连续介质力学中体积黏性的定义[44-45]

其中,β为无量纲体积黏性系数;cd为膨胀波速;Le为单元的特征尺寸;为体积应变率

其中,u为位移场.定义体积黏性力态的标量形式为

其中,b1为无量纲体积黏性系数;v0为式(20)定义的弹性理论波速;VH为邻域的体积;ζ为应变率.在键基模型中,键的相对伸长量s为一种应变测量,可用于应变率的计算

如图8 所示,∥ ξ∥ 为键的初始长度(蓝线),η ˙·n为键的速度在变形后键的方向上的投影(红线).

图8 应变率的运动学表示Fig.8 The kinematic expression of strain rate

基于上述定义,可得体积黏性力态的标量表达式

当b1αp=γc时,式(54)与式(39)相同,证明了3.1 中黏性效应分析的正确性.

由3.1 中黏性效应分析可知,式(53)本质上是一种线性黏性,无法实现对不同频率成分的选择性抑制.考虑到高频成分波色散的严重性,有必要对其进行非线性改进.这里,将引入二次项,以放大对高频成分的抑制作用

其中,b2为无量纲二次黏性系数,其余变量的定义与式(51)一致.为充分发挥一次项和二次项的黏性效应,将式(53)和式(55)进行多项式组合

其中,b1和b2分别为一次项和二次项的黏性系数,得到的黏性力态称为多项式黏性.

式(48)给出了引入黏性的理论框架,式(56)则从体积变形出发,考虑了对高频成分的选择性抑制,构造了多项式黏性力态,完成了对近场动力学方法的改进.

4 数值算例

本章基于波传播的数值算例,开展数值模拟研究.旨在对理论分析进行验证的同时,考察波间断性的影响.在近场动力学中,通常采用EMU 方法[46]进行数值计算,该方法对控制方程进行空间离散

其中,N为邻域内的粒子数目,∆V为离散粒子的体积.由于波的传播为动力问题,采用显式中心差分方法进行求解

其中,u,和分别为加速度,速度和位移;∆t为满足稳定条件的时间步长.

4.1 一维波传播问题

对一维算例,首先考虑带理论解的波传播问题,进行对比研究,以验证理论分析的正确性.再在此基础上,考虑极端载荷下的激波传播问题,探究波间断性的影响.

4.1.1 带理论解的波传播问题

为了验证理论分析的正确性,并说明对不同间断性波的适用性,这里取具有解析表达式的光滑正弦波和间断方波为研究对象,对一维波的传播问题进行了研究.将1 m 长的弹性杆一端固定,另一端输入目标波,提取t=77.3 μs 时沿杆长的位移分布,将键基模型、线性黏性和多项式黏性的计算结果与理论解进行对比分析.图9 和图10 分别给出的是正弦波和方波的计算结果.显然,两种波的键基模型计算结果均在波后出现了数值振荡.结合理论分析结果,一维近场动力学系统中存在波色散,且大部分频率成分的波速小于理论波速,表现为拖尾波的形式,第二章中谱分析的正确性得到了验证.

图9 带理论解的一维正弦波传播结果对比Fig.9 Comparative study of 1D sine wave propagation problem with theoretical result

图10 带理论解的一维方波传播结果对比Fig.10 Comparative study of 1D square wave propagation problem with theoretical result

提取波后振荡进行进一步研究,发现在两种波形的计算结果中,黏性的引入均可以有效地抑制数值振荡.且多项式黏性能够在抑制波后振荡的同时,更好地保持原始波峰和波形,得到了更为合理的波传播结果.这验证了第3 章所提方法的正确性.需要注意的是,相比于正弦波,方波的波后振荡更为显著,说明近场动力学中的色散行为与波的间断性有关,这将在后续算例中进行研究和讨论.

4.1.2 激波传播问题

考虑1 m 长弹性杆,在两端受爆炸载荷作用的波传播问题.如图11 所示,两端炸药距杆端均为0.5 m,且同时发生爆炸.为探究波间断性的影响,设置小当量炸药(W=1 kg)和大当量炸药(W=80 kg)两类工况进行计算.杆的材料参数为 ρ=E=193 GPa.将一维杆离散成1000 个粒子进行计算,邻域半径取 δ=3.015h,采用压力时程模拟爆炸载荷,并忽略负压区的影响.

图11 一维波传播算例示意图 (单位: mm)Fig.11 The schematic of numerical example for 1D wave propagation (unit: mm)

(1)小当量工况(W=1 kg)

图12 和图13 分别为5 个典型时刻(t=50,100,150,200,250 μs)沿杆长方向的位移和应力分布.根据杆长和材料特性,波的传播周期为

计算结果表明,爆炸载荷产生压力波,压力波由杆的两端向中间传播.t=100 μs 时在中点(A点)相遇后,继续向前进行传播.t=200 μs 时达到杆端后,反射产生拉伸波,向相反方向进行传播.模拟得到的波的传播周期与式(59)相符,验证了数值计算的正确性.相比于图12,图13 中的应力波并不光滑,在波后出现了数值振荡.结合第2 章的谱分析结果,一维近场动力学系统中存在波色散,且大部分频率成分的波速小于理论波速,表现为拖尾波的形式,谱分析的正确性得到了验证.

图12 沿一维杆长方向的位移分布Fig.12 Displacement distribution along the 1D bar

图13 沿一维杆长方向的应力分布Fig.13 Stress distribution along the 1D bar

进一步,提取A点的应力波时程,对黏性效应进行研究.如图14 所示,两端的压力波在t=100 μs 时相遇叠加,反射产生的拉伸波在t=300 μs 时再次相遇叠加,波的传播周期与理论解相符.对波峰处放大进行观察,键基系统的计算结果在波后出现了数值振荡,这与图13 的计算结果相一致.施加黏性后,数值振荡得到了有效地抑制,证明了黏性引入的有效性.

对比两种黏性力态,多项式黏性能够在抑制波后振荡的同时,更好地保持原始波峰和波形,得到更为合理的波传播结果.这是因为数值耗散的引入会消耗波的能量,从而影响波峰和波形.而多项式黏性中含有二次项,可以放大对高频成分的抑制作用,在同等的抑制效果下减小了能量耗散,降低了对原始波峰和波形的影响.为了更直观地从能量角度进行说明,对图14 进行频域分析,得到应力波的能量谱密度曲线.如图15 所示,相比于体积黏性,多项式黏性能够降低能量耗散,且二者的差异主要集中在105～106的频率区间.

图14 A 点的应力波时程Fig.14 Time history of stress wave at point A

图15 A 点应力波的能量谱密度曲线Fig.15 Energy spectrum density of stress wave at point A

(2)大当量工况(W=80 kg)

下面,纳入对波间断性的考虑,对大当量工况进行计算.此时,t=50,100,150,200,250 μs 时刻沿杆长方向的位移和应力分布分别如图16 和图17 所示.显然,相比于小当量工况,此时的位移和应力的波形变得更为尖锐,即波的间断性变强,对应于激波的产生.与小当量的计算结果相比,位移也出现了波后振荡,应力的数值振荡也更加明显.说明强间断性波的色散行为更为显著,呈现出Gibbs 不稳定性.因此,色散在对波传播特性产生不利影响的同时,还将引起数值不稳定性.

图16 沿一维杆长方向的位移分布Fig.16 Displacement distribution along the 1D bar

图17 沿一维杆长方向的应力分布Fig.17 Stress distribution along the 1D bar

仍提取点A的应力波时程进行进一步分析,应力波时程和能量谱密度曲线分别如图18 和图19 所示.与图17 的计算结果相一致,图18 中波的间断性变强,波后的振荡加剧.与小当量工况相同,黏性能够有效地抑制色散造成的数值振荡,且多项式黏性能够更好地保持原始波峰和波形,说明所提方法对激波同样适用.图19 则从能量角度说明了多项式黏性的优势,即减小了能量耗散,从而降低了对波峰和波形的影响.

图18 A 点的应力波时程Fig.18 Time history of stress wave at point A

图19 A 点应力波的能量谱密度曲线Fig.19 Energy spectrum density of stress wave at point A

综合一维计算结果,近场动力学系统中的波色散在一维波传播问题中表现为波后振荡,强间断性波的色散行为更为显著,呈现出Gibbs 不稳定性.黏性能够有效地抑制波后振荡,且多项式黏性能够更好地保持原始波峰和波形.

4.2 二维波传播问题

对二维算例,考虑内径0.15 m,外径0.3 m 的圆环在爆炸载荷作用下的波传播问题.如图20 所示,炸药在圆环的中心位置处爆炸,距圆环内壁0.15 m.圆环的材料参数为,E=30 GPa.这里,取炸药质量W=70 kg,对大当量爆炸下的激波传播进行模拟.在数值计算中,粒子间距取h=10 mm,邻域半径取 δ=3.015h,爆炸载荷仍采用压力时程进行模拟,并忽略负压区的影响.

图20 二维波传播算例示意图 (单位: mm)Fig.20 The schematic of numerical example for 2D wave propagation (unit: mm)

中心爆炸下圆环中波的传播本质上是球形腔加压问题,存在理论解[47].这里,关注压力云图,并采用对比研究的方式,以考察色散对波的空间传播行为的影响.由图21 (a1)～ 图21(c1) 可知,爆炸理论上会产生压力波,压力波沿周向均匀分布,并逐渐沿径向向外传播.而在键基系统的计算结果中,如图21 (a2)～图21(c2)所示,压力波的分布并不均匀,呈现出沿45°的对称性.结合第2 章的谱分析结果,二维系统中高频成分振荡的色散关系具有沿45°方向的对称性,数值计算结果验证了谱分析的正确性.

图21 二维圆环的压力云图Fig.21 The pressure cloud map of 2D ring

此外,相比于一维问题,二维情况下的拖尾波现象更为显著,拖尾波的峰值(3 GPa)甚至超过了波前峰值(1.5 GPa).拖尾波的存在不仅会在波后引起非物理变形,还会消耗波前峰值,不利于对波传播过程的正确捕捉.因此,在二维情况下,波色散会通过影响波的空间分布和传播过程,对波的空间传播行为产生不利影响.

进一步,观察多项式黏性的计算结果,以探究黏性在二维问题中的表现.如图21 (a3)～ 图21(c3)所示,此时的拖尾波得到了有效的抑制,降低了对波前峰值的消耗,波的传播过程得到了较好地捕捉.同时,压力波的分布更为均匀,与理论解更为一致.因此,在黏性的作用下,近场动力学系统中波的空间传播行为更为合理,证明了黏性引入的有效性.此时,压力波的分布宽度相比理论解有所增加,这应当与特征长度和数值耗散有关.近场动力学引入了特征长度,即邻域半径 δ,在特征长度范围内考虑粒子间的相互作用.特征长度的存在会在一定程度上扩大影响范围,使波的分布宽度变大.另一方面,由一维数值模拟结果可知,数值耗散会消耗波的能量,在一定程度上增大波宽.但总体来说,黏性的引入,使得波前峰值和波的空间分布均得到了较好地捕捉,证明了所提方法同样适用于多维问题.

5 结论

本文采用谱分析方法,对近场动力学系统中不同维度、不同频率成分的色散特性进行了研究,并提出了相应的改进方法.同时进行数值计算,验证理论分析结果,并探究波间断性的影响,主要结论如下.

(1) 相比于低频成分,高频成分的色散关系还呈现出振荡趋势和零能模式,色散问题更为突出.

(2) 高频成分的色散行为会随波的传播方向发生变化,呈现出沿45°的对称性.

(3) 色散在影响波传播特性的同时,还会引起Gibbs不稳定性,表现为强间断性波的色散现象更为显著.

(4) 黏性能够引入数值耗散,抑制色散带来的不利影响.二次项的引入可以放大对高频成分的抑制作用,在同等的抑制效果下,减小对原始波峰和波形的影响.

(5) 色散在一维波传播中引起波后振荡,在二维问题中将引起更为严重的拖尾波和波空间分布的不均匀性.

本文研究可为减小近场动力学系统中波色散的不利影响,获得正确的波传播过程和固体行为提供重要参考.同时,尚存在一些值得进一步研究的问题,如在冲击问题中的应用、对零能模式的抑制等.事实上,若想纳入针对零能模式的抑制作用,需要先识别出近场动力学中的零能变形模式,再针对其进行抑制.零能模式是一种不引起外力的虚假变形模式,可以从零能模式的特点出发,对近场动力学中的零能变形模式进行研究,确定其变形特征.再在本文提出的数值黏性的基础上,根据研究出的近场动力学中零能模式的变形特征,对变形相关项进行修改,以实现对零能模式的抑制效果.