混凝土高壓狀態方程實驗與數值模擬研究*

孫玉祥，王 杰,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 無氧銅靶片中均存在著復雜的沖擊波反射與透射。當飛片低速撞擊靶片時，由于前驅壓縮波造成材料的壓實，在后繼拉伸波快速追趕卸載作用下，壓力降低；而高速撞擊條件下，前驅沖擊波始終保持了高于后繼拉伸波的波速，未能發生沖擊波的追趕卸載作用。