混凝土高压状态方程实验与数值模拟研究*

孙玉祥，王 杰,2，武海军，周婕群，李金柱，皮爱国，黄风雷

（1. 北京理工大学爆炸科学与技术国家重点实验室，北京 100081；2. 中国工程物理研究院化工材料研究所，四川 绵阳 621999；3. 中国工程物理研究院总体工程研究所，四川 绵阳 621999）

混凝土是一种脆性多孔材料，因其出色的抗压性能、低廉的造价，广泛应用于防护工事的建造中。在爆炸或高速侵彻、碰撞过程中，高压力荷载作用于混凝土局部区域，混凝土材料处于高静水压力状态。在这种状态下，材料的抗畸变能力（即剪切强度）可近似地忽略不计，此时，其本构关系简化成为静水压力p与体积应变μ之间的关系，混凝土可以视为无黏性的可压缩流体[1]。

Gebbeken 等[4-5]、Riedel 等[6-7]、王永刚等[8]对不同强度混凝土材料开展平板撞击实验，得到了用于构建混凝土高压状态方程的冲击Hugoniot 数据，通过对比，可发现不同强度混凝土其压力与密度关系不一致，材料强度是其状态方程的影响因素之一。此外，Grady 等[9]开展了压力范围3～25 GPa 的平板撞击实验，分析了Hugoniot 应力与残余应变的关系，认为混凝土冲击压缩下的压力-体积关系与混凝土孔隙和强度有关；Kipp 等[10]获得了两种骨料尺寸的混凝土在0.08～0.55 GPa 低压力段的平板撞击实验数据，并拟合得到波速与粒子速度关系；Hall 等[11]探究了密度接近、骨料尺寸ASTM 编号分别为57 和7 的两种混凝土高压状态差异，发现骨料尺寸对实验中混凝土加载与卸载响应影响不大；Tsembelis 等[12-13]通过平板撞击实验验证了不含骨料的水泥砂浆与含骨料的混凝土状态方程无明显差别，后将撞击压力提高至18 GPa，研究了水泥砂浆的剪应力强度和失效破坏。综上所述，静水压力为2 GPa 以上的平板撞击实验数据较少；其次，多项实验结果验证了骨料对混凝土高压状态方程无明显影响，而不同抗压强度的混凝土状态方程存在明显差异，表明了混凝土高压下的体积压缩特性与其基体抵抗压缩的能力有关。当前，同时开展多个强度混凝土状态方程的实验研究极少，缺乏对混凝土状态方程强度相关性的作用机理的深入认识。

因此，本文中基于火炮发射技术，对实测抗压强度为26.5 MPa 和42.1 MPa 的混凝土开展一维应变条件下的反向平板撞击实验，压力范围为2～11 GPa，拟合出混凝土材料在该压力范围的波速与粒子速度（us-up）、压力-体积应变（p-μ）关系。对比以往实验研究结果，分析混凝土抗压强度对其一维应变条件下压缩特性的作用机理。依据孔隙率确定混凝土压实密度，确定了HJC 模型高压段状态方程参数，通过反向平板撞击数值模拟与实验结果对比，验证模型参数的准确性，并讨论冲击波在混凝土材料中的传播特点。

1 平板撞击实验

1.1 实验装置与试样

利用 ∅ 58 mm 的火炮驱动混凝土飞片以v0的初速度撞击TU1 无氧铜靶板，撞击产生的平面冲击波传播至靶板自由面，造成自由面表面粒子速度变化，由安装在自由面后的多普勒探针系统（Doppler probe system，DPS）的探针探测并记录粒子速度变化，实验装置如图1 所示，图2 为DPS 测速装置与靶室布置情况。撞击过后，由于飞片与靶板横截面尺寸远大于其厚度，质点微团相互约束，不能侧向移动，处于三向应力状态，近似为一维应变条件。

测速装置DPS 的时间分辨率为50 s-1，位移分辨率为0.5 μm，测速范围为0.1 m/s～9 km/s，测速误差不超过1%，可保证测速精准。根据实验需求，共布置7 路DPS 探针，探针1 和探针2 位于靶板中心；探针3～6 沿以靶板中心为圆心、半径10 mm 的圆周均匀分布，其测量数据可与探针1 和探针2 的数据对比，排除飞片骨料与基体不均匀性的影响，并且可根据探针3～6 粒子速度变化的先后时间，计算飞片的碰靶角度；探针7 距离靶板中心26 mm，用于测量飞片的碰靶速度。

图1 飞片撞击实验装置示意图Fig. 1 Flyer-impact experimental device

