基于孔隙演化的砂土沖擊絕熱關系研究

高 飛， 邱艷宇， 王明洋， 張先鋒， 程怡豪

(1. 南京理工大學 智能彈藥國防重點學科實驗室， 南京 210094; 2. 解放軍理工大學 爆炸沖擊防災減災國家重點實驗室， 南京 210007)

基于孔隙演化的砂土沖擊絕熱關系研究

高 飛1，2， 邱艷宇1，2， 王明洋1，2， 張先鋒1， 程怡豪2

(1. 南京理工大學 智能彈藥國防重點學科實驗室， 南京 210094; 2. 解放軍理工大學 爆炸沖擊防災減災國家重點實驗室， 南京 210007)

基于砂土多孔、含水的特性，對動載荷作用下孔隙演化對砂土沖擊絕熱關系的影響進行了研究。基于砂土基體的不可壓縮假設，采用單個球形氣孔等效模型和廣義Mises屈服準則推導了砂土考慮孔隙壓實的演化方程；根據沖擊Hugoniot突躍條件和Grüneisen型方程，獲得了砂土考慮孔隙演化的狀態方程；根據混合物沖擊波關系和孔隙演化狀態方程，獲得了砂土基體材料沖擊絕熱關系、干砂土和濕砂土沖擊絕熱關系，并與已有實驗結果進行了對比。結果表明：計算得到的絕熱線與實驗數據吻合較好，該孔隙演化狀態方程能夠較為準確的反映動載作用下砂土的動態響應過程。

砂土； 孔隙演化； Grüneisen型狀態方程； 沖擊絕熱關系

砂土是自然界普遍存在的地質材料，是由砂土顆粒、水和空氣構成的多組分、多相、非均勻混合物。砂土及其它多孔顆粒介質動態力學行為的研究對許多工程應用，如地震工程、侵徹研究、陶瓷裝甲性能研究等都有重要的參考價值[1]。目前，國內對于砂土的研究主要集中在強度理論、剪脹理論及本構模型方面[2-6]，其動態力學行為的研究較少涉及。砂土動態響應過程很大程度上受顆粒尺寸、孔隙率及含水量的影響，尤其是含水量及孔隙率的不同，對砂土中沖擊波的傳播規律有著顯著的影響。

研究表明[7]，含水量(濕度為3%～20%)對砂土的動態應力應變關系影響很大，濕砂比干砂具有更好的可壓縮性，且濕度為7%的濕砂表現的最軟(可壓縮性最好)，并不是含水量越高的濕砂可壓縮性越好。Proud等[8]研究在不同應變率范圍內含水量對砂力學響應的影響規律，結果表明含水量在不同應變率區間的作用機理十分復雜，在低應變率時水起到潤滑作用，而在高應變率時水的介入可能會出現砂越來越難壓的現象。Brown等[1]通過飛片沖擊實驗，得到了砂的Hugoniot關系，修正了經典P-α模型。Borg等[9]研究沖擊壓縮條件下壓實模型對于描述不同初始孔隙率砂的動力學響應時，發現P-α型壓實模型與Grüneisen型狀態方程的組合所得到的壓力時程曲線能夠更好地再現實驗結果。

本文將微觀不均勻多孔介質等效為均勻介質處理，將砂土內部分散的孔隙等效為僅包含單個氣孔的球形顆粒介質，基于廣義Mises準則推導了另一種形式的孔隙壓實P-α方程，建立了考慮孔隙壓實效應的多孔材料狀態方程，計算了砂土基體材料和多孔材料的沖擊絕熱關系，與已有的Hugoniot實驗數據對比，證明了孔隙演化方程和狀態方程的合理性。該簡化模型與研究思路對深入研究砂土及其它多孔介質動態力學行為具有參考意義。

1 基本模型與孔隙演化方程

1.1基本假設與等效模型

對砂土這種非均勻多孔介質，把它看作是水和石英組成的雙組分混合物，基體是水和石英，夾雜物是氣孔。孔隙率變化對砂土力學性能的影響不容忽視，為此有必要建立考慮孔隙動態演化的方程。在推導這個方程時建立了如下模型，認為含大量初始孔隙的砂土與內部僅包含單個氣孔的球形顆粒的表現是相似的[10]，如圖1所示，氣孔周圍是砂土基體材料，模型內、外半徑分別為a,b，相應內、外半徑初始值分別為a0、b0，且孔隙體積、基體體積均與原物質相同，這就保證了初始孔隙率的不變。距離平面中心b處作用有徑向壓力pm。在推導孔隙動態演化方程過程中假設如下：

