冰雹動態本構建模與驗證

王計真

摘 要：飛機飛行過程中，機身外表容易遭受冰雹撞擊，造成機體破壞，危害乘客安全。為研究結構抗冰雹沖擊性能，本文基于拉格朗日（Lagrange）方法、耦合歐拉-拉格朗日方法（CEL）和光滑粒子動力學（SPH）方法，建立冰雹動態本構模型，并根據冰撞剛性靶試驗數據，修正冰雹模型參數，驗證三種方法的有效性。數值計算結果表明，三種方法均能較好地描述冰雹動態破碎行為，且在中低速冰撞分析時，基于應變率強化效應和張力失效準則的Lagrange方法與試驗有更好的一致性。

關鍵詞：冰雹撞擊； 動態本構； 冰撞剛性靶； 應變率強化； 張力失效

中圖分類號：TG146.23 文獻標識碼：A DOI：10.19452/j.issn1007-5453.2023.08.007

基金項目： 航空科學基金（2016ZA23005）

民機在飛行或起降過程中，容易遭受冰雹的沖擊威脅。冰雹沖擊過程屬于小質量高速沖擊，對飛機結構尤其是復合材料結構安全產生重大影響，嚴重威脅乘員的生命安全，帶來嚴重的經濟損失。飛機結構抗冰雹撞擊研究具有重要理論與工程意義。

冰有20余種結構形式（晶體結構或非晶體狀態），因此，構建一種普適材料模型描述其動態行為極為困難，且冰雹撞擊分析僅在航空、軌道交通和風電等少數行業存在需求，因而針對冰雹動態本構研究有限。20世紀70年代，Haynes[1]首次實現了冰的靜態及準靜態壓縮性能測試。隨后，Schulson等[2-3]和Dempsey等[4]研究了不同溫度和晶體結構下冰的靜態力學性能和裂紋擴展屬性。研究結果表明，冰雹的拉伸強度遠小于壓縮強度，且靜態壓縮過程表現出韌性破壞特性。近年來，冰雹動態本構的研究逐步開展，表明冰為典型的應變率敏感材料。Botto[5]開展了低應變率下（10？8 ～10？3s？1）冰結構的動態力學性能測試，觀察到隨應變率增加材料特性由韌性向脆性轉變。Kim等[6]和Shazly等[7]采用霍普金森桿測試高應變率動態壓縮力學性能，結果表明冰的壓縮強度隨應變率的提高而增大。宋振華等[8-9]修正前人試驗中的不合理因素，獲取更高應變率下壓縮強度。相關文獻中很少有關于冰雹動態本構建模的研究，彈塑性本構和彈塑性水動力本構[10]都曾被用于冰雹撞擊分析，但由于未考慮高速沖擊過程的應變率效應，模擬精度不高。Tippmann等[11]也采用彈塑性本構，但考慮了屈服強度的應變率強化效應和張力失效特性，建立了一種冰雹模型，該模型得到的撞擊載荷峰值較為準確，但峰值出現時間不夠準確，且網格依賴性太高，求解效率有待提高。最近Sain等[12]較為深入地研究了冰雹撞擊破壞機理，并在總結前人研究成果的基礎上，建立了一種唯象本構模型，該模型將彈塑性模型和狀態方程進行結合，考慮流變應力的效應以更真實模擬材料失效后的承載特性。張笑宇等[13]針對碳纖維復合材料蜂窩夾芯板，進行了冰雹沖擊及多次沖擊數值研究，針對多次沖擊工況，分析了不同沖擊能量、沖擊角度以及沖擊位置的多次沖擊損傷，以及在冰雹多次沖擊下的損傷累積和疊加。

近十年來，美國聯邦航空條例（FAR）、歐洲聯合航空要求（JAR）和中國民用航空規章（CCAR）都對近代飛機結構提出專門要求，以保證結構受到大量冰雹撞擊后仍能保持其安全飛行和著陸的功能。1984年美國聯邦航空局（FAA）在其發布的具有指導作用的咨詢通報《復合材料飛機結構》（AC 20-107A）中，對包括冰雹在內的離散源撞擊飛機結構的沖擊安全問題有專門規定，所有即將投入商業運行的飛機的設計必須嚴格滿足相應的離散源撞擊指標。中國民用航空規章第25部《運輸類飛機適航標準》（CCAR-25-R4）也規定了結構元件不能因包括冰雹在內的離散源損傷而失效。

