剛性彈體對素混凝土厚靶侵徹響應(yīng)的LDPM數(shù)值模擬研究

馮 君， 李文彬， 徐 磊， 陳 宇， 陳 曦

(1. 南京理工大學(xué) 瞬態(tài)物理國家重點實驗室，南京 210094； 2. 南京理工大學(xué) 機械工程學(xué)院，南京 210094； 3. 河海大學(xué) 水利水電學(xué)院，南京 210098； 4. 南京炮兵學(xué)院 偵測系，南京 211132)

隨著鉆地武器發(fā)展及防護結(jié)構(gòu)評估需要，混凝土侵徹問題成為當(dāng)前沖擊工程領(lǐng)域研究熱點問題。彈體對混凝土侵徹是一個伴隨著靶體材料大變形、高應(yīng)變率和高壓的過程，涉及斷裂、擠壓、破碎等復(fù)雜響應(yīng)[1]。采用空腔膨脹分析等理論方法獲得控制方程組的求解非常復(fù)雜且侵徹過程的工程模型存在太多簡化和假設(shè)[2-4]，因此基于理論的工程分析模型在研究復(fù)雜工況的侵徹問題存在局限性。近年來隨著多尺度模擬等計算力學(xué)技術(shù)的發(fā)展,出現(xiàn)了越來越豐富的數(shù)值計算方法(SPH、RKPM、物質(zhì)點法、近場動力學(xué)法等[5-6])，仿真技術(shù)成為研究混凝土侵徹問題的重要手段。研究表明，離散元法在模擬分析混凝土破壞失效等不連續(xù)問題中具有巨大優(yōu)勢[7-10]。

基于離散單元基本假設(shè)，Cusatis等[11-12]開發(fā)出了一種模擬水泥基體內(nèi)部骨料相互作用的細觀離散元模型(Lattice Discrete Particle Model，LDPM)。根據(jù)混凝土骨料級配信息在幾何體內(nèi)部隨機生成骨料顆粒，并在其周圍構(gòu)建潛在裂紋微面以及細觀本構(gòu)關(guān)系。LDPM模型能夠模擬混凝土材料斷裂、多軸響應(yīng)、尺寸效應(yīng)等，并在模擬結(jié)構(gòu)受到?jīng)_擊與爆炸等動態(tài)載荷工況中獲得較好結(jié)果[13-14]。對比有限元模型中采用刪除單元來模擬侵徹、裂紋擴展等現(xiàn)象，離散元模型避免了單元刪除帶來的質(zhì)量、能量不守恒。針對23 MPa強度素混凝土三軸壓縮實驗結(jié)果，本文標定了LDPM參數(shù)，通過模擬混凝土侵徹驗證了模型適用性。并對彈頭形狀和彈體尺寸等因素對侵徹響應(yīng)的影響進一步數(shù)值模擬研究。

1 LDPM簡介

1.1 LDPM離散元基本假設(shè)

首先，根據(jù)Fuller級配曲線將骨料顆粒等效成不同直徑的圓球隨機投放在混凝土體內(nèi)部。由于水泥基體材料以及骨料與基體交界面相比于骨料，其強度和剛度都要弱很多，LDPM假設(shè)混凝土裂紋都發(fā)生在基體和基體-骨料接觸面，通過Delaunay技術(shù)等生成三角形潛在裂紋微面包絡(luò)骨料顆粒圓球的LDPM單元[15]，如圖1所示。

圖1 LDPM單元及其內(nèi)部骨料顆粒[15] Fig.1 LDPM element and aggregate particle[15]

采用剛體動力學(xué)描述細觀模型系統(tǒng)的變形，將細觀應(yīng)變用潛在裂紋微面中心位置的位移[uc]來表示：

(1)

式中：l是相鄰骨料間距，nT，lT和mT是定義潛在裂紋微面上局部坐標系的單位向量。在文獻[16-17]中，Cusatis證明了式(1)的細觀應(yīng)變定義與連續(xù)介質(zhì)力學(xué)中應(yīng)變在多面體微面上投影分量一致。

1.2 LDPM本構(gòu)關(guān)系

