第六届“钱令希计算力学青年奖”特邀论文

爆炸和冲击载荷下金属材料及结构的动态失效仿真*

柳占立, 初东阳, 王 涛, 王毅刚

(清华大学 航天航空学院 工程力学系, 北京 100084)

摘要: 通过数值模拟研究爆炸冲击载荷下金属材料和结构的动态失效规律,对表征爆炸冲击毁伤效应及设计新型抗冲击结构有重要意义强动载下金属材料失效涉及材料大变形、热力耦合、材料状态变化等多个复杂物理过程,给数值仿真带来了极大挑战,其中包括裂纹、剪切带等复杂失效模式的几何描述、动态失效准则的确定、塑性与损伤耦合演化的描述等问题针对这些挑战性问题,基于能量变分建立描述金属动态失效过程的热弹塑性相场理论和计算模型,实现了断裂与剪切带失效模式的统一描述,并提出了其显式有限元高效求解策略进一步将该模型应用于爆炸冲击载荷下金属脆韧失效模式转变、绝热剪切带(ASBs)自组织及冲击波作用下薄壁圆盘失效形式转变三个典型金属动态失效问题,验证了理论模型的准确性及计算模型的稳健性该工作为后续开展基于仿真的爆炸冲击毁伤评估及防护结构设计研究奠定了基础

关 键 词: 爆炸和冲击; 动态失效; 断裂; 剪切带; 数值模拟

引 言

爆炸和冲击载荷下,金属材料和结构呈现出复杂的失效形式,其中包括韧脆断裂、绝热剪切带(ASBs)等多种失效模式,如图1所示,其形成与材料局部应力状态及变形速率密切相关研究者们针对爆炸冲击载荷下金属材料和结构的动态失效开展了大量实验表征工作结果表明,对梁、板壳型结构,在动态冲击下不同的局部应力状态导致了拉伸或剪切型的韧性断裂如固支梁[1]和周边固支的低碳钢方板[2-3]在不同脉冲载荷的作用下发生了拉伸断裂和剪切破坏,而对于周边固支圆板则出现了盘状和花瓣状失效模式[4-5]当冲击载荷导致局部高应变率、低应力三轴度状态时,金属材料局部储耗能能力显著提高,塑性变形高度局部化,最终导致剪切带产生研究者们通过红外线温度探测器验证了由局部高储耗能导致剪切带内的高温,并观测到其由一系列瞬态和周期性的“热点”组成[6-8]而在真实爆炸冲击事件中,往往多个ASBs同时发生并相互作用,如在典型的厚壁圆筒动态坍塌(TWC)实验中,可以观察到ASBs的间距和模式具有自组织性[9-11]在经典的Kalthoff平板冲击实验[12]中,还发现了随冲击速率的改变试件会发生韧性断裂与ASBs失效两种失效形式之间的转变由此可见爆炸冲击载荷下金属材料和结构的失效形式及背后物理力学机制的复杂性

图1 爆炸冲击下金属材料发生弹塑性动态变形失效
Fig. 1 The diagram of thermal elastic-plastic dynamic failure issues under blast and impact loading

通过数值仿真研究爆炸冲击载荷下金属材料和结构的复杂失效形式,对表征金属材料及结构的动态失效规律及设计新型抗冲击结构有重要意义但强动载下金属材料失效涉及材料大变形、热力耦合、材料状态变化等多个复杂物理过程,给数值仿真带来了极大挑战其中难点包括复杂失效形式的几何描述、动态失效准则的确定、材料塑性与损伤耦合演化的描述等近年来,扩展有限元方法、传统损伤方法、近场动力学等数值方法在该领域的研究中取得了不同程度的进展,其中基于能量变分的相场法因其在失效准则及数值实现方面的优势得到广泛关注Francfort和Marigo[13]结合Griffith脆性断裂理论和变分方法,提出了通过全局势能最小求解断裂问题之后经过Ambrosio和Tortorelli[14]、Bourdin等[15]、Miehe等[16-17]的不断发展,模拟材料脆性断裂过程的相场模型逐渐趋于成熟,并在动态裂纹扩展问题模拟中也显示出了较好的效果[18-22]近年来研究者们又将上述方法用于金属材料韧性断裂模拟[23-28],主要在原有脆性断裂相场模型基础上从刚度退化函数、相场演化驱动能、塑性屈服面软化等方面进行修正 [23-25]McAuliffe和Waisman [27]及Miehe等[25]分别从热塑性温度软化、应变率效应考虑,建立了描述断裂与ASBs失效的相场模型但是,目前相场模型在描述金属韧性断裂或ASBs失效时的系统能量表达、塑性与损伤耦合关系等方面仍不统一,将两种失效形式统一描述的理论模型较少,将该方法应用于爆炸冲击载荷下的复杂失效模式仿真仍需进一步发展

