基于混凝土細觀骨料模型的SHPB仿真模擬研究

郭瑞奇， 任輝啟， 張 磊， 龍志林， 吳祥云， 李澤斌

(1. 湘潭大學 土木工程與力學學院，湖南 湘潭 411105； 2. 軍事科學院 國防工程研究院，河南 洛陽 471023)

混凝土是常見的工程材料，對防護工程而言需要考慮其動態力學性能。分離式Hopkinson壓桿(Split Hopkinson Pressure Bar，SHPB)實驗系統是廣泛應用于測試材料動態力學性能的裝置，可實測材料在高應變率(102～104s-1)下的應力應變曲線，也是研究混凝土材料動態力學性能的最基本的實驗手段。現今已有諸多學者利用不同直徑的SHPB裝置對混凝土材料進行了實驗研究[1-4]。

而隨著數值分析技術以及高速計算機的快速發展，有部分學者開始利用有限元分析技術來取代部分實驗，并利用數值仿真的方法來模擬SHPB實驗過程，從而定性地研究混凝土材料的動態力學性能，如巫緒濤等[5-6]模擬了用于混凝土沖擊試驗的O/100 mm的SHPB裝置，并討論了應力波在該裝置中的傳播特性以及罰函數算法中罰因子合理數值的選取問題；賈彬等[7]模擬了中國工程物理研究院結構力學研究所O/50 mm的SHPB實驗系統，確定了一套混凝土HJC(Holmquist-Johnson-Cook)模型參數確定方法；余道興等[8]分別運用了K & C(Karagozian and Case)，HJC和CSC(Continuous Surface Cap)三種材料模型對混凝土的SHPB試驗進行模擬。然而大多數研究都是從宏觀水平上將混凝土試件視為均勻材料，典型的如圖1所示。

圖1 將混凝土試件視為均勻材料的數值模擬[9]Fig.1 Uniform numerical model of concrete

由于混凝土是一種非均勻材料，其內部含有大量隨機分布的粗骨料顆粒，用SHPB對混凝土進行實驗時不可避免地要面對非均勻性問題。這種非均勻性問題主要來自于兩個方面：一方面是由夾雜顆粒導致的材料非均勻性；另一方面是應力波效應，SHPB實驗不考慮試件中波的傳播，也就是兩個基本假定之一——“短試件的應力/應變沿其長度均勻分布(動態平衡)[10]”，這對于常規的金屬材料是易于滿足的，而對混凝土材料而言，波在“骨料——砂漿基體”界面傳播時會存在復雜的反射和透射，因此使用均勻的混凝土模型來模擬SHPB實驗或許與真實實驗結果有些偏差。

本文將混凝土視為由砂漿基體和粗骨料組成的雙相復合材料，將骨料簡化為球形，建立了圓柱狀三維混凝土隨機骨料模型。以O/100 mm的SHPB實驗裝置為例，建立了整體實驗裝置的有限元模型，并通過數值模擬的方法從細觀層次上研究非均勻模型在SHPB仿真模擬中的有效性。通過五種不同的梯形波加載以及五種不同骨料含量模型的計算研究了應變率及骨料含量對混凝土動態強度和應力均勻性的影響。最后建立了四種不同骨料尺寸的混凝土模型，并分析了骨料尺寸對模擬結果的影響。

1 混凝土細觀骨料模型的建立

近十幾年來，從Wittmann等[11]提出的二維多邊形骨料模型到Sheng等[12-13]提出的三維凸多面體隨機骨料模型，逼近真實混凝土材料的高含量細觀骨料模型的建模方法已經較為成熟，作者也對混凝土細觀骨料模型的建立、求解算法和求解效率進行了相關研究[14-16]。近年來，已有學者開始利用細觀模型對鋼筋混凝土梁的抗沖擊性能[17]、混凝土的抗侵徹性能[18]以及鋼纖維混凝土的動態強度[19]進行了數值模擬研究。

本文采用郭瑞奇等研究中的方法，建立了二級配圓柱狀三維混凝土骨料模型，其中小石粒徑為5～20 mm，中石20～40 mm，體分比為55∶45。為了保證所生成的骨料顆粒不超過給定的圓柱區域，需要對球心坐標公式進行修正

(1)