(1) 基體材料是均勻和各向同性的、不可壓且滿足廣義Mises屈服準則；

(2) 基體中氣孔的分布在各向均勻；

(3) 忽略氣孔的表面能和氣孔中氣體的壓力，假定基體材料與多孔材料的內能是恒等的(亦即多孔材料的內能是由基體的應變決定的)；

圖1 砂土的孔隙演化等效模型

(4) 在分析氣孔周圍材料處于塑性狀態時，忽略了氣孔的彈性階段微小壓縮量。

1.2孔隙演化動力學方程

對多孔材料用參數α表示其密實程度：

(1)

式中：v為多孔介質的總體積；vm為基體的體積；ρ為多孔介質的密度；ρm為基體的密度；當α=1時材料達到完全密實狀態。

如果忽略氣孔的表面能和氣孔中氣體的壓力，則多孔材料的狀態方程與基體材料具有相同的形式，而多孔材料中的壓力P與基體材料中的壓力Pm有如下關系式[11]：

(2)

引入中心在平面上的球坐標系(r,θ,φ)，其即時半徑用a表示。設距離平面中心b處作用有應力：

σr=-Pm=-αP

(3)

式中：

(4)

假定基體材料是不可壓縮的，則球形坐標下的控制方程組可表示為

(5)

(6)

σr-σθ=Y0+kPm=Y0-k(σr+2σθ)/3

(7)

式中，ur為徑向速度;σr、σθ分別為徑向應力和環向應力;Y0、k分別為黏結力和摩擦因數。

含氣孔的球形顆粒內外邊界的邊界條件是：

σr(b,t)=-αP,σr(a,t)=0

(8)

對式(6)積分，可得：

(9)

式中：

(10)

利用邊界點(a,a0)、(b,b0)，可得：

(11)

由式(4)，可知：

(12)

即可確定函數：

(13)

并得到加速度：

(14)

這樣，考慮到式(7)，方程(5)可表示為

(15)

求解該方程，可得解：

(16)

分別利用式(8)的兩個邊界條件并考慮到下式：

(17)

最終可得：

(18)

式中：

(19a)

(19b)

(19c)

式(18)在下列條件時(氣孔周圍材料已進入塑性狀態，此時孔隙的變化可觀察到)可用于確定α：

(20)

(21)

圖2是根據方程式(21)得到的壓力-孔隙率變化時砂土的加載-卸載曲線圖，是另一種形式的P-α模型。在1-2初始加載階段，砂土壓縮時孔隙率α0沒有變化，因為氣孔收縮的彈性和彈塑性分量很小，忽略不計。當壓力按照方程式(21)的變化規律繼續增大時(2-3段)，砂土顆粒材料產生了不可逆地壓實 (α逐漸減小逼近1)。卸載階段，材料在孔隙率恒定的條件下，由點3卸載至點4 (P=0)。

對于孔隙介質沖擊壓縮曲線，也可以用Herrmann[12]提出的P-α模型描述壓力與孔隙率的關系，分彈性和塑性階段，如圖3所示。當壓力小于彈性極限壓力Pe時，材料處于彈性階段，且孔隙率的變化是可逆的，壓力與孔隙率有如下關系[13]：

圖2 砂土壓力與孔隙率的加卸載演化曲線

圖3 典型的孔隙介質壓力與孔隙率變化曲線[12]

(22a)

(22b)

(22c)

αe=α(P=Pe)

(22d)

利用式(22)可以計算由初始孔隙率α0到彈性極限孔隙率αe的微小變化，由于孔隙率變化很小，一般假設在彈性階段滿足：dα/dP=0，即認為ce=cm0。

在P>Pe的塑性階段，通常用多項式的形式表示壓力與孔隙率的關系[14]：

(23)

式中：Pe和Ps均為材料常數，分別表示砂土剛進入塑性階段和壓密時對應的壓力值。

