基于混凝土黏弹性变形机制的近场动力学方法

路德春，宋志强，王国盛

（北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124）

在实际工程中，混凝土结构通常是带裂缝工作的，并且破坏过程具有强非连续特性。传统连续介质力学方法为数值计算提供了丰富的理论基础，基于连续介质力学的部分数值方法为混凝土非连续变形的模拟提供了方法和可能性，例如扩展有限元方法［1］、边界元法［2］、有限差分法［3］、无网格法［4-5］等，成功解决了混凝土结构的一系列工程问题。然而，上述方法的控制方程仍然采用微分方程，在处理混凝土裂缝开展带来的非连续问题时，仍然会导致对裂纹扩展的不合理预测。Silling［6］提出的近场动力学方法将物体离散为具有体积和质量的材料点。所有材料点自身的位移共同决定了物体的变形，材料点之间的相互作用是由材料点的相对位移产生，这种相互作用（也被称为键）的极限状态决定了材料的强度。描述材料点运动的控制方程采用空间积分格式，避免了在裂纹位移不连续处计算空间导数的无效性，从而为模拟不连续变形提供了新的思路和手段。目前，近场动力学已经成为研究材料裂纹形成、发展和分叉等非连续问题的重要方法［7-14］。混凝土在加载过程中会出现开裂现象，其开裂模式、变形和强度特性在很大程度上取决于加载速率及其梯度［15-19］，这也表现为黏性效应、热活化效应和惯性效应［20］。近场动力学处理非连续变形问题因其优势，已经被广泛用于模拟混凝土的动态开裂过程［21-26］。然而，原始近场动力学［6］材料点的运动方程仅包含与惯性效应对应的加速度项，缺少与黏性效应对应的速度项。从而带来了对动态荷载作用下混凝土变形预测不合理的问题。Silling［27-28］发现，在原始近场动力学运动方程中引入与阻尼相关的速度项可以解决波传播的不稳定性问题，相关的工作在最近的研究中取得了长足进展［14，29-30］。由于阻尼力是一种非保守力，因此使用仅包含保守力所做功的虚功原理推导近场动力学运动方程的方案不再适用。要合理模拟材料的变形，必须基于能量守恒原理引入非保守力做功，发展一种新的严格遵循热力学原理的近场动力学方法。另一方面，原始近场动力学的断裂准则也无法合理反映混凝土强度的速率相关特性［22，31-32］。混凝土的强度本质上是材料本构层面的内生特性，应通过强度或断裂准则来表征［19，33-36］。因此，在键中引入率相关断裂准则也是十分必要的。

本文目的是发展包含非保守力做功的近场动力学方法。基于混凝土的变形机理，确定材料点的黏弹性相互作用模式，发展黏弹性近场动力学材料点的运动方程和键的断裂准则。对于材料点的运动，通过包含非保守力所做功的能量方程获得控制方程。对于键的破坏，利用混凝土动态单轴S 强度准则发展的率相关断裂准则来描述。并通过数值算例验证和应用所提出的近场动力学方法。

1 黏弹性近场动力学的理论框架

近场动力学的基本思想是将连续材料离散为具有体积V和质量m的有限个材料点，任意材料点i与一定范围H（i）内的其他材料点都具有相互作用。材料点i的域内所有材料点相互作用的集合控制着点i的运动，该运动产生的惯性力与总的相互作用力保持平衡，其运动方程为［6］

式中：k为微弹性模量；s为键拉伸，其等于键变形后长度与初始长度的差与初始长度的比值；μ为历史函数。

现有研究表明［20］，混凝土的动态承载力包括惯性承载力和真实动态强度。真实动态强度受热活化效应和黏性效应控制，这2 种效应可以分别用弹簧和牛顿黏壶模型来表征。受混凝土动态物理机制的启发，建立了键型近场动力学材料点黏弹性相互作用模型，如图1 所示。模型中的阻尼力是一种非保守力，用于实现系统能量耗散的功能，所发展的模型反映了由非保守力所驱动的能量耗散。

图1 近场动力学黏弹性相互作用模型Fig.1 Viscoelastic interaction model of peridynamics

1.1 黏弹性运动方程

基于发展的黏弹性相互作用模型，通过包括非保守力做功的哈密尔顿原理［39-40］推导近场动力学黏弹性运动方程，其表达式为

式中：δ为变分符号；t为时间；Wnc为非保守力所做的功；T为总动能；U为所有保守力产生的总势能。T、U表达式如式（4）：

新建立的运动方程既包括由弹簧力和体力组成的保守力项，也包括由阻尼力组成的非保守力项。当非保守力项等于零时，黏弹性运动方程式退化为原始近场动力学运动方程，它是发展的运动方程一个特例。

