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

張靜靜 董龍雷

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

關(guān)鍵詞：物質(zhì)點法；塑性變形；沖擊溫度；超高速沖擊

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

文獻標志碼：B

0 引 言

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

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

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

1 物質(zhì)點法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 塑性功

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

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

2.3 沖擊溫度

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

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

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

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

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

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

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

3 結(jié) 論

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

參考文獻：

[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ì)點法[M]. 北京： 清華大學出版社， 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]. 力學進展， 1998， 28（4）： 479-487.

（編輯 武曉英）