在彈性響應(yīng)階段，法向和切向細觀應(yīng)力與對應(yīng)的細觀應(yīng)變呈線性關(guān)系：tN=ENeN；tM=ETeM；tL=ETeL， 其中EN=E0，ET=αE0,E0是等效法向彈性模量，α是切向-法向耦合系數(shù)。對于非線性響應(yīng)階段，LDPM主要考慮以下三個方面的細觀力學(xué)響應(yīng)：斷裂、空隙坍塌和壓實以及壓剪耦合內(nèi)摩擦。

法向壓縮強度滿足-σbc(eD,eV)≤tN≤0, 其中σbc是體應(yīng)變相關(guān)的細觀壓縮邊界條件，偏應(yīng)變eD=eN-eV。彈性響應(yīng)過后，由于基體內(nèi)部孔隙坍塌，導(dǎo)致法向應(yīng)力以較小的剛度Hc隨著體應(yīng)變增大而增大，σbc=σc0+〈-eV-ec0〉Hc(rDV)，其中：Hc(rDV)=Hc0/(1+κc2〈rDV-κc1〉)，σc0是細觀壓縮強度；rDV=eD/eV表示偏應(yīng)變與體應(yīng)變的比,κc1和κc2是材料參數(shù)。當(dāng)-eV≥ec1時，密實硬化發(fā)生，σbc=σc1(rDV)exp[(-eV-ec1)Hc(rDV)/σc1(rDV) ]，其中σc1(rDV)=σc0(rDV)+(ec1-ec0)Hc(rDV)。

基于黏聚斷裂理論[18]，LDPM中裂紋開裂速率相關(guān)的動態(tài)屈服準則為：

(2)

2 LDPM模擬23 MPa強度混凝土侵徹

2.1 23 MPa強度混凝土LDPM參數(shù)標定

Warren等[19]根據(jù)美國ASTM標準對111.0 mm高，50.3 mm直徑的混凝土圓柱體進行三軸壓縮實驗，裝置示意圖如圖2所示。仿真中通過增大試件表面壓力，獲得混凝土試件靜水壓力響應(yīng)；而三軸壓縮加載則分為兩步，首先增大試件表面壓力至所需的圍壓值并保持壓力不變，再通過控制軸向位移進行軸向壓縮。

圖2 三軸壓縮實驗示意圖[19] Fig.2 Triaxial compression test setup[19]

LDPM模擬混凝土材料需要其對應(yīng)的配比信息來生成離散的骨料顆粒，根據(jù)文獻[19]中給出的23 MPa混凝土配比信息，將1 m3單位體積混凝土中所含的各組分質(zhì)量列在表1中。

表1 單位體積混凝土各組分質(zhì)量Tab.1 Mixture weight of unit volume concrete

Cusatis等研究中LDPM參數(shù)取值，模擬混凝土圓柱體靜水壓力實驗，標定參數(shù)如下：等效法向彈性模量E0=16 500 MPa，切向-法向耦合系數(shù)α=0.25，細觀壓縮強度σc0=42 MPa，初始硬化模量Hc0=0.45E0，漸

近轉(zhuǎn)化應(yīng)變比κc0=6。實驗與標定結(jié)果仿真對比如圖3所示。

通過三軸壓縮仿真(如圖4)，標定LDPM參數(shù)如下：漸近轉(zhuǎn)化應(yīng)力σN0=400 MPa，初始內(nèi)摩擦因數(shù)μ0=0.4，漸近內(nèi)摩擦因數(shù)μ∞=0，偏應(yīng)變閾值比κc1=2，偏量損傷系數(shù)κc2=2。其余參數(shù)見Cusatis等的研究取值為：細觀拉伸強度σt=2.7 MPa，細觀剪切強度σs=2.0σt，拉伸特征長度lt=100 mm。

圖3 靜水壓力實驗仿真 Fig.3 Hydrostatic pressure test simulation

圖4 三軸圍壓實驗仿真 Fig.4 Triaxial compression test simulation

圖5 不同圍壓混凝土三軸壓縮仿真試件破壞形態(tài) Fig.5 Numerical modelling of concrete damage mode under various confined pressures

