STUDY OF WAVE DISPERSION AND PROPAGATION IN PERIDYNAMICS
-
摘要: 近场动力学是一类基于非局部思想的新固体力学方法, 其采用积分形式的控制方程, 自然地适用于极端载荷下材料破碎和裂纹发展的模拟, 被广泛用于国防安全等领域的研究. 但是, 非局部性会引入色散效应, 对波的传播产生不利影响, 制约其对断裂等固体行为的捕捉能力. 为此, 采用谱分析方法, 对近场动力学系统的色散行为进行了全面的研究. 发现相比于低频成分, 高频成分的色散关系呈现出振荡趋势和零能模式, 色散问题更为严重. 高频域的色散行为还随波的传播方向发生改变, 呈现出沿45°的对称性. 而近场动力学系统本身缺乏数值耗散, 无法抑制色散问题带来的不利影响. 因此, 从引入数值耗散的角度出发, 在合理保留传统近场动力学理论框架的基础上, 建立了黏性引入的控制方程. 并考虑固体中常见的体积变形和对高频成分的选择性抑制, 构造了相应的黏性力态. 最后, 在数值研究中模拟了极端载荷下激波的产生, 以探究波的间断性对色散行为的影响. 发现间断性强的波表现出更为显著的色散行为, 呈现出Gibbs不稳定性. 这些均能有效地被黏性力态所抑制, 验证了所提方法的正确性. 这为在近场动力学系统中实现对波传播过程的正确捕捉, 获得正确的固体行为提供了重要参考, 从而为国防安全领域研究提供了技术支撑和借鉴.Abstract: Periydnamics (PD) is a new nonlocal method reformulated from solid mechanics. It adopts the integral form of governing equation and is naturally suitable to model fragments and cracks under extreme events, thus widely applied in the field of national defense security. However, the nonlocality in PD introduces the dispersion effect and imposes adverse effect on wave propagation, which will greatly restrict its capability in capturing solid behaviors, especially the fractures. For this purpose, we employ the spectral analysis method to study the dispersion behavior of PD system comprehensively. It is found that compared to the low frequencies, the dispersion relation of high frequency components shows an oscillation trend and zero-energy modes, leading to more serious dispersion problems. The dispersion behavior of high frequencies changes with the wave propagation direction and shows 45° symmetry in the spatial wave propagation. As the PD system itself is non-dissipative, the adverse effect of the dispersion problem can not be suppressed. As a result, the simulation accuracy may be greatly influenced. To introduce the numerical dissipation for dispersion effect suppression, the governing equation of viscosity introduction is proposed as a minimum variation of conventional PD. Both the typical deformation in solids and the selective suppression on high frequencies are considered then the corresponding viscous force state is constructed. Finally, a numerical study is conducted to model the shock waves under extreme events and investigate the influence of wave discontinuity. It is indicated that the wave discontinuity aggravates the dispersion problem and shows Gibbs instability in the wave propagation. These can be effectively suppressed by the viscous force state, which verifies the proposed method. This provides an important reference to reproduce the correct wave propagation process and obtain the reasonable solid behavior in PD, thus helps to support and guide the research of national defense security field.
-
Keywords:
- peridynamics /
- wave dispersion /
- high frequency /
- wave discontinuity /
- viscosity
-
引 言
在载荷的作用下, 固体的变形并不是瞬间完成的[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 基本思想
近场动力学通过引入特征长度
$ \delta $ , 定义了邻域范围$$ {{{H}}_x} = {{{H}}_x}\left( {{\boldsymbol{x}},\delta } \right): = \left\{ {{\boldsymbol{x'}} \in R:\left\| {{\boldsymbol{x'}} - {\boldsymbol{x}}} \right\| \leqslant \delta } \right\} $$ (1) 如图1所示, 在邻域范围内, 粒子间由键相连, 键随着粒子运动而发生变形. 对某一空间域
$\varOmega$ , 假设在某一时刻$ t $ , 域内任一物质点${{\boldsymbol{x}}}$ 只与其邻域范围内的其他物质点${{\boldsymbol{x}}'}$ 之间存在相互作用力, 则根据牛顿第二定律可得近场动力学的控制方程为[38]$$ \rho \left( {\boldsymbol{x}} \right){\boldsymbol{\ddot u}}\left( {{\boldsymbol{x}},t} \right) = \int_{{{{H}}_x}} {{\boldsymbol{f}}\left( {{\boldsymbol{x}},{\boldsymbol{x'}}} \right){\text{d}}V'} + {\boldsymbol{b}}\left( {{\boldsymbol{x}},t} \right) $$ (2) 其中,
$ \rho $ 为物质密度;${{\boldsymbol{u}}}$ 为物质点的位移;${{\boldsymbol{b}}}$ 为体力密度;${{\boldsymbol{f}}}$ 为力密度函数, 包含了材料的本构关系.根据本构关系不同, 可将近场动力学分为键基近场动力学和态基近场动力学. 键基近场动力学由于模型简单、计算方便, 在断裂问题中的应用更为广泛. 因此, 本文的研究将基于键基模型进行展开. 这里, 对键基近场动力学作简要介绍.
1.2 键基近场动力学
相比于态基近场动力学, 键基近场动力学的模型较为简化. 该模型认为键与键之间相互独立, 力密度函数只与单一键有关[13]. 若定义粒子间的相对位置和相对位移为
$$ \left. \begin{gathered} {{\boldsymbol{\xi}} }{\text{ = }}{{\boldsymbol{x}}'} - {{\boldsymbol{x}}} \\ {{\boldsymbol{\eta}} } = {{\boldsymbol{u}}}\left( {{{\boldsymbol{x}}'}} \right) - {{\boldsymbol{u}}}\left( {{\boldsymbol{x}}} \right) \\ \end{gathered} \right\} $$ (3) 则未变形的键和变形后的键可以分别用
${{\boldsymbol{\xi}} }$ 和${{\boldsymbol{\eta}} } + {{\boldsymbol{\xi}} }$ 表示. 由于采用了键间独立性的假定, 粒子间的相互作用可视为中心弹簧, 将微观势函数表示为$$ w\left( {{{\boldsymbol{\eta}} },{{\boldsymbol{\xi}} }} \right) = \frac{{\mathfrak{C}\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right){s^2}\left\| {{\boldsymbol{\xi}} } \right\|}}{2},{\text{ }}s = \frac{{\left\| {{{\boldsymbol{\eta}} } + {{\boldsymbol{\xi}} }} \right\| - \left\| {{\boldsymbol{\xi}} } \right\|}}{{\left\| {{\boldsymbol{\xi}} } \right\|}} $$ (4) 其中,
$ s $ 为键的相对伸长量,$\mathfrak{C}\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)$ 为微弹性模量. 对势函数求导可得力密度函数的表达式$$ {{\boldsymbol{f}}}\left( {{{\boldsymbol{\eta}} },{{\boldsymbol{\xi}} }} \right) = \frac{{\partial w\left( {{{\boldsymbol{\eta}} },{{\boldsymbol{\xi}} }} \right)}}{{\partial {{\boldsymbol{\eta}} }}} = \mathfrak{C}\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)s\left\| {{\boldsymbol{\xi}} } \right\|\frac{{\partial s}}{{\partial {{\boldsymbol{\eta}} }}} = \mathfrak{C}\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)s{{\boldsymbol{n}}} $$ (5) 其中
${{\boldsymbol{n}}}$ 为变形后键的方向$$ {{\boldsymbol{n}}} = \frac{{{{\boldsymbol{\eta}} } + {{\boldsymbol{\xi}} }}}{{\left\| {{{\boldsymbol{\eta}} } + {{\boldsymbol{\xi}} }} \right\|}} $$ (6) 微弹性模量有多种函数形式, 当取为常数
$ c $ 时, 键基近场动力学就变成了微观弹脆性模型, 将被用于后面的理论分析和数值计算. 该模型根据应变能等效, 可得不同维度下的微弹性模量$$ c = \left\{\begin{array}{l}2E/(A{\delta }^{2}),\quad\;\,\,\, 1{\rm{D}}\\ 9E/(B\text{π} {\delta }^{3}),\quad\, 2{\rm{D}}\\ 12E/(\text{π} {\delta }^{4}),\quad\,\, 3{\rm{D}}\end{array}\right. $$ (7) 其中,
$ E $ 为弹性模量,$ A $ 为一维杆的横截面积,$ B $ 为二维板的厚度, 二维情况只选取了平面应力问题作为代表. 此时, 力密度函数可写为$$ {{\boldsymbol{f}}}\left( {{{\boldsymbol{\eta}} },{{\boldsymbol{\xi}} }} \right) = cs{{\boldsymbol{n}}} $$ (8) 2. 波的传播特性分析
这里, 将从微观弹脆性模型入手, 兼顾低频和高频成分, 对近场动力学系统中波的传播特性展开全面的分析. 采用小变形假设, 微观弹脆性模型的控制方程可以写为[28]
$$ \rho \left( {\boldsymbol{x}} \right){\boldsymbol{\ddot u}}\left( {{\boldsymbol{x}},t} \right) = \int_{{{{H}}_x}} {C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\left[ {{{\boldsymbol{u}}}\left( {{\boldsymbol{x'}},t} \right) - {{\boldsymbol{u}}}\left( {{\boldsymbol{x}},t} \right)} \right]{\text{d}}V'} $$ (9) 其中,
$C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right) = \dfrac{c}{{\left\| {{\boldsymbol{\xi}} } \right\|}}$ , 体力项为0. 基于式(9), 采用谱分析方法[39], 先对一维情况进行研究, 再扩展到二维, 探究波的空间传播特性.2.1 一维谱分析
对于一维问题, 式(9)可以写为
$$ \rho \left( x \right)\ddot u\left( {x,t} \right) = \int_{{{ - }}\delta }^\delta {C\left( {\left\| \xi \right\|} \right)\left[ {u\left( {x{\text{ + }}\xi ,t} \right) - u\left( {x,t} \right)} \right]A{\text{d}}\xi } $$ (10) 其中,
$ A $ 为一维模型的横截面积. 为了研究波的传播特性, 有必要在频域内展开分析. 这里关注线弹性问题, 因此, 可对位移项进行傅里叶展开$$ u\left( {x,t} \right) = a\exp \left[ {{\rm{i}}\kappa \left( {x - vt} \right)} \right] $$ (11) 其中,
${\rm{i}}$ 为虚数单位, 满足${{\rm{i}}^2}{\text{ = }} - 1$ ;$ a $ 为幅值;$ v $ 为相速度;$ \kappa $ 为$2\text{π}$ 长度内的波数. 将式(11)代入式(10), 可约去时间项得$$ - \rho {\kappa ^2}{v^2} = A\int_{{{ - }}\delta }^\delta {C\left( {\left\| \xi \right\|} \right)\left[ {\exp \left( {{\rm{i}}\kappa \xi } \right) - 1} \right]{\text{d}}\xi } $$ (12) 进一步, 考虑邻域的对称性, 可约去奇数项得
$$ - \rho {\kappa ^2}{v^2} = A\int_{{{ - }}\delta }^\delta {C\left( {\left\| \xi \right\|} \right)\left[ {\cos \left( {\kappa \xi } \right) - 1} \right]{\text{d}}\xi } $$ (13) 若定义与波数有关的项为
$$ {M_\kappa } = A\int_{{{ - }}\delta }^\delta {C\left( {\left\| \xi \right\|} \right)\left[ {1 - \cos \left( {\kappa \xi } \right)} \right]{\text{d}}\xi } $$ (14) 则式(13)可视为相速度
$ v $ 的方程$$ {v^2} - \frac{{{M_\kappa }}}{{\rho {\kappa ^2}}} = 0 $$ (15) 求解可得相速度为
$$ v{\text{ = }}\sqrt {\frac{{{M_\kappa }}}{{\rho {\kappa ^2}}}} $$ (16) 显然, 式(16)中的虚数项为0, 故近场动力学系统中没有数值耗散. 当波色散产生不利影响时, 如数值振荡、非物理变形等, 系统本身无法对其进行抑制. 随着不利因素的逐渐累积, 数值计算的精度将会受到严重的影响, 无法保证对固体行为的正确捕捉. 这一方面说明了进行近场动力学系统中波色散研究的必要性, 另一方面为抑制色散现象提供了方向.
进一步, 考虑数值计算中的空间离散, 对近场动力学系统的色散特性进行分析. 根据邻域的对称性, 经过数值离散, 式(14)变为
$$ {M_\kappa } = \sum\limits_{j = 1}^{\left( {m - 1} \right)/2} {\frac{{2c}}{{\left\| {{\xi _j}} \right\|}}\left[ {1 - \cos \left( {\kappa hj} \right)} \right]Ah} $$ (17) 其中,
$ m $ 为一维邻域内的粒子数,$ h $ 为粒子间距,$ c $ 为微弹性模量. 采用式(7)中的微弹性模量, 式(17)可写为$$ {M_\kappa } = \frac{{4E}}{{{\delta ^2}}}\sum\limits_{j = 1}^{\left( {m - 1} \right)/2} {\frac{{\left[ {1 - \cos \left( {\kappa hj} \right)} \right]}}{j}} $$ (18) 将式(18)代入式(16), 可得一维近场动力学系统的色散关系
$$ \left| {\frac{v}{{{v_0}}}} \right| = \frac{2}{{\kappa \delta }}\sqrt {\sum\limits_{j = 1}^{\left( {m - 1} \right)/2} {\frac{{\left[ {1 - \cos \left( {\kappa hj} \right)} \right]}}{j}} } $$ (19) 其中,
$ {v_0} $ 为弹性理论波速$$ {v_0} = \sqrt {\frac{E}{\rho }} $$ (20) 其中,
$ E $ 为弹性模量,$ \rho $ 为密度. 由于$ \kappa \delta $ 和$ \kappa {\xi _i} $ 均与$ \kappa h $ 有关, 式(19)可视为${N_h} = {{\kappa h} \mathord{\left/ {\vphantom {{\kappa h} ({2\text{π} })}} \right. } ({2\text{π} })}$ 的函数, 表示粒子间距内的波数. 若取邻域半径为$ \delta {\text{ = }}3.015 h $ , 式(19)中的色散关系将如图2所示.由图2可知, 在一维近场动力学系统中, 相速度与理论波速的比值偏离局部理论的结果, 表现出色散行为. 局部理论各个频率成分的波的传播速度相同, 在传播时会同时达到. 而近场动力学中各个频率成分的波的传播速度不同, 在传播时会分先后达到. 对于传播速度小于理论波速的频率成分, 在传播中会落后于其他波, 表现为拖尾波. 因此, 定义临界点为与局部理论色散关系的交点, 即
$\dfrac{v}{{{v_0}}} = 1$ 的点. 可以看出, 大部分频率成分位于临界点以下, 在波的传播中表现为拖尾波.为区分色散行为和研究方便, 将
${N_h} \geqslant 1$ 定义为高频域,$0 \leqslant {N_h} < 1$ 定义为低频域. 即高频域与低频域的划分与波长有关, 高频域的波长范围为$\lambda \leqslant h$ , 低频域的波长范围为$\lambda > h$ , 其中$h$ 为粒子间距. 显然, 低频成分(虚线段)和高频成分(实线段)分别呈现出不同的色散特性. 在高频域, 波速不仅存在振荡的趋势, 还在某些频率处等于0. 根据文献[28-29], 这些点的相速度为0, 对应于零能模式点. 另有学者将零能模式定义为群速度为0的频散带, 在这种定义下, 首先有必要给出一维模型的频率表达式. 基于式(19), 有$$ \begin{split} &\left| {\frac{\omega }{{{\omega _0}}}} \right| = \frac{{2\kappa h}}{{\kappa \delta }}\sqrt {\sum\limits_{j = 1}^{\left( {m - 1} \right)/2} {\frac{{\left[ {1 - \cos \left( {\kappa hj} \right)} \right]}}{j}} } =\\ &\qquad\frac{{2h}}{\delta }\sqrt {\sum\limits_{j = 1}^{\left( {m - 1} \right)/2} {\frac{{\left[ {1 - \cos \left( {2\text{π} {N_h}j} \right)} \right]}}{j}} } \end{split}$$ (21) 显然, 由于余弦函数具有周期性, 频率将呈现在
${N_h}$ 取整区间内的周期性变化. 为了更加直观地说明此周期性, 将频率$\left| {\dfrac{\omega }{{{\omega _0}}}} \right|$ 与${N_h}$ 的关系绘制成图3.如图3中蓝色线所示, 频率随波数呈现周期性变化, 并在波数取整处为0. 与图2不同的是, 这种周期性在低频域(虚线)和高频域(实线)呈现相同的变化规律. 由于频率呈现从0到0的周期性变化, 对应的斜率应当呈现从正到负或从负到正的周期性变化. 而频率的斜率正比于群速度
$\dfrac{{\partial \omega }}{{\partial \kappa }}$ , 说明存在群速度为0的点. 为了进一步验证上述结论, 基于频率的表达式(21), 求导得到群速度与波数的关系, 绘制图像进行表示. 如图3中红色线所示, 群速度会在整波数区间内从正变为负, 并在从正变为负的过程中出现零点. 在这种情况下, 零能模式点对应于图3中红色线与y = 0水平线的交点.综上所述, 无论是采用何种形式的零能模式定义, 一维近场动力学模型中均存在零能模式点. 振荡趋势和零能模式均是非物理响应, 会加重色散行为带来的不利影响. 因此, 相比于低频成分, 高频成分的色散问题更为严重, 应当给予关注.
2.2 二维谱分析
下面, 考虑波的空间传播, 对二维系统展开分析. 此时, 除方向
$ {x_1} $ 外, 还需引入方向$ {x_2} $ , 控制方程(9)变为$$ \begin{split} &\rho \left( {{\boldsymbol{x}}} \right){\ddot u}\left( {{{\boldsymbol{x}}},t} \right) = \int_{{{ - }}\delta }^\delta \int_{ - \sqrt {{\delta ^2} - {r^2}} }^{\sqrt {{\delta ^2} - {r^2}} } C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\left[ {{\boldsymbol{u}}}\left( {{{\boldsymbol{x}}} + {{\boldsymbol{\xi}} },t} \right) -\right.\\ &\qquad \left.{{\boldsymbol{u}}}\left( {{{\boldsymbol{x}}},t} \right) \right]B {\text{d}}w{\text{d}}r\end{split} $$ (22) 其中,
$ B $ 为二维模型的厚度,${{\boldsymbol{\xi}} }$ 为相对位置$$ {{\boldsymbol{\xi}} }{\text{ = }}{{\boldsymbol{x}}'} - {{\boldsymbol{x}}} = r{{{\boldsymbol{x}}}_1} + w{{{\boldsymbol{x}}}_2} $$ (23) 同时, 位移项的傅里叶展开变为
$$ {{\boldsymbol{u}}}\left( {{{\boldsymbol{x}}},t} \right) = {{\boldsymbol{a}}}\exp \left[ {{\rm{i}}\kappa \left( {{{\boldsymbol{x}}} \cdot {{\boldsymbol{b}}} - vt} \right)} \right] $$ (24) 其中,
${{\boldsymbol{a}}}$ 为与位移${{\boldsymbol{u}}}$ 共线的幅度向量;${{\boldsymbol{b}}}$ 为波传播方向的单位矢量; 其他变量的定义同式(11). 将式(24)代入式(22), 可同时约去矢量项和时间项, 有$$ - \rho {\kappa ^2}{v^2} = \int_{{{ - }}\delta }^\delta {\int_{ - \sqrt {{\delta ^2} - {r^2}} }^{\sqrt {{\delta ^2} - {r^2}} } {C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\left[ {\exp \left( {{\rm{i}}\kappa {{\boldsymbol{\xi}} } \cdot {{\boldsymbol{b}}}} \right) - 1} \right]B} {\text{d}}w{\text{d}}r} $$ (25) 同样地, 根据邻域的对称性, 可约去奇数项得
$$ - \rho {\kappa ^2}{v^2} = \int_{{{ - }}\delta }^\delta {\int_{ - \sqrt {{\delta ^2} - {r^2}} }^{\sqrt {{\delta ^2} - {r^2}} } {C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\left[ {\cos \left( {\kappa {{\boldsymbol{\xi}} } \cdot {{\boldsymbol{b}}}} \right) - 1} \right]B} {\text{d}}w{\text{d}}r} $$ (26) 仍可将与波数有关的项用一个变量进行表示
$$ {N_\kappa } = \int_{{{ - }}\delta }^\delta {\int_{ - \sqrt {{\delta ^2} - {r^2}} }^{\sqrt {{\delta ^2} - {r^2}} } {C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\left[ {1 - \cos \left( {\kappa {{\boldsymbol{\xi}} } \cdot {{\boldsymbol{b}}}} \right)} \right]B} {\text{d}}w{\text{d}}r} $$ (27) 此时, 关于波速的方程可写为
$$ {v^2} - \frac{{{N_\kappa }}}{{\rho {\kappa ^2}}} = 0 $$ (28) 解得相速度为
$$ v{\text{ = }}\sqrt {\frac{{{N_\kappa }}}{{\rho {\kappa ^2}}}} $$ (29) 与一维的谱分析结果相同, 相速度的虚数项为0, 没有数值耗散, 无法抑制色散带来的不利影响; 实数项与波数有关, 表现出色散行为.
进一步, 考虑数值计算, 对
$ {N_\kappa } $ 进行空间离散, 有$$ {N_\kappa } = \sum\limits_{j = 1}^n {\frac{c}{{\left\| {{{{\boldsymbol{\xi}} }_j}} \right\|}}} \left[ {1 - \cos \left( {\kappa {{{\boldsymbol{\xi}} }_j} \cdot {{\boldsymbol{b}}}} \right)} \right]B{h^2} $$ (30) 其中,
$ n $ 为二维邻域内的粒子数,$ h $ 为粒子间距,${{{\boldsymbol{\xi}} }_j} \cdot {{\boldsymbol{b}}}$ 为相对位置在波传播方向上的投影. 仍采用式(7)中的微弹性模量, 式(30)可写为$$ {N_\kappa } = \frac{{9E}}{{\text{π} {\delta ^3}}}\sum\limits_{j = 1}^n {\frac{{\left[ {1 - \cos \left( {\kappa {{{\boldsymbol{\xi}} }_j} \cdot {{\boldsymbol{b}}}} \right)} \right]}}{{\left\| {{{{\boldsymbol{\xi}} }_j}} \right\|}}} {h^2} $$ (31) 将式(31)代入式(29), 可以得到二维近场动力学系统的色散关系
$$ \left| {\frac{v}{{{v_0}}}} \right| = \frac{3}{{\kappa \delta }}\sqrt {\frac{h}{{\text{π} \delta }}\sum\limits_{j = 1}^n {\frac{h}{{\left\| {{{{\boldsymbol{\xi}} }_j}} \right\|}}\left[ {1 - \cos \left( {\kappa {{{\boldsymbol{\xi}} }_j} \cdot {{\boldsymbol{b}}}} \right)} \right]} } $$ (32) 其中,
$ {v_0} $ 为式(20)定义的弹性理论波速. 显然, 在二维问题中, 波的色散行为不仅与频率有关, 还与波的传播方向有关.如图4所示, 对于对称邻域, 对关于
$ {x_2} $ 轴对称的传播方向${{\boldsymbol{b}}'}$ , 一定可以同样对称地找到关于$ {x_2} $ 轴对称的相对位置向量${{{\boldsymbol{\xi}} }_j}^\prime$ , 使式(30)中的投影相等$$ {{{\boldsymbol{\xi}} }_j}^\prime \cdot {{\boldsymbol{b}}' = }\left( { - \xi _j^{{x_1}}} \right)\cos \left( {\text{π} - \theta } \right) + \xi _j^{{x_2}}\sin \theta = {{{\boldsymbol{\xi}} }_j} \cdot {{\boldsymbol{b}}} $$ (33) 由于邻域内每个相对位置向量
${{{\boldsymbol{\xi}} }_j}$ , 都存在关于$ {x_2} $ 轴对称的相对位置向量${{{\boldsymbol{\xi}} }_j}^\prime$ 满足投影相等关系, 故最终求和得到的式(30)也相等. 同样地, 对于$ {x_1} $ 轴也满足此对称性$$ {{{\boldsymbol{\xi}} }_j}^\prime \cdot {{\boldsymbol{b}}' = }\xi _j^{{x_1}}\cos \theta + \left( { - \xi _j^{{x_2}}} \right)\sin \left( { - \theta } \right) = {{{\boldsymbol{\xi}} }_j} \cdot {{\boldsymbol{b}}} $$ (34) 因此, 可将波的传播方向缩小至
$ \theta = 0 \sim {90^ \circ } $ 范围内进行研究. 显然, 在$ {0^ \circ } \sim {90^ \circ } $ 范围内, 关于$ {45^ \circ } $ 方向存在对称性$$ \xi _j^{{x_2}}\cos \left( {\frac{\text{π} }{2} - \theta } \right) + \xi _j^{{x_1}}\sin \left( {\frac{\text{π} }{2} - \theta } \right) = \xi _j^{{x_1}}\cos \theta + \xi _j^{{x_2}}\sin \theta $$ (35) 因此, 二维键基模型的色散行为将表现出沿
$ {45^ \circ } $ 的对称性.基于此对称性, 对波传播方向与
$ {x_1} $ 方向的夹角在$ \theta = {0^ \circ } \sim {45^ \circ } $ 之间的色散行为进行研究. 这里, 仍取${N_h} = {{\kappa h} \mathord{\left/ {\vphantom {{\kappa h} ({2\text{π} })}} \right. } ({2\text{π} })}$ 为横坐标, 并取波的传播方向分别为$ \theta = {0^ \circ },{30^ \circ },{45^ \circ } $ , 式(32)中的色散关系将如图5所示. 由图5可知, 当波沿$ {x_1} $ 方向传播时, 二维近场动力学系统中色散关系的变化趋势与一维系统较为一致, 即在低频成分和高频成分表现出不同的色散行为, 并在高频域呈现出振荡趋势和零能模式. 但随着传播方向改变, 高频成分色散关系的振荡程度发生变化, 零能模式点的位置也在移动.为了更加直观地表征不同方向的频散行为, 提取高频成分色散关系变化率的最大值作为振荡程度的衡量指标, 将
$0^\circ\sim{360^\circ }$ 范围内高频成分色散关系的振荡程度绘制雷达图进行对比分析. 如图6所示, 对于对称邻域, 高频成分色散关系的变化趋势呈现出关于$ {45^ \circ } $ 方向的对称性, 验证了前述分析的正确性. 在$0^\circ\sim{45^\circ }$ 方向范围内, 随着波的传播方向改变, 振荡指标的大小发生变化, 与图5的分析结果相一致.综合上述分析, 相比于低频成分, 高频成分的色散行为更为严重, 表现出振荡趋势和零能模式. 高频域的色散关系会随波的传播方向变化, 呈现出沿45°的对称性. 而近场动力学系统本身没有数值耗散, 无法抑制色散带来的不利影响, 会严重影响对固体行为的正确捕捉.
3. 近场动力学方法的改进
事实上, 作为一种非局部理论, 近场动力学在高频段呈现频散行为是非局部效应的合理体现. 对于非均质材料, 频散行为是真实存在的, 可以通过近场动力学来模拟非均质材料中的非线性色散关系, 且近场作用半径与非均质材料的特征尺寸具有强相关性[29]. 但是, 对于变异性不大的材料, 通常可以认为是均质材料, 近场动力学的频散行为就会带来不利影响. 根据第2章的研究, 近场动力学的频散行为在高频域表现出振荡趋势. 也就是说, 理论上应该同时达到的不同频率成分的波, 在近场动力学的计算中传播速度会出现差异, 这种差异会在数值模拟中引起非物理变形和损伤.
对于极端载荷作用下的问题, 一方面这种不利影响会加重, 严重影响数值模拟的精度. 另一方面, 极端载荷作用往往伴随着裂纹的产生, 而裂纹的扩展与波的传播紧密相关, 异常的波传播过程必定会对裂纹扩展的预测产生不利影响. 而极端载荷作用下工程结构的安全性问题是国防安全领域的研究重点, 近场动力学方法的优势正是对于裂纹扩展捕捉的能力. 因此, 有必要抑制近场动力学中的频散行为, 提高其对波传播的模拟精度, 进而保证其对极端载荷作用下计算的可靠性, 为国防安全邻域研究提供理论指导和技术支撑. 本章将从引入数值耗散的角度出发, 对近场动力学体系的色散行为进行抑制. 并考虑高频域色散更为严重的特点, 对高频成分进行选择性抑制.
3.1 黏性效应的研究
在固体力学中, 通常认为黏性与数值耗散有关, 为近场动力学中耗散机制的引入提供了方向. 首先, 有必要对黏性效应进行研究, 证明黏性引入的有效性. 这里, 仍从微观弹脆性模型入手. 考虑到黏性与应变率的相关性, 可将黏性项取为力密度函数的率形式
$$ {{\boldsymbol{f}}_{{\text{vs}}}}\left( {{\boldsymbol{x}},{\boldsymbol{x'}}} \right) = \gamma {\boldsymbol{\dot f}}\left( {{\boldsymbol{x}},{\boldsymbol{x'}}} \right) $$ (36) 其中,
$ \gamma $ 为无量纲系数. 在小变形下, 有$$ \begin{split} & {{\boldsymbol{f}}_{{\text{vs}}}}\left( {{\boldsymbol{x}},{\boldsymbol{x'}}} \right) = \gamma C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\frac{{{\text{d}}\left[ {{{\boldsymbol{u}}}\left( {{{\boldsymbol{x}}}{\text{ + }}{{\boldsymbol{\xi}} },t} \right) - {{\boldsymbol{u}}}\left( {{{\boldsymbol{x}}},t} \right)} \right]}}{{{\text{d}}t}}= \\ &\qquad D\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\left[ {{\dot {\boldsymbol{u}}}\left( {{{\boldsymbol{x}}}{\text{ + }}{{\boldsymbol{\xi}} },t} \right) - {\dot {\boldsymbol{u}}}\left( {{{\boldsymbol{x}}},t} \right)} \right] \end{split} $$ (37) 其中,
$D\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right){\text{ = }}\gamma C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right){\text{ = }}\dfrac{{\gamma c}}{{\left\| {{\boldsymbol{\xi}} } \right\|}}$ . 引入黏性效应, 控制方程(9)变为$$ \begin{split} & \rho \left( {\boldsymbol{x}} \right){{\ddot {\boldsymbol{u}}}}\left( {{\boldsymbol{x}},t} \right) = \int_{{{{H}}_x}} {C\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\left[ {{{\boldsymbol{u}}}\left( {{{\boldsymbol{x}}}{\text{ + }}{{\boldsymbol{\xi}} },t} \right) - {{\boldsymbol{u}}}\left( {{\boldsymbol{x}},t} \right)} \right]}{\text{ + }} \\ &\qquad D\left( {\left\| {{\boldsymbol{\xi}} } \right\|} \right)\left[ {{\dot {\boldsymbol{u}}}\left( {{{\boldsymbol{x}}}{\text{ + }}{{\boldsymbol{\xi}} },t} \right) - {\dot {\boldsymbol{u}}}\left( {{{\boldsymbol{x}}},t} \right)} \right]{\text{d}}V{'} \end{split} $$ (38) 以一维问题为例, 式(38)可化简为
$$ \begin{split} & \rho \left( x \right)\ddot u\left( {x,t} \right) = A\int_{{{ - }}\delta }^\delta {C\left( {\left\| \xi \right\|} \right)\left[ {u\left( {x{\text{ + }}\xi ,t} \right) - u\left( {x,t} \right)} \right]}{\text{ + }} \\ &\qquad D\left( {\left\| \xi \right\|} \right)\left[ {\dot u\left( {x{\text{ + }}\xi ,t} \right) - \dot u\left( {x,t} \right)} \right]{\text{d}}\xi \end{split} $$ (39) 其中,
$ A $ 为一维模型的横截面积. 沿用第2章的谱分析方法, 采用式(11)中的傅里叶展开, 式(39)可写为$$ \begin{split} & - \rho {\kappa ^2}{{v'}^2} = A\int_{{{ - }}\delta }^\delta {\left[ {C\left( {\left\| \xi \right\|} \right) - {\rm{i}}\kappa v'D\left( {\left\| \xi \right\|} \right)} \right]\left[ {\exp \left( {{\rm{i}}\kappa \xi } \right) - 1} \right]{\text{d}}\xi }= \\ &\qquad A\int_{ - \delta }^\delta {\left[ {C\left( {\left\| \xi \right\|} \right) - {\rm{i}}\kappa v'D\left( {\left\| \xi \right\|} \right)} \right]\left[ {\cos \left( {\kappa \xi } \right) - 1} \right]{\text{d}}\xi }\\[-12pt] \end{split} $$ (40) 其中,
$ v' $ 为黏性作用下的相速度. 仍采用式(14)中对波数相关项的定义, 式(40)可表示为相速度的方程$$ {v'^2} + \frac{{\gamma {M_\kappa }}}{{\rho \kappa }}{\rm{i}}v' - \frac{{{M_\kappa }}}{{\rho {\kappa ^2}}} = 0 $$ (41) 求解可得相速度
$$ v' = {v_i}{\rm{i}} + {v_r} = - \frac{{\gamma {M_\kappa }}}{{2\rho \kappa }}{\rm{i}} + \sqrt {\frac{{{M_\kappa }}}{{\rho {\kappa ^2}}} - {{\left( {\frac{{\gamma {M_\kappa }}}{{2\rho \kappa }}} \right)}^2}} $$ (42) 此时, 虚数项不再为0
$$ {v_i} = - \frac{{\gamma {M_\kappa }}}{{2\rho \kappa }} $$ (43) 说明黏性项能够引入数值耗散, 这与文献[40-43]的研究结论相一致, 证明了黏性引入的有效性.
进一步, 考虑高频域色散更为严重的特点, 对黏性项在不同频率成分的表现进行研究. 由于数值耗散的引入是为了抑制色散行为带来的不利影响, 这里关注虚数项(43)与原始相速度(16)的比值. 借助式(18), 比值可以表示为
$$ \left| {\frac{{{v_i}}}{v}} \right| = \frac{{\gamma {M_\kappa }}}{{2\rho \kappa }}\sqrt {\frac{{\rho {\kappa ^2}}}{{{M_\kappa }}}} = \frac{{\gamma {v_0}}}{\delta }\sqrt {\sum\limits_{i = 1}^{\left( {m - 1} \right)/2} {\frac{{\left[ {1 - \cos \left( {\kappa h{{i}}} \right)} \right]}}{{{i}}}} } $$ (44) 其中,
$ {v_0} $ 为式(20)定义的弹性理论波速,$ \gamma $ 为无量纲系数,$ \delta $ 为邻域半径,$ h $ 为粒子间距. 仍取${N_h} = {{\kappa h} \mathord{\left/ {\vphantom {{\kappa h} ({2\text{π} })}} \right. } ({2\text{π} })}$ 为横坐标, 归一化后式(44)中的比值关系将如图7所示. 由图7可知, 数值耗散对色散行为的抑制作用在低频域(虚线段)和高频域(实线段)呈现相同的周期变化, 即黏性效应对低频和高频成分的表现相同. 实际上, 式(37)本质上是一种线性黏性, 无法实现对不同频率成分的选择性抑制.3.2 黏性效应的引入
因此, 可将黏性效应引入近场动力学, 抑制波色散带来的不利影响. 并有必要在黏性项的构造中, 考虑对高频成分的选择性抑制. 基于近场动力学的非局部性和键基模型的键间独立性假设, 将在键层次以非局部作用的方式引入黏性效应. 相应地, 将控制方程定义为
$$ \rho \left( {\boldsymbol{x}} \right){\boldsymbol{\ddot u}}\left( {{\boldsymbol{x}},t} \right) = \int_{{{{H}}_x}} {{\boldsymbol{f}}\left( {{\boldsymbol{x}},{\boldsymbol{x'}}} \right){\text{ + }}{{\boldsymbol{f}}_{{\text{vs}}}}\left( {{\boldsymbol{x}},{\boldsymbol{x'}}} \right){\text{d}}V'} + {\boldsymbol{b}}\left( {{\boldsymbol{x}},t} \right) $$ (45) 其中,
$ {\boldsymbol{f}} $ 为力密度函数, 是传统力态;$ {{\boldsymbol{f}}_{{\text{vs}}}} $ 为黏性力态. 同时, 考虑到在键基模型中, 力密度函数的方向为变形后键的方向, 故满足线动量守恒$$ - {\boldsymbol{f}}\left( {{{\boldsymbol{\eta}} },{{\boldsymbol{\xi}} }} \right) = {\boldsymbol{f}}\left( { - {{\boldsymbol{\eta}} }, - {{\boldsymbol{\xi}} }} \right) $$ (46) 和角动量守恒
$$ {\boldsymbol{f}}\left( {{{\boldsymbol{\eta}} },{{\boldsymbol{\xi}} }} \right) \otimes \left( {{{\boldsymbol{\eta}} }{\text{ + }}{{\boldsymbol{\xi}} }} \right) = 0 $$ (47) 为了保持上述守恒关系, 取黏性力态的方向与力密度函数的方向相同. 此时, 黏性力态可以分解为标量项和方向项
$$ {{\boldsymbol{f}}_{{\text{vs}}}}\left( {{\boldsymbol{x}},{\boldsymbol{x'}}} \right) = {f_{{\text{vs}}}}\left( {{\boldsymbol{x}},{\boldsymbol{x'}}} \right){{\boldsymbol{n}}} $$ (48) 其中,
$ {f_{{\text{vs}}}} $ 为黏性力态的标量形式,${{\boldsymbol{n}}}$ 为变形后键的方向. 显然, 这种引入方式不仅保留了传统近场动力学模型的理论框架, 而且只需要给出黏性力态标量形式的定义.下面, 将针对体积变形, 这一常见的固体变形模式, 构造黏性力态的标量表达式. 参照连续介质力学中体积黏性的定义[44-45]
$$ p = \beta \rho {c_{\text{d}}}{L_{\text{e}}}{\dot \varepsilon _{{\text{vol}}}} $$ (49) 其中,
$ \beta $ 为无量纲体积黏性系数;$ {c_{\text{d}}} $ 为膨胀波速;$ {L_{\text{e}}} $ 为单元的特征尺寸;$ {\dot \varepsilon _{{\text{vol}}}} $ 为体积应变率$$ {\dot \varepsilon _{{\text{vol}}}}{\text{ = }}{\dot \varepsilon _{ii}} = \nabla \cdot {\dot {\boldsymbol{u}}} $$ (50) 其中,
${{\boldsymbol{u}}}$ 为位移场. 定义体积黏性力态的标量形式为$$ {f_1} = {b_1}\frac{{\rho {v_0}}}{{{V_{\text{H}}}}}\zeta $$ (51) 其中,
$ {b_1} $ 为无量纲体积黏性系数;$ {v_0} $ 为式(20)定义的弹性理论波速;$ {V_{\text{H}}} $ 为邻域的体积;$ \zeta $ 为应变率. 在键基模型中, 键的相对伸长量$ s $ 为一种应变测量, 可用于应变率的计算$$ \zeta {\text{ = }}\dot s = \frac{{{\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}}}{{\left\| {{\boldsymbol{\xi}} } \right\|}} $$ (52) 如图8所示,
$\left\| {{\boldsymbol{\xi}} } \right\|$ 为键的初始长度(蓝线),${\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}$ 为键的速度在变形后键的方向上的投影(红线).基于上述定义, 可得体积黏性力态的标量表达式
$$ {f_1} = {b_1}\frac{{\rho {v_0}}}{{{V_{\text{H}}}}}\frac{{{\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}}}{{\left\| {{\boldsymbol{\xi}} } \right\|}}{\text{ = }}{b_1}{\alpha _{\text{p}}}\frac{{{\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}}}{{\left\| {{\boldsymbol{\xi}} } \right\|}} $$ (53) 其中,
${\alpha _{\text{p}}}{\text{ = }}\dfrac{{\rho {v_0}}}{{{V_{\text{H}}}}}$ . 将式(53)代入式(48), 当发生小变形时, 体积黏性作用下的近场动力学的控制方程为$$ \begin{split} & \rho \left( x \right)\ddot u\left( {x,t} \right) = A\int_{{{ - }}\delta }^\delta {C\left( {\left\| \xi \right\|} \right)\left[ {u\left( {x{\text{ + }}\xi ,t} \right) - u\left( {x,t} \right)} \right]}{\text{ + }} \\ &\qquad \frac{{{b_1}{\alpha _{\text{p}}}}}{{\left\| {\xi } \right\|}}\left[ {\dot u\left( {x{\text{ + }}\xi ,t} \right) - \dot u\left( {x,t} \right)} \right]{\text{d}}\xi \end{split} $$ (54) 当
$ {b_1}{\alpha _{\text{p}}}{\text{ = }}\gamma c $ 时, 式(54)与式(39)相同, 证明了3.1中黏性效应分析的正确性.由3.1中黏性效应分析可知, 式(53)本质上是一种线性黏性, 无法实现对不同频率成分的选择性抑制. 考虑到高频成分波色散的严重性, 有必要对其进行非线性改进. 这里, 将引入二次项, 以放大对高频成分的抑制作用
$$ {f_2} = {\text{sign}}\left( \zeta \right)\frac{{\rho \delta }}{{{V_{\text{H}}}}}{\left( {{b_2}\zeta } \right)^2}{\text{ = sign}}\left( {\frac{{{\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}}}{{\left\| {{\boldsymbol{\xi}} } \right\|}}} \right){\alpha _{\text{q}}}{\left( {{b_2}\frac{{{\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}}}{{\left\| {{\boldsymbol{\xi}} } \right\|}}} \right)^2} $$ (55) 其中,
$ {b_2} $ 为无量纲二次黏性系数,${\alpha _{\text{q}}}{\text{ = }}\dfrac{{\rho \delta }}{{{V_{\text{H}}}}}$ , 其余变量的定义与式(51)一致. 为充分发挥一次项和二次项的黏性效应, 将式(53)和式(55)进行多项式组合$$ {f_{{\text{vs}}}}{\text{ = }}{f_1}{\text{ + }}{f_2}{\text{ = }}{b_1}{\alpha _{\text{p}}}\frac{{{\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}}}{{\left\| {{\boldsymbol{\xi}} } \right\|}}{\text{ + sign}}\left( {\frac{{{\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}}}{{\left\| {{\boldsymbol{\xi}} } \right\|}}} \right){\alpha _{\text{q}}}{\left( {{b_2}\frac{{{\dot {\boldsymbol{\eta}} } \cdot {{\boldsymbol{n}}}}}{{\left\| {{\boldsymbol{\xi}} } \right\|}}} \right)^2} $$ (56) 其中,
$ {b_1} $ 和$ {b_2} $ 分别为一次项和二次项的黏性系数, 得到的黏性力态称为多项式黏性.式(48)给出了引入黏性的理论框架, 式(56)则从体积变形出发, 考虑了对高频成分的选择性抑制, 构造了多项式黏性力态, 完成了对近场动力学方法的改进.
4. 数值算例
本章基于波传播的数值算例, 开展数值模拟研究. 旨在对理论分析进行验证的同时, 考察波间断性的影响. 在近场动力学中, 通常采用EMU方法[46]进行数值计算, 该方法对控制方程进行空间离散
$$ \rho \left( {\boldsymbol{x}} \right){\boldsymbol{\ddot u}}\left( {{\boldsymbol{x}},t} \right) = \sum\limits_{i = 1}^N {{\boldsymbol{f}}\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_i}} \right)} \Delta {V_i} + {\boldsymbol{b}}\left( {{\boldsymbol{x}},t} \right) $$ (57) 其中,
$ N $ 为邻域内的粒子数目,$ \Delta V $ 为离散粒子的体积. 由于波的传播为动力问题, 采用显式中心差分方法进行求解$$ \left. \begin{gathered} {{{\dot {\boldsymbol{u}}}}_{t + 1}} = {{{\dot {\boldsymbol{u}}}}_t} + \frac{{{{{\ddot {\boldsymbol{u}}}}_t} + {{{\ddot {\boldsymbol{u}}}}_{t + 1}}}}{2}\Delta t \\ {{{\boldsymbol{u}}}_{t + 1}} = {{{\boldsymbol{u}}}_t} + {{{\dot {\boldsymbol{u}}}}_t}\Delta t + \frac{{{{{\ddot {\boldsymbol{u}}}}_t} + {{{\ddot {\boldsymbol{u}}}}_{t + 1}}}}{4}\Delta {t^2} \\ \end{gathered} \right\} $$ (58) 其中,
${{\boldsymbol{u}}}$ ,${\dot {\boldsymbol{u}}}$ 和${\ddot {\boldsymbol{u}}}$ 分别为加速度, 速度和位移;$ \Delta t $ 为满足稳定条件的时间步长.4.1 一维波传播问题
对一维算例, 首先考虑带理论解的波传播问题, 进行对比研究, 以验证理论分析的正确性. 再在此基础上, 考虑极端载荷下的激波传播问题, 探究波间断性的影响.
4.1.1 带理论解的波传播问题
为了验证理论分析的正确性, 并说明对不同间断性波的适用性, 这里取具有解析表达式的光滑正弦波和间断方波为研究对象, 对一维波的传播问题进行了研究. 将1 m长的弹性杆一端固定, 另一端输入目标波, 提取
$t = 77.3$ μs时沿杆长的位移分布, 将键基模型、线性黏性和多项式黏性的计算结果与理论解进行对比分析. 图9和图10分别给出的是正弦波和方波的计算结果. 显然, 两种波的键基模型计算结果均在波后出现了数值振荡. 结合理论分析结果, 一维近场动力学系统中存在波色散, 且大部分频率成分的波速小于理论波速, 表现为拖尾波的形式, 第二章中谱分析的正确性得到了验证.![]()
图 9 带理论解的一维正弦波传播结果对比Figure 9. Comparative study of 1D sine wave propagation problem with theoretical result![]()
图 10 带理论解的一维方波传播结果对比Figure 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)两类工况进行计算. 杆的材料参数为
$\rho = 8027\;{{{\text{kg}}}\mathord{\left/ {\vphantom {{{\text{kg}}} {{{\text{m}}^{\text{3}}}}}} \right. } {{{\text{m}}^{\text{3}}}}}$ ,$E = 193\;{\text{GPa}}$ . 将一维杆离散成1000个粒子进行计算, 邻域半径取$ \delta = 3.015 h $ , 采用压力时程模拟爆炸载荷, 并忽略负压区的影响.![]()
图 11 一维波传播算例示意图 (单位: mm)Figure 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 = \frac{L}{v} = \frac{L}{{\sqrt {{E \mathord{\left/ {\vphantom {E \rho }} \right. } \rho }} }} = \frac{1}{{\sqrt {{{1.93 \times {{10}^{11}}} \mathord{\left/ {\vphantom {{} {8027}}} \right. } {8027}}} }} = 2.04 \times {10^{ - 4}}\;{\text{s}} $$ (59) 计算结果表明, 爆炸载荷产生压力波, 压力波由杆的两端向中间传播. t=100 μs时在中点(A点)相遇后, 继续向前进行传播. t=200 μs时达到杆端后, 反射产生拉伸波, 向相反方向进行传播. 模拟得到的波的传播周期与式(59)相符, 验证了数值计算的正确性. 相比于图12, 图13中的应力波并不光滑, 在波后出现了数值振荡. 结合第2章的谱分析结果, 一维近场动力学系统中存在波色散, 且大部分频率成分的波速小于理论波速, 表现为拖尾波的形式, 谱分析的正确性得到了验证.
进一步, 提取A点的应力波时程, 对黏性效应进行研究. 如图14所示, 两端的压力波在
$ t = 100 $ μs时相遇叠加, 反射产生的拉伸波在$ t = 300 $ μs时再次相遇叠加, 波的传播周期与理论解相符. 对波峰处放大进行观察, 键基系统的计算结果在波后出现了数值振荡, 这与图13的计算结果相一致. 施加黏性后, 数值振荡得到了有效地抑制, 证明了黏性引入的有效性.对比两种黏性力态, 多项式黏性能够在抑制波后振荡的同时, 更好地保持原始波峰和波形, 得到更为合理的波传播结果. 这是因为数值耗散的引入会消耗波的能量, 从而影响波峰和波形. 而多项式黏性中含有二次项, 可以放大对高频成分的抑制作用, 在同等的抑制效果下减小了能量耗散, 降低了对原始波峰和波形的影响. 为了更直观地从能量角度进行说明, 对图14进行频域分析, 得到应力波的能量谱密度曲线. 如图15所示, 相比于体积黏性, 多项式黏性能够降低能量耗散, 且二者的差异主要集中在105 ~ 106的频率区间.
(2)大当量工况(W = 80 kg)
下面, 纳入对波间断性的考虑, 对大当量工况进行计算. 此时,
$ t = 50,100,150,200,250 $ μs时刻沿杆长方向的位移和应力分布分别如图16和图17所示. 显然, 相比于小当量工况, 此时的位移和应力的波形变得更为尖锐, 即波的间断性变强, 对应于激波的产生. 与小当量的计算结果相比, 位移也出现了波后振荡, 应力的数值振荡也更加明显. 说明强间断性波的色散行为更为显著, 呈现出Gibbs不稳定性. 因此, 色散在对波传播特性产生不利影响的同时, 还将引起数值不稳定性.仍提取点A的应力波时程进行进一步分析, 应力波时程和能量谱密度曲线分别如图18和图19所示. 与图17的计算结果相一致, 图18中波的间断性变强, 波后的振荡加剧. 与小当量工况相同, 黏性能够有效地抑制色散造成的数值振荡, 且多项式黏性能够更好地保持原始波峰和波形, 说明所提方法对激波同样适用. 图19则从能量角度说明了多项式黏性的优势, 即减小了能量耗散, 从而降低了对波峰和波形的影响.
综合一维计算结果, 近场动力学系统中的波色散在一维波传播问题中表现为波后振荡, 强间断性波的色散行为更为显著, 呈现出Gibbs不稳定性. 黏性能够有效地抑制波后振荡, 且多项式黏性能够更好地保持原始波峰和波形.
4.2 二维波传播问题
对二维算例, 考虑内径0.15 m, 外径0.3 m的圆环在爆炸载荷作用下的波传播问题. 如图20所示, 炸药在圆环的中心位置处爆炸, 距圆环内壁0.15 m. 圆环的材料参数为
$\rho = 2000\;{{{\text{kg}}} \mathord{\left/ {\vphantom {{{\text{kg}}} {{{\text{m}}^{\text{3}}}}}} \right. } {{{\text{m}}^{\text{3}}}}}$ ,$E = 30\;{\text{GPa}}$ . 这里, 取炸药质量W = 70 kg, 对大当量爆炸下的激波传播进行模拟. 在数值计算中, 粒子间距取$ h = 10\;{\text{mm}} $ , 邻域半径取$ \delta = 3.015 h $ , 爆炸载荷仍采用压力时程进行模拟, 并忽略负压区的影响.![]()
图 20 二维波传播算例示意图 (单位: mm)Figure 20. The schematic of numerical example for 2D wave propagation (unit: mm)中心爆炸下圆环中波的传播本质上是球形腔加压问题, 存在理论解[47]. 这里, 关注压力云图, 并采用对比研究的方式, 以考察色散对波的空间传播行为的影响. 由图21 (a1) ~ 图21(c1) 可知, 爆炸理论上会产生压力波, 压力波沿周向均匀分布, 并逐渐沿径向向外传播. 而在键基系统的计算结果中, 如图21 (a2) ~ 图21(c2)所示, 压力波的分布并不均匀, 呈现出沿45°的对称性. 结合第2章的谱分析结果, 二维系统中高频成分振荡的色散关系具有沿45°方向的对称性, 数值计算结果验证了谱分析的正确性.
此外, 相比于一维问题, 二维情况下的拖尾波现象更为显著, 拖尾波的峰值(3 GPa)甚至超过了波前峰值(1.5 GPa). 拖尾波的存在不仅会在波后引起非物理变形, 还会消耗波前峰值, 不利于对波传播过程的正确捕捉. 因此, 在二维情况下, 波色散会通过影响波的空间分布和传播过程, 对波的空间传播行为产生不利影响.
进一步, 观察多项式黏性的计算结果, 以探究黏性在二维问题中的表现. 如图21 (a3) ~ 图21(c3)所示, 此时的拖尾波得到了有效的抑制, 降低了对波前峰值的消耗, 波的传播过程得到了较好地捕捉. 同时, 压力波的分布更为均匀, 与理论解更为一致. 因此, 在黏性的作用下, 近场动力学系统中波的空间传播行为更为合理, 证明了黏性引入的有效性. 此时, 压力波的分布宽度相比理论解有所增加, 这应当与特征长度和数值耗散有关. 近场动力学引入了特征长度, 即邻域半径
$ \delta $ , 在特征长度范围内考虑粒子间的相互作用. 特征长度的存在会在一定程度上扩大影响范围, 使波的分布宽度变大. 另一方面, 由一维数值模拟结果可知, 数值耗散会消耗波的能量, 在一定程度上增大波宽. 但总体来说, 黏性的引入, 使得波前峰值和波的空间分布均得到了较好地捕捉, 证明了所提方法同样适用于多维问题.5. 结 论
本文采用谱分析方法, 对近场动力学系统中不同维度、不同频率成分的色散特性进行了研究, 并提出了相应的改进方法. 同时进行数值计算, 验证理论分析结果, 并探究波间断性的影响, 主要结论如下.
(1) 相比于低频成分, 高频成分的色散关系还呈现出振荡趋势和零能模式, 色散问题更为突出.
(2) 高频成分的色散行为会随波的传播方向发生变化, 呈现出沿45°的对称性.
(3) 色散在影响波传播特性的同时, 还会引起Gibbs不稳定性, 表现为强间断性波的色散现象更为显著.
(4) 黏性能够引入数值耗散, 抑制色散带来的不利影响. 二次项的引入可以放大对高频成分的抑制作用, 在同等的抑制效果下, 减小对原始波峰和波形的影响.
(5) 色散在一维波传播中引起波后振荡, 在二维问题中将引起更为严重的拖尾波和波空间分布的不均匀性.
本文研究可为减小近场动力学系统中波色散的不利影响, 获得正确的波传播过程和固体行为提供重要参考. 同时, 尚存在一些值得进一步研究的问题, 如在冲击问题中的应用、对零能模式的抑制等. 事实上, 若想纳入针对零能模式的抑制作用, 需要先识别出近场动力学中的零能变形模式, 再针对其进行抑制. 零能模式是一种不引起外力的虚假变形模式, 可以从零能模式的特点出发, 对近场动力学中的零能变形模式进行研究, 确定其变形特征. 再在本文提出的数值黏性的基础上, 根据研究出的近场动力学中零能模式的变形特征, 对变形相关项进行修改, 以实现对零能模式的抑制效果.
-
![]()
图 9 带理论解的一维正弦波传播结果对比
Figure 9. Comparative study of 1D sine wave propagation problem with theoretical result
![]()
图 10 带理论解的一维方波传播结果对比
Figure 10. Comparative study of 1D square wave propagation problem with theoretical result
![]()
图 11 一维波传播算例示意图 (单位: mm)
Figure 11. The schematic of numerical example for 1D wave propagation (unit: mm)
![]()
图 20 二维波传播算例示意图 (单位: mm)
Figure 20. The schematic of numerical example for 2D wave propagation (unit: mm)
-
[1] 高光发. 波动力学基础. 北京: 科学出版社, 2019 Gao Guangfa. Fundamentals of Wave Dynamics. Beijing: National Scientific Press, 2019 (in Chinese)
[2] 李永池. 波动力学. 合肥: 中国科学技术大学出版社, 2015 Li Yongchi. Wave Dynamics. Hefei: University of Science and Technology of China Press, 2015 (in Chinese)
[3] Cox BN, Gao H, Gross D, et al. Modern topics and challenges in dynamic fracture. Journal of the Mechanics and Physics of Solids, 2005, 53(3): 565-596 doi: 10.1016/j.jmps.2004.09.002
[4] Freund LB. Dynamic Fracture Mechanics. Cambridge: Cambridge University Press, 1998
[5] Ravi-Chandar K. Dynamic Fracture. New York: Elsevier, 2004
[6] Liang Y, Zhang X, Liu Y. Extended material point method for the three-dimensional crack problems. International Journal for Numerical Methods in Engineering, 2021, 122(12): 3044-3069 doi: 10.1002/nme.6653
[7] Wu J, Wang D, Lin Z, et al. An efficient gradient smoothing meshfree formulation for the fourth-order phase field modeling of brittle fracture. Computational Particle Mechanics, 2020, 7(2): 193-207 doi: 10.1007/s40571-019-00240-5
[8] Imtiaz H, Liu B. An efficient and accurate framework to determine the failure surface/envelop in composite lamina. Composites Science and Technology, 2021, 201: 108475 doi: 10.1016/j.compscitech.2020.108475
[9] Bobaru F, Foster JT, Geubelle PH, et al. Handbook of Peridynamic Modeling. CRC Press, 2016
[10] Hu Y, Feng G, Li S, et al. Numerical modelling of ductile fracture in steel plates with non-ordinary state-based peridynamics. Engineering Fracture Mechanics, 2020, 225: 106446 doi: 10.1016/j.engfracmech.2019.04.020
[11] Yang D, He X, Liu X, et al. A peridynamics-based cohesive zone model (PD-CZM) for predicting cohesive crack propagation. International Journal of Mechanical Sciences, 2020, 184: 105830 doi: 10.1016/j.ijmecsci.2020.105830
[12] Xia Y, Meng X, Shen G, et al. Isogeometric analysis of cracks with peridynamics. Computer Methods in Applied Mechanics and Engineering, 2021, 377: 113700 doi: 10.1016/j.cma.2021.113700
[13] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175-209 doi: 10.1016/S0022-5096(99)00029-0
[14] Yu H, Chen X, Sun Y. A generalized bond-based peridynamic model for quasi-brittle materials enriched with bond tension–rotation–shear coupling effects. Computer Methods in Applied Mechanics and Engineering, 2020, 372: 113405 doi: 10.1016/j.cma.2020.113405
[15] Silling SA, Epton M, Weckner O, et al. Peridynamic states and constitutive modeling. Journal of Elasticity, 2007, 88(2): 151-184 doi: 10.1007/s10659-007-9125-1
[16] Behzadinasab M, Foster JT. A semi-lagrangian constitutive correspondence framework for peridynamics. Journal of the Mechanics and Physics of Solids, 2020, 137: 103862 doi: 10.1016/j.jmps.2019.103862
[17] Zhang H, Qiao P. A two-dimensional ordinary state-based peridynamic model for elastic and fracture analysis. Engineering Fracture Mechanics, 2020, 232: 107040 doi: 10.1016/j.engfracmech.2020.107040
[18] Zhou XP, Wang YT, Shou YD, et al. A novel conjugated bond linear elastic model in bond-based peridynamics for fracture problems under dynamic loads. Engineering Fracture Mechanics, 2018, 188: 151-183 doi: 10.1016/j.engfracmech.2017.07.031
[19] Zhu QZ, Ni T. Peridynamic formulations enriched with bond rotation effects. International Journal of Engineering Science, 2017, 121: 118-129 doi: 10.1016/j.ijengsci.2017.09.004
[20] Gu X, Zhang Q, Madenci E. Non-ordinary state-based peridynamic simulation of elastoplastic deformation and dynamic cracking of polycrystal. Engineering Fracture Mechanics, 2019, 218: 106568 doi: 10.1016/j.engfracmech.2019.106568
[21] Liu Z, Bie Y, Cui Z, et al. Ordinary state-based peridynamics for nonlinear hardening plastic materials' deformation and its fracture process. Engineering Fracture Mechanics, 2019, 223: 106782
[22] 王涵, 黄丹, 徐业鹏等. 非常规态型近场动力学热黏塑性模型及其应用. 力学学报, 2018, 50(4): 810-819 (Wang Han, Huang Dan, Xu Yepeng, et al. Non-ordinary state-based peridynamic thermal-viscoplastic model and its application. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(4): 810-819 (in Chinese) doi: 10.6052/0459-1879-18-113 [23] Wang H, Xu YP, Hang D. A non-ordinary state-based peridynamic formulation for thermo-visco-plastic deformation and impact fracture. International Journal of Mechanical Sciences, 2019, 159: 336-344 doi: 10.1016/j.ijmecsci.2019.06.008
[24] 胡祎乐, 余音, 汪海. 基于近场动力学理论的层压板损伤分析方法. 力学学报, 2013, 45(4): 624-628 (Hu Yile, Yu Yin, Wang Hai. Damage analysis method for laminates based on peridynamic theory. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(4): 624-628 (in Chinese) doi: 10.6052/0459-1879-12-368 [25] 章青, 顾鑫, 郁杨天. 冲击载荷作用下颗粒材料动态力学响应的近场动力学模拟. 力学学报, 2016, 48(1): 56-63 (Zhang Qing, Gu Xin, Yu Yangtian. Peridynamics simulation for dynamic response of granular materials under impact loading. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1): 56-63 (in Chinese) doi: 10.6052/0459-1879-15-291 [26] Wu P, Yang F, Chen Z, et al. Stochastically homogenized peridynamic model for dynamic fracture analysis of concrete. Engineering Fracture Mechanics, 2021, 253: 107863 doi: 10.1016/j.engfracmech.2021.107863
[27] Huang XP, Kong XZ, Chen ZY, et al. Peridynamics modelling of dynamic tensile failure in concrete. International Journal of Impact Engineering, 2021, 155: 103918 doi: 10.1016/j.ijimpeng.2021.103918
[28] Bažant ZP, Luo W, Chau VT, et al. Wave dispersion and basic concepts of peridynamics compared to classical nonlocal damage models. Journal of Applied Mechanics, 2016, 83(11): 111004 doi: 10.1115/1.4034319
[29] Butt SN, Timothy JJ, Meschke G. Wave dispersion and propagation in state-based peridynamics. Computational Mechanics, 2017, 60(5): 725-738 doi: 10.1007/s00466-017-1439-7
[30] Wildman RA. Discrete micromodulus functions for reducing wave dispersion in linearized peridynamics. Journal of Peridynamics and Nonlocal Modeling, 2019, 1(1): 56-73 doi: 10.1007/s42102-018-0001-0
[31] Gu X, Zhang Q, Huang D, et al. Wave dispersion analysis and simulation method for concrete SHPB test in peridynamics. Engineering Fracture Mechanics, 2016, 160: 124-137 doi: 10.1016/j.engfracmech.2016.04.005
[32] Li S, Jin Y, Lu H, et al. Wave dispersion and quantitative accuracy analysis of bond-based peridynamic models with different attenuation functions. Computational Materials Science, 2021, 197: 110667 doi: 10.1016/j.commatsci.2021.110667
[33] Wildman RA, Gazonas GA. A finite difference-augmented peridynamics method for reducing wave dispersion. International Journal of Fracture, 2014, 190(1): 39-52
[34] Alebrahim R, Packo P, Zaccariotto M, et al. Improved wave dispersion properties in 1D and 2D bond-based peridynamic media. Computational Particle Mechanics, 2021, 9: 597-614
[35] Zhou G, Hillman M. A non-ordinary state-based Godunov-peridynamics formulation for strong shocks in solids. Computational Particle Mechanics, 2020, 7(2): 365-375 doi: 10.1007/s40571-019-00254-z
[36] Ren XD, Zhu JG. Temporally stabilized peridynamics methods for shocks in solids. Computational Mechanics, 2021, 69: 489-504
[37] Zhu F, Zhao J. Peridynamic modelling of blasting induced rock fractures. Journal of the Mechanics and Physics of Solids, 2021, 153: 104469 doi: 10.1016/j.jmps.2021.104469
[38] 黄丹, 章青, 乔丕忠等. 近场动力学方法及其应用. 力学进展, 2010, 40(4): 448-459 (Huang Dan, Zhang Qing, Qiao Pizhong, et al. A review of peridynamics (PD) method and its application. Advances in mechanics, 2010, 40(4): 448-459 (in Chinese) doi: 10.6052/1000-0992-2010-4-J2010-002 [39] Jia F, Gao Z, Don WS. A spectral study on the dissipation and dispersion of the WENO schemes. Journal of Scientific Computing, 2015, 63(1): 49-77 doi: 10.1007/s10915-014-9886-1
[40] Ren B, Fan H, Bergel GL, et al. A peridynamics–SPH coupling approach to simulate soil fragmentation induced by shock waves. Computational Mechanics, 2015, 55(2): 287-302 doi: 10.1007/s00466-014-1101-6
[41] Monaghan JJ. On the problem of penetration in particle methods. Journal of Computational physics, 1989, 82(1): 1-15 doi: 10.1016/0021-9991(89)90032-6
[42] Silling SA, Parks ML, Kamm JR, et al. Modeling shockwaves and impact phenomena with Eulerian peridynamics. International Journal of Impact Engineering, 2017, 107: 47-57 doi: 10.1016/j.ijimpeng.2017.04.022
[43] Lai X, Liu L, Li S, et al. A non-ordinary state-based peridynamics modeling of fractures in quasi-brittle materials. International Journal of Impact Engineering, 2018, 111: 130-146 doi: 10.1016/j.ijimpeng.2017.08.008
[44] Caramana EJ, Shashkov MJ, Whalen PP. Formulations of artificial viscosity for multi-dimensional shock wave computations. Journal of Computational Physics, 1998, 144(1): 70-97 doi: 10.1006/jcph.1998.5989
[45] Landshoff R. A numerical method for treating fluid flow in the presence of shocks. Los Alamos National Lab NM, 1955
[46] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Computers & Structures, 2005, 83(17-18): 1526-1535
[47] Achenbach J. Wave Propagation in Elastic Solids. New York: Elsevier, 2012
-
期刊类型引用(1)
1. 梁绍敏，冯云田，赵婷婷，王志华. 颗粒材料破碎行为数值分析方法研究综述. 力学学报. 2024(01): 1-22 . 
本站查看
其他类型引用(1)
下载:
计量
- 文章访问数: 1090
- HTML全文浏览量: 348
- PDF下载量: 165
- 被引次数: 2