式中：R為圓柱底面半徑；Ri為第i顆骨料的半徑；ZB和ZT分別為圓柱體底部和頂部邊界坐標值；rdm1，rdm2，rdm3分別為位于(0,1)之間的三個獨立隨機數。

首先在Fortran程序中隨機生成ANSYS參數化建模所需要的命令流，建立幾何模型以后使用Hypermesh進行有限元網格剖分，最后使用LS-DYNA進行有限元計算，整個模型建立的流程如圖2所示。

圖2 混凝土細觀骨料模型建立流程圖Fig.2 Flow chart of the concrete aggregate model

以簡單的球形骨料為例，建立了骨料含量為46.06%的混凝土模型，整個試件厚度為50 mm，底面半徑為50 mm，由于細觀骨料模型的高度幾何復雜性，選用solid168十節點四面體單元進行網格剖分。入射桿長5 m，透射桿長3 m，直徑均為100 mm，對于這種規則形狀采用solid164六面體單元映射網格劃分。有限元模型建立的各個階段如圖3所示，依次是幾何模型、骨料顆粒網格剖分、砂漿基體網格剖分和整體有限元模型x=0 mm(中心軸線處)截面示意圖。

圖3 數值模型建立過程Fig.3 Modeling procedure of the numerical model

2 混凝土細觀骨料模型有效性驗證

為了驗證混凝土細觀骨料模型用于SHPB仿真實驗的有效性，以上文建立的骨料含量46.06%的混凝土試件為例，進行動力學模擬分析，并將結果與實際實驗數據進行比較。

2.1 材料本構模型及參數選擇

HJC本構模型常用于混凝土材料的動態力學仿真中，使用特征化等效應力對HJC強度模型的描述為

(2)

HJC模型狀態方程分為彈性壓縮(OA段)、壓實變形(AB段)和密實后變形(BC段)三個階段，如圖4所示。

圖4 HJC材料模型狀態方程曲線Fig.4 Equation of state curve for HJC model

方秦等[20]從塑性屈服面理論出發，合理確定了巖石材料的HJC本構模型參數，因此本文首先由文獻[21]中的數據分別確定骨料和砂漿基體的HJC材料參數，然后根據方秦等研究中對花崗巖材料的研究對骨料參數進行修正，砂漿的模型參數和骨料的修正數據分別如表1和表2所示，計算時統一換算成cm-g-μs單位制。

表1 砂漿HJC模型參數Tab.1 HJC model parameters for mortar

表2 花崗巖HJC模型參數修正值Tab.2 HJC model parameters revised for moorstone

本文并未引入額外的失效準則來控制單元失效，僅使用HJC模型中自帶的失效類型參數作為單元的破壞準則，即FS=0.002。 入射桿和透射桿采用線彈性模型，取彈性模量E=200 GPa， 密度ρ0=7.8×103kg/m3， 泊松比ν=0.3。在三種材料之間分別設置三對接觸，其中骨料與砂漿之間采用面面固連接觸，接觸剛度罰因子取默認值。另外兩種采用自動面面接觸，罰因子參照巫緒濤等的研究取為“2.0”。

2.2 加載計算及試件應力均勻性分析

由于混凝土細觀骨料模型的幾何復雜性，無法采用對稱模型。在SHPB仿真模擬中為保證求解精度，整體模型單元數量較多，計算規模較大，因此本文采用直接在入射桿端部施加應力脈沖的方式加載。

在桿端z=0處分別作用圖5所示的五種不同梯形脈沖，分作兩組對比實驗。第一組為加載時間相同但應力幅值不同的三種梯形荷載，分別為σ1=100 MPa，σ2=150 MPa，σ3=200 MPa， 其上升沿為60 μs，總加載歷時為240 μs。第二組為應力幅值均為100 MPa但加載時間不同的三種梯形載荷σ1，σ′1，σ″1中上升沿分別為60 μs，80 μs，120 μs，總加載歷時分別為240 μs，360 μs，480 μs，程序總計算時間2 400 μs。

圖5 五種不同的梯形載荷Fig.5 Five kinds of trapezoidal loads