在忽略彈性階段孔隙率的微小變化后，式(21)表示的P-α關系曲線與Herrmann提出的經典P-α模型在砂土變形的塑性階段相似，僅僅是函數形式和模型參數不同。

2 砂土的狀態方程與沖擊絕熱關系

本節首先建立孔隙演化砂土Grüneisen型狀態方程；然后利用混合物沖擊波關系，求砂土基體的沖擊Hugoniot線；最后根據砂土狀態方程求干砂和濕砂的沖擊Hugoniot線，并與實驗數據對比。

2.1基于孔隙演化的砂土狀態方程

因為砂土中夾雜空氣的質量與砂土顆粒相比小幾個數量級，故氣孔內空氣的影響可以忽略不計，因此干砂就可以看作以石英為基體的多孔材料，夾雜物是氣孔；濕砂則看作以石英和水為基體的簡單二組分混合物，夾雜物是氣孔。利用基體的Hugoniot參數和Grüneisen型狀態方程，便可得到砂土的狀態方程。若利用基體材料沖擊波速Dm與質點速度Um的簡單線性關系：

Dm=cm0+qm0Um

(24)

根據沖擊波波陣面的突躍條件和Grüneisen型方程，基體材料的狀態方程可表示為

Pm=PH+γmρm(Em-EH)

(25)

(26)

式中:PH為沖擊絕熱線上的壓力;γm為與ρm對應的Grüneisen系數;ρm0為基體初始密度;EH為沖擊波波陣面上基體的比內能;Em為與ρm、PH對應的基體比內能。

根據沖擊波波陣面上的守恒條件：

(27)

對砂土而言，假設條件與2.1節相同，忽略氣孔的表面能和氣孔中氣體的壓力，假定基體材料與多孔材料的內能相等(即Em=E(ρ，P)，E為砂土的比內能)，將式(27)、式(26)、式(1)、式(2)代入式(25)，并利用Grüneisen假定關系：γmρm=γm0ρm0，則壓縮時砂土多孔材料的狀態方程為

(28)

(29)

2.2砂土基體材料沖擊絕熱關系

根據混合物沖擊波關系式：

(30)

將式(30)代入式(24)，可得：

(31)

對于砂土混合物的基體(即暫不考慮孔隙)而言，以(νm，Pm)為變量的沖擊絕熱關系可表示為

(32)

式中：m1=ω0/(1+ω0)，m2=1/(1+ω0)分別為砂土中水和石英的質量分數(m1+m2= 1)，ω0為砂土的質量濕度；νi0、ρi0(i=1, 2)分別為水、石英的初始比容和密度。計算用到的水、石英及其混合物線性關系Hugoniot參數見表1。

表1 水、石英和兩者混合物的線性關系Hugoniot參數[17]

圖4是以(νm，Pm)和(Um，Dm)為變量的砂土基體材料沖擊絕熱曲線，水的絕熱曲線是1，石英的絕熱曲線是2，密度ρ0=1.8 g/cm3，濕度為28%砂土基體混合物(水和石英)的絕熱曲線是3。由圖4 (a)的結果得出，在100 GPa壓力范圍內以混合物沖擊波關系和線性關系表示的Pm-νm線，水由初始比容1 cm3/g壓縮到了0.45 cm3/g，壓縮量為55%；石英由初始比容0.376 cm3/g壓縮到了0.25 cm3/g，壓縮量為33.5%；水和石英的基體混合物由初始比容0.513 cm3/g壓縮到了0.3 cm3/g，壓縮量為41.5%。基體混合物的Hugoniot線界于水和石英的之間，這也是符合混合物理論基本規律的。

(a) 比容與壓力的關系

(b) 質點速度與波速的關系

由圖4 (b)可以看出，由式(30)表示基體沖擊波速度與質點速度關系時，就某單個組分而言，水或石英的Dm-Um線的線性關系十分好，且在質點速度為0 km/s時對應的沖擊波速度(即體積聲速c0)也與表1的值吻合很好。水和石英混合物的Dm-Um線可以分為兩段：在低壓段(當物體具有彈性力學行為時)，Dm-Um線在Um=1.15 km/s 時向下偏折，主要原因是，在建立固體高壓狀態方程時一般不計剪切強度的影響，這種假設在高壓階段較為合理，但低壓段剪切強度已逐漸恢復，若仍不考慮剪切強度將會低估壓力值(或波速)的大小，發生曲線向下偏折的現象，通常理論計算得到的波速-質點速度關系都是如此[15]；在高壓段(當物體具有流體力學行為時)，即Um>1.15 km/s 時，Dm-Um線的線性規律還是比較好的。在不考慮局部線性偏離的情況下，在相當寬廣的壓力范圍內，Dm-Um線基本能夠表征砂土基體材料的沖擊壓縮特性。