圖5給出了LDPM模擬0 MPa、50 MPa、100 MPa、200 MPa和400 MPa五種圍壓下，混凝土圓柱試件三軸壓縮的裂紋分布及其破壞形態(tài)，五種工況下軸向應(yīng)變εa分別為0.005、0.082、0.107、0.123和0.139。當(dāng)圍壓為0時，無側(cè)限單軸壓縮載荷下圓柱試件破壞模式為劈裂破壞；隨著側(cè)向圍壓的增大，試件破壞模式變成剪切和粉碎破壞。

2.2 混凝土侵徹數(shù)值模型

針對23 MPa強度素混凝土厚靶，F(xiàn)orrestal等[20]進行了彈體發(fā)射速度在139～379 m/s之間的系列侵徹試驗。試驗中所用的尖卵形頭部彈體直徑為80.01 mm，CRH為3和6的兩種彈體具體尺寸如圖6所示。由于本文研究的侵徹試驗后彈體頭部幾乎沒有形狀變化和質(zhì)量損失，為了提高計算效率，本文選擇剛性彈體進行下文侵徹仿真分析。

圖6 彈體示意圖[20] Fig.6 Diagram of projectiles[20]

彈體與混凝土之間接觸設(shè)置為罰函數(shù)接觸，其接觸剛度為時間步長相關(guān)的函數(shù)，其具體表達式如下

(3)

式中：mmin是發(fā)生接觸的顆?；蛘邌卧匈|(zhì)量最小值，Δt是顯示計算步中的時間步長。罰函數(shù)接觸力Fn=Ksif·p，p是接觸體相互侵入量。

2.3 LDPM模擬混凝土侵徹結(jié)果分析

文獻[20]中8發(fā)混凝土侵徹試驗及其LDPM仿真侵徹深度結(jié)果如表2所列，其中1-6號試驗中采用CRH

為3的彈體，7號和8號試驗采用CRH為6的彈體。LDPM模擬剛性彈體以139～379 m/s撞擊速度侵徹23 MPa強度混凝土厚靶，仿真侵徹深度與試驗值基本吻合，其誤差在8.3%以內(nèi)。而CRH為6的7號、8號侵徹試驗的預(yù)測值與試驗更加接近，誤差分別為3.4%和3.2%。對比試驗6和8，在相同彈體撞擊速度情況下，仿真和試驗均表明CRH為6的彈體侵徹深度值更大。

表2 混凝土厚靶侵徹深度仿真與實驗對比Tab.2 Numerical and experimentalpenetrationdepth comparison

圖7給出了不同侵徹工況下仿真和試驗的彈體減速度隨著時間變化的對比，LDPM仿真得到的彈體過載峰值、曲線變化趨勢以及過載持續(xù)時間基本和試驗曲線相符。LDPM模擬彈體以139 m/s撞擊速度侵徹23 MPa強度混凝土，在0.1 ms、0.2 ms、0.3 ms和0.4 ms時刻彈靶作用情況如圖8所示。結(jié)合圖7(a)彈體過載曲線，在0.4 ms時刻，彈體基本完成侵徹，此時可以看出靶體迎彈面附近產(chǎn)生明顯的開坑區(qū)域。

圖7 彈體過載時程變化仿真與實驗對比 Fig.7 Experimental and numerical comparison of projectile deceleration vs. time

圖8 彈體以139 m/s撞擊速度侵徹厚靶的裂紋變化過程仿真模擬 Fig.8 Concrete target cracksrevolutionsimulation during penetration by 139 m/s striking velocity projectile

3 侵徹阻應(yīng)力LDPM模擬分析

結(jié)合前文標定的23 MPa強度混凝土參數(shù)，本節(jié)通過彈體恒定速度侵徹混凝土厚靶進一步探討彈體形狀和直徑對于侵徹阻應(yīng)力的影響規(guī)律。

3.1 CRH影響分析

