基于非常规态近场动力学对混凝土动态拉伸断裂的数值模拟研究

刘振华; 孔祥振; 洪建; 方秦

doi:10.11883/bzycj-2024-0485

基于非常规态近场动力学对混凝土动态拉伸断裂的数值模拟研究

doi: 10.11883/bzycj-2024-0485

陆军工程大学爆炸冲击防灾减灾国家重点试验室，江苏南京 210007

基金项目: 国家自然科学基金（52178515，52208513）

详细信息

作者简介:
刘振华（1999－　），男，硕士研究生，liuzhenhua1025@163.com

通讯作者:
孔祥振（1988－　），男，工学博士，副教授，ouckxz@163.com

中图分类号: O383
计量
- 文章访问数: 212
- HTML全文浏览量: 38
- PDF下载量: 75
- 被引次数: 0
出版历程
- 收稿日期: 2024-12-11
- 修回日期: 2025-02-17
- 网络出版日期: 2025-02-18

Numerical investigation on dynamic tensile fracture in concrete material by non-ordinary state-based peridynamics

State Key Laboratory of Disaster Prevention & Mitigation of Explosion & Impact, Army Engineering University of PLA, Nanjing 210007, Jiangsu, China

摘要

摘要: 为了准确预测冲击爆炸荷载作用下混凝土材料的动态拉伸断裂破坏，基于非常规态近场动力学理论框架，本文首先建立了修正的Monaghan人工体积黏性计算方法用于消除数值振荡，然后将前期建立的等效计算应变率方法植入到课题组前期研发的Kong-Fang混凝土材料模型中，用以准确计算应变率突变时的应变率效应。在此基础上，开展了一维杆中的弹性波传播的数值模拟，发现在力矢量状态上额外附加修正的Monaghan人工体积黏性力矢量状态，可有效地抑制由变形梯度近似导致的非物理数值振荡现象，进而讨论分析了人工体积黏性参数的影响并给出参数建议值。最后将模型用于混凝土试件层裂的数值模拟，对比分析了人工体积黏性、不同应变率效应计算方法对动态拉伸断裂预测结果的影响规律，数值模拟结果表明，准确预测混凝土材料动态拉伸断裂破坏需同时考虑修正的Monaghan人工体积黏性和等效计算应变率，建立的考虑修正的Monaghan人工体积黏性和等效计算应变率的非常规态近场动力学模型可以较好地预测裂缝位置和数量，为冲击爆炸荷载作用下混凝土材料动态拉伸断裂破坏的数值模拟提供了新思路。
- 动态断裂 /
- 修正的Monaghan人工体积黏性 /
- 等效计算应变率 /
- 非常规态近场动力学
Abstract: To accurately predict the dynamic tensile fracture in concrete materials subjected to impact and blast loadings, this study first establishes a modified Monaghan artificial bulk viscosity computation method within the framework of a non-ordinary state-based peridynamics (NOSB-PD) theory to eliminate numerical oscillations. Subsequently, the corrected strain-rate computation method, previously developed, is integrated into the Kong-Fang concrete material model, which was proposed earlier by the research group to calculate accurately the strain-rate effect during sudden changes. Based on the two methods above, numerical simulations of elastic wave propagation in a one-dimensional rod are conducted, and the results demonstrate that the additional inclusion of the modified Monaghan artificial bulk viscosity force vector state into the original force vector state can effectively suppress the non-physical numerical oscillations caused by the deformation gradient approximation. The superiority of the modified Monaghan artificial bulk viscosity is validated through comparative analysis with the original Monaghan artificial bulk viscosity. Furthermore, the influence of the modified Monaghan artificial bulk viscosity parameters is investigated, and recommended values for these parameters are provided. Finally, the aforementioned model is used to numerically simulate the spall test in concrete specimens, where the effects of including or excluding the modified Monaghan artificial bulk viscosity and different strain-rate computation methods on the prediction results of dynamic tensile fracture are compared and analyzed. The numerical simulation results demonstrate that accurately predicting the dynamic tensile fracture in concrete materials requires simultaneous consideration of the modified Monaghan artificial bulk viscosity and corrected strain-rate computation. The established non-ordinary state-based peridynamics model that accounts for both the modified Monaghan artificial bulk viscosity and corrected strain-rate computation demonstrates strong capabilities in predicting crack locations and quantities based on both qualitative and quantitative analysis metrics. This work provides new insights into the numerical simulation of dynamic tensile fracture in concrete materials under impact and blast loadings.
- Dynamic tensile fracture /
- Modified Monaghan artificial bulk viscosity /
- Corrected strain-rate /
- Non-ordinary state-based peridynamics

HTML全文

混凝土材料广泛应用于军民用防护工程中，在服役周期内可能会受到蓄意或偶然的冲击爆炸作用^[1]。冲击爆炸荷载作用近区，混凝土材料受高压高应变率作用并产生体积压碎、剪切流动等大变形动态破坏^[2-3]；而在冲击爆炸荷载中远区，应力波幅值随距离逐渐衰减，在混凝土结构自由表面反射为拉伸波并可能产生动态拉伸断裂破坏^[2-3]，宏观上表现为裂缝和震塌等破坏现象。对动态拉伸断裂的准确预测是阐明防护结构在冲击爆炸荷载作用下动态响应的重要基础^[4-5]，而随着计算能力的提升和数值算法的发展，数值模拟已成为预测动态拉伸断裂的重要技术手段。

目前对混凝土动态拉伸断裂破坏的数值模拟主要基于连续介质力学的有限元方法，通过在材料模型中引入损伤来预测拉伸断裂，即认为断裂破坏通过损伤带描述^[6]。但损伤和应变软化会导致偏微分动量守恒方程失去正定性（双曲性），导致动态断裂的模拟结果强烈依赖于网格尺寸，不具备预测能力^[6-7]。为解决上述问题，Kong等^[7]提出了基于损伤的非局部计算方法，可较好解决网格尺寸的敏感性问题。也有部分研究采用非连续方法，即当计算某一物理量（如等效塑性应变、最大主应变等）达到给定阈值时，删除单元来模拟裂缝^[8-9]，但该方法强烈依赖于单元删除准则和网格尺寸，也不具备预测能力。考虑到断裂破坏的大变形特性，有学者采用光滑粒子流体动力学（smoothed particle hydrodynamics, SPH）^[10]方法和光滑粒子伽辽金（smoothed particle Galerkin, SPG）^[11]算法对冲击爆炸荷载作用下混凝土材料动态拉伸断裂破坏开展了数值模拟研究，但也存在破坏准则（类似上述单元删除准则）和粒子间距（类似单元尺寸）的强敏感性问题。

混凝土动态拉伸断裂破坏是一个典型由连续到不连续的过程，而上述方法控制方程均为偏微分的动量守恒方程，依赖于连续性假设，与动态断裂的不连续特性相悖，这也是预测结果出现网格不收敛、破坏准则强敏感性问题的本质。注意到Silling等^[12]提出的近场动力学方法采用积分型控制方程，有效解决了上述偏微分控制方程在不连续（裂缝）处空间导数不存在的奇异性问题，为预测混凝土材料的动态拉伸断裂破坏提供了具有前景的方法。近场动力学方法包括键型近场动力学、常规态型近场动力学以及非常规态型近场动力学^[12-13]，其中非常规态近场动力学（non-ordinary state-based peridynamics, NOSB-PD）通过对非局部变形梯度的近似，可以将连续介质力学框架下的材料模型引入，非常便于描述冲击爆炸荷载作用下混凝土材料复杂的动态力学行为^[14]。然而，采用NOSB-PD预测混凝土材料动态拉伸断裂破坏，仍需要解决两个关键问题。