本文基于能量变分建立针对金属材料动态失效的相场理论模型和显式计算格式,发展统一断裂与剪切带失效的描述及模型参数标定方法,将模型应用于爆炸冲击载荷下金属脆韧失效模式转变、ASBs自组织性及冲击波作用下薄壁圆盘失效形式转变三个典型金属动态失效问题分析第1节介绍了理论及计算模型的建立;第2节介绍了三个典型金属动态失效问题的仿真;第3节给出了主要结论

1 基于相场的金属动态失效理论及计算模型

韧性断裂及ASBs失效是金属动态失效的两种主要形式,其均在连续体内产生不连续面,相场法可对两种失效形式产生的不连续面进行描述如图2所示,各向同性热弹塑性变形体Ω内同时存在位移场u、温度场θ、韧性裂纹ΓT及剪切带ΓS的演化,∂Ωv,∂Ωt,∂Ωθ,∂ΩJ分别为速度、牵引力、温度及热流边界条件,且∂Ωv∪∂Ωt=∂Ω,∂Ωv∩∂Ωt=∅及∂Ωθ∪∂ΩJ=∂Ω,∂Ωθ∩∂ΩJ=∅韧性裂纹ΓT及剪切带ΓS可通过相场变量d∈[0,1]描述,d=0,1分别代表材料未失效与完全失效的状态相场演化方程通过使系统自由能密度取极小值进行推导,系统自由能密度由固体中弹塑性存储能密度、断裂失效能密度及热存储能密度三部分组成:

ψ(u,θ,d)=ψs(u,θ,d)+ψd(d)+ψθ(θ,d)

(1)

图2 热弹塑性动态失效问题示意图
Fig. 2 The sketch of thermal elastic-plastic dynamic failure issues

1.1 考虑相场的热弹塑性本构模型

根据有限变形框架,总变形率可分解为弹性、塑性及热膨胀变形率:

D=De+Dp+Dθ

(2)

材料本构为基于Cauchy应力的Jaumann客观率定义的次弹性本构,这里考虑材料失效演化后应力σ在材料固有应力σ0的基础上进行折减:

(3)

其中g(d)=(1-d)2为典型的退化函数,为材料固有弹性刚度

塑性部分采用适用于金属的各向同性强化J2流动模型,塑性变形率可表达为

(4)

其中S=trace(σ)I/3为偏应力,为等效应力,为等效塑性应变率采用考虑应变硬化、应变率硬化和热软化效应的Johnson-Cook屈服条件,为模拟剪切带中塑性演化,材料失效后对固有屈服应力σy0进行折减:

(5)

其中C,m为材料参数,θ0θm为参考温度及熔化温度,为参考应变率,为应变硬化部分,具体形式将会因材料而异

考虑各向同性热膨胀,膨胀系数为α,则热变形速率为

(6)

1.2 系统自由能密度形式

相场控制方程由式(1)的系统自由能密度具体形式最终决定,其中材料存储能包括弹性储存能密度和塑性存储能密度两部分,与次弹性本构对应的材料固有弹性存储能密度为

ψe0=σ0Dedτ,

(7)

一般将其分解为对材料失效贡献的及无贡献的两部分,模型采用膨胀体积应变能与偏斜应变能对材料失效贡献的分解方法,这与静水压不导致金属塑性演化相一致:

(8)

塑性存储能密度为塑性功中的一部分,对金属材料而言,一般认为χ∈(0,1)的塑性功转化为热,这将体现在热传导方程的体热源项中,χ通常取为常数0.9,故塑性存储能密度表达为

ψp0=(1-χ)σ0Dpdτ

(9)

为准确捕捉材料失效起始,引入能量密度阈值w0,其可理解为材料损伤起始前必要的微结构演化(如再结晶)所需的总能量对金属材料,塑性功远大于弹性储存能,因此假设只有弹塑性存储能密度之和高于(1-χ)w0时,材料失效演化起始,故材料存储能密度表达为

(10)

根据相场方法对不连续面的描述,失效能密度可通过相场变量d表达为

(11)

其中Gc为从失效演化到最终产生间断面对应的有效能量释放率,lc是表征间断面宽度的特征长度参数

对密度为ρ、比热为c的材料,存储热能密度为

(12)

1.3 控制方程

考虑相场、温度场演化的金属弹塑性动态失效问题的求解域如下:

(13)

