張靜靜 董龍雷

摘要：針對超高速?zèng)_擊使結(jié)構(gòu)瞬間產(chǎn)生極高的溫度并造成結(jié)構(gòu)軟化現(xiàn)象的問題，提出考慮沖擊溫度影響的改進(jìn)物質(zhì)點(diǎn)法算法，采用Taylor桿的撞擊過程模擬分析這一現(xiàn)象。將模擬結(jié)果與經(jīng)典物質(zhì)點(diǎn)法算法對同一算例的計(jì)算結(jié)果進(jìn)行對比，結(jié)果顯示改進(jìn)物質(zhì)點(diǎn)法算法計(jì)算結(jié)果精度可提高0.77%。

關(guān)鍵詞：物質(zhì)點(diǎn)法；塑性變形；沖擊溫度；超高速?zèng)_擊

中圖分類號：TB302.3；TB115.1

文獻(xiàn)標(biāo)志碼：B

0 引 言

越來越多的學(xué)者研究有限元模擬大變形過程遇到的網(wǎng)格畸變問題。物質(zhì)點(diǎn)法是1994年由SULSKY等[1-2]提出的，該算法采用拉格朗日質(zhì)點(diǎn)和歐拉網(wǎng)格雙重描述，是一種完全的拉格朗日質(zhì)點(diǎn)類法。該方法將連續(xù)體離散為一組離散物質(zhì)點(diǎn)，材料的物理信息，如速度、應(yīng)力和位置等，存儲(chǔ)在每一塊材料區(qū)域?qū)?yīng)的物質(zhì)點(diǎn)上；歐拉網(wǎng)格僅用于動(dòng)量方程的求解和空間導(dǎo)數(shù)的計(jì)算，不攜帶任何物質(zhì)信息；在每一個(gè)計(jì)算時(shí)間步中，物質(zhì)點(diǎn)和計(jì)算網(wǎng)格完全固連，不存在相對運(yùn)動(dòng)，避免歐拉法因非線性對流項(xiàng)產(chǎn)生的數(shù)值計(jì)算困難，且易于跟蹤物質(zhì)界面。另外，由于物質(zhì)點(diǎn)攜帶物質(zhì)的所有信息，因此每一步計(jì)算結(jié)束后，將已變形的背景網(wǎng)格丟棄，在新一步的計(jì)算過程中采用新的背景網(wǎng)格，因此物質(zhì)點(diǎn)法可避免大變形問題中網(wǎng)格畸變產(chǎn)生的數(shù)值計(jì)算困難，在涉及大變形（如爆炸、沖擊等）問題上具有明顯優(yōu)勢。

物質(zhì)點(diǎn)法最初采用顯式積分法，適用于載荷作用時(shí)間短的瞬態(tài)問題，如爆炸、沖擊等。很多學(xué)者根據(jù)各自需求對算法進(jìn)行改進(jìn)。目前，物質(zhì)點(diǎn)法已廣泛應(yīng)用于各類工程模擬研究中。SULSKY等[1-2]模擬Taylor桿和鋼球侵徹鋁靶問題；HUANG等[3]和MA等[4]采用物質(zhì)點(diǎn)接觸算法模擬中低速?zèng)_擊侵徹問題。在模擬材料失效方面，CHEN等[5-7]采用物質(zhì)點(diǎn)法模擬沖擊載荷下脆性材料的動(dòng)態(tài)失效問題和材料在局部加熱情況下的失效問題。

溫度是影響材料特性的重要因素之一。在已有文獻(xiàn)中，溫度模擬過程大多只考慮大變形對結(jié)構(gòu)造成溫度升高的影響。在超高速?zèng)_擊問題中，如Taylor桿撞擊試驗(yàn)等，在沖擊絕熱壓縮和沖擊載荷作用下，材料的沖擊波耗散效應(yīng)會(huì)引起結(jié)構(gòu)發(fā)生劇烈的溫度變化，稱為沖擊溫度。然而，經(jīng)典的物質(zhì)點(diǎn)法并沒有將沖擊溫度考慮在內(nèi)，引起模擬結(jié)果誤差，因此本文在物質(zhì)點(diǎn)法中考慮沖擊溫度的影響。