選擇CRH分別為3、4和6的60 mm直徑剛性尖卵形頭部彈體，分別以400 m/s、500 m/s、600 m/s、700 m/s和800 m/s恒定速度侵徹LDPM描述的混凝土厚靶。彈體在混凝土厚靶內(nèi)部侵徹時，LDPM仿真得到彈體周圍骨料顆粒的速度場如圖9所示。其中左側(cè)圖顯示：只有彈頭附近混凝土骨料顆粒速度較大，且混凝土顆粒擾動區(qū)域與未擾動區(qū)域交接面呈現(xiàn)比彈頭形狀更加鈍的錐形；而右側(cè)放大圖顯示速度較大骨料集中在彈體頭部幾乎等厚度的一層區(qū)域。

從仿真結(jié)果提取彈體受力，圖10(a)～圖10(c)分別給出了CRH為3、4、6時，400 m/s、600 m/s、和800 m/s恒定速度侵徹受到的阻力隨著侵徹行程的變化曲線。由于彈體與混凝土骨料隔離之間罰函數(shù)接觸關(guān)系的設(shè)定，侵徹過程中彈體與隨機分布骨料的接觸通過判別關(guān)系不斷更新，因此彈體受到阻力的行程曲線呈現(xiàn)很大的跳躍性。侵徹開始階段，由于彈體頭部沒有完全進入靶體內(nèi)部，侵徹阻力呈現(xiàn)上升階段。等彈體頭部完全進入靶體后，阻力趨于跳躍式穩(wěn)定階段。數(shù)值模擬表明：侵徹阻力隨著彈體侵徹速度增大而增大，隨著CRH的增大而減?。磺褻RH=6時，侵徹阻力跳躍性明顯地降低，這是因為彈體頭部更尖導(dǎo)致接觸面積增大，彈體受力更加均勻。

圖9 深侵徹靶體內(nèi)部骨料顆粒速度分布 Fig.9 Velocity distribution of particles inside the concrete target

圖10 不同彈頭CRH彈體深侵徹阻力情況 Fig.10 Penetration velocity dependent resistive force for projectiles with different CRH

Forrestal等[20]認為彈體在穩(wěn)定侵徹階段侵徹阻力可以表示為

F=πa2(R+NρV2)

(4)

本文采用公式(4)中Forrestal侵徹阻力形式對LDPM仿真獲得的侵徹阻力穩(wěn)定階段進行分析，對比不同彈頭形狀彈體以不同速度侵徹工況下，靶體靜態(tài)阻應(yīng)力R的變化情況。如圖11所示，60 mm直徑彈體侵徹混凝土厚靶分析發(fā)現(xiàn)，靶體靜態(tài)阻應(yīng)力R不隨侵徹速度或尖卵形彈體頭部形狀變化而改變，這與空腔膨脹分析理論[22]相符。

圖11 不同彈頭CRH的60 mm直徑彈體侵徹 靶體靜態(tài)阻應(yīng)力項 Fig.11 Target static resistancefor various CRHs with 60 mm projectile diameter

LDPM仿真得到的靶體內(nèi)一部分侵徹隧道如圖12所示，其左右分別為CRH為3和4的彈體侵徹混凝土厚靶在形成的裂紋分布。侵徹隧道直徑基本與彈體直徑一致，這與深侵徹實驗觀察到的現(xiàn)象一致。

圖12 靶體內(nèi)部侵徹隧道 Fig.12 Terminalballistic tunnel inside concrete target

3.2 彈徑尺寸影響分析

類似地，結(jié)合式(4)分析穩(wěn)定侵徹段阻力，得到靶體靜態(tài)阻應(yīng)力R變化規(guī)律如圖15所示，靶體靜態(tài)阻應(yīng)力R明顯地隨著彈體直徑增大而減小。根據(jù)侵徹試驗中描述，本研究中混凝土靶體采用的骨料最大直徑約為10 mm。彈體直徑為最大骨料直徑3倍、6倍和8倍時，LDPM仿真得到侵徹23 MPa強度混凝土的靶體靜態(tài)阻應(yīng)力值分別為260 MPa、175 MPa、163 MPa。其中80 mm直徑彈體侵徹對應(yīng)的靶體靜態(tài)阻應(yīng)力R與Forrestal論文中分析結(jié)果165 MPa非常接近。最近，劉志林等[7]通過二維離散單元模擬認為彈徑為骨料最大直徑10倍以上時，彈徑尺寸效應(yīng)幾乎不變，這與本文得到的阻應(yīng)力彈徑尺寸效應(yīng)規(guī)律趨勢相接近。