本文在冰雹撞擊剛性靶試驗的基礎上，基于有限元軟件ABAQUS，采用拉格朗日（Lagrange）、CEL和SPH三種方法建立了冰雹本構模型，并對比研究了三種模型的適用性。

1 冰雹撞擊剛性靶試驗

1.1 冰彈制備

用注射器將去離子純凈水注入模具；將模具放置在冰箱中，調節冰箱為-18℃±3℃溫度環境，保持溫差±2℃以內冰凍6h以上；取出冰彈及模具，浸入25℃±5℃的水中直至冰彈可在模具中滑動，將冰彈取出放置于-10℃±3℃的環境中待用。本次試驗根據ASTM F320[14]要求選擇最大尺寸（直徑50mm）的球形冰彈，質量約為60g。

1.2 試驗方法

剛性靶板冰撞試驗系統如圖1所示，包括D80氣炮系統、控制系統、激光瞄準系統、數據采集系統、力傳感器、剛性靶板和高速攝像機測速系統等。試驗過程中通過力傳感器測量冰撞載荷數據，高速攝像機對冰撞過程進行全程拍攝。力傳感器安裝如圖2所示，單個力傳感器垂向量程為150kN，根據載荷大小選擇單載荷傳感器或雙載荷傳感器并聯。

1.3 試驗結果分析

自然環境中的冰雹凝結溫度在一定范圍內變化，為研究凍結溫度對冰雹力學性能的影響，選擇溫度范圍-5～-30℃制作冰彈，以200m/s速度開展剛性靶板冰撞試驗，獲取載荷峰值與凍結溫度的關系，見表1。可以看出，冰撞載荷峰值變化小于1.25%，表明冰彈力學性能對凍結溫度不敏感，后續冰彈制作溫度均選為-10℃。

為研究冰雹撞擊過程，在50～300m/s速度范圍內開展純凈冰彈剛性靶板沖擊試驗。圖3是沖擊速度為50m/s時冰彈沖擊剛性靶板過程。由圖3可以看出，在t=0時，冰彈開始與靶板垂直接觸，在0～0.3ms時，冰彈開始碰撞剛性靶板至完全破碎，冰彈碎屑開始沿靶板平面橫向飛散。同時通過觀測高速攝像可以發現，冰雹沖擊破壞形式主要以脆性斷裂為主，斷裂后，冰雹開始沿靶板平面橫向擴散，呈現明顯的塑性流動特性。

圖4給出了50m/s、150m/s和300m/s三組不同速度下，冰彈的沖擊載荷曲線，圖5給出了冰彈沖擊載荷峰值及達到載荷峰值時刻隨沖擊速度的變化。由圖4、圖5可以看出，隨著沖擊速度的增加，沖擊載荷峰值逐漸增大，擬合載荷峰值與沖擊速度關系曲線，得到沖擊載荷峰值與沖擊速度的平方近似呈線性關系；隨著沖擊速度的提高，達到載荷峰值的時刻提前。

2 冰雹動態本構模型

2.1 材料模型

冰從力學屬性上可分為三類：單晶冰、多晶冰和弱冰（如雪）。本文中冰從制作過程及其微觀檢測結果均表明，在結構上和性能上屬于多晶冰范疇。

結合前人分析，總結冰雹的材料屬性包括：冰雹與應變率有很強的依賴關系，當從低應變率增大到高應變率時，冰的力學行為從韌性轉變為脆性；冰雹的壓縮強度隨應變率的增大而增大，拉伸強度比壓縮強度要小得多；冰雹高速撞擊時，在相當大的變形下完全破裂，之后表現出流體特性。結合實驗室制備的冰球模型和前人的研究，給出冰球基本力學參數，見表2。

2.2 Lagrange方法