图2 DPS 测速装置与靶室布置Fig. 2 DPS and target room arrangement

为防止飞片与靶板侧向稀疏波对中心区域冲击波的干扰，保证冲击波是均匀的平面波，飞片与靶板宽厚比须大于2[14]。因此，混凝土飞片设计厚度为12 mm，直径50 mm，端面平行度在0.05 mm 以内，表面平整度在0.02 mm 以内，骨料的直径2～5 mm；混凝土飞片试件，如图3 所示，飞片和弹托，如图4 所示。靶板材料使用密度为8 930 kg/m3的TU1 无氧铜，直径40 mm，厚度4 mm，实验前打磨去除金属氧化层，TU1 无氧铜靶片见图5。

实验设计两种强度混凝土，由42.5R 水泥、煤灰、中砂、细骨料按一定配比浇注而成，养护28 d。C25 混凝土实测抗压强度为26.5 MPa，密度为2 201.7 kg/m3，孔隙率为14.59%；C40 混凝土实测抗压强度为42.1 MPa，密度为2 176.8 kg/m3，孔隙率为15.73%。

图3 混凝土飞片试件Fig. 3 Concrete flyer specimen

图4 飞片与弹托Fig. 4 Flyer and sabot

图5 TU1 无氧铜靶片试件Fig. 5 TU1 copper specimen

为获得压力不低于2 GPa 的混凝土状态数据，拟开展C25 和C40 混凝土飞片实验各4 发，C25 飞片拟发射速度分别为500、1 000、1 200、1 500 m/s，C40 飞片拟发射速度分别为800、1 000、1 200、1 500 m/s。

1.2 实验原理

当飞片以v0的初速度撞击靶板后，压缩波从飞片与靶板撞击界面向两侧传播，波传播的X-t如图6所示。

由于冲击波在TU1 无氧铜中的传播速度快，压缩波在t1时刻到达靶板自由面，引起表面粒子速度跳跃。同时，自由面向介质内部反射拉伸波，而拉伸波经撞击界面反射为压缩波于t2时刻再次传播至靶板自由面，粒子速度再次跳跃，DPS 探针可记录到图7 中2 次压缩波到达自由面时间间隔（经估算约2.03 μs）内的粒子速度平台，其幅值大小为up1。

图6 冲击波传播X-t 图Fig. 6 Shock wave propagation

图7 TU1 靶板自由面粒子速度历史示意图Fig. 7 Particle velocity history curve on TU1 target

已知TU1 无氧铜波速us1与粒子速度up1关系如下：

式中：TU1 无氧铜体波波速cb=3 940 m/s，其冲击波波速与粒子速度（us-up）关系斜率S1=1.489。

冲击波在TU1 无氧铜靶板中的传播满足Hugoniot 跳跃条件，可求得靶板中的冲击Hugoniot 应力σH：

由于飞片与靶板撞击界面的连续性，σH同为飞片的Hugoniot 应力。一维应变条件下，Hugoniot 应力σH与静水压力p满足[15]：

式中：τmax表示混凝土最大剪应力，远小于σH，可忽略不计。因此：

在混凝土飞片中应用Hugoniot 跳跃条件，可分别求得飞片中的波速us2与粒子速度up2：

Grady[9]、Gebbeken 等[4-5]的实验结果表明，混凝土材料的冲击绝热关系可拟合为us=Sup+c形式，得到系数S和截距c。根据质量守恒和动量守恒方程，由下式：

即可计算实验中各个静水压力值p对应的密度ρ、体积应变μ。

1.3 实验结果与分析

每发实验均可获得6 路DPS 探针数据，图8 为C40 混凝土飞片以989 m/s 的速度撞击TU1 无氧铜靶板自由面粒子速度历史影像图，6 路探针原始影响经图像算法处理得到图9 所示靶板自由面粒子速度历史数据。为消除混凝土飞片不均性导致的数据误差，对各实验工况数据取6 路探针数据的平均值如图10～11 所示。取粒子速度第一平台段平均值为up1，依据式(1)～(6)计算Hugoniot 应力σH、混凝土波速us和粒子速度up，表1 为混凝土平板撞击实验数据处理结果。经测算，8 组平板撞击实验最大碰靶角度为0.31°，可近似认为混凝土飞片与无氧铜靶板发生正撞击。

图8 6 路DPS 探针原始影像图Fig. 8 Images by six DPS probes

图9 速度为989 m/s 的C40 混凝土飞片撞击TU1 靶板的自由面粒子速度历史Fig. 9 Particle velocity of TU1 free surface in C40 concrete flyer impact test at 989 m/s