在Hopkinson桿實驗中，試件內部的應力均勻與否是影響實驗結果有效性的一個重要因素。應力均勻性假定對于金屬類均質材料往往是易于實現的，而對于混凝土材料而言，其試件厚度比一般金屬類材料大了一個量級，且具有破壞應變小、含有粗骨料夾雜顆粒、波速較低等特點，因此必須對其在實驗過程中的應力均勻性進行分析，以驗證實驗結果的有效性。通過比較試件兩端的應力歷史即可檢測其應力平衡過程，因此研究人員先后提出了類似的應力差之比的公式作為試件應力平衡的準則[22-23]

(3)

式中：σ1(t)和σ2(t)分別為試件前后端面的應力歷史， 應力平衡因子δ(t)已經被一些學者用于評估SHPB實驗中試件的應力均勻性問題[24]， 通常認為δ(t)<0.05時試件達到了較好的應力平衡狀態。

相對于均勻的混凝土模型而言，細觀骨料模型更能逼近真實的試件結構。本文在對混凝土細觀骨料模型用于SHPB仿真實驗的結果分析之前，首先對試件在應力波作用過程中的應力均勻性進行研究。

在試件前后端面取三對單元，它們與中軸線的距離分別為0 mm，25 mm，50 mm，三種載荷作用下的應力平衡因子—時間曲線如圖6所示。可以看到在100 MPa的梯形載荷作用下試件的整體應力均勻性相對較好，在1 070 μs左右達到應力平衡最低點以后，δ(t)處于小幅震蕩狀態。同時從“3.1”節的相關結果中可以看到，在1 200 μs左右時由于試件表面開始出現宏觀裂紋，最外側單元的δ(t)出現了波動。除此之外，在應力波作用的大部分時間里，試件內外單元的應力平衡因子變化趨勢都相對穩定。而在另外兩種較高幅值的梯形載荷作用下，試件的應力均勻性相對較差，由于宏觀裂紋的過早出現，整個試件難以保持穩定的應力平衡狀態，內外單元的應力平衡因子到達最低點以后便開始劇烈震蕩，而且試件破壞越嚴重，這種震蕩就越劇烈。

此外，無論哪種載荷作用下，試件達到應力平衡狀態(δ(t)≤0.05)或應力平衡因子到達最低點都需要一段較長的時間(不低于100 μs)。這表明在實際實驗中，應該在保證達到目標應變率的前提下，適當延長應力脈沖升時，從而使混凝土試件有充足的時間達到應力平衡狀態。

圖6 不同幅值載荷作用下試件的應力平衡因子Fig.6 Balance factor of the concrete specimen under different loads

2.3 模型有效性驗證

取入射桿和透射桿中部單元作為數據采集點，由LS-DYNA計算出的五種軸向應力波形如圖7所示。

圖7 五種不同的梯形載荷作用下入射桿和透射桿 中部單元軸向應力波形Fig.7 Axial stress wave of the middle elements in incident bar and transmitted bar under five kinds of trapezoidal loads

(4)

(5)

(6)

式中：A0為桿的橫截面積；As和ls為試件的橫截面積和長度；E0和c0分別為入射桿的彈性模量和波速。首先通過式(5)求得試件在五種載荷作用下的應變率—時間曲線，如圖8所示。

圖8 混凝土試件的應變率-時間曲線Fig.8 Strain rate-time curve of concrete specimen

從圖7和“3.1”節內容中可以看到，隨著應變率的提高，混凝土破壞后反射波信號會快速增加，直接處理得到的應變率曲線將不正確。因此，我們在圖8中標注了混凝土細觀骨料模型細小損傷出現的時間點，在混凝土破壞之前的應變率才是混凝土試樣的主要加載應變率。使用幅值為100 MPa的三種梯形脈沖加載時，試件在應力波作用階段其應變率主要位于20～50 s-1，隨著載荷作用時間的增加，試件的應變率相對而言更加穩定。而當加載脈沖幅值為150 MPa及200 MPa時，試件的應變率主要位于50～100 s-1，雖然這兩種載荷升時均為60 μs，但應力幅值越高，加載脈沖上升沿越陡，應變率曲線的平臺段也越短，難以實現恒應變率加載。通過式(4)和式(6)分別求出試件的應力—時間曲線和應變—時間曲線后，消去時間軸即可得到試件的平均應力—應變曲線。巫緒濤等和呂太洪等[26]給出了混凝土材料在SHPB實驗中得到的4個應變率下的應力—應變曲線，劉傳雄等[27]也對骨料尺寸為15～20 mm的混凝土材料試樣進行了應變率范圍30～180 s-1的動態壓縮實驗。本文首先將加載時間為240 μs而應力幅值不同的三種波形作用下混凝土試件的動態應力應變曲線進行比較，然后參考文獻中對應應變率范圍內的應力應變曲線，將實驗數據與數值模擬的結果進行比較分析，如圖9所示。