圖13 30 mm彈徑恒定速度 侵徹阻力行程曲線 Fig.13 Perforation resistance curves for 30 mm diameter projectile

圖14 80 mm彈徑恒定速度 侵徹阻力行程曲線 Fig.14 Perforation resistance curves for 80 mm diameter projectile

圖15 不同彈徑時靶體靜態(tài)阻應(yīng)力 項與侵徹速度關(guān)系 Fig.15 Target static resistance vs. penetration velocity for different diameter projectile

4 結(jié) 論

基于LDPM細觀離散元模型，本研究建立了剛性彈侵徹素混凝土厚靶的數(shù)值計算模型，驗證了數(shù)值模型模擬混凝土侵徹問題的準確性。對不同工況混凝土侵徹仿真分析發(fā)現(xiàn)：

(1) 通過三軸壓縮響應(yīng)實驗標定LDPM參數(shù)，LDPM能夠較準確地預(yù)測侵徹23 MPa強度混凝土厚靶問題中侵徹深度以及彈體過載時程曲線，并能模擬得到靶前開坑區(qū)域。

(2) 在400～800 m/s撞擊速度范圍內(nèi)，混凝土厚靶侵徹分析中靶體靜態(tài)阻應(yīng)力與彈體侵徹速度和尖卵形彈體頭部形狀均無關(guān)。彈體頭部越尖，侵徹阻力越均勻。

(3) 彈體直徑在8倍最大骨料直徑范圍內(nèi)，靶體靜態(tài)阻應(yīng)力隨著彈體直徑的增大而減小。本文研究的混凝土厚靶受到彈體400～800 m/s速度范圍內(nèi)撞擊，彈徑為最大骨料直徑3倍、6倍和8倍的彈體受到的靶體靜態(tài)阻應(yīng)力分別為260 MPa、175 MPa和163 MPa。

[ 1 ] 張雄, 廉艷平, 劉巖, 等. 物質(zhì)點法[M]. 北京: 清華大學(xué)出版社, 2013.

[ 2 ] DURBAN D, MASRI R. Dynamic spherical cavity expansion in a pressure sensitive elastoplastic medium[J]. International Journal of Solids & Structures, 2004, 41(41): 5697-5716.

[ 3 ] 劉志林，孫巍巍，王曉鳴，等. 基于蓋帽模型的混凝土動態(tài)球形空腔膨脹模型和侵徹阻力分析[J]. 兵工學(xué)報，2015, 36(12): 2209-2216.

LIU Zhilin, SUN Weiwei, WANG Xiaoming, et al. Spherical cavity-expansion model for concrete targets based on Cap model and penetration resistance analysis [J]. Acta Armamentarii, 2015, 36(12): 2209-2216.

[ 4 ] FENG J, LI W B, PAN G W, et al. Cavitation analysis of spherical shock wave evolution in concrete medium [J]. Acta Mechanica, 2017, 228: 401-414.

[ 5 ] KIRANE K, SU Y, BA?ANT Z P. Strain-rate-dependent microplane model for high-rate comminution of concrete under impact based on kinetic energy release theory [C]//Proc. R. Soc. A. The Royal Society, 2015, 471(2182): 20150535.

[ 6 ] 顧鑫, 章青, 黃丹. 基于近場動力學(xué)方法的混凝土板侵徹問題研究 [J]. 振動與沖擊, 2016, 35(6): 52-58.

GU Xin, ZHANG Qing, HUANG Dan. Peridynamics used in solving penetration problem of concrete slabs [J]. Journal of Vibration and Shock, 2016, 35(6): 52-58.

[ 7 ] 劉志林，孫巍巍，王曉鳴. 基于顆粒流離散元模型的彈丸侵徹細觀混凝土數(shù)值模擬方法研究[J]. 振動與沖擊，2016，35(4): 162-169.

