張靜靜 董龍雷

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

關鍵詞：物質點法；塑性變形；沖擊溫度；超高速沖擊

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

文獻標志碼：B

0 引 言

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

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

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

1 物質點法

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

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

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

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

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

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

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

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

2 考慮沖擊溫度影響的物質點法算法

采用考慮材料的溫度軟化效應的Johnson-Cook本構模型進行驗證，其對應公式為

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

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

2.1 溫度軟化效應驗證

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

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

2.2 塑性功

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

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

2.3 沖擊溫度

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

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

由熱力學關系可知

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

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

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

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

3 結 論

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

參考文獻：

[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] 張雄， 廉艷平， 劉巖， 等. 物質點法[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.

（編輯 武曉英）