一是非局部变形梯度的近似会引起数值振荡和零能模式^[15]，通常表现为位移场的不规则振荡以及有位移但无应变能的非物理现象，这会导致对波传播预测的不准确，进而导致对动态拉伸断裂预测的不准确。已有研究通常采取附加额外力矢量状态的方法来抑制零能模式和数值振荡。Littlewood^[15]建立了当前构型中物质点（PD点）实际位置与预测位置的差值成比例的罚函数，并将其额外施加在力密度矢量状态上，有效地抑制了零能模式；Breitenfeld等^[16]通过在原始力密度矢量状态中额外附加三种（弹簧键力、平均位移状态和罚函数）不同的力密度矢量状态函数来抑制在裂缝尖端高应变梯度区的数值振荡，开展了在拉伸作用下一维杆的数值模拟并分析不同方法的优缺点并给出最佳控制参数值；Gu等^[17]提出了包含沙漏力的罚函数方法来抑制数值不稳定问题，通过简支弹性梁变形和二维杆弹性波传播问题的模拟验证了方法的可靠性。也有研究采取变形梯度的高阶修正方法，如Yaghoobi和Chorzepa^[18]将影响函数的泰勒级数展开引入了非局部变形梯度张量的高阶近似，较好地抑制了一维和二维问题下的零能模式，然而对于物体表面附近粒子仍然存在非物理振荡；Chen^[19]提出了一种键关联的变形梯度表达式，通过每个单独键附近物质点的变形状态来计算键关联的变形梯度，该方法能够较好抑制二维和三维问题下的数值不稳定；Gu等^[20]借鉴近场动力学微分算子的概念引入了变形梯度和力密度矢量状态的高阶表达式，提出了键关联的高阶NOSB-PD，有效地提高了数值精度，同时数值振荡也得到较好抑制。上述研究表明由于非局部变形梯度的近似计算方法引起的零能模式和数值不稳定问题，会导致应力波传播以及断裂破坏预测的不准确，因此合理的方法及其相应的参数至关重要。

二是对材料动态拉伸力学行为尤其是高应变率效应的准确描述。混凝土是应变率敏感材料，其动态拉伸强度会随着应变率的增加而增加。通常认为动态拉伸试验中观测的应变率效应由3种不同因素决定：（1）Stefan效应^[21]，在动态加载时混凝土中的自由水运动产生了黏性力，从而使混凝土的裂纹扩展阻力变大；（2）微裂缝的惯性效应^[22-23]，即微裂缝扩展存在速度上限，这会抑制裂缝发展；（3）结构惯性效应^[24-27]，即在动态加载混凝土试件时由泊松效应引起横向变形导致的惯性效应。Stefan效应和微裂缝惯性使得混凝土材料的应变率效应与时间相关的，特别表现在当应变率突然变化时，应力响应存在明显的迟滞效应，试验数据^[28-29]也证实迟滞效应的存在。而目前常用混凝土材料模型，如Holmquist-Johnson-Cook (HJC) 模型^[30]、Karagozian and Case (K&C) 模型^[31]和Riedel-Hiermaier-Thoma (RHT) 模型^[32]等，应变率效应多采用动态强度增强因子（dynamic increase factor, DIF）放大屈服面的简化方法，无法描述上述物理和试验观测的迟滞效应。而动态拉伸断裂过程中应变率急剧变化，忽略上述迟滞效应可能会带来对拉伸断裂预测的不准确。最近Kong等^[5]提出了一个新的等效计算应变率模型，较好解决了迟滞效应的描述，但该研究未能解决拉伸断裂不连续的描述问题。

为解决上述问题，并能够可靠预测混凝土材料的动态拉伸断裂破坏，本文中首先将修正的Monaghan人工体积黏性^[33]植入到NOSB-PD框架中，然后采用团队研发的Kong-Fang模型描述混凝土材料的动态力学行为，并对Kong-Fang模型中的应变率效应计算方法进行改进。进一步采用上述模型对一维杆中的弹性波传播进行了数值模拟，分析了人工体积黏性参数的影响并给出了参数建议值；然后将模型用于混凝土试件层裂的数值模拟，对比分析了不同应变率效应计算方法对预测结果的影响。

1. 考虑人工体积黏性和等效计算应变率的非常规态近场动力学模型

首先对非常规态近场动力学理论框架进行简要介绍，然后提出修正的Monaghan人工体积黏性用以消除应力波传播模拟中的数值振荡；进而简要介绍课题组提出的等效计算应变率方法并将其植入到课题组研发的混凝土材料Kong-Fang本构模型中，用以准确描述应变率急剧变化情况下的应力迟滞效应；最后介绍数值实现方法。

1.1 非常规态近场动力学模型

非常规态近场动力学理论通常将物体离散为一系列包含所有物性信息的物质点，其中每个物质点与特定邻域内的其他物质点之间存在相互作用，特定邻域称为近场范围，δ为近场范围半径。如图1所示，每个物质点由其坐标x_i (i = 1, 2, …… ∞) 表示，记物质点x_i的近场范围内物质点x_j的集合为 ${{\boldsymbol{H}}_{{{{\boldsymbol{x}}}_i}}} = \left\{ {{{{\boldsymbol{x}}}_i} \in R:\left| {{{{\boldsymbol{x}}}_j} - {{{\boldsymbol{x}}}_i}} \right| {{＜}} \delta } \right\}$ 。位置矢量状态 $\underline{X} \left( {{{\xi }_{ij}}} \right)$ 和位移矢量状态 $\underline{U} \left( {{{\xi }_{ij}}} \right)$ 分别定义为：

图 1 NOSB-PD理论中质点的运动

Figure 1. Kinematics of material points in NOSB-PD theory

下载: 全尺寸图片幻灯片

$\underline{X}\left(\xi_{ij}\right)={\boldsymbol{x}}_j-{\boldsymbol{x}}_i=\xi_{ij},\quad\quad\quad \underline{U}\left(\xi_{ij}\right)={\boldsymbol{u}}\left({\boldsymbol{x}}_j,t\right)-{\boldsymbol{u}}\left({\boldsymbol{x}}_i,t\right)={\boldsymbol{\eta}}_{ij}$

(1)

式中：ξ_ij为参考构型中物质点x_i和x_j之间的相对位置矢量，η_ij为当前构型中的相对位移矢量，u(x_i, t)和u(x_j, t)分别为t时刻物质点x_i和x_j的位移矢量状态。

在当前构型中物质点x_i和x_j在t时刻的位置矢量分别由y(x_i, t) = x_i + u(x_i, t)和y(x_j, t) = x_j + u(x_j, t)求得，其变形矢量状态 $\underline{Y} \left( {{{\xi }_{ij}}} \right)$ 定义为：

$\underline{Y} \left( {{{\xi }_{ij}}} \right) = {{\xi }_{ij}} + {{\eta }_{ij}} = {y}\left( {{{x}_j},t} \right) - {y}\left( {{{x}_i},t} \right) = \left( {{{x}_j} + {u}\left( {{{x}_j},t} \right)} \right) - \left( {{{x}_i} + {u}\left( {{{x}_i},t} \right)} \right)$

(2)

在NOSB-PD理论中^[13]，近场范围 ${{{\boldsymbol{H}}}_{{{x}_i}}}$ 中各物质点的积分型运动方程表达式如下：