1 物質(zhì)點(diǎn)法

物質(zhì)點(diǎn)法將一個(gè)連續(xù)體離散為一系列物質(zhì)點(diǎn)，見圖1。連續(xù)體的密度可近似表達(dá)為

在網(wǎng)格的計(jì)算方式上，物質(zhì)點(diǎn)法與有限元法非常相近，兩者的區(qū)別[8]在于：

（1）有限元法采用高斯積分，將積分轉(zhuǎn)化為被積函數(shù)在各個(gè)高斯點(diǎn)處的值與該高斯點(diǎn)所代表的體積之積的和。物質(zhì)點(diǎn)法采用物質(zhì)點(diǎn)積分，將積分轉(zhuǎn)化為被積函數(shù)在各物質(zhì)點(diǎn)處的值與該物質(zhì)點(diǎn)所代表的體積之積的和。

（2）有限元法的計(jì)算網(wǎng)格始終與物體固連；物質(zhì)點(diǎn)法的背景網(wǎng)格只在每個(gè)時(shí)間步內(nèi)與物體固連，在每個(gè)時(shí)間步結(jié)束時(shí)，將已變形的背景網(wǎng)格丟棄，并在下一個(gè)時(shí)間步內(nèi)采用新的規(guī)則的背景網(wǎng)格進(jìn)行計(jì)算。由于物質(zhì)點(diǎn)已經(jīng)攜帶物體的所有物質(zhì)信息，因此，在下一個(gè)計(jì)算時(shí)間步內(nèi)，可將物質(zhì)點(diǎn)信息映射到新的背景網(wǎng)格上得到網(wǎng)格信息。因此，物質(zhì)點(diǎn)法在計(jì)算大變形的過程中不會(huì)出現(xiàn)網(wǎng)格畸變的現(xiàn)象。

物質(zhì)點(diǎn)法與無網(wǎng)格方法（以光滑粒子流體動(dòng)力學(xué)方法為例）相比，兩者都需要對物質(zhì)點(diǎn)進(jìn)行離散并且將信息攜帶在各個(gè)物質(zhì)點(diǎn)上，兩者的區(qū)別在于：

（1）光滑粒子流體動(dòng)力學(xué)方法中物質(zhì)點(diǎn)與物質(zhì)點(diǎn)的聯(lián)系通過搜索物質(zhì)點(diǎn)的鄰域進(jìn)行。物質(zhì)點(diǎn)法在將物質(zhì)點(diǎn)信息映射到背景網(wǎng)格時(shí)，已使各物質(zhì)點(diǎn)之間產(chǎn)生聯(lián)系，因此物質(zhì)點(diǎn)法可避免物質(zhì)點(diǎn)鄰域搜索步驟，提高計(jì)算效率。

（2）光滑粒子流體動(dòng)力學(xué)方法通過求解物質(zhì)點(diǎn)組的動(dòng)力學(xué)方程和跟蹤每個(gè)物質(zhì)點(diǎn)的運(yùn)動(dòng)軌道求得整個(gè)系統(tǒng)的力學(xué)行為，即其動(dòng)量方程的求解也是在物質(zhì)點(diǎn)組上實(shí)現(xiàn)的。物質(zhì)點(diǎn)法將物質(zhì)點(diǎn)的質(zhì)量、速度等參數(shù)映射到背景網(wǎng)格上進(jìn)行動(dòng)量方程的求解。

開源的物質(zhì)點(diǎn)法分析結(jié)構(gòu)溫度造成的影響因素僅考慮大變形情況，本文針對刨削產(chǎn)生和發(fā)展過程中產(chǎn)生高溫現(xiàn)象的影響因素，如沖擊、摩擦等，提出對標(biāo)準(zhǔn)物質(zhì)點(diǎn)法的改進(jìn)，并通過經(jīng)典的Taylor桿撞擊試驗(yàn)予以驗(yàn)證。

2 考慮沖擊溫度影響的物質(zhì)點(diǎn)法算法