相场求解为在Td内找到使系统自由能极小的相场分布,即由此可得相场演化方程:

(14)

考虑到不可逆损伤失效条件,引入局部历史场变量记录相场驱动能量密度历史最大值:

(15)

为采用显式时间积分,将相场控制方程改为与时间相关的形式:

(16)

其中ζ是控制相场演化的黏性参数

结合位移场的动量守恒及温度场的热传导方程,相场、温度场演化的金属弹塑性动态失效问题的控制方程及边界条件为

(17)

其中b为体积力,J是假定与温度梯度成比例的内部热流,J=-kθθkθ为热导率,γ为单位体积的体热源,由塑性功中产热部分提供,故γ=χσDp

1.4 数值实现

为在ABAQUS/Explicit中实现金属弹塑性动态失效相场模型,首先用有限元方法对控制方程进行时空离散,对求解域Ω进行离散后,(u,θ,d)可用一阶有限元形函数近似:

(18)

式中ue,θe,de分别是单元e的位移、温度和相场;分别是单元e中节点I的位移、温度和相场值;n是单元节点数;是有限元形函数利用Galerkin近似得到金属弹塑性动态失效相场模型的空间离散方程组,如下:

(19)

式中u={ue},θ={θe},d={de}是包含整个求解域中uθd的时变节点自由度的位移、温度和损伤相场矢量式(19)中矩阵的显式表达参见文献[29]

在时间积分上采用显式中心差分求解位移场,用显式前向差分求解温度场和相场速度矢量在每步(tktk+1)的中间时刻tk+1/2计算,分别被标记为位移场积分的中心差分格式为

(20)

温度场积分的前向差分格式为

(21)

相场积分的前向差分格式为

(22)

以上数值形式通过ABAQUS/Explicit用户单元子程序实现,具体方案可参见文献[29]

2 典型动态失效问题研究

本节应用上面发展的基于变分原理的相场计算模型研究爆炸冲击载荷下金属脆韧失效模式转变、ASBs自组织及冲击波作用下薄壁圆盘失效形式转变三个典型金属动态失效问题,用来验证理论模型的准确性及计算模型的稳健性

2.1 Kalthoff实验脆韧失效模式转变

Kalthoff实验中出现了脆韧失效模式转变,是评估失效模型是否具有同时模拟拉伸断裂与剪切带能力的经典案例[12]根据实验,子弹以速度v0冲击两侧对称预制缺口的高强度钢试样,在冲击速度相对较低时,如20 m/s,试样从缺口处发生拉伸断裂,裂纹与缺口成70°;而在高速冲击时,如39 m/s,试样大致沿缺口方向发生ASBs失效该实验模型及实验结果如图3所示

(a) Kalthoff实验几何及边界条件[12](b) 实验中冲击速度为20 m/s及39 m/s试样中出现的拉伸断裂和剪切带
(a) The geometry and boundary conditions for the Kalthoff test [12] (b) The tensile crack and shear band appearing in the experiments conducted by Kalthoff for impact velocities of 20 m/s and 39 m/s
图3 Kalthoff实验模型及实验结果
Fig. 3 The experimetal model and results for the Kalthoff experiment

ASBs发生条件通常为低应力三轴度状态,而拉伸断裂为高应力三轴度状态,为了用弹塑性动态失效相场模型同时描述拉伸断裂与ASBs两种失效形式,引入与应力三轴度相关的失效参数有效临界能量释放率Gc、能量密度阈值w0,其中应力三轴度定义为η=σm/σeq,σm=tr(σ)/3如此应力三轴度可作为区分两种失效形式的关键因子,Gcw0与应力三轴度的关系通过双曲正切函数表示:

(23)

其中下标S代表ASBs对应的失效参数,下标T代表拉伸断裂对应的失效参数,ηcr为决定两种失效形式转变的临界应力三轴度,参数δ决定了转变区域的宽度图4曲线表示了失效参数随应力三轴度的变化,通常ASBs对应的能量密度阈值与有效临界能量释放率高于拉伸断裂所对应的值

图4 通过Gc,w0随应力三轴度的转变描述剪切带与拉伸断裂失效参数间的转变
Fig. 4 The failure properties transition between the shear banding and the tensile fracture represented by the transition of Gc and w0 with stress triaxiality η

算例所用高强度钢材料屈服函数应变硬化部分为计算模型材料参数参考文献[30],其中由于模型表述的差异GcS=22.4 N/mm,w0S=201.6 N/mm2计算结果如图5所示,最终断裂面及剪切带分布与实验结果相吻合若考虑失效阻力,相场演化起始前能量密度阈值可视为其中的一部分,根据相场对间断面的描述理论,该部分失效阻力可通过2w0lc近似表征,相场演化起始后失效阻力可用Gc表征,故通过Gc+2w0lc可综合评价材料失效演化阻力,其云图同样在图5中给出显然,拉伸断裂的阻力要明显低于剪切带失效