圖9 仿真模擬結果有效性驗證Fig.9 Validation of analog simulation results effectiveness

由于應力波的傳播速度遠高于裂紋擴展的速度，會使混凝土試件表現出“應變滯后效應[28]”，在圖9(a)所示的三種不同應變率下混凝土試件的應力應變曲線上表現為應力相同時試樣的應變隨著應變率的提高而降低，也就是材料動態彈性模量隨應變率的提高而逐漸增大(應變率效應)。幅值為150 MPa的梯形波作用下試件的應力應變曲線與100 MPa的結果相比，應力峰值沒有明顯變化，但是試件的動態彈性模量有所增強，且卸載段由黏彈性回滯現象變為應變軟化現象。在沖擊荷載作用下，混凝土材料內部同時出現應變率硬化和損傷軟化效應，從這兩種載荷作用下試件應力應變曲線的比較可以看到，應變率的提高使得軟化效應更為明顯，損傷效應和應變率硬化效應達到平衡，這種飽和現象在這里表現為應力峰值并沒有明顯地隨應變率的提高而增加，王道榮等[29]對骨料顆粒較大的混凝土材料進行SHPB實驗時也發現了同樣的現象。幅值為200 MPa的梯形波作用下試件的應力應變曲線與150 MPa的結果相比，試件的動態彈性模量和應力峰值都有所提高。

通過實驗數據與仿真結果的比較可以看到，由于不同研究人員所采用的實驗裝置各不相同，而且制作的混凝土試件其骨料類型、骨料含量等都具有一定的隨機性，所以實際的試驗數據本身也有一定的區別。與劉傳雄等研究中的結果相比，本文通過仿真模擬所得的應力應變曲線在彈性階段和卸載段都相對較陡，我們認為前者的主要原因是混凝土試件在實際加工制作中必然會存在一定的平行度公差，這樣會導致實驗測得的彈性模量較低，而模擬中試樣兩端是理想狀態下絕對平行的，因此曲線在彈性階段相對較陡。后者是由于在仿真模擬中，當單元達到破壞應變后被刪除，所形成的孔洞和裂紋不會被填補，而實際實驗中混凝土存在一個壓緊密實的過程，因此其曲線在卸載段相對緩和。除此之外，盡管數值模擬中的加載波形、材料參數及實驗裝置與實際實驗有所不同，但使用混凝土細觀骨料模型仿真所得到的應力—應變曲線，其體現的動態力學行為與實驗結果是十分相近的，包括應變率較低時出現的黏彈性回滯現象(見圖9(b))和應變率較高時的損傷軟化現象(見圖9(c))，因此使用混凝土細觀骨料模型可有效模擬SHPB實驗并反應材料的動態力學性能。

3 不同應變率及骨料含量對模擬結果的影響與分析

3.1 不同應變率下混凝土細觀骨料模型的破壞描述

圖10和圖11給出了幅值為100 MPa的梯形載荷作用下應力脈沖由入射桿末端穿過短試件進入透射桿的整個過程：t=980 μs時，應力波到達入射桿末端，即將進入混凝土試件；t=1 020 μs時，應力波剛剛穿過混凝土試件，由于材料的不均勻性，可以看到整個試件軸向應力的不均勻分布，且透射波波陣面已不再保持為平截面；t=1 060 μs時，試件兩端面開始出現應力集中現象，且端面有部分損傷出現；t=1 160 μs時，試件表面開始出現應力集中現象，試件與入射桿和透射桿接觸面出現明顯損傷；t=1 200 μs時由于骨料與砂漿之間的相互作用，試件表面出現應力集中現象的部位開始出現條狀鋸齒樣裂紋；t=1 380 μs時，應力脈沖已經穿過試件。