基于Lagrange方法的有限元模型，采用彈塑性本構關系、應變率強化效應和張力失效準則定義冰的本構模型。在冰雹的高速撞擊條件下，碰撞物體會發生硬化現象，高應變率的硬化作用改變了原定的冰雹力學過程，屈服強度發生動態上升，因為屈服強度并不唯一，所以本研究對塑性階段屈服強度值的設定，采用了基于應變率的強度設置，本研究的應變率強化效應見表3。

ABAQUS 的軟件動態失效模型，其中常見的有考慮應變率效應的屈服準則、Johnson-Cook失效準則、剪切失效準則、張力失效準則等，其中剪切失效準則與張力失效準則均適用于高應變率的動態問題。張力失效準則在宏觀上描述了物質因為張力加載作用而被破壞失效的過程，在微觀層面上，它用靜水壓強（壓力）應力來衡量模型的動態碎裂和切斷力，同時該準則也可以作為中間準則與其他準則相連接，共同定義破壞失效?；趶埩κУ膭討B失效模型，可較好地模擬冰雹沖擊破壞過程。張力失效準則受預設場變量值的作用，當材料單元節點靜水壓力應力達到準則設定值時，該點就發生失效。

冰雹材料進行的設置是將偏應力失效類型設置為脆性，壓應力失效類型設置為韌性。在這種設置下，一旦壓應力達到失效臨界值，偏應力將會被設置為零，并長期保持為零，根據準則此時在單元上的應力受到切斷應力臨界值的限制，當靜水張力或靜水壓力超過了切斷應力，應力就只表現靜水壓力和靜水張力，失效單元表現出類似于流體的行為。當應力狀態將要超過截斷壓力，應力狀態進入后失效狀態時，切應力減小為零，單元將不再承受切應力只承受靜水壓應力。通過與試驗對比，此模型可較好地模擬冰雹的動態碎裂過程。

2.3 SPH方法和CEL方法

該模型可正確地代表冰雹在撞擊早期階段的冰雹行為，開始時是用大剛度來表征其特點，而在隨后階段，即撞擊后冰雹開裂，其行為特性表現為流體特性。材料模型是用拉伸張力失效判據來表征的，當達到拉伸失效應力時，偏應力分量置零，并且材料只能保持壓縮應力。在計算中，只采用彈塑性水動力材料模型的失效判據，而不用塑性失效應變判據，設置張力失效準則，同Lagrange方法一致。

3 冰雹撞擊傳感器的仿真分析

3.1 有限元模型

圖6所示為冰雹沖擊傳感器的數值仿真模型和試驗模型。為簡化分析，認為傳感器為剛體，且底部固定。為提高計算效率，在有限元仿真時采用1/4模型。

3.2 仿真結果

圖7為通過Lagrange方法、CEL方法和SPH方法進行數值仿真給出的冰雹撞擊變形過程。從變形來看，三者冰雹變形形式較為相似，均能體現出冰雹碎裂之后的塑性流動過程，和試驗現象較為一致。

圖8為當沖擊速度為49.75m/s時，采用不同方法的直徑55mm冰彈撞擊力-時間歷程的比較，仿真中的撞擊力采用的冰雹與傳感器的法向撞擊力。從圖8中可以看出，三種分析方法均能模擬出冰雹的材料硬化性質和脆性斷裂后的塑性流動過程，在中低速時SPH方法收斂性較差，不能較好地模擬出載荷歷程，而CEL方法的載荷峰值與試驗結果誤差偏小，Lagrange方法的峰值較為接近。但仿真中的峰值時間較試驗值提前，主要原因是本文未考慮測力傳感器的剛度，仿真過程中對接觸過程的處理也較為理想化，還需要進一步完善。

4 結論

本研究根據冰雹粒子的物理特性，分析冰雹沖擊動態本構，研究復合材料動態本構和損傷模型，并基于ABAQUS軟件對冰雹撞擊進行數值模擬?？傻玫饺缦陆Y論：

（1） 開展50～300m/s的沖擊試驗，從試驗數據得出沖擊載荷峰值與沖擊速度的平方近似呈線性關系，隨著沖擊速度的提高載荷峰值時刻提前。

（2） 冰雹沖擊破壞形式主要以脆性斷裂為主，脆性斷裂后，呈現明顯的塑性流動特性，破壞前無明顯的屈服階段，沖擊壓縮載荷作用下的冰雹呈現明顯應變率強化效應。