(a) v0=20 m/s (b) v0=39 m/s
图5 不同冲击速度下对应的相场及相场失效演化阻力分布
Fig. 5 The distributions of phase fields and the evolutionary resistances with different impact velocities

为了解释图中的颜色,读者可以参考本文的电子网页版本,后同

为了进一步研究拉伸断裂与ASBs间转变机制,对断裂与剪切带失效起始点处的相场损伤驱动因子进行了输出,如图6所示,其中驱动因子的定义为2Hlc/Gc,即相场演化方程右端驱动能量密度与左端阻力能量密度的比值,该比值越大,相场演化越剧烈结果显示,当压缩波到达缺口位置,应力集中导致低的应力三轴度状态,当冲击速度高时,材料存储能足以使ASBs开动;而当冲击速度较低时,存储能无法使剪切带开动,应力波继续传播,在试样右界面反射形成拉伸波后传到缺口处,由应力集中导致较高三轴度状态,材料发生拉伸断裂

2.2 自组织ASBs的数值模拟

为将热弹塑性相场失效模型真正用于工程中ASBs的模拟,需根据基本实验给出失效模型参数的标定方法热弹塑性相场失效模型是局部材料点(输入)的应力-应变响应,而整个结构的响应是在实验(输出)中获得的由于剪切带的高度局部化,剪切带内材料的响应与整个结构的响应有很大的不同,因此不能直接用结构的整体响应来表示材料的局部响应合理的模型输入参数,特别是w0GcS,可以通过它们之间的一对一映射来确定下面以HY100钢为例,对Marchand和Duffy[6]纯剪试验的全剪应力-剪应变曲线进行了标定,得到了局部材料点的合理参数纯剪切模型是通过薄壁圆筒扭转来实现的,模型的几何图和荷载图如图7所示圆柱体中心添加一个缺陷,以指定剪切带开始的位置

图6 不同冲击速度下剪切带失效与拉伸断裂失效起始点相场驱动因子随时间的演化
Fig. 6 The curves of the damage driving factors at the tensile fracture initiation point and the shear banding initiation point varying with time for different impact velocities

模拟使用最小网格尺寸为0.01 mm,远小于长度尺度参数lc=0.1 mm,材料黏塑性变形参数参见文献[31]计算得到的损伤和温度的总体分布如图8(a)和(b)所示,标称应变率为1 600 s-1从图中可以看出,在模型的中间,整个环上形成了一个明显的剪切带剪切带的最高温度为244.45 ℃图8(c)~(g)所示为若干名义剪切应变下剪切带的局部变形最初,整个模型的剪切变形是均匀的当名义剪切应变达到一定值(如图8(e)中的γNOM=47%)时,剪切局部化发生在缺陷附近,并迅速扩展成圆形剪切带图9显示了整个结构的响应曲线和剪切带中的局部响应, 并将它们与实验结果进行了比较通过选择合适的w0GcS(w0=0.33 J/mm3GcS=8.0 N/mm),整个结构的应力-应变曲线和剪切带的局部响应与实验结果吻合较好特别是在整体应力-应变曲线的软化段,模型仍能与实验值很好地吻合,而在这个阶段剪切带发展迅速,这说明将w0GcS作为剪切带演化的参数是合理的

ASBs是厚壁圆筒动态坍塌(TWC)实验中的典型失效形式, 其特征间距和模式具有自组织特征[11]对这些自组织剪切带的位置、 间距和宽度演化规律的基本和定量认识, 对控制材料失效和利用剪切带具有指导意义

考虑一个厚壁圆筒, 材料为304不锈钢(SS304L), 外径为38.0 mm, 内径为22.0 mm, 夹在两个铜环(外部驱动环和内部止动环)之间外铜环的直径为40.0 mm, 内铜环的直径为20.0 mm(内外铜环厚度均为2.0 mm)初始状态下,三个圆柱环之间有0.1 mm的间隙在动态压缩载荷作用下,它们可以相互接触,也可以分离从以往的实验结果可知,在ASBs的形成和演化过程中,系统的对称性将被破坏,因此,这里采用完整模型外铜柱的外边界受到随时间呈指数衰减的动态边界压力,即

p(t)=p0exp(-t/tf),

(24)