目前對于固體在高壓作用下的動力性質研究較多，而對于低應力區與高應力區之間的過渡區的研究工作相對較少，這一過渡區的范圍對于強度大的物體來講(金屬、巖石及混凝土)大約從數百兆帕到數十吉帕，在這一壓力范圍內實現了從彈性狀態到內摩擦狀態，甚至到流體動力學狀態的轉換，十分復雜[16]。由公式P=ρDmUm，計算得到偏折點處對應的壓力值約10 GPa，該點恰好處于這一應力過渡區，這種應力狀態的轉換也解釋了發生偏折現象的物理本質。

為分析初始孔隙率和濕度對砂土Hugoniot線的影響，根據方程(21)和表2參數計算了干砂土(α0=2.062)和濕砂土(α0=1.083，ω0=28%)的P-α曲線，如圖5。圖5(a)為干砂土的P-α曲線，發現干砂土達到密實態(α=1)所需的壓力值為20 GPa；而圖5(b)濕砂土達到密實態的壓力值為5 GPa，可能是因為水起到了潤滑作用，砂粒間摩擦力減小，使原來松散多孔介質更容易向密實態轉變。可見，同樣是砂土，當初始孔隙率和濕度不同時，材料沖擊壓縮響應會有較大差異，進一步明晰了孔隙率和濕度在研究砂土類多孔介質時的重要性。

(a) 干砂土孔隙演化曲線

(b) 濕砂土孔隙演化曲線

材料ρ0/(g·cm-3)ω0/%ν空氣/%ρm0/(g·cm-3)cm0/(m·s-1)qm0γm0Y0/GPak干砂土1.29051.52.6636802.121.00.10.75濕砂土1.80287.751.9524501.861.280.0250.5

2.3砂土材料沖擊絕熱關系

根據2.1節的砂土多孔材料狀態方程和表2參數計算了以(ν，P)和(U，D)為變量的干砂土和濕砂土的Hugoniot線，如圖6所示，實線和虛線代表計算值，實心黑點和空心圓圈代表實驗值[17]。

圖6 (a)給出了12 GPa壓力范圍內計算得到的干砂土和濕砂土P-νHugoniot線，與實驗結果吻合較好，在高壓段的計算值與實驗值的偏差比低壓段稍大。干砂土由初始比容0.775 cm3/g壓縮到了0.334 cm3/g，壓縮量為57%；濕砂土由初始比容0.555 cm3/g壓縮到了0.382 cm3/g，壓縮量為31.2%。值得注意的是，干砂土和濕砂土的P-ν線在ν=0.5 cm3/g時交叉，在ν=0.5-0.775 cm3/g 時干砂土的P-νHugoniot線幾乎為水平線，主要是因為干砂土的初始孔隙率(α0=2.062)大，而壓縮氣孔所需的壓力值很小，故此階段壓力變化不明顯，當ν=0.5 cm3/g (對應α=1.33)時，氣孔幾乎完全崩塌而逐漸開始壓縮固體顆粒，故壓力值陡增；對于濕度為28%的濕砂土，由于水的存在填充了大部分孔隙，使得材料在初始階段(α0=1.083)就已經接近密實態，故在一開始比容的微弱變化就使壓力明顯增大。對比圖6 (a)干砂土和濕砂土P-νHugoniot線上聯結終態點與初態點的弦線(即Rayleigh線)斜率，發現濕砂Rayleigh線的斜率絕對值明顯大于干砂Rayleigh線的斜率絕對值，表明相同沖擊條件下濕砂中沖擊波物質波速大于干砂中沖擊波物質波速，主要原因是在干砂土中，存在相當數量的氣孔，氣孔的存在使得沖擊波的能量很大程度上耗散在氣孔變形與壓潰，進而影響了沖擊波的傳播。而水的介入，一方面增加了濕度，另一方面也減小了孔隙率，起到填充孔隙均衡材料的作用，使得沖擊波速度有著顯著的提高。