采用考慮材料的溫度軟化效應(yīng)的Johnson-Cook本構(gòu)模型進(jìn)行驗(yàn)證，其對應(yīng)公式為

式（4）中：第一項(xiàng)A+Bεnp為應(yīng)變的表達(dá)式，A為材料名義屈服強(qiáng)度，B和n反映應(yīng)力硬化的影響程度；第二項(xiàng)和第三項(xiàng)分別反映應(yīng)變率和溫度的影響程度。在應(yīng)變和應(yīng)變率均非常小且材料溫度為室溫的情況下，式（4）可簡化為σy=A。當(dāng)應(yīng)力或應(yīng)變率增大時(shí)，σy增加；當(dāng)材料溫度接近熔融溫度時(shí)，溫度影響項(xiàng)1-（T*）m接近于0，該現(xiàn)象導(dǎo)致材料屈服強(qiáng)度趨于0，結(jié)構(gòu)極易發(fā)生變形。因此，當(dāng)材料溫度趨于熔融溫度時(shí)，結(jié)構(gòu)出現(xiàn)溫度軟化現(xiàn)象，結(jié)構(gòu)屈服強(qiáng)度非常小，容易發(fā)生變形。以Taylor桿撞擊試驗(yàn)[9]為例對該現(xiàn)象進(jìn)行驗(yàn)證。

Taylor桿為銅桿，以190 m/s的速度向剛性墻撞擊，銅桿初始長度L0=25.40 mm。在物質(zhì)點(diǎn)法模擬過程中，物質(zhì)點(diǎn)間距為0.38 mm，Taylor桿共離散為21 172個(gè)物質(zhì)點(diǎn)。模擬時(shí)間80 μs，此時(shí)動(dòng)能為0，基本不再變化，不考慮任何能量產(chǎn)生或損耗因素。Taylor桿的材料參數(shù)見表1。在該試驗(yàn)中，Taylor桿的最終長度L=16.20 mm，底端直徑D=13.50 mm，距離底部0.2L0處直徑W=10.10 mm。

2.1 溫度軟化效應(yīng)驗(yàn)證

模擬材料溫度為30、300、1 000和1 500 K的Taylor桿撞擊，在80 μs時(shí)Taylor桿沿Y軸中心位置的切面見圖2。在不同初始溫度條件下，Taylor桿的長度L、底部直徑D、距離底部0.2L0位置處的直徑W的計(jì)算結(jié)果見表2。

圖2和表2表明，在相同的初始條件和計(jì)算參數(shù)下，隨著材料溫度的升高，結(jié)構(gòu)的變形增大，說明隨著溫度的升高，材料硬度降低，結(jié)構(gòu)越容易發(fā)生變形。刨削的產(chǎn)生和發(fā)展過程伴隨著大量的熱量產(chǎn)生，因此在對其過程進(jìn)行數(shù)值分析時(shí)必須考慮溫度，主要可以從塑性變形和沖擊2個(gè)方面考慮物質(zhì)點(diǎn)法中溫度的影響。

2.2 塑性功

低應(yīng)變率和高應(yīng)變率變形過程的處理往往是不同的：低應(yīng)變率現(xiàn)象常常處理為等溫過程；高應(yīng)變率現(xiàn)象常常處理為絕熱過程。絕熱過程中結(jié)構(gòu)變形所做的功常轉(zhuǎn)換為熱能進(jìn)而導(dǎo)致結(jié)構(gòu)溫度升高。因此，當(dāng)材料發(fā)生高應(yīng)變率變形時(shí)，大部分塑性功會(huì)轉(zhuǎn)換為熱能。

以Taylor桿試驗(yàn)為例，分析材料初始溫度為30 K、考慮塑性變形做功的影響下物質(zhì)點(diǎn)法算法的計(jì)算結(jié)果。Taylor桿沿Y軸中心位置的切面在80 μs時(shí)的結(jié)果見圖3，最終桿長L=16.09 mm，底部直徑D=12.89 mm，距離底部0.2L0位置處的直徑W=9.28 mm。與試驗(yàn)結(jié)果相比，計(jì)算結(jié)果誤差為4.44%。

2.3 沖擊溫度