其中,峰值压力p0为0.833 GPa,持续时间tf为60 μs特征长度尺度参数lc=0.1 mm,网格尺寸h=lc/2,黏性参数ζ取5.0×10-8 kN·s/mm2

图7 六角法兰扭转试样的几何图和载荷图[6](单位: mm)
Fig. 7 The geometric and loading sketch of the hexagonal flange torsion specimen[6](unit: mm)

(a) 损伤分布(γNOM=67%) (b) 温度分布(γNOM=67%)
(a) Damage distribution for γNOM=67% (b) Temperature distribution for γNOM=67%

(c) γNOM=25%(d) γNOM=36% (e) γNOM=47%(f) γNOM=55% (g) γNOM=67%
图8 纯剪切情况的模拟结果
Fig. 8 The numerical results of ASBs

计算得到不同时刻的相场分布(即ASBs的分布)如图10所示当外载荷产生的应力波传递到内筒时,内筒发生初始损伤随着损伤的演化,损伤局部化现象逐渐显现,形成了ASBs,最终的ASBs显示出螺旋轨迹和周期分布的特点具体演化过程为:初始加载时,结构中的应力均匀演化;随着扰动的累积,逐渐产生了小的非对称波动,并迅速放大,形成大量小的局部变形带(图10(a)),即ASBs随着时间的推移,相邻剪切带之间存在竞争,只有其中的一部分能够形成更大的ASBs(图10(b))直到几个或几十个大的ASBs延伸到外部边界,结构完全坍塌(图10(c))在ASBs起始时,ASBs和加载方向(径向)之间的角度为0°在ASBs演化的整个过程中,ASBs延伸方向与局部主应力的方向的夹角始终保持约为45°以上演化行为也被称为ASBs的自组织行为

图11显示了几个典型时刻整个场的温度(温升)分布,在t=60 μs时,最大温升达到了67.61 ℃从图中可知,ASBs的演化是伴随着温升的,因为在ASBs内,材料在极短的时间内发生了较为严重的塑性变形,产生的热量无法及时扩散出去,形成了这些温升较大的自组织带状结构,而温升导致的材料软化又会促进剪切的局部化此外,在ASBs形成和演化的早期阶段,会出现一些不连续的“热点”(图11(a)),这通常被认为是ASBs快速演化的前兆,最终的温度分布(图11(c))与最终的ASBs分布(图10(c))是一致也就是说,温升主要发生在ASBs的内部,这与前人的实验观察和模拟结果是一致的

(a) 整体应力-应变曲线与实验结果的比较 (b) 剪切带的局部剪切应变与整体剪切应变曲线
(a) The global stress-strain curves (b) The local shear strain-global shear strain curves
图9 参数w0GcS的校准过程
Fig. 9 The calibration of parameters w0 and GcS

(a) t=40 μs (b) t=50 μs (c) t=60 μs
图10 不同时刻的相场分布(即ASBs分布)
Fig. 10 The distributions of phase fields representing ASBs

(a) t=40 μs (b) t=50 μs (c) t=60 μs
图11 不同时刻的温升分布
Fig. 11 The distributions of temperature increases

根据研究者们对ASBs的研究,ASBs演化过程中出现的材料软化主要由热软化和损伤软化两种机制引起为了定量评估这两种软化机制的影响,图12给出了损伤软化因子f(d)=(1-d)2、热软化因子f(θ)=1-((θ-θ0)/(θm-θ0))m的变化,以及ASBs中一个单元(内环附近)对应的Mises应力随时间的变化通过分析这两个软化因子可知: 热软化发生在损伤软化之前; 与热软化相比,损伤软化对SS304L厚壁圆筒的坍塌和ASBs的演化贡献更大,这表明热软化不是SS304L厚壁圆筒剪切局部化失效的主导因素

图12 ASBs中单元的Mises应力、热软化系数和损伤软化系数随时间的变化
Fig. 12 The Mises stress, the thermal softening factor and the damage softening factor varying with time

虽然热软化在SS304L厚壁圆筒坍塌中贡献较小,但并不意味着热软化在ASBs的形成和演化中不重要事实上,对于SS304L,ASBs演化初始阶段的热软化系数大于损伤软化系数换言之,热软化在ASBs的早期演化中起主导作用由于热软化,大量初始ASBs被诱导出来并最终导致损伤软化主导的自组织ASBs另外,对于其他一些材料,如钨合金,其ASBs演化过程中的温升可达上千摄氏度,此时的热软化主导了ASBs演化的全过程

2.3 圆形金属板在冲击波作用下的动态失效分析