图10 C25 混凝土飞片撞击TU1 靶板自由面粒子速度历史Fig. 10 Particle velocity of TU1 free surface in C25 concrete flyer impact tests

分别对C25 和C40 混凝土波速与粒子速度关系作线性拟合，如图12 所示，C25 拟合关系式为us=2.39(±0.27)up+1 251.7(±254.7)，C40 拟合关系式为us=1.82(±0.27)up+2 194.3(±257.9)。由式(7)～(9)计算得到各个静水压力值p对应的密度ρ、体积应变μ，如图13～14 所示。

平板撞击实验结果表明，对于同一强度的混凝土试件，其密度、体积应变随着静水压力的升高而增大。微观上，在高静水压力作用下，材料分子或原子间距离缩小，造成微观尺度上极小的体积压缩；宏观上，混凝土是一种多孔脆性材料，由于孔隙分布的随机性，质点微团在复杂的应力状态下发生剪切断裂，孔隙塌陷，表现出密度、体积应变的显著变化。

图11 C40 混凝土飞片撞击TU1 靶板自由面粒子速度历史Fig. 11 Particle velocity of TU1 free surface in C40 concrete flyer impact tests

图12 混凝土波速-粒子速度（us-up）关系拟合Fig. 12 Shock velocity and particle velocity linear fitting

表1 混凝土平板撞击实验数据Table 1 Flyer-impact test data of concrete

图13 混凝土压力-密度Fig. 13 Pressure-density data for concrete

图14 混凝土压力-体积应变Fig. 14 Pressure-volumetric strain data for concrete

此外，对比C25 和C40 混凝土压力体积应变关系发现，初始密度、孔隙率相近，强度不同混凝土压力与体积应变数据曲线是接近平行、区别明显的两条线，相同的静水压力p作用下，混凝土强度越高，体积应变越小。Grady[9]、Kipp 等[10]、Hall 等[11]已通过实验排除骨料对混凝土状态方程的影响，在初始密度、孔隙率接近的前提下，基体强度越高，质点微团发生剪切断裂需要达到的应力状态越高，孔隙塌陷越不完全，体积更难被压缩减小，因此体积应变更小。

2 HJC 状态方程参数确定

图15 HJC 模型状态方程[2]Fig. 15 EOS of the HJC model

如图15 所示，混凝土状态方程反映了该材料在受压缩过程的压力与体积关系的变化，HJC 模型中，混凝土状态方程分为3 段，分别为弹性段、过渡段和压实段。根据平板撞击实验结果，重新确定HJC 模型压实段参数pL、μL、K1、K2、K3。

HJC 本构模型中，为获得不含常数项的压力与体积应变的三次多项式关系，压实段压力-体积应变（p-μ）关系修正为压力-压实的体积应变（p- μ¯ ）关系。通过确定混凝的压实密度：

分别由C25 和C40 混凝土拟合的p-μ关系与式(10)～(14)确定p- μ¯ 的三次多项式系数，p-μ曲线与原始HJC 模型过渡段直线交点分别确定C25 和C40 混凝土对应的pL为0.9 GPa 和1.2 GPa。至此，HJC 模型中混凝土高压段状态方程参数确定完成，见表2。

表2 混凝土HJC 模型状态方程参数Table 2 EOS parameters of the HJC model for concrete

3 平板撞击实验数值模拟

基于LS-DYNA 动力学有限元分析软件，开展平板撞击实验数值模拟，检验基于该实验标定的HJC 模型状态方程参数的准确性，并进一步探究高压下混凝土状态响应。

3.1 有限元模型与材料参数

平板撞击实验数值模拟采用Lagrange 算法，建立四分之一模型，对称面设置对称边界，按照实验试件及其安装构件实际尺寸建立。模型分为3 个部分，分别是混凝土飞片、TU1 无氧铜靶片和靶片固定装置，各部件接触均为自动面-面接触。靶片固定装置边缘设置固定边界条件和非反射边界条件。设置5 个观测点如图16 所示，O点位于TU1 无氧铜靶板中心，记录O点速度变化，在混凝土飞片轴线上，从撞击界面开始，间隔3 mm，分布有A、B、C、D等4 个观测点，记录4 个位置点压力变化。

混凝土飞片材料模型为HJC 模型，除状态方程参数外，其他参数按照HJC 模型默认方法给出，设置单元不失效，主要材料参数见表3。TU1 无氧铜靶片与铝制靶片固定结构均采用Johnson-Cook 材料模型，采用Grüneisen 状态方程。