圖10 應力波作用于混凝土細觀骨料模型過程的應力云圖Fig.10 Stress nephogram of concrete model with mesoscopic aggregates subjected to stress wave

隨著加載幅值和應變率的升高，混凝土細觀骨料模型表面裂紋的演化更加明顯，細觀骨料模型的破壞過程主要表現為試件兩端和表面出現的細小損傷沿粗骨料間隙的擴展與延伸。幅值150 MPa和幅值200 MPa梯形載荷作用下1 200 μs時試件與粗骨料的軸向應力云圖分別如圖12和圖13所示，可以看到，加載幅值越高，損傷演化越迅速，同一時間下高幅值加載時試件表面的裂紋已經開始融匯貫通，破壞更嚴重。同時從圖7中也可以看到，試件的迅速破壞會導致波阻抗降低，透射波隨之減小，相應的反射波會快速增加。

圖11 不同時間下混凝土細觀骨料模型的應力云圖Fig.11 Stress nephogram of concrete model with mesoscopic aggregates at different times

圖12 幅值150 MPa梯形載荷作用下混凝土試件應力云圖Fig.12 Stress nephogram of concrete specimen subjected to trapezoidal load with the amplitude of 150 MPa

圖13 幅值200 MPa梯形載荷作用下混凝土試件應力云圖Fig.13 Stress nephogram of concrete specimen subjected to trapezoidal load with the amplitude of 200 MPa

3.2 不同骨料含量對試件應力均勻性及動態強度的影響

為了研究粗骨料含量對混凝土試件的應力均勻性及動態強度的影響，建立了四種不同骨料含量的混凝土細觀骨科模型(見圖14)，模擬過程中統一施加σ1=100 MPa的梯形載荷，按照上文所述的方式計算出每種試件的應力平衡因子，結果如圖15所示。

圖14 四種不同骨料含量的混凝土模型Fig.14 Four concrete models of different aggregates content

圖15 骨料含量對試件應力平衡性的影響Fig.15 Effects of aggregates content on stress uniformity of concrete specimen

由于每種模型中的骨料為隨機分布，在試件兩端面所取單元應力時程計算出的應力平衡因子只能反映出局部的平衡狀態，從圖15中可以看到試件內外單元的應力平衡性好壞并沒有明顯的規律性。但是結合圖6(a)可以發現，隨著骨料含量的增加，混凝土試件越來越密實，試件內外的應力平衡變化逐漸趨于一致，試件整體的應力均勻性更好。

除了一些輕骨料混凝土和珊瑚骨料混凝土，通常粗骨料的強度要高于水泥石的強度以及混凝土的系統強度，這種情況下適當增加骨料含量能夠起到提高混凝土整體強度的效果。取入射桿和透射桿中部單元的應力波形做出每種試件的應力應變曲線，如圖16所示。結果表明，隨著骨料含量的增加，試件的動態彈性模量和峰值應力都有所提高。隨著試件整體動態強度的提高，產生單位變形所需要的能量增加，在圖16中表現為同一載荷作用下，骨料含量越高則試件所能達到的最大應變逐漸減小。

圖16 骨料含量對混凝土試件動態強度的影響Fig.16 Effects of aggregates content for dynamic strength of concrete specimen

3.3 骨料尺寸對試件應力均勻性及動態強度的影響

為了研究粗骨料尺寸對混凝土試件的應力均勻性及動態強度的影響，建立了四種不同骨料尺寸的混凝土細觀骨料模型(見圖17)。為了避免骨料含量對結果的影響，將骨料體積率控制在30%左右，四種試件的相關參數見表3。

圖17 四種不同骨料尺寸的混凝土模型Fig.17 Four concrete models with aggregates at different size

試樣編號骨料粒徑范圍/mm最小粒徑最大粒徑骨料含量/%骨料個數a51030.02472b101530.09119c152030.1245d202530.7321

從上面可以看到骨料尺寸越小，在體積率相同的情況下其顆粒數量大大增加，在試件內部的分布就越近似于均勻狀態。四種試件的應力平衡因子—時間曲線如圖18所示，隨著骨料尺寸的減小，試件的應力平衡因子在到達低點以后變化更趨于穩定。而且相對于大尺寸骨料的混凝土模型而言，小尺寸骨料混凝土模型內外應力平衡因子趨勢的一致性更好，表明其整體應力均勻性更好。

