冰雹动态本构建模与验证

王计真

摘 要：飞机飞行过程中，机身外表容易遭受冰雹撞击，造成机体破坏，危害乘客安全。为研究结构抗冰雹冲击性能，本文基于拉格朗日（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失效准则、剪切失效准则、张力失效准则等，其中剪切失效准则与张力失效准则均适用于高应变率的动态问题。张力失效准则在宏观上描述了物质因为张力加载作用而被破坏失效的过程，在微观层面上，它用静水压强（压力）应力来衡量模型的动态碎裂和切断力，同时该准则也可以作为中间准则与其他准则相连接，共同定义破坏失效。基于张力失效的动态失效模型，可较好地模拟冰雹冲击破坏过程。张力失效准则受预设场变量值的作用，当材料单元节点静水压力应力达到准则设定值时，该点就发生失效。

冰雹材料进行的设置是将偏应力失效类型设置为脆性，压应力失效类型设置为韧性。在这种设置下，一旦压应力达到失效临界值，偏应力将会被设置为零，并长期保持为零，根据准则此时在单元上的应力受到切断应力临界值的限制，当静水张力或静水压力超过了切断应力，应力就只表现静水压力和静水张力，失效单元表现出类似于流体的行为。当应力状态将要超过截断压力，应力状态进入后失效状态时，切应力减小为零，单元将不再承受切应力只承受静水压应力。通过与试验对比，此模型可较好地模拟冰雹的动态碎裂过程。

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软件对冰雹撞击进行数值模拟。可得到如下结论：

（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