图16 数值模拟模型Fig. 16 Numerical simulation model

表3 混凝土HJC 模型主要参数Table 3 Parameters of HJC model for concrete

3.2 数值模拟结果与分析

提取O点对应单元的速度作为该点的粒子速度，图17 和图18 分别为该点在C25 和C40 混凝土反向平板撞击实验数值模拟粒子速度与实验粒子速度对比。数值模拟结果在粒子速度走势、平台段脉宽和幅值等方面与实验吻合度较高，说明数值模拟参数合理。对比同一发射速度下的前后两次粒子速度跳跃，后者增大速度变缓，反映了在TU1 靶板中冲击波扰动由强扰动变为弱扰动。

图17 C25 混凝土粒子速度数值模拟与实验结果对比Fig. 17 Comparison of simulation and test results for C25 concrete

图18 C40 混凝土粒子速度数值模拟与实验结果对比Fig. 18 Comparison of simulation and test results for C40 concrete

提取A、B、C、D点的压力曲线如图19，数值计算工况为C25 混凝土飞片以479 m/s 发射速度撞击靶板。飞片与TU1 靶片发生碰撞后，压缩波使得飞片中A点压力迅速升高，压力波形震荡剧烈，约2.03 μs 后，由TU1 无氧铜自由面反射的拉伸波到达A点，压力下降，后经过多次拉伸波卸载，最终压力卸载为零。在压缩波加载下，A点首先出现压力的升高，随后，B、C、D依次升高。由于压缩波作用过后，混凝土密实程度增大，波的传导能力增强，拉伸波将以更高的速度到达B、C、D点，因此A点首个矩形压力脉冲宽度最大，B、C、D依次变小。图20 为混凝土飞片压力云图。

图19 速度为479 m/s C25 混凝土飞片压力时程曲线Fig. 19 Pressure-time curves of C25 flyer at 479 m/s

图20 速度为479 m/s C25 混凝土飞片压力云图Fig. 20 Pressure cloud in C25 concrete flyer at 479 m/s

对比图19 和图21 不同撞击速度各点压力变化，较低的撞击速度下，前驱压缩波波速慢，靶板反射拉伸波追赶并造成前驱波压力卸载，因此，撞击速度为479 m/s 的混凝土中压力峰值由A向D依次减小，C、D两点峰值卸载尤为明显；而当飞片发射速度较高，混凝土前驱压缩波始终保持高于拉伸波的波速，拉伸波未能造成前驱波压力降低。其他实验数值模拟中，混凝土飞片中压力变化也呈现相似规律。

图21 速度为1 017 m/s C25 混凝土飞片压力时程曲线Fig. 21 Pressure-time curve of C25 flyer at 1 017 m/s

4 结 论

为探究高压下混凝土动态压缩特性，明确强度对混凝土状态方程影响，本文对C25 和C40 混凝土进行了一维应变动态力学行为的实验研究，通过反向平板撞击实验，获得了压力范围在2～11 GPa 的冲击Hugoniot 数据，由实验数据拟合了混凝土材料波速与粒子速度（us-up）关系，进一步推导出用于构建混凝土高压状态方程的压力-密度（p-ρ）、压力-体积应变（p-μ）数据，并依据实验结果讨论了混凝土基体强度造成混凝土状态方程差异的作用机理。通过混凝土孔隙率计算了混凝土压实密度，重新确定了两种强度混凝土HJC 状态方程参数，并通过反向平板撞击实验的数值模拟，分析了混凝土在高静水压下冲击波传播特征，主要结论如下：

（1）实验结果表明，混凝土材料在2～11 GPa 的高静水压力作用下持续被压缩，密度与体积应变随压力增大而增大，造成这一现象的主要原因是混凝土孔隙塌陷。

（2）对比C25 和C40 两种初始密度、孔隙率相近的混凝土实验结果，发现两者压力-密度或体积应变数据分布趋势相同，但不重合，呈平行关系。同一静水压力下，混凝土抗压强度越高，基体孔隙越难塌陷，体积应变就越小。但由于不同强度混凝土实验数据较少，尚不能给出量化的混凝土抗压强度与其状态方程的关系。

（3）平板撞击实验数值模拟结果显示，混凝土与TU1 无氧铜靶片中均存在着复杂的冲击波反射与透射。当飞片低速撞击靶片时，由于前驱压缩波造成材料的压实，在后继拉伸波快速追赶卸载作用下，压力降低；而高速撞击条件下，前驱冲击波始终保持了高于后继拉伸波的波速，未能发生冲击波的追赶卸载作用。