沖擊溫度是材料在沖擊壓縮下的重要物理參數(shù)。[10]從原理上講，沖擊溫度可通過材料的熱力學(xué)函數(shù)、守恒方程和沖擊絕熱線計(jì)算獲得。

根據(jù)文獻(xiàn)統(tǒng)計(jì)，沖擊溫度的計(jì)算大致可分為3種近似方法：（1）利用Mie-Grüneisen狀態(tài)方程和Hugoniot曲線計(jì)算；（2）利用三項(xiàng)式狀態(tài)方程計(jì)算；（3）利用以等熵線為參考線的Mie-Grüneisen狀態(tài)方程計(jì)算。[11-16]湯文輝等[17]對沖擊溫度近似算法的研究結(jié)果表明：方法（1）涉及的參數(shù)少，沒有考慮電子熱運(yùn)動(dòng)和固液相變的影響；方法（2）考慮的因素最全面，既考慮電子熱運(yùn)動(dòng)的貢獻(xiàn)，又考慮固液相變的影響，因此涉及的參數(shù)較多；方法（3）與等熵線有關(guān)，不考慮電子熱運(yùn)動(dòng)和固液相變的影響，其可靠性取決于等熵方程的選取。從計(jì)算結(jié)果與試驗(yàn)結(jié)果對比分析看，方法（2）在很寬的壓強(qiáng)范圍內(nèi)與試驗(yàn)結(jié)果最符合，而方法（1）僅僅在固相區(qū)與試驗(yàn)結(jié)果符合。這是由于方法（1）未考慮固液相變的影響，因此不適用于材料發(fā)生熔融后的狀態(tài)區(qū)域。本文在研究過程中不考慮熔融現(xiàn)象，數(shù)值仿真過程僅考慮固相結(jié)構(gòu)，由于方法（1）涉及的參數(shù)少，因此將Mie-Grüneisen狀態(tài)方程和Hugoniot曲線計(jì)算沖擊溫度的方法添加到物質(zhì)點(diǎn)法算法中。

由熱力學(xué)關(guān)系可知

同樣以Taylor桿試驗(yàn)為例，分析材料初始溫度為30 K、考慮沖擊溫度下物質(zhì)點(diǎn)法算法的計(jì)算結(jié)果。

Taylor桿沿y軸中心位置的切面在80 μs時(shí)的結(jié)果見圖4，最終桿長L=15.88 mm，底部直徑D=13.56 mm，距離底部0.2L0位置處的直徑W=9.25 mm。與試驗(yàn)結(jié)果相比，計(jì)算結(jié)果誤差為3.61%。

綜合考慮塑性功和沖擊溫度對結(jié)構(gòu)計(jì)算結(jié)果的影響，Taylor桿沿y軸中心位置的切面在80 μs時(shí)的計(jì)算結(jié)果見圖5，最終桿長L=16.00 mm，底部直徑D=13.40 mm，距離底部0.2L0位置處的直徑W=9.26 mm。與試驗(yàn)結(jié)果相比，計(jì)算結(jié)果誤差為3.31%。

綜合以上結(jié)果，采用物質(zhì)點(diǎn)法對高速?zèng)_擊過程進(jìn)行模擬須考慮沖擊溫度對結(jié)構(gòu)溫度分布的影響。

3 結(jié) 論

針對物質(zhì)點(diǎn)法算法中溫度的影響展開研究，在物質(zhì)點(diǎn)法中同時(shí)考慮塑性功、沖擊溫度和摩擦對溫度的影響，對算法進(jìn)行改進(jìn)。以Taylor桿撞擊試驗(yàn)為例，溫度對結(jié)構(gòu)計(jì)算影響很大，隨著溫度的升高，金屬出現(xiàn)溫度軟化效應(yīng)，且沖擊溫度造成的溫升影響遠(yuǎn)遠(yuǎn)高于大變形造成的溫升影響，驗(yàn)證典型物質(zhì)點(diǎn)法算法中考慮沖擊溫度等因素的必要性。這一改進(jìn)可為刨削機(jī)理的研究打下基礎(chǔ)。

參考文獻(xiàn)：