本小节应用相场损伤模型对四周固定的圆形金属板在冲击波载荷下的动态失效形式进行数值分析圆板尺寸如图13所示,内径为32 mm,外径为40 mm,厚度为0.5 mm,其中圆板外环深色部分为固定边界,冲击波载荷作用在圆板的浅色部分冲击波载荷为随时间指数衰减的压力曲线p=p0e-t/ξp0为压力峰值,ξ为衰减常数p0=20 MPa,ξ=50 μs时的压力时程曲线如图14所示

图13 薄壁金属板示意图 图14 压力载荷时程曲线(p0=20 MPa,ξ=50 μs)
Fig. 13 The sketch of the thin-walled metal plateFig. 14 The pressure loading time histories (p0=20 MPa, ξ=50 μs)

圆板在该冲击载荷下通常发生3种失效形式[4-5]:在较低冲击压力下,圆板内仅发生塑性弯曲和拉伸变形,但不发生破裂;在中等冲击压力下,塑性铰从圆板的周围移动到中心,随后在板的中心或周围处发生拉伸断裂;在较大冲击压力下,塑性铰还未移动到中心,圆板即从周围发生断裂

本算例使用不同压力峰值p0、相同衰减常数ξ=50 μs的冲击载荷,对材料为HY100钢的圆板动态失效形式进行分析分别使用Johnson-Cook塑性模型和相场断裂模型描述材料的塑性屈服和损伤演化行为,相关材料参数与文献[31]相同特征长度参数lc=0.1 mm,黏性参数ζ=5×10-8 kN·s/mm2

(a) p0=15 MPa (b) p0=25 MPa

(c) p0=35 MPa (d) p0=42 MPa

(e) p0=45 MPa (f) p0=50 MPa
图15 金属圆盘在不同爆炸冲击载荷下的失效形式
Fig. 15 The failure modes of metal plates with different blast loadings

图15为金属板在不同峰值压力载荷下的失效形式如图15(a)所示,当p0=15 MPa时,塑性铰从圆板的周围移动到中心,圆板中心处出现损伤,但并未破裂,失效形式为模式1图15(b)是当p0=25 MPa时圆板的失效形式,在塑性铰到达中心后,圆板会继续向前运动,在拉伸应力的作用下即从中心处出现花瓣状裂纹,失效形式为模式2p0=35 MPa时,板的中心和周围均出现拉伸断裂,如图15(c)所示而从图15(d)所示的圆板在p0=42 MPa时的失效形式中可见,塑性铰移动到了板的中心,但中心处还未发生破裂,周围即发生拉伸断裂p0≥45 MPa时,塑性铰还未移动到中心,圆板直接在周围被撕裂,如图15(e)和图15(f)所示,该失效形式为模式3

通过使用峰值压力p0在15 MPa到50 MPa范围内的冲击载荷进行计算发现,当p0<17 MPa时,圆板失效形式为模式1;当17 MPa≤p0<45 MPa时,圆板的失效形式为模式2;当p0≥45 MPa时,圆板的失效形式为模式3由此可见,金属弹塑性动态失效相场模型能定性捕捉到爆炸冲击下四周固定圆板的不同失效形式,能够为金属结构的抗爆设计研究提供支持

3 结 论

本文通过基于能量变分建立描述金属动态失效的热弹塑性相场理论模型,解决了金属动态失效数值仿真中复杂间断面描述、动态失效准则及塑性损伤耦合演化描述等诸多困难,并提出了动态冲击下断裂与剪切带失效的统一描述及模型参数标定方法结合ABAQUS用户子程序,实现了理论和计算模型的显式有限元求解通过将模型应用于爆炸冲击载荷下金属脆韧失效模式转变、ASBs自组织性及薄壁圆盘失效形式转变三个典型金属动态失效问题,验证了理论模型的合理性、准确性及计算模型的稳健性该工作为后续开展基于仿真的爆炸冲击毁伤评估及防护结构设计研究奠定了基础

参考文献References):

[1] MENKES S B, OPAT H J. Tearing and shear failures in explosively loaded broken beams[J]. Explosion Mechanics, 1973, 13: 480-486.

[2] TEELING-SMITH R G, NURICK G N. The deformation and tearing of thin circular plates subjected to impulsive loads[J]. International Journal of Impact Engineering, 1991, 11(1): 77-91.

[3] NURICK G N, SHAVE G C. The deformation and tearing of thin square plates subjected to impulsive loads: an experimental study[J]. International Journal of Impact Engineering, 1996, 18(1): 99-116.

[4] LEE Y W, WIERZBICKI T. Fracture prediction of thin plates under localized impulsive loading, part Ⅰ: dishing[J]. International Journal of Impact Engineering, 2005, 31(10): 1253-1276.