1.2 混凝土率相关断裂准则

在黏弹性近场动力学中，同样采用应变率相关的最大拉应变准则作为材料点的断裂准则，结合断裂准则的键力密度矢量可表示为

式中：FS为弹簧的弹力；FD为黏壶的阻尼力；μ被重新定义为

f(s0，ṡ)=s0d是断裂准则，s0d为键拉伸率为ṡ条件下的动态极限键拉伸。

对于黏弹性近场动力学方法的率相关断裂准则，通过定义动态极限键拉伸来描述极限键拉伸随加载速率的增加，键拉伸的动态增长因子定义为：DIFt=s0d/s0。对于强度的DIFs，目前已经提出了许多率相关公式［20，33，41-43］。在这些公式中，Lu等［20］提出的动态单轴S准则能够合理描述混凝土的真实动态强度。通过类比S 准则给出了DIFt的率相关公式，如式（12）：

式中：ϖ为对数键拉伸率，ϖ=lgṡ；λ和χmax是2 个参数，其确定方法与Lu等［20］提出的方法相同。

2 材料参数的标定方法

发展的黏弹性近场动力学中包含2 类材料参数：微弹性参数和微黏性参数，通过将黏弹性近场动力学和连续介质力学之间的能量密度进行等效，即可获得近场动力学的材料参数。

在二维条件下，连续介质力学的黏弹性本构关系被用来计算耗散能密度，其方程如式（13）：

任何变形模式都可以分解为2种基本的变形模式：体积变形和纯剪切变形。对于体积变形模式，只有正应变ε在板中产生，可得到体积变形模式下微弹性参数κ与微黏性参数τ。

同理，对于纯剪切变形模式，当物体变形时，只有工程剪应变γ在板中产生，可以得到纯剪切变形模式下的微弹性参数κ与微黏性参数τ。

由于2个材料点之间的弹簧力和阻尼力始终沿其轴向，是一个一维分量，因此，在黏弹性近场动力学中应该只有一个微弹性参数和一个微黏性参数。联立式（19）和式（20），得到K=2G、η1=2η2。根据K和G以及η1和η2之间的转换关系，可以得出弹性与黏性泊松比在二维情况下始终都为常数1/3。因此，微弹性参数和微黏性参数可以统一表示为

3 数值算例

采用FORTRAN 90 语言编写了黏弹性近场动力学程序，给出2个算例：动荷载作用下裂纹的分叉模拟和含圆孔各向同性平板的拉伸破坏。

3.1 拉伸动荷载作用下裂纹的分叉模拟

图2为带初始裂缝的混凝土平板在不同冲击应力条件下的开裂模式，对应的冲击应力大小分别为σ1= 0.2 MPa，σ2= 0.5 MPa，σ3= 0.6 MPa，σ4=1.4 MPa［44］。与试验相同工况的数值模型如图3 所示，离散材料点的直径d=0.25 mm，每个材料点的域半径选择为Δ=3d=0.75mm。混凝土的材料参数为杨氏模量E=19.2GPa，黏性参数η=8.28 x 10-22N·s·m-2，ρ= 2 700 kg·m-3，S 准则的参数λ=1.4，χmax=4。边界的体积校正和表面效应的处理与Silling等［37］和Madenci等［38］的处理方法相同。荷载通过应力边界条件施加在板的两端边界区域上，荷载条件与试验相同。

图2 试验中不同应力状态下板的破坏模式Fig.2 Failure modes of plates in different stress states in experiment

图3 模型几何尺寸及荷载条件Fig.3 Geometric dimensions and load conditions of the model

在σ2=0.5MPa 的应力状态下，黏弹性近场动力学模拟的裂纹扩展过程材料点的损伤如图4b 所示，在加载的初期，裂纹从初始裂纹的尖端起裂并且新裂纹沿着初始裂纹方向扩展。当扩展达到一定长度之后，裂纹出现了分叉现象，由1 条裂纹变成了2条裂纹，2 条分叉裂纹分别沿着与初始裂纹呈大约30°的方向向右边界扩展，最终达到右边界，整个试件破坏。模拟的结果与其他学者模拟结果［45］吻合较好，进而验证了发展的黏弹性近场动力学方法的有效性。在不同应力条件下，裂纹的扩展过程有很大区别。当σ1=0.2MPa 时，混凝土板的新裂纹完全沿着初始裂纹的方向发展，近似一条直线并且没有分叉，这种情况下的破坏模式与准静态破坏相同。当应力增加σ2=0.5MPa 时，出现了一次分叉，有2条主裂纹。当应力增加σ3= 0.6MPa 时，出现了2次分叉，有3条主裂纹。其中一次分叉形成的2条裂纹中有1 条产生二次分叉。当应力增加σ4=1.4MPa时，出现了3次分叉，产生4条主裂纹以及更多的小裂纹。混凝土是一种率敏感的准脆性材料，动态荷载会影响其裂纹的扩展模式，并且随着加载速率的增加，裂纹的分叉次数会逐渐增多，这是因为加载速率越大，平衡输入的能量所需要的裂纹区域更多。黏弹性近场动力学方法能很好地捕获混凝土的率相关行为，并且对混凝土开裂等非连续问题能进行合理描述。