[1] SULSKY D， CHEN Z， SCHREYER H L. A particle method for history-dependent materials[J]. Computer Methods in Applied Mechanics and Engineering， 1994， 118（1-2）： 179-196. DOI： 10.1016/0045-7825（94）90112-0.

[2] SULSKY D， ZHOU S J， SCHREYER H L. Application of a particle-in-cell method to solid mechanics[J]. Computer Physics Communications， 1995， 87（1-2）： 236-252. DOI： 10.1016/0010-4655（94）00170-7.

[3] HUANG P， ZHANG X， MA S， et al. Contact algorithms for the material point method in impact and penetration simulation[J]. International Journal for Numerical Methods in Engineering， 2011， 85（4）： 498-517. DOI： 10.1002/nme.2981.

[4] MA Z， ZHANG X， HUANG P. An object-oriented MPM framework for simulation of large deformation and contact of numerous grains[J]. Computer Modeling in Engineering and Sciences， 2010， 55（1）： 61-87. DOI： 10.3970/cmes.2010.055.061.

[5] CHEN Z， HU W， SHEN L， et al. An evaluation of the MPM for simulating dynamic failure with damage diffusion[J]. Engineering Fracture Mechanics， 2002， 69（17）： 1873-1890. DOI： 10.1016/S0013-7944（02）00066-8.

[6] CHEN Z， FENG R， XIN X， et al. A computational model for impact failure with shear-induced dilatancy[J]. International Journal for Numerical Methods in Engineering， 2003， 56（14）： 1979-1997. DOI： 10.1002/nme.651.

[7] CHEN Z， GAN Y， CHEN J. A coupled thermo-mechanical model for simulating the material failure evolution due to localized heating[J]. Computer Modeling in Engineering and Sciences， 2008， 26（2）： 123-137. DOI：10.3970/cmes.2008.026.123.

[8] 張雄， 廉艷平， 劉巖， 等. 物質(zhì)點(diǎn)法[M]. 北京： 清華大學(xué)出版社， 2013： 41-42.

[9] JOHNSON G R， HOLMQUIST T J. Evaluation of cylinder： Impact test data for constitutive model constants[J]. Journal of Applied Physics， 1988， 64（8）： 3901-3910. DOI： 10.1063/1.341344.

[10] WILLIAMS Q， JEANLOZ R， BASS J， et al. Melting curve of iron to 250 Gigapascals： A constraint on temperature at Earths center[J]. Science， 1987， 236（4798）： 181-183. DOI： 10.1126/science.236.4798.181.

[11] BROWN J M， MCQUEEN R G. Phase transitions， Grüneisen parameter， and elasticity for shocked iron between 77 GPa and 400 GPa[J]. Journal of Geophysical Research： Solid Earth， 1986， 91（B7）： 7485-7494. DOI： 10.1029/JB091iB07p07485.

[12] RICE M， MCQUEEN R G， WALSH J. Compression of solids by strong shock waves[J]. Solid State Physics， 1958， 6： 1-63. DOI： 10.1016/S0081-1947（08）60724-9.

[13] MCQUEEN RG， MARSH S P. The equation of state of solids from shock wave studies[C]∥Proceedings of High Velocity Impact Phenomena. New York： Elsevier， 1970： 293.

[14] LYZENGA G， AHRENS T J. Multiwavelength optical pyrometer for shock compression experiments[J]. Review of Scientific Instruments， 1979， 50（11）： 1421-1424. DOI： 10.1063/1.1135731.

[15] MILLER G H， AHRENS T J， STOLPER E M. The equation of state of molybdenum at 1 400 ℃[J]. Journal of Applied Physics， 1988， 63（9）： 4469-4475. DOI： 10.1063/1.341124.

[16] AHRENS T J. Shock-induced optical radiation from solids[C]// Proceedings of 2nd International Symposium on Intense Dynamic Loading and its Effects. Chengdu： Sichuan University Press， 1994.

[17] 湯文輝， 張若棋， 胡金彪， 等. 沖擊溫度的近似計(jì)算方法[J]. 力學(xué)進(jìn)展， 1998， 28（4）： 479-487.

（編輯 武曉英）