[5] YOUNG-WOONG L, WIERZBICKI T. Fracture prediction of thin plates under localized impulsive loading, part Ⅱ: discing and petalling[J]. International Journal of Impact Engineering, 2005, 31(10):1277-1308.

[6] MARCHAND A, DUFFY J. An experimental study of the formation process of adiabatic shear bands in a structural steel[J]. Journal of the Mechanics and Physics of Solids,1988, 36(3): 251-283.

[7] ZHOU M, ROSAKIS A J, RAVICHANDRAN G. Dynamically propagating shear bands in impact-loaded prenotched plates, Ⅰ: experimental investigations of temperature signatures and propagation speed[J]. Journal of the Mechanics and Physics of Solids,1996, 44(6): 981-1006.

[8] GUDURU P R, RAVICHANDRAN G, ROSAKIS A J. Observations of transient high temperature vortical microstructures in solids during adiabatic shear banding[J]. Physical Review E, 2001, 64(3): 036128.

[9] NESTERENKO V F, MEYERS M A, WRIGHT T W. Self-organization in the initiation of adiabatic shear bands[J]. Acta Materialia, 1998, 46(1): 327-340.

[10] XUE Q, MEYERS M A, NESTERENKO V F. Self organization of shear bands in stainless steel[J]. Materials Science and Engineering: A, 2004, 384(1): 35-46.

[11] LOVINGER Z, RIKANATI A, ROSENBERG Z, et al. Electromagnetic collapse of thick-walled cylinders to investigate spontaneous shear localization[J]. International Journal of Impact Engineering, 2011, 38(11): 918-929.

[12] KALTHOFF J F. Failure methodology of mode-Ⅱ loaded cracks[C]//International Conference on Complexity and Frontiers in Strength and Fracture. Sendai, Japan, 2001.

[13] FRANCFORT G A, MARIGO J J. Revisiting brittle fracture as an energy minimization problem[J]. Journal of the Mechanics and Physics of Solids, 1998, 46(8): 1319-1342.

[14] AMBROSIO L, TORTORELLI V M. Approximation of functional depending on jumps by elliptic functional via Γ-convergence[J]. Communications on Pure and Applied Mathematics, 1990, 43(8): 999-1036.