图4 不同应力状态下板的破坏模式Fig.4 Failure modes of plates in different stress states

3.2 含圆孔平板的动荷载拉伸

数值试验研究了在板中心具有圆孔的二维混凝土板在动态拉伸荷载下的模拟，圆孔的直径是10mm，施加在板上的几何参数和荷载条件如图5所示。黏弹性近场动力学材料参数与3.1节算例的材料参数保持一致，材料点的边长和域半径也与3.1节的算例相同。荷载通过位移边界条件施加在板的两端边界区域上，3 种加载速率分别为：v1= 5 m·s-1，v2=20 m·s-1，v3=50 m·s-1。在图5中边界的体积校正和表面效应的处理与Silling 等［37］和Madenci等［38］的处理方法相同。

图5 模型几何尺寸及荷载条件Fig.5 Geometric dimensions and load conditions of the model

原始近场动力学和黏弹性近场动力学模拟的混凝土板在不同加载速率下的荷载-位移曲线与损伤破坏模式分别如图6 和图7 所示。对于黏弹性近场动力学，荷载-位移曲线的刚度受加载速率的影响十分明显，随着加载速率的增加而增加，并且极限承载能力也随着加载速率而提高。不同加载速率下损伤范围的区别也十分显著，随着加载速率增加，损伤破坏区域由一条非常细的直线逐渐变宽直至变成V形。这种现象的产生是因为混凝土中涉及变形的材料点的数量随着加载速率的增加而增加，更多的材料点参与能量耗散过程。然而在不同加载速率下，原始近场动力学模拟得到的荷载-位移曲线的斜率基本没有变化，极限承载力有部分提升，并不显著，同时损伤破坏区域受加载速率影响的现象也十分微弱。通过对比得出，黏弹性近场动力学能够合理捕获混凝土刚度和强度都受加载速率影响的特性，并且能够反应混凝土非连续破坏的加载速率相关现象。

图6 不同加载速度下原始近场动力学（PD）与黏弹性近场动力学（VPD）的荷载-位移曲线对比Fig.6 Comparison of load-displacement curves between peridynamics and viscoelastic peridynamics at different loading velocities

图7 原始近场动力学与黏弹性近场动力学破坏模式对比Fig.7 Comparison of failure modes between original peridynamics and viscoelastic peridynamics

4 结论

在键型近场动力学方法基础上，发展了黏弹性近场动力学材料点相互作用模型，根据材料点相互作用模型推导了黏弹性键型近场动力学运动方程，通过与连续介质力学能量密度等效获得了材料参数的确定方法。利用动荷载下混凝土板的数值算例，验证了所提黏弹性近场动力学方法的有效性。主要结论如下：

（1）基于混凝土的变形机制，建立了材料点的黏弹性相互作用模型，基于哈密尔顿原理，推导了黏弹性运动方程。该方程合理反映了混凝土在变形过程中由非弹性变形所产生的能量耗散以及混凝土材料的动态强度与变形行为。通过近场动力学和连续介质力学之间的能量密度等效的方法，建立了近场动力学微弹性和微黏性参数的标定方法。该方法所得到的微弹性参数与原始近场动力学微弹性参数一致。

（2）利用黏弹性近场动力学方法模拟了不同加载速率下的混凝土构件的破坏模式，验证了该方法在模拟动荷载下混凝土强度、变形和裂纹扩展模式方面的合理性。通过模拟混凝土带有圆孔板在动态拉伸荷载下的力学性能，得到了受加载速率影响的破坏模式。

所提的黏弹性近场动力学方法能够合理反映混凝土在非弹性变形下的能量耗散以及不同加载速率和应力状态下的裂纹扩展和分叉，为混凝土动态力学性能的描述提供了一种有效的数值方法和工具。对于混凝土在更加真实的三维情况下的力学性能和破坏状态的模拟还有待进一步研究，并且对于复杂应力状态下的力学性能研究也亟待发展。

作者贡献声明：

路德春：项目负责人，提供研究平台、构思论文、指导模型构建及分析数据、修改论文。

宋志强：构建模型、分析数据、实现程序、呈现结果及撰写论文。

王国盛：项目负责人，提供理论及创新思路、提供编程帮助、构思与修改论文。