四種模型的應力應變曲線如圖19所示，可以看到骨料粒徑在10～25 mm內的三個試件動態強度沒有明顯區別，而且與圖16中骨料含量30.25%的二級配混凝土模型的動態應力應變曲線重合度較高。但是當骨料尺寸在5～10 mm內時，含量為30.02%的混凝土動態彈性模量已經相當于骨料含量40.78%的二級配混凝土，所能達到的應力峰值相當于骨料含量為46.06%的二級配混凝土。

圖18 骨料尺寸對試件應力平衡性的影響Fig.18 Effects of aggregate size on stress uniformity of concrete specimen

圖19 骨料尺寸對混凝土試件動態強度的影響Fig.19 Effects of aggregate size for dynamic strength of concrete specimen

我們對骨料尺寸為5～10 mm內的混凝土模型以及20～25 mm內的混凝土模型進行比較分析，做出1 200 μs時沿試件中心軸線處的截面圖以及對應的軸向應力云圖，如圖20所示。由于本文采用了較為堅硬的花崗巖骨料，混凝土的破壞形式表現為裂紋全部繞骨料與砂漿的界面進行擴展，而骨料本身沒有破壞。這種情況下骨料尺寸的增加會引起骨料周圍較為嚴重的應力集中現象，在應力波作用過程中會產生更多的微裂紋，從而降低混凝土的整體強度。而骨料尺寸較小的情況下，數量較多的骨料在試件內部的分布更均勻，這樣能夠使骨料充分受力從而更好地發揮其強度以及骨架作用。

除此之外，實際實驗中混凝土材料在養護時由于骨料和砂漿的收縮不同會產生大量的微裂紋。隨著骨料尺寸的增加會導致其與砂漿模量的差別更加明顯，產生更多的初始裂紋，從而進一步降低混凝土材料的整體強度[30-31]。

圖20 1 200 μs時不同骨料尺寸的混凝土模型應力云圖Fig.20 Stress nephogram of concrete models with different aggregate size at 1 200 μs

4 結 論

本文建立了圓柱狀三維混凝土細觀骨料模型，并將其應用于SHPB的仿真實驗中。通過五種不同梯形載荷的加載計算，得出了試件不同應變率下的應力—應變曲線，并將模擬結果分別與相應應變率下的試驗結果進行對比分析，驗證了本文所建立的混凝土細觀骨料模型用于SHPB仿真實驗的有效性。

通過對比試件在三種不同幅值載荷作用下的應力平衡因子發現，在較低幅值梯形載荷作用下試件的整體應力均勻性相對較好，在應力波作用的大部分時間里，試件內外單元的應力平衡因子變化趨勢都相對穩定。而在另外兩種較高幅值的梯形載荷作用下，試件的應力均勻性相對較差，由于宏觀裂紋的過早出現，整個試件難以達到應力平衡狀態，內外單元的應力平衡因子到達最低點以后便開始劇烈震蕩，而且試件破壞越嚴重，這種震蕩就越劇烈。在實際實驗中，應該在保證達到目標應變率的前提下，適當延長應力脈沖升時，從而使混凝土試件有充足的時間達到應力平衡狀態。

對五種不同骨料含量混凝土模型的計算分析結果表明，在一定的范圍內隨著骨料含量的增加，試件內外的應力平衡變化逐漸趨于一致，試件整體的應力均勻性更好；骨料含量越高，混凝土試件整體的動態強度越高，但同時所能達到的最大應變也逐漸減小。

對四種不同骨料尺寸混凝土模型的計算分析結果表明，當骨料尺寸減小到一定程度時，相同體積率下骨料數量大大增加，分布更加均勻。試件的應力均勻性更好，強度更高。

此外，本文所建立的基于細觀骨料模型的SHPB整體試驗裝置有限元模型，只需改變相應的材料類型和相關參數，其同樣適用于泡沫鋁、珊瑚砂、礁灰巖等多孔材料和含夾雜顆粒的非均質材料的仿真實驗研究，為多種材料的SHPB細觀數值模擬提供方便。