[15] BOURDIN B, FRANCFORT G A, MARIGO J J. Numerical experiments in revisited brittle fracture[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(4): 797-826.

[16] MIEHE C, WELSCHINGER F, HOFACKER M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[J]. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1273-1311.

[17] MIEHE C, HOFACKER M, WELSCHINGER F. A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits[J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45/48): 2765-2778.

[18] LARSEN C J. Models for dynamic fracture based on Griffith’s criterion[C]//IUTAM Symposium on Variational Concepts With Applications to the Mechanics of Materials. Dordrecht, 2010: 131-140.

[19] LARSEN C J, ORTNER C, SÜLI E. Existence of solutions to a regularized model of dynamic fracture[J]. Mathematical Models and Methods in Applied Sciences, 2010, 20(7): 1021-1048.

[20] BOURDIN B, LARSEN C J, RICHARDSON C L. A time-discrete model for dynamic fracture based on crack regularization[J]. International Journal of Fracture, 2011, 168(2): 133-143.

[21] BORDEN M J, VERHOOSEL C V, SCOTT M A, et al. A phase-field description of dynamic brittle fracture[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 217: 77-95.

[22] HOFACKER M, MIEHE C. A phase field model of dynamic fracture: robust field updates for the analysis of complex crack patterns[J]. International Journal for Numerical Methods in Engineering, 2013, 93(3): 276-301.

[23] AMBATI M, GERASIMOV T, DE LORENZIS L. Phase-field modeling of ductile fracture[J]. Computational Mechanics, 2015, 55(5): 1017-1040.

[24] AMBATI M, KRUSE R, DE LORENZIS L. A phase-field model for ductile fracture at finite strains and its experimental verification[J]. Computational Mechanics, 2016, 57(1): 149-167.

[25] MIEHE C, HOFACKER M, SCHNZEL L M, et al. Phase field modeling of fracture in multi-physics problems, part Ⅱ: coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[J]. Computer Methods in Applied Mechanics and Engineering, 2015, 294: 486-522.

[26] MIEHE C, ALDAKHEEL F, RAINA A. Phase field modeling of ductile fracture at finite strains: a variational gradient-extended plasticity-damage theory[J]. International Journal of Plasticity, 2016, 84: 1-32.

[27] MCAULIFFE C, WAISMAN H. A unified model for metal failure capturing shear banding and fracture[J]. International Journal of Plasticity, 2015, 65: 131-151.

[28] BORDEN M J, HUGHES T J R, LANDIS C M, et al. A phase-field formulation for fracture in ductile materials: finite deformation balance law derivation, plastic degradation, and stress triaxialityeffects[J]. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 130-166.

[29] WANG T, YE X, LIU Z L, et al. A phase-field model of thermo-elastic coupled brittle fracture with explicit time integration[J]. Computational Mechanics, 2020, 65: 1305-1321.

[30] CHU D Y, LI X, LIU Z L, et al. A unified phase field damage model for modeling the brittle-ductile dynamic failure mode transition in metals[J]. Engineering Fracture Mechanics, 2019, 212: 197-209.

[31] WANG T, LIU Z L, CUI Y N, et al. A thermo-elastic-plastic phase-field model for simulating the evolution and transition of adiabatic shear band, part Ⅰ: theory and model calibration[J]. Engineering Fracture Mechanics, 2020, 232: 107028. DOI: 10.1016/j.engfracmech.2020.107028.

Dynamic Failure Simulation of Metal Materials and Structures Under Blast and Impact Loading

LIU Zhanli, CHU Dongyang, WANG Tao, WANG Yigang

(Department of Engineering MechanicsSchool of Aerospace Engineering, Tsinghua UniversityBeijing 100084, P.R.China)

Abstract: Studying the dynamic failure laws of metal materials and structures under blast and impact loading through numerical simulation is of great significance for characterizing the damage effects of explosive shock and designing novel impact-resistant structures. The metals’ failure under strong dynamic loading involves multiple complex physical processes, such as the large deformation, the thermo-mechanical coupling and the material state changes. These complex physical processes bring great challenges to numerical simulation, including the geometric description of complex dynamic failure modes such as cracks and shear bands, the determination of failure criteria and the description of plasticity-damage coupled evolution, etc. In response to these challenging issues, a theoretical and computational thermal elasto-plastic phase field model was established based on the energy variational principle to describe metals’ dynamic failure. For the model, a unified description of the crack and shear band was realized, and an efficient solving strategy of explicit finite elements was proposed. The model was further applied to 3 typical metals’ dynamic failure issues under blast and impact loading: the transition between brittle and ductile failure modes of metals, the self-organization of adiabatic shear bands
(aSBs) and the transition between failure modes of thin-walled disks under shock waves. The results verify the accuracy of the theoretical model and the robustness of the computational model. This work lays a foundation for the subsequent development of damage assessment and protective structure design against blast and impact loading based on simulation.

Key words: blast and impact; dynamic failure; fracture; shear band; numerical simulation

ⓒ 应用数学和力学编委会,ISSN 1000-0887

http://www.applmathmech.cn

*收稿日期: 2020-09-07;

修订日期:2020-12-07

基金项目: 科学挑战专题基金(TZ2018002);民用飞机专项科研项目(MJ-2017-F-20)

作者简介: 柳占立(1981—),男,副教授,博士,博士生导师(通讯作者. E-mail: liuzhanli@tsinghua.edu.cn).

引用格式: 柳占立, 初东阳, 王涛, 王毅刚. 爆炸和冲击载荷下金属材料及结构的动态失效仿真[J]. 应用数学和力学, 2021, 42(1): 1-14.

(我刊编委柳占立来稿)

中图分类号: O232

文献标志码:A

DOI: 10.21656/1000-0887.410262

(Contributed by LIU Zhanli, M. AMM Editorial Board)

特邀作者简介: 柳占立,清华大学航天航空学院长聘副教授、博导2004、2009年在清华大学工程力学系获学士和博士学位,2009年至2012年在美国西北大学机械工程系从事博士后研究现任清华大学航天航空学院工程力学系副系主任主要围绕固体强度与断裂的数值仿真和工程设计开展研究,包括材料及结构的动态失效分析、爆震毁伤和防护、基于机器学习的计算力学及反向设计等研究成果为爆炸冲击波人体防护、页岩水力压裂、飞行器穿盖弹射救生等国家重大工程需求提供力学计算指导,作为课题负责人承担国防973、工程物理研究院科学挑战计划等国家重大项目在固体力学和计算力学重要期刊JMPSIJSSCMAMEIJNME等发表学术论文100余篇,出版中英文专著3部2011年教育部全国百篇优秀博士论文获得者,2015年获中国力学青年科技奖,2017年获国家自然科学基金优秀青年基金支持,2018年获教育部自然科学奖一等奖(排名2),2020年获钱令希计算力学青年奖