(a) 比容與壓力曲線

(b) 質點速度與波速曲線

Fig.6 Comparison of computation and experiment (adapted from Belov et al.[17]) shock adiabatic curves for dry sand and moist sand

圖6 (b)給出了計算得到的干砂土和濕砂土的D-U線，計算值與實驗結果吻合較好，需要注意的是，在低壓區干砂土和濕砂土的D-U線發生了向上偏折的現象，但整體而言還是符合基本規律的。Herrmann[12]和Белов[18]在利用Grüneisen型狀態方程計算多孔金屬鐵及Al2O3多孔陶瓷的D-U線時，低壓段都發生了向上偏折的現象，但均未給出合理的解釋。狀態方程式(28)在描述物質熱力學特性時準確性是有限的，這主要受限于參數γ取常數，實際上Grüneisen系數取決于密度和溫度。此外，這個狀態方程可以成功的描述在高密度和高溫度下物質的熱力屬性，不能很好的描述密度很低時的壓力特性，這也是低壓區發生偏折的原因。總體而言，本文基于孔隙演化推導的砂土狀態方程在描述其沖擊Hugoniot線時是可行的，無論是P-ν線還是D-U線，計算值與實驗值吻合都較好。

3 結 論

基于單個球形氣孔顆粒模型和廣義Mises屈服準則，推導了孔隙壓縮演化方程，建立了砂土考慮孔隙演化的狀態方程，計算了砂土基體、干砂土、濕砂土的Hugoniot關系，與實驗結果吻合較好。主要結論如下：

(1) 基于相關假設推導的孔隙演化動力學方程在忽略氣孔彈性階段的微小變化后，與經典P-α模型在變形的塑性階段相似，僅僅是函數形式和模型參數不同。

(2) 由孔隙演化動力學方程的計算結果可知，干砂土和濕砂土達到密實態所需的壓力值差別很大，在研究砂土類多孔介質沖擊壓縮特性時要重點關注孔隙率和濕度這兩個特征量。

(3) 根據孔隙演化狀態方程計算的干砂土和濕砂土Hugoniot線與實驗結果吻合較好，證明了簡化模型與假設的合理性，表明狀態方程可用于砂土的沖擊壓縮計算；由于數據的不足，該模型是否適用于其他多孔介質的計算還需要進一步實驗驗證。

(4) 建立的狀態方程在計算D-U線時，在低壓過渡段發生了向上偏折現象，并給出了初步的解釋；低應力區與高應力區之間的過渡區的研究工作相對較少，有必要在未來的研究中加以重視。

[1] BROWN J L, VOGLER T J, GRADY D E, et al. Dynamic compaction of sand[C]// American Institute of Physics Conference Proceedings. American Institute of Physics, 2007, 1363-1366.

[2] 周健, 池永. 砂土力學性質的細觀模擬[J]. 巖土力學, 2003, 24(6): 901-906.

ZHOU Jian, CHI Yong. Mesomechanical simulation of sand mechanical properties[J]. Rock and Soil Mechanics, 2003, 24(6): 901-906.

[3] 蔡正銀, 李相菘. 砂土的剪脹理論及其本構模型的發展[J]. 巖土工程學報, 2007, 29(8): 1122-1128.