（3） CEL方法、SPH方法和Lagrange方法均能描述冰雹動態破壞行為，中低速時基于應變率強化效應和張力失效準則的Largrange方法和試驗一致性更好。

參考文獻

[1]Haynes F D. Effect of temperature on the strength of snow-ice[D]. Hanover： Cold Regions Research and Engineering Lab， 1978.

[2]Currier J H， Schulson E M. The tensile strength of ice as a func‐tion of grain size[J]. Acta Metallurgica， 1982， 30（8）： 1511-1514.

[3]Nixon W A， Schulson E M. A micromechanical view of the fracture toughness of ice[J]. Le Journal de Physique Colloques， 1987， 48（C1）： 313-319.

[4]Dempsey J P， DeFranco S J， Adamson R M， et al. Scale effects on the in-situ tensile strength and fracture of ice. Part I： Large grained freshwater ice at Spray Lakes Reservoir， Alberta[M]. Fracture Scaling： Springer Netherlands， 1999.

[5]Batto R A， Schulson E M. On the ductile-to-brittle transition in ice under compression[J]. Acta Metallurgica et Materialia， 1993， 41（7）： 2219-2225.

[6]Kim H， Keune J N. Compressive strength of ice at impact strain rates[J]. Journal of Materials Science， 2007， 42（8）： 2802-2806.

[7]Shazly M， Prakash V， Lerch B A. High strain-rate behavior of ice under uniaxial compression[J]. International Journal of Solids and Structures， 2009， 46（6）： 1499-1515.

[8]宋振華. 冰載荷作用下碳纖維復合材料桁條加筋曲面板的沖擊動力響應研究[D]. 廣州：暨南大學， 2014. Song Zhenhua. The dynamic response of stringer stiffened curved composite panels under the hail ice impact[D]. Guang‐zhou： Ji’nan University，2014.（in Chinese）

[9]Song Z H， Wang Z， Kim H， et al. Pulse shaper design and dynamic mechanical property tests on ice using the split hopkinson pressure bar[J]. International Journal of Impact Engineering， 2016， 12（2）：1-5.

[10]Anghileri M， Castelletti L M L， Invernizzi F， et al. A survey of numerical models for hail impact analysis using explicit finite element codes[J]. International Journal of Impact Engineering， 2005， 31（8）： 929-944.

[11]Tippmann J D， Kim H， Rhymer J D. Experimentally validated strain rate dependent material model for spherical ice impact simulation[J]. International Journal of Impact Engineering， 2013， 57：43-54.

[12]Sain T， Narasimhan R. Constitutive modeling of ice in the high strain rate regime[J]. International Journal of Solids and Structures， 2011， 48（5）： 817-827.

[13]張笑宇，徐緋，張玉林，等.復合材料蜂窩夾芯板冰雹沖擊及多次沖擊數值研究[J].航空科學技術，2021，32（12）：74-83. Zhang Xiaoyu， Xu Fei， Zhang Yulin， et al. Numerical study on single and repeated impact of hail on composite honeycomb sandwich panels[J].Aeronautical Science & Technology， 2021， 32（12）：74-83.（in Chinese）

[14]ASTM International. Standard test method for hail impact resistance of aerospace transparent enclosures[R]. Report No. F320-10， 2010.

Modeling and Verification of the Dynamic Constitutive of the Hailstone

Wang Jizhen

National Key Laboratory of Strength and Structural Integrity， Aircraft Strength Research Institute of China， Xi’an 710065， China

Abstract： During the flight， the fuselage may be hit by the hailstone， causing structures damage and endangering passenger safety. In order to study the structural hail impact resistance， a hailstone dynamic constitutive model based on Lagrange method， Coupled Euler-Lagrange （CEL） method and Smooth Particle Dynamics（SPH ） method were established. Then the model was verificated by comparing with test data of ice impact rigid panel. The numerical results show that the three methods can describe the dynamic breaking behavior of hailstone. During the low and medium velocity impact， the Lagrange method， based on strain rate enhancement and tension failure criterion， has better consistency with test data.

Key Words： hailstone impact； dynamic constitutive； ice impact rigid panel； strain rate enhancement； tension failure