$\rho \ddot {\boldsymbol{u}}\left( {{{{\boldsymbol{x}}}_i},t} \right) = \int\limits_{{{{\boldsymbol{H}}}_{{{x}_i}}}} {\left\{ {\underline{T} \left( {{{{\boldsymbol{x}}}_i},t} \right)\left( {{{\xi }_{ij}}} \right) - \underline{T}\left( {{{{\boldsymbol{x}}}_j},t} \right)\left( {{{\xi }_{ij}}} \right)} \right\}} {\text{d}}{V_{{{{\boldsymbol{x}}}_j}}} + {\boldsymbol{b}}\left( {{x_i},t} \right)$

(3)

式中：ρ为当前构型的物质质量密度， ${\ddot {\boldsymbol{u}}}\left( {{{{\boldsymbol{x}}}_i},t} \right)$ 为t时刻物质点x_i的加速度， $\underline{T}\left( {{{{\boldsymbol{x}}}_i},t} \right)\left( {{{\xi }_{ij}}} \right)$ 为物质点x_j施加在物质点x_i的力矢量状态，同样有 $\underline{T}\left( {{{{\boldsymbol{x}}}_j},t} \right)\left( {{{\xi }_{ji}}} \right)$ ， ${\text{d}}{V_{{{{\boldsymbol{x}}}_j}}}$ 为物质点x_j的体积微元，b(x_i, t)为体力密度矢量。

式（3）中的力矢量状态 $\underline{T}\left( {{{{\boldsymbol{x}}}_i},t} \right)\left( {{{\xi }_{ij}}} \right)$ 定义为：

(4)

式中： $\omega \left( {{{\xi }_{ij}}} \right) = {{\text{e}}^{{{ - {{\left| {{\xi _{ij}}} \right|}^2}} / {{\delta ^2}}}}}$ 为高斯影响函数。

式（4）中 ${{\boldsymbol{K}}}_{{{{\boldsymbol{x}}}_i}}^{ - 1}$ 为形状张量 ${{{\boldsymbol{K}}}_{{{{\boldsymbol{x}}}_i}}}$ 的逆矩阵，即：

${{{\boldsymbol{K}}}_{{{{\boldsymbol{x}}}_i}}} = \int\limits_{{{{\boldsymbol{H}}}_{{{x}_i}}}} {\omega \left( {{{\xi }_{ij}}} \right)} {{\xi }_{ij}} \otimes {{\xi }_{ij}}{\text{d}}{V_{{{x}_j}}}$

(5)

式中： ${{{\boldsymbol{P}}}_{{{{\boldsymbol{x}}}_i}}}$ 是与连续介质力学中柯西应力σ建立联系的第一类Piola-Kirchhoff应力张量，其表达式为：

${{{\boldsymbol{P}}}_{{{{\boldsymbol{x}}}_i}}} = \det \left( {{{{\boldsymbol{F}}}_{{{{\boldsymbol{x}}}_i}}}} \right){{\boldsymbol{\sigma F}}}_{{{{\boldsymbol{x}}}_i}}^{ - {\text{T}}}$

(6)

式中： ${{{\boldsymbol{F}}}_{{{{\boldsymbol{x}}}_i}}}$ 为非局部变形梯度， ${{\boldsymbol{F}}}_{{{{\boldsymbol{x}}}_i}}^{{\rm{ - T}}}$ 为 ${{{\boldsymbol{F}}}_{{{{\boldsymbol{x}}}_i}}}$ 转置的逆矩阵， ${{{\boldsymbol{F}}}_{{{{\boldsymbol{x}}}_i}}}$ 表示为：

${{{\boldsymbol{F}}}_{{{{\boldsymbol{x}}}_i}}} = {{{\boldsymbol{N}}}_{{{{\boldsymbol{x}}}_i}}}{{\boldsymbol{K}}}_{{{{\boldsymbol{x}}}_i}}^{ - 1}$

(7)

式中： ${{\boldsymbol{K}}}_{{{{\boldsymbol{x}}}_i}}^{ - 1}$ 为形状张量的逆矩阵， ${{{\boldsymbol{N}}}_{{{{\boldsymbol{x}}}_i}}}$ = $\displaystyle \int\limits_{{{{\boldsymbol{H}}}_{{{{\boldsymbol{x}}}_i}}}} {\omega \left( {{{\xi }_{ij}}} \right)} \underline{Y}\left( {{{\xi }_{ij}}} \right) \otimes {{\xi }_{ij}}{\text{d}}{V_{{{{\boldsymbol{x}}}_j}}}$ 。

1.2 修正的Monaghan人工体积黏性

注意到非局部变形梯度（式(7)）是近似得到，会引起明显的数值振荡，而在模拟应力波传播时的数值振荡会引起对拉伸断裂破坏模拟的不准确，提出修正的Monaghan人工体积黏性用以解决该问题。

Monaghan人工体积黏性^[33]广泛应用于SPH方法中，主要用于解决强间断冲击波不连续问题的模拟，其在更新的压力中额外附加黏性项，黏性项定义为：