CAI Zhengyin, LI Xiangsong. Development of dilatancy theory and constitutive model of sand[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(8): 1122-1128.

[4] 常在, 楊軍, 程曉輝. 砂土強度和剪脹性的顆粒力學分析[J]. 工程力學, 2010, 27(4): 95-104.

CHANG Zai, YANG Jun, CHENG Xiaohui. Granular mechanical analysis of the strength and dilatancy of sands[J]. Engineering Mechanics, 2010, 27(4): 95-104.

[5] 李學豐, 黃茂松, 錢建固. 宏細觀結合的砂土各向異性破壞準則[J]. 巖石力學與工程學報, 2010, 29(9): 1885-1892.

LI Xuefeng, HUANG Maosong, QIAN Jiangu. Failure criterion of anisotropic sand with method of macro-meso incorporation[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(9): 1885-1892.

[6] 張建民. 砂土動力學若干基本理論探究[J]. 巖土工程學報, 2012, 34(1): 1-50.

ZHANG Jianmin. New advances in basic theories of sand dynamics[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(1): 1-50.

[7] MARTIN B E, CHEN Weinong, SONG Bo, et al. Moisture effects on the high strain-rate behavior of sand[J]. Mechanics of Materials, 2009, 41: 786-798.

[8] PROUD W G, CHAPMAN D J, WILLIAMSON D M, et al. The dynamic compaction of sand and related porous systems[C]// American Institute of Physics Conference Proceedings. American Institute of Physics, 2007, 1403-1408.

[9] BORG J P, COGAR J R, LLOYD A, et al. Computational simulations of the dynamic compaction of porous media[J]. International Journal of Impact Engineering, 2006, 33: 109-118.

[10] CARROLL M M, HOLT A C. Static and dynamic porecollapse relations for ductile porous materials[J]. Journal of Applied Physics, 1972, 43(4): 1626-1636.

[11] CARROLL M M, HOLT A C. Suggested modification of the Pα model for porous materials[J]. Journal of Applied Physics, 1972, 43(2): 759-761.

[12] HERRMANN W. Constitutive equation for the dynamic compaction of ductile porous materials[J]. Journal of Applied Physics, 1969, 40(6): 2490-2499.

[13] JUTZI M, BENZ W, MICHEL P. Numerical simulations of impacts involving porous bodies: I. Implementing sub-resolution porosity in a 3D SPH Hydrocode[J]. ICARUS, 2008, 198: 242-255.

[14] BORG J P, MORRISSEY M P, PERICH C A, et al. In situ velocity and stress characterization of a projectile penetrating a sand target: experimental measurements and continuum simulations[J]. International Journal of Impact Engineering, 2013, 51: 23-35.

[15] 經福謙. 實驗物態方程導引[M]. 北京: 科學出版社, 1986.

[16] 戚承志, 錢七虎. 巖體動力變形與破壞的基本問題[M]. 北京: 科學出版社, 2009.

[17] BELOV N N, KOPANITSA D G, KUMPYAK O G, et al. Calculation of reinforced-concrete structures for explosive and impact loads[M]. (in Russian), STT, Northampton, Tomsk, 2004.

[18] БЕЛОВ Н Н, ЮГОВ Н Т, АФАНАСЬЕВА С А, et al. Исследование процессов деформирования и разрушения хрупких материалов[J]. Механика Композиционных Материалов и Конструкций, 2001, 7(2): 131-142.

Shock-hugoniotrelationshipsofsandconsideringporosityevolution

GAO Fei1, 2, QIU Yanyu1, 2, WANG Mingyang1, 2, ZHANG Xianfeng1, CHENG Yihao2

(1. Ministerial Key Laboratory of ZNDY，Nanjing University of Science and Technology, Nanjing 210094, China;2. State Key Lab for Disaster Prevention & Mitigation of Explosion & Impact, PLA University of Science and Technology, Nanjing 210007, China)

Effects of porosity evolution on sand shock-hugoniot relationships under dynamic loading were investigated based on porous and hydrous characteristics of sand. Based on the incompressible hypothesis of sand matrix, using a single spherical pore equivalent model and the generalized Mises strength criterion, the porosity evolution equation of sand was derived. According to Hugoniot jump condition and Grüneisen-type state equation, the state equation of sand considering porosity dynamic compaction was deduced. The shock-hugonoit curves for sand matrix, dry sand and moist sand were gained by using the mixture-shock wave relationship and porosity evolution state equation. Results showed that the calculated shock-hugonoit curves for sand agree well with the available test data in literature; the porosity evolution state equation can more correctly reflect the dynamic response process of sand under dynamic loading.

sand; porosity evolution; Grüneisen-type equation of state; shock-hugoniot relationships

國家自然科學基金面上項目(51478466)；國家自然科學基金青年項目(51508568)；爆炸科學與技術國家重點實驗室開放基金項目(KFJJ15-07M)

2016-01-26 修改稿收到日期：2016-07-28

高飛 男，博士生，1990年生

王明洋 男，教授，博士生導師，1966年生

E-mail:wmyrf@163.com

O383

： A

10.13465/j.cnki.jvs.2017.17.021