LIU Zhilin, SUN Weiwei, WANG Xiaoming. Numerical simulation for projectile penetrating meso-scale concrete based on particle flow discrete element model [J]. Journal of Vibration and Shock, 2016, 35(4): 162-169.

[ 8 ] SHIU W, DONZé F V, DAUDEVILLE L. Penetration prediction of missiles with different nose shapes by the discrete element numerical approach[J]. Computers & Structures, 2008, 86(21): 2079-2086.

[ 9 ] 陳俊, 黃曉明. 基于三維離散元法的瀝青混合料斷裂過程模擬[J]. 華南理工大學(xué)學(xué)報(自然科學(xué)版), 2012, 40(7): 21-26.

CHEN Jun, HUANG Xiaoming. Simulation of fracture process of asphalt mixture using three-dimension discrete element method[J]. Journal of South China University of Technology (Nature Science Edition), 2012, 40(7): 21-26.

[10] 楊利福, 常曉林, 周偉,等. 基于變形離散元法的重力壩地震開裂分析[J]. 振動與沖擊, 2016, 35(7): 49-55.

YANG Lifu, CHANG Xiaolin, ZHOU Wei, et al. Seismic cracking analysis of a gravity dam based on deformable discrete element method[J]. Journal of Vibration and Shock, 2016, 35(7): 49-55.

[11] CUSATIS G, PELESSONE D, MENCARELLI A. Lattice discrete particle model (LDPM) for failure behavior of concrete. I: Theory [J]. Cement and Concrete Composites, 2011, 33(9): 881-890.

[12] CUSATIS G, MENCARELLI A, PELESSONE D, et al. Lattice discrete particle model (LDPM) for failure behavior of concrete. II: Calibration and validation [J]. Cement and Concrete Composites, 2011, 33(9): 891-905.

[13] SMITH J, CUSATIS G, PELESSONE D, et al. Discrete modeling of ultra-high-performance concrete with application to projectile penetration [J]. International Journal of Impact Engineering, 2014, 65: 13-32.

[14] FENG J, SUN W, LIU Z, et al. An armour-piercing projectile penetration in a double-layered target of ultra-high-performance fiber reinforced concrete and armour steel: Experimental and numerical analyses [J]. Materials & Design, 2016, 102: 131-141.

[15] WAN L, WENDNER R, LIANG B, et al. Analysis of the behavior of ultra high performance concrete at early age [J]. Cement and Concrete Composites, 2016, 74: 120-135.

[16] CUSATIS G, SCHAUFFERT E A. Discontinuous cell method (DCM) for cohesive fracture propagation [C]//Proceedings of the 7th International Conference on Fracture Mechanics of Concrete and Concrete Structures (FraMCos 7), Korea Concrete Institute, Jeju, South Korea, 2010: 23-28.

[17] REZAKHANI R, CUSATIS G. Asymptotic expansion homogenization of discrete fine-scale models with rotational degrees of freedom for the simulation of quasi-brittle materials [J]. Journal of the Mechanics and Physics of Solids, 2016, 88: 320-34.

[18] CUSATIS G. Strain-rate effects on concrete behavior [J]. International Journal of Impact Engineering, 2011, 38(4): 162-170.

[19] WARREN T L, FOSSUM A F, FREW D J. Penetration into low-strength (23 MPa) concrete: target characterization and simulations [J]. International Journal of Impact Engineering, 2004, 30(5): 477-503.

[20] FORRESTAL M J, FREW D J, HICKERSON J P, et al. Penetration of concrete targets with deceleration-time measurements [J]. International Journal of Impact Engineering, 2003, 28(5): 479-497.

[21] FREW D J, HANCHAK S J, GREEN M L, et al. Penetration of concrete targets with ogive-nose steel rods [J]. International Journal of Impact Engineering, 1998, 21(6): 489-497.

[22] FENG J, LI W, WANG X, et al. Dynamic spherical cavity expansion analysis of rate-dependent concrete material with scale effect [J]. International Journal of Impact Engineering, 2015, 84: 24-37.