${\varPi }_{ij}=\left\{\begin{array}{ll}{\overline{\rho }}_{ij}\left(-\alpha {\overline{c}}_{ij}{\phi }_{ij}+\beta {\phi }_{ij}^{2}\right)&\quad\quad {\boldsymbol{x}}_{ij}\cdot {\boldsymbol{v}}_{ij} {{＜}} 0\\ 0&\quad\quad {\boldsymbol{x}}_{ij}\cdot {\boldsymbol{v}}_{ij}{{≥}} 0\end{array}\right.$

(8)

式中：α和β为常数， ${\overline c_{ij}},\;{\overline \rho _{ij}}$ 分别为材料声速和密度，x_ij = x_i − x_j为相对位置矢量，v_ij = v_i − v_j为相对速度矢量。

大量试算表明，在NOSB-PD框架中，采用式(8)仍无法很好消除拖尾波；此外考虑到本文中的人工体积黏性主要用于消除应力波的数值振荡，将上述黏性项修改为：

${\varPi _{ij}} = {\overline \rho _{ij}}\left[ {\alpha {{\bar c}_{ij}}{\phi _{ij}} + {\text{sign}}\left( {{\phi _{ij}}} \right)\beta \phi _{ij}^2} \right]$

(9)

式中： ${\text{sign}}$ 为符号函数。考虑到NOSB-PD的非局部特点，式(9)中各参数定义如下：

${\phi _{ij}} = \frac{{{{\overline \delta }_{ij}}{{{\boldsymbol{v}}}_{ij}} \cdot {{{\boldsymbol{x}}}_{ij}}}}{{{{\left| {{{{\boldsymbol{x}}}_{ij}}} \right|}^2} + {{\left( {\varphi {{\overline \delta }_{ij}}} \right)}^2}}},\quad\quad{\overline c_{ij}} = \frac{1}{2}\left( {{c_i} + {c_j}} \right),\quad\quad{\overline \rho _{ij}} = \frac{1}{2}\left( {{\rho _i} + {\rho _j}} \right),\quad\quad{\overline \delta _{ij}} = \frac{1}{2}\left( {{\delta _i} + {\delta _j}} \right)$

(10)

式中： ${\overline \delta _{ij}}$ 为近场范围的均值， $\left| {{{x}_{ij}}} \right|$ 为两物质点之间的距离，φ = 0.01为防止分母变为零的小量。将式(9)变为柯西应力形式为：

${\sigma _{\text{M}}} = {\varPi _{ij}}\kappa$

(11)

式中：κ为二阶单位张量。类比1.1节式(4)中力矢量状态 $\underline{T}\left( {{{{\boldsymbol{x}}}_i},t} \right)\left( {{{\xi }_{ij}}} \right)$ 的确定方法，人工体积黏性的力矢量状态 $\underline {T} _{{\rm{M}}}\left( {{{{\boldsymbol{x}}}_i},t} \right)\left( {{{\xi }_{ij}}} \right)$ 表达式如下：

$\underline T_{{\rm{M}}}\left({\boldsymbol{x}}_{i},t\right)\left({\xi}_{ij}\right)=\omega \left({\xi}_{ij}\right){ \boldsymbol{P}}_{{\boldsymbol{x}}_{i},{\text{M}}}{\boldsymbol{K}}_{{\boldsymbol{x}}_{i}}^{-1}\left({\xi}_{ij}\right),{\boldsymbol{P}}_{{\boldsymbol{x}}_{i},{\text{M}}}=\mathrm{det}\left({\boldsymbol{F}}_{{\boldsymbol{x}}_{i}}\right){\sigma}_{\rm{M}}{\boldsymbol{F}}_{{\boldsymbol{x}}_{i}}^{-\text{T}}$

(12)

因此，将式（12）与式（4）相加即为考虑人工体积黏性的力矢量状态：

$\underline {T} \left( {{{{\boldsymbol{x}}}_i},t} \right)\left( {{{\xi }_{ij}}} \right) = \omega \left( {{{\xi }_{ij}}} \right){{{\boldsymbol{P}}}_{{{{\boldsymbol{x}}}_i}}}{{\boldsymbol{K}}}_{{{{\boldsymbol{x}}}_i}}^{ - 1}\left( {{{\xi }_{ij}}} \right) + \omega \left( {{{\xi }_{ij}}} \right){{{\boldsymbol{P}}}_{{{{\boldsymbol{x}}}_i},\;{\rm{M}}}}{K}_{{{x}_i}}^{ - 1}\left( {{{\xi }_{ij}}} \right)$

(13)

此外，注意到式（13）将经典连续介质力学的一点应力状态与非常规态近场动力学的运动方程建立联系，通过考虑修正的Monaghan人工体积黏性的柯西应力 ${{\sigma }_{\text{M}}}$ 来计算Piola-Kirchhoff应力张量 ${P_{{{x}_i},\;{\text{M}}}}$ ，进而得到修正的Monaghan人工体积黏性的力矢量状态 $\underline {T}_{{\rm{M}}}\left( {{{{\mathrm{x}}}_i},t} \right)\left( {{{\xi }_{ij}}} \right)$ ，并增加在原始力矢量状态 $\underline{T} \left( {{{{\mathrm{x}}}_i},t} \right)\left( {{{\xi }_{ij}}} \right)$ 上，从而抑制强间断弹性波传播的数值振荡。相比之下，已有研究^[16]仅在运动方程上单独附加黏性项，其应力张量无法反映人工黏性。

1.3 考虑等效计算应变率的Kong-Fang模型

为了更好描述混凝土材料在冲击爆炸作用下的力学行为，课题组近年来研发了Kong-Fang材料模型，该模型由Kong等^[4]针对普通混凝土材料提出，近年来逐渐应用于岩石^[34]、超高性能混凝土（ultra-high performance concrete, UHPC）^[35]、活性粉末混凝土（Reactive powder concrete, RPC）^[36]和地聚合物混凝土（geopolymer based high performance concrete, G-HPC）^[37]等材料，得到了广泛应用与认可。本文中采用Kong-Fang模型描述混凝土动态力学行为，模型介绍可参见^[4]，这里不再赘述。

应变率效应的准确描述，尤其是动态拉伸断裂时应变率突变下高应变率效应的准确描述，对准确模拟混凝土材料的动态拉伸断裂破坏至关重要。注意到Kong-Fang模型采用瞬时应变率计算动态强度增强因子（dynamic intensity factor, DIF），并通过DIF径向放大强度面的方法考虑应变率效应，已有研究指出该方法无法描述应变率突变时应力响应的迟滞效应^[5]。为解决该问题，Kong等^[5]提出了等效计算应变率方法，其表达式如下：

${\overline {\dot \varepsilon} _{n + 1}} = \left\{ \begin{array}{ll} \dfrac{{{{\overline {\dot \varepsilon }}_n} + \left( {{{\Delta t}/ \eta }} \right){{\dot \varepsilon }_{n + 1}}}}{{1 + {{\Delta t} / \eta }}}&\quad\quad D {{＞}} 0 \\ \dot \varepsilon &\quad\quad D {{≤}} 0 \end{array} \right.$

(14)

式中： ${\overline {\dot \varepsilon} _n}$ 和 ${\overline {\dot \varepsilon} _{n + 1}}$ 分别为当前和下一时步的等效应变率， $\dot \varepsilon$ 为瞬时应变率，Δt为时间步长，η为时间松弛因子，D为总损伤。本文中将上述等效计算应变率方法植入到Kong-Fang模型中，用以准确模拟动态拉伸断裂时的应变率效应。

1.4 数值实现

以下首先将NOSB-PD中积分型控制方程进行空间离散化处理，然后给出考虑等效计算应变率的Kong-Fang模型应力更新方法，最后通过预测矫正的Newmark方法对离散化的控制方程进行数值求解。

将物体离散为有限物质点，则积分型运动方程（式(3)）变为：

$\rho \ddot {\boldsymbol{u}}\left( {{{\boldsymbol{x}}_i},t} \right) = \sum\limits_{j = 1}^N {\left\{ {\underline T \left( {{{\boldsymbol{x}}_i},t} \right)\left( {{\xi _{ij}}} \right) - \underline T \left( {{{{\boldsymbol{x}}}_j},t} \right)\left( {{\xi _{ij}}} \right)} \right\}} {V_{{{\boldsymbol{x}}_j}}} + {{\boldsymbol{b}}}\left( {{{\boldsymbol{x}}_i},t} \right)$

(15)

式中：N为物质点x_i近场范围内邻近物质点x_j的总数， ${V_{{x_j}}}$ 为物质点x_j的体积。式(5)中的形状张量 ${{\boldsymbol{K}}_{{{\boldsymbol{x}}_i}}}$ 和式(7)中的非局部变形梯度 ${{\boldsymbol{F}}_{{{\boldsymbol{x}}_i}}}$ 的离散形式也可相应得到。

考虑等效计算应变率的Kong-Fang模型用于在给定应变增量Δε下，求解式(6)中的柯西应力σ。在NOSB-PD中，Δε可通过非局部变形梯度 ${{\boldsymbol{F}}_{{{\boldsymbol{x}}_i}}}$ 得到，即：

$\Delta \varepsilon = {\boldsymbol{D}}\Delta t ,\quad\quad\quad {\boldsymbol{D}} = \frac{1}{2}\left( {{\boldsymbol{L}} + {{\boldsymbol{L}}^{\text{T}}}} \right) , \quad\quad\quad{\boldsymbol{L}} = \dot {{\boldsymbol{F}}}_{{{\boldsymbol{x}}_i}} {\boldsymbol{F}}_{{{\boldsymbol{x}}_i}}^{ - 1}$

(16)

式中：Δt为时间步长，D为变形率张量，L为空间速度梯度， $\dot {\boldsymbol{F}}_{{{\boldsymbol{x}}_i}}$ 为非局部变形梯度 ${{\boldsymbol{F}}_{{{\boldsymbol{x}}_i}}}$ 的导数。然后根据Kong-Fang模型的应力更新算法（详见文献[4]）可得到柯西应力σ，进而更新式(4)中的力矢量状态。

求解物质点x_i在t时刻的速度矢量 ${\dot {\boldsymbol{u}}_{{{\boldsymbol{x}}_i},\;t}}$ 和位移增量矢量 $\Delta {{\boldsymbol{u}}_{{{\boldsymbol{x}}_i},\;t}}$ 的Newmark方法为：

${\dot {\boldsymbol{u}}_{{{\boldsymbol{x}}_i},\;t}} = {\dot {\boldsymbol{u}}_{{{\boldsymbol{x}}_i},\;\left( {t - \Delta t} \right)}} + \left[ {\left( {1 - \gamma } \right){{\ddot {\boldsymbol{u}}}_{{{\boldsymbol{x}}_i},\;\left( {t - \Delta t} \right)}} + \gamma {{\ddot {\boldsymbol{u}}}_{{{\boldsymbol{x}}_i},\;t}}} \right]\Delta t$

(17)

$\Delta {{{\boldsymbol{u}}}_{{{\boldsymbol{x}}_i},\;t}} = {\dot {\boldsymbol{u}}_{{{\boldsymbol{x}}_i},\;\left( {t - \Delta t} \right)}}\Delta t + \left[ {\left( {\frac{1}{2} - \psi } \right){{\ddot {\boldsymbol{u}}}_{{{\boldsymbol{x}}_i},\;\left( {t - \Delta t} \right)}} + \psi {{\ddot {\boldsymbol{u}}}_{{{\boldsymbol{x}}_i},\;t}}} \right]\Delta {t^2}$

(18)

式中：γ和ψ为Newmark常数。式(15)和式(17)～(18)可采用预测校正方法求解，得到物质点x_i的位移、速度以及加速度（详见文献[14]）。其中时间步长Δt需要满足Courant-Friedrichs-Levy (CFL) 条件：

$c\frac{{\Delta t}}{{\Delta x}} {{≤}} {C_{\max }}$

(19)

式中：c为声速，C_max≤1.0。

2. 一维杆中的弹性波传播

为验证修正的Monaghan人工体积黏性对于消除应力波传播模拟中数值振荡的可靠性，基于考虑人工体积黏性和等效计算应变率的非常规态近场动力学模型，对一维杆中的弹性波传播进行数值模拟，在开展收敛性分析（δ收敛和m收敛）基础上，验证了提出的人工体积黏性方法的可靠性，进一步讨论了人工体积黏性参数（α和β）的影响并给出参数建议值。

2.1 收敛性分析

如图2所示，考虑一细长一维杆，杆长为L=1 m，为避免横向惯性影响，横截面积A应取充分小量，这里取为1.6×10⁻⁵m²；弹性材料假定为钢，弹性模量为E=210 GPa，密度为ρ=7800 kg/m³，弹性波速为c=5188 m/s。杆左侧施加矩形波，压力幅值为p=100 MPa，作用时间为t=42 μs，通过设置一层虚拟粒子实现，杆右侧为自由边界。

图 2 一维杆中的弹性波传播数值模型示意图

Figure 2. Numerical model for wave propagation in one-dimensional rod

下载: 全尺寸图片幻灯片

PD作为非局部理论，其中粒子间距Δx和邻域半径δ均会对模拟结果产生影响，为准确预测一维杆中的弹性波传播首先进行收敛性分析^[38]，即δ收敛性和m（= δ/Δx）收敛性。其中δ收敛性通过固定m=3，改变δ（0.001 2、0.003 和0.006 m）取值研究其收敛性；m收敛性通过固定δ=0.001 2 m，改变m（1.5、3.0和6.0）取值研究其收敛性。此外，采用均匀离散方式，时间步长根据最小粒子间距下（Δx=0.000 2 m）的稳定时间步长确定。

图3给出了不考虑人工体积黏性（α = 0、β =0）时一维杆中0.25 m处测点的应力时程曲线，可以看出，模拟结果的δ收敛性和m收敛性均较好，而对于波峰处的δ收敛性存在差异，主要由于不同粒子间距Δx所导致，当Δx过大造成升压时间相对滞后。此外，由于未施加人工体积黏性，弹性波传播过程中发生明显的数值振荡，尤其是相比理论解出现了拉伸拖尾波，该非物理的拖尾波可能会导致材料出现非物理的动态拉伸断裂。

图 3 不考虑人工体积黏性时一维杆中0.25 m处测点的应力时程曲线

Figure 3. Stress-time history of gauge located at 0.25 m without consideration of artificial bulk viscosity

下载: 全尺寸图片幻灯片

图4给出了考虑人工体积黏性（α = 0.1、β =3.0）时一维杆中0.25 m处测点的应力时程曲线，同样可以看出，模拟结果的δ收敛性和m收敛性均较好。此外，由于引入了修正的Monaghan人工体积黏性，数值振荡尤其是拉伸拖尾波得到有效抑制，且随着邻域半径δ的减小，模拟结果逐渐趋近于理论解。

图 4 考虑人工体积黏性时一维杆中0.25 m处测点的应力时程曲线

Figure 4. Stress-time history of gauge located at 0.25 m with consideration of artificial bulk viscosity

下载: 全尺寸图片幻灯片

2.2 黏性参数的影响

上述分析表明，未引入修正的Monaghan人工体积黏性时，弹性波传播出现明显的数值振荡和拉伸拖尾波，然而当引入修正的Monaghan人工体积黏性时，有效地抑制了数值振荡和拉伸拖尾波。考虑到人工体积黏性参数α和β会对上述抑制效果有明显影响，因此需对两参数进行敏感性分析，并给出合适的建议值。考虑到计算成本和计算精度的影响，基于上述收敛性分析，进行两参数敏感性分析时，粒子间距Δx取为0.0004 m、邻域半径δ取为0.001 2 m，则m为3。

为了满足黏性参数的相互独立条件，通过固定β=3.0改变α取值研究其影响；固定α=0.1改变β取值研究其影响。图5~7分别给出了不同α和β取值下数值模拟预测的0.25和0.50 m处测点的应力时程曲线。可以看出，当α=0.03, 0.06时，最大峰值应力处数值振荡抑制效果较差，同时产生幅值较大的拖尾波（约13 MPa）；当α=0.2, 0.3时，虽然拖尾波得到较好抑制，但随着弹性波传播的数值解逐渐偏离理论解，形成升压时间增加和持续时间减小现象；而当α=0.1时，不但最大峰值应力的数值振荡和拉伸拖尾波得到有效抑制，且矩形波持续时间与理论解吻合较好，同时随着波的传播拉伸拖尾波现象逐渐削弱；此外，α对波传播影响较大而β影响较小。综合上述分析，修正的Monaghan人工体积黏性参数建议值为α=0.1、β=3.0。

图 5 不同黏性参数下0.25 m处测点的应力时程曲线

Figure 5. Stress-time history of gauge located at 0.25 m with different viscosity parameters

下载: 全尺寸图片幻灯片

图 6 不同黏性参数下0.50 m处测点的应力时程曲线

Figure 6. Stress-time history of gauge located at 0.50 m with different viscosity parameters

下载: 全尺寸图片幻灯片

图 7 不同黏性参数下不同测点的应力时程曲线

Figure 7. Stress-time history of different gauges with different viscosity parameters

下载: 全尺寸图片幻灯片

为了验证修正的Monaghan人工体积黏性的优越性，图8给出了未修正和修正的Monaghan人工体积黏性数值模拟预测的0.25和0.50 m处测点的应力时程曲线。可以看出相同黏性参数值（α=0.1、β=3.0）时，修正的Monaghan人工体积黏性更好地抑制了数值振荡和拉伸拖尾波，同时更加趋于理论解。

图 8 不同黏性方式下不同测点的应力时程曲线

Figure 8. Stress-time history of different gauges with different viscosity modes

下载: 全尺寸图片幻灯片

3. 混凝土层裂试验验证

上述分析表明了引入修正的Monaghan人工体积黏性的必要性并给出参数建议取值，基于上述参数建议值，并采用考虑等效计算应变率的Kong-Fang模型，对Schuler等^[39]开展的混凝土层裂试验进行数值模拟，首先进行δ收敛性和m收敛性分析，然后对比分析了人工体积黏性和不同应变率计算方法对预测混凝土动态拉伸断裂破坏的影响。

3.1 收敛性分析

Schuler等^[39]开展了混凝土试件层裂试验，试验中采用直径为74.2 mm、长度为60 mm的子弹撞击直径为74.2 mm、长度为5 000 mm的入射杆产生压缩波，然后作用于直径为74.2 mm、长度为250 mm、抗压强度为35 MPa的混凝土试件上。压缩波传播至混凝土自由端反射形成拉伸波使试件发生动态拉伸断裂破坏，且随着子弹撞击速度（4.1、7.6和11.1 m/s）的增加，混凝土的动态拉伸断裂从单条裂缝变成复杂的多条裂缝。

建立相应的数值模型如图9所示，将数值模型均匀离散为包含所有物性信息的物质点，考虑到物质点间距Δx和邻域半径δ均会对混凝土动态拉伸断裂的模拟结果产生一定影响，因此以撞击速度4.1 m/s为例，首先进行了δ收敛性和m（m=δ/Δx）收敛性分析，其中δ的收敛性通过固定m=3.0，改变δ值（0.009、0.012和0.015 m）来研究；m的收敛性通过固定δ=0.012 m，改变m值（2.0、3.0和4.0）来研究。

图 9 混凝土层裂数值模型示意图

Figure 9. Numerical model for concrete spall test

下载: 全尺寸图片幻灯片

子弹和入射杆均采用弹性材料描述，其中弹性模量和泊松比分别为200 GPa和0.3；混凝土材料采用1.3节考虑等效计算应变率的Kong-Fang模型描述，模型参数根据抗压强度35 MPa自动生成^[4]，前期研究表明当时间松弛因子（式(14)中的η）取为0.1 ms时能够较好地预测混凝土材料在应变率变化时的迟滞效应^[5]，主要参数见表1。此外，子弹与入射杆、入射杆与混凝土试件之间的相互作用通过Liu等^[14]提出的点对体接触算法实现。

表 1 层裂试验PD模拟参数

Table 1. Parameters required in PD simulations for spall tests

参数	取值	参数	取值
粒子间距Δx	4 × 10⁻³ m	密度ρ	2 400 kg/m³
邻域半径δ	3Δx	强度面参数a₁	0.5876
时间步Δt	4 × 10⁻⁷ s	强度面参数a₂	0.025/f_c
Newmark常数γ	0.6	损伤参数d₁	0.04
Newmark常数ψ	0.6	损伤参数d₂	1.5
抗压强度f_c	35 MPa	流动法则参数ω	0.5
抗拉强度f_t	3.22 MPa	断裂应变ε_frac	0.01
弹性模量E	28 GPa	时间迟滞因子η	1 × 10⁻⁴ s
泊松比v	0.3	状态方程	文献[4]

下载: 导出CSV

| 显示表格

图10和图11分别给出了子弹冲击速度4.1 m/s工况下，不同δ和m取值预测的拉伸断裂情况，可以看出不同参数取值下虽然裂缝带宽度预测不同（主要由于粒子间距Δx的不同导致），但裂缝位置基本一致，表现出较好的收敛性。为平衡计算效率和计算精度，后续模拟中选取Δx=0.004 m, δ=3Δx = 0.012 m。

图 10 子弹冲击速度4.1 m/s工况下的δ收敛性

Figure 10. Convergence of δ under the condition of projectile impact velocity with 4.1 m/s

下载: 全尺寸图片幻灯片

图 11 子弹冲击速度4.1 m/s工况下的m收敛性

Figure 11. Convergence of m under the condition of projectile impact velocity with 4.1 m/s

下载: 全尺寸图片幻灯片

3.2 对比分析

基于上述收敛性分析确定的Δx和δ以及表1给出的主要材料参数，其中选取合适的粒子间距Δx=0.004 m，邻域半径δ=3Δx，由Δx确定稳定时间步长Δt=4 × 10⁻⁷s，预测矫正Newmark方法中常数γ和ψ均为0.6。此外，混凝土材料属性参数与文献[39]一致，模型参数由抗压强度35 MPa自动生成^[4]。本节开展了子弹不同冲击速度下混凝土层裂的数值模拟，并对比分析了是否施加人工黏性和不同应变率计算方法对混凝土动态拉伸断裂破坏预测结果的影响。施加人工体积黏性时α=0.1、β=3.0；不施加时α=0、β=0。考虑等效计算应变率时，η=0.1 ms；考虑常用的瞬时应变率时，η趋近于0（见式14），数值计算时可取一非零小量代替。

图12分别给出了试验后试件的断裂情况（图12(a)），以及不考虑黏性且考虑等效计算应变率（图12(b)）、考虑黏性且考虑等效计算应变率（图12(c)）、考虑黏性且考虑常用的瞬时应变率（图12(d)）预测的试件断裂及破坏情况。可以看出，随着冲击速度的增加，混凝土试件由单条裂缝变成多条裂缝破坏现象；当不考虑黏性且考虑等效计算应变率时，模拟结果虽然再现了混凝土裂缝的产生和数量，但裂缝位置与试验不符，且试件左侧出现非物理的压缩破坏，这是由于未考虑人工体积黏性使得预测的应力波出现严重的数值振荡且峰值过大导致的；当考虑黏性且考虑常用的瞬时应变率时，模拟结果中裂缝位置和数量均与试验不符，这是由于常用的瞬时应变率计算方法忽略了动态拉伸断裂时应变率突变引起的应力迟滞效应导致的；当同时考虑黏性和等效计算应变率计算方法时，数值模拟预测的裂缝位置和数量均与试验吻合较好。

图 12 不同子弹冲击速度下不同计算方法预测的混凝土试件断裂破坏及与试验的对比

Figure 12. Experimental and numerical fracture and failure in concrete specimens with different projectile impact velocity

下载: 全尺寸图片幻灯片

图13进一步给出了子弹不同冲击速度下数值模拟预测的裂缝位置和试验结果^[39]的对比曲线，需要指出的是文献[39]中未提供11.1 m/s冲击下的裂缝位置信息，可以看出数值模拟预测的裂缝位置与试验吻合很好。

图 13 不同冲击速度下试验与PD模拟的裂缝位置对比

Figure 13. Comparisons of crack locations between test and PD simulations under different impact velocities

下载: 全尺寸图片幻灯片

经上述对比分析可知，基于NOSB-PD预测混凝土材料动态拉伸断裂破坏时存在两个关键问题，即数值振荡和应变率效应，通过不同速度下混凝土层裂的模拟验证了修正的人工体积黏性抑制数值振荡和采用等效计算应变率方法准确计算应变率效应的可靠性。

4. 结　论

为了准确预测冲击爆炸荷载作用下混凝土材料的动态拉伸断裂破坏，首先建立了修正的Monaghan人工体积黏性计算方法用于消除数值振荡，然后采用前期建立的等效计算应变率方法用以准确计算应变率突变时的应变率效应，进而将两种方法植入到NOSB-PD框架中，开展了一维杆中的弹性波传播的数值模拟、分析了人工体积黏性参数的影响并给出参数建议值，最后将模型用于混凝土试件层裂的数值模拟，对比分析了人工体积黏性、不同应变率效应计算方法对动态拉伸断裂预测结果的影响规律。主要结论如下：

(1)在力矢量状态上额外附加修正的Monaghan人工体积黏性力矢量状态，可以有效地抑制NOSB-PD中由变形梯度的近似导致的非物理数值振荡现象。

(2)通过讨论修正的Monaghan人工体积黏性参数（α和β）对一维杆中的弹性波传播影响规律，给出适合研究强间断波传播问题的人工体积黏性参数建议值：α = 0.1、β =3.0。

(3)通过对比分析人工体积黏性和不同应变率计算方法对混凝土层裂预测结果的影响，表明预测动态拉伸断裂破坏，需同时考虑人工体积黏性和等效计算应变率，建立的考虑人工体积黏性和等效计算应变率的非常规态近场动力学模型可以很好预测裂缝位置和数量，后期研究会将模型进一步用于弹体侵彻和装药爆炸下混凝土材料动态破坏的数值模拟。

图 1 NOSB-PD理论中质点的运动

Figure 1. Kinematics of material points in NOSB-PD theory

下载: 全尺寸图片幻灯片

图 2 一维杆中的弹性波传播数值模型示意图

Figure 2. Numerical model for wave propagation in one-dimensional rod

下载: 全尺寸图片幻灯片

图 3 不考虑人工体积黏性时一维杆中0.25 m处测点的应力时程曲线

Figure 3. Stress-time history of gauge located at 0.25 m without consideration of artificial bulk viscosity

下载: 全尺寸图片幻灯片

图 4 考虑人工体积黏性时一维杆中0.25 m处测点的应力时程曲线

Figure 4. Stress-time history of gauge located at 0.25 m with consideration of artificial bulk viscosity

下载: 全尺寸图片幻灯片

图 5 不同黏性参数下0.25 m处测点的应力时程曲线

Figure 5. Stress-time history of gauge located at 0.25 m with different viscosity parameters

下载: 全尺寸图片幻灯片

图 6 不同黏性参数下0.50 m处测点的应力时程曲线

Figure 6. Stress-time history of gauge located at 0.50 m with different viscosity parameters

下载: 全尺寸图片幻灯片

图 7 不同黏性参数下不同测点的应力时程曲线

Figure 7. Stress-time history of different gauges with different viscosity parameters

下载: 全尺寸图片幻灯片

图 8 不同黏性方式下不同测点的应力时程曲线

Figure 8. Stress-time history of different gauges with different viscosity modes

下载: 全尺寸图片幻灯片

图 9 混凝土层裂数值模型示意图

Figure 9. Numerical model for concrete spall test

下载: 全尺寸图片幻灯片

图 10 子弹冲击速度4.1 m/s工况下的δ收敛性

Figure 10. Convergence of δ under the condition of projectile impact velocity with 4.1 m/s

下载: 全尺寸图片幻灯片

图 11 子弹冲击速度4.1 m/s工况下的m收敛性

Figure 11. Convergence of m under the condition of projectile impact velocity with 4.1 m/s

下载: 全尺寸图片幻灯片

图 12 不同子弹冲击速度下不同计算方法预测的混凝土试件断裂破坏及与试验的对比

Figure 12. Experimental and numerical fracture and failure in concrete specimens with different projectile impact velocity

下载: 全尺寸图片幻灯片

图 13 不同冲击速度下试验与PD模拟的裂缝位置对比

Figure 13. Comparisons of crack locations between test and PD simulations under different impact velocities

下载: 全尺寸图片幻灯片

表 1 层裂试验PD模拟参数

Table 1. Parameters required in PD simulations for spall tests

参数	取值	参数	取值
粒子间距Δx	4 × 10⁻³ m	密度ρ	2 400 kg/m³
邻域半径δ	3Δx	强度面参数a₁	0.5876
时间步Δt	4 × 10⁻⁷ s	强度面参数a₂	0.025/f_c
Newmark常数γ	0.6	损伤参数d₁	0.04
Newmark常数ψ	0.6	损伤参数d₂	1.5
抗压强度f_c	35 MPa	流动法则参数ω	0.5
抗拉强度f_t	3.22 MPa	断裂应变ε_frac	0.01
弹性模量E	28 GPa	时间迟滞因子η	1 × 10⁻⁴ s
泊松比v	0.3	状态方程	文献[4]

下载: 导出CSV

参考文献(39)

[1]	方秦, 陈小伟. 冲击爆炸效应与工程防护专辑·编者按 [J]. 中国科学: 物理学力学天文学, 2020, 50(2): 024601. DOI: 10.1360/SSPMA-2019-0404. FANG Q, CHEN X W. Special issue on impact and explosion effect and engineering protection [J]. Scientia Sinica Physica, Mechanica & Astronomica, 2020, 50(2): 024601. DOI: 10.1360/SSPMA-2019-0404.
[2]	KONG X Z, FANG Q, WU H, et al. Numerical predictions of cratering and scabbing in concrete slabs subjected to projectile impact using a modified version of HJC material model [J]. International Journal of Impact Engineering, 2016, 95: 61–71. DOI: 10.1016/j.ijimpeng.2016.04.014.
[3]	KONG X Z, FANG Q, LI Q M, et al. Modified K&C model for cratering and scabbing of concrete slabs under projectile impact [J]. International Journal of Impact Engineering, 2017, 108: 217–228. DOI: 10.1016/j.ijimpeng.2017.02.016.
[4]	KONG X Z, FANG Q, CHEN L, et al. A new material model for concrete subjected to intense dynamic loadings [J]. International Journal of Impact Engineering, 2018, 120: 60–78. DOI: 10.1016/j.ijimpeng.2018.05.006.
[5]	KONG X Z, FANG Q, ZHANG J H, et al. Numerical prediction of dynamic tensile failure in concrete by a corrected strain-rate dependent nonlocal material model [J]. International Journal of Impact Engineering, 2020, 137: 103445. DOI: 10.1016/j. ijimpeng.2019.103445. DOI: 10.1016/j.ijimpeng.2019.103445.
[6]	KHOE Y S, WEERHEIJM J. Limitations of smeared crack models for dynamic analysis of concrete [C]//The 12th International LS-DYNA Users Conference. Constitutive Models, 2012.
[7]	KONG X Z, FANG Q, HONG J. A new damage-based nonlocal model for dynamic tensile failure of concrete material [J]. International Journal of Impact Engineering, 2019, 132: 103336. DOI: 10.1016/j.ijimpeng.2019.103336.
[8]	XU K, LU Y. Numerical simulation study of spallation in reinforced concrete plates subjected to blast loading [J]. Computers & Structures, 2006, 84(5/6): 431–438. DOI: 10.1016/j.compstruc.2005.09.029.
[9]	LI J, HAO H. Numerical study of concrete spall damage to blast loads [J]. International Journal of Impact Engineering, 2014, 68: 41–55. DOI: 10.1016/j.ijimpeng.2014.02.001.
[10]	孔祥振, 方秦. 基于SPH方法对强动载下混凝土结构损伤破坏的数值模拟研究 [J]. 中国科学: 物理学力学天文学, 2020, 50(2): 25–33. DOI: 10.1360/SSPMA-2019-0186. KONG X Z, FANG Q. Numerical predictions of failures in concrete structures subjected to intense dynamic loadings using the smooth particle hydrodynamics method [J]. Scientia Sinica Physica, Mechanica & Astronomica, 2020, 50(2): 25–33. DOI: 10.1360/SSPMA-2019-0186.
[11]	YANG Y Z, FANG Q, KONG X Z. Failure mode and stress wave propagation in concrete target subjected to a projectile penetration followed by charge explosion: experimental and numerical investigation [J]. International Journal of Impact Engineering, 2023, 177: 104595. DOI: 10.1016/j.ijimpeng.2023.104595.
[12]	SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces [J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209. DOI: 10.1016/S0022-5096(99)00029-0.
[13]	SILLING S A, EPTON M, WECKNER O, et al. Peridynamic states and constitutive modeling [J]. Journal of Elasticity, 2007, 88(2): 151–184. DOI: 10.1007/s10659-007-9125-1.
[14]	LIU X, KONG X Z, FANG Q, et al. Peridynamics modelling of projectile penetration into concrete targets [J]. International Journal of Impact Engineering, 2025, 195: 105110. DOI: 10.1016/j.ijimpeng.2024.105110.
[15]	LITTLEWOOD D J. Simulation of dynamic fracture using peridynamics, finite element modeling, and contact [C]//ASME International Mechanical Engineering Congress and Exposition. American Society of Mechanical Engineers, 2010, 44465: 209–217.
[16]	BREITENFELD M S, GEUBELLE P H, WECKNER O, et al. Non-ordinary state-based peridynamic analysis of stationary crack problems [J]. Computer Methods in Applied Mechanics and Engineering, 2014, 272: 233–250. DOI: 10.1016/j.cma.2014.01.002.
[17]	GU X, ZHANG Q, YU Y T. An effective way to control numerical instability of a nonordinary state-based peridynamic elastic model [J]. Mathematical Problems in Engineering Theory Methods and Applications, 2017(1): 1750876. DOI: 10.1155/2017/1750876.
[18]	YAGHOOBI A, CHORZEPA M G. Higher-order approximation to suppress the zero-energy mode in non-ordinary state-based peridynamics [J]. Computers & Structures, 2017, 188: 63–79. DOI: 10.1016/j.compstruc.2017.03.019.
[19]	CHEN H L. Bond-associated deformation gradients for peridynamic correspondence model [J]. Mechanics Research Communications, 2018, 90: 34–41. DOI: 10.1016/j.mechrescom.2018.04.004.
[20]	GU X, ZHANG Q, MADENCI E, et al. Possible causes of numerical oscillations in non-ordinary state-based peridynamics and a bond-associated higher-order stabilized model [J]. Computer Methods in Applied Mechanics and Engineering, 2019, 357: 112592. DOI: 10.1016/j.cma.2019.112592.
[21]	ROSSI P. A physical phenomenon which can explain the mechanical behaviour of concrete under high strain rates [J]. Materials and Structures, 1991, 24: 422–424. DOI: 10.1007/BF02472015.
[22]	REINHARDT H W, WEERHEIJM J. Tensile fracture of concrete at high loading rates taking account of inertia and crack velocity effects [J]. International Journal of Fracture, 1991, 51(1): 31–42. DOI: 10.1007/BF00020851.
[23]	CURBACH M, EIBL J. Crack velocity in concrete [J]. Engineering Fracture Mechanics, 1990, 35(1/3): 321–326. DOI: 10.1016/0013-7944(90)90210-8.
[24]	LI Q M, MENG H. About the dynamic strength enhancement of concrete-like materials in a split Hopkinson pressure bar test [J]. International Journal of Solids and Structures, 2003, 40(2): 343–360. DOI: 10.1016/S0020-7683(02)00526-7.
[25]	ZHANG M, WU H J, LI Q M, et al. Further investigation on the dynamic compressive strength enhancement of concrete-like materials based on split Hopkinson pressure bar tests: part I: experiments [J]. International Journal of Impact Engineering, 2009, 36(12): 1327–1334. DOI: 10.1016/j.ijimpeng.2009.04.009.
[26]	LI Q M, LU Y B, MENG H. Further investigation on the dynamic compressive strength enhancement of concrete-like materials based on split Hopkinson pressure bar tests: part II: numerical simulations [J]. International Journal of Impact Engineering, 2009, 36(12): 1335–1345. DOI: 10.1016/j.ijimpeng.2009.04.010.
[27]	LU Y B. , LI Q M. About the dynamic uniaxial tensile strength of concrete-like materials [J]. International Journal of Impact Engineering, 2011, 38(4): 171–180. DOI: 10.1016/j.ijimpeng.2010.10.028.
[28]	SCHMINDT M J. High pressure and high strain rate behavior of cementitious materials: experiments and elastic/viscoplastic modeling [M]. University of Florida, 2003.
[29]	TANDON S, FABER K T, BAZANT Z P, et al.Cohesive crack modeling of influence of sudden changes in loading rate on concrete fracture [J]. Engineering Fracture Mechanics, 1995, 52(6): 987–997. DOI: 10.1016/0013-7944(95)00080-F.
[30]	HOLMQUIST T J, JOHNSON G R. A computational constitutive model for concrete subjected to larger strains, high strain rates and high pressure [C]//The 14th International Symposium Ballistics. USA: American Defense Preparedness Association, 1995, 591–600.
[31]	MALVAR L J, CRAWFORD J E, WESEVICH J W, et al. A plasticity concrete material model for DYNA3D [J]. International Journal of Impact Engineering, 1997, 19: 847–873. DOI: 10.1016/S0734-743X(97)00023-7.
[32]	RIEDEL W, THOMA K, HIERMAIER S, et al. Penetration of reinforced concrete by BETA-B-500 numerical analysis using a new macroscopic concrete model for hydrocodes [C]//Proceedings of the 9th International Symposium on the Effects of Munitions with Structures. Berlin-Strausberg Germany, 1999; 315.
[33]	MONAGHAN J J, GINGOID R A. Shock simulation by the particle method SPH [J]. Journal of Computational Physics, 1983, 52(2): 374–389. DOI: 10.1016/0021-9991(83)90036-0.
[34]	HUANG X P, KONG X Z, CHEN Z Y, et al. A computational constitutive model for rock in hydrocode [J]. International Journal of Impact Engineering, 2020, 145: 103687. DOI: 10.1016/j.ijimpeng.2020.103687.
[35]	YANG S B, KONG X Z, WU H, et al. Constitutive modelling of UHPCC material under impact and blast loadings [J]. International Journal of Impact Engineering, 2021, 153: 103860. DOI: 10.1016/j.ijimpeng.2021.103860.
[36]	XU S, WU P, LI Q, et al. Experimental investigation and numerical simulation on the blast resistance of reactive powder concrete subjected to blast by embedded explosive [J]. Cement and Concrete Composites, 2021, 119: 103989. DOI: 10.1016/j.cemconcomp.2021.103989.
[37]	YUAN P C, XU S C, LIU J, et al. Experimental and numerical study of blast resistance of geopolymer based high performance concrete sandwich walls incorporated with metallic tube core [J]. Engineering Structures, 2023, 278: 115505. DOI: 10.1016/j.engstruct.2022.115505.
[38]	HA Y D, BOBARU F. Studies of dynamic crack propagation and crack branching with peridynamics [J]. International Journal of Fracture, 2010, 162(1/2): 229–244. DOI: 10.1007/s10704-010-9442-4.
[39]	SCHULER H, MAYRHOFER C, THOMA K. Spall experiments for the measurement of the tensile strength and fracture energy of concrete at high strain rates [J]. International Journal of Impact Engineering, 2006, 32(10): 1635–1650. DOI: 10.1016/j.ijimpeng.2005.01.010.