Dynamic deformation model of thin-walled ellipsoidal shells under impact loading
-
摘要: 为探究薄壁椭球壳在局部冲击载荷作用下的变形特征,开展了实验和数值模拟研究。利用轻气枪开展弹丸冲击实验,采用三维数字图像相关技术记录变形过程,得到了圆柱弹丸以不同速度冲击作用下椭球壳的全局变形形貌以及中心凹陷深度和凹陷边界。在弹丸冲击椭球壳的数值模拟中,重点研究了3种曲率半径变化对椭球壳凹陷深度、凹陷长短轴的影响,采用量纲分析方法得到了无量纲变形特征所依赖的主要无量纲自变量,通过参数敏感性分析消减影响较小的参数,在保持相同材料、弹体尺寸与壳体厚度同一缩比时,得到了无量纲变形特征与3种曲率半径和速度之间的响应面函数表达式,并提出了根据凹陷深度和凹陷边界预测全局变形的公式,所得表达式尺寸效应良好,预测精度较高,可为工程中大尺寸曲面薄壳冲击载荷防护设计提供参考。
-
关键词:
- 薄壁椭球壳 /
- 冲击载荷 /
- 三维数字图像相关技术 /
- 量纲分析 /
- 响应面模型
Abstract: In order to study the deformation characteristics of thin-walled ellipsoidal shells under localized impact loading, experimental investigations and numerical simulations were conducted. The global deformation characteristics, central dent depth and dent boundary of the recovery ellipsoidal shell impacted by cylindrical projectiles at different velocities were obtained by projectile impact tests on a light gas gun apparatus and three-dimensional digital image correlation (DIC) technology for deformation process record. The simulation analysis focused on the effects of three different curvature radii on the depression depth and the lengths of the major and minor axes of the ellipsoidal shell. The primary dimensionless independent variables on which the dimensionless deformation characteristics depend were determined by means of dimensional analysis. The influence of less significant parameters was reduced through parameter sensitivity analysis. Under the condition of maintaining consistent scaling ratios for material properties, projectile dimensions, and shell thickness, specific response surface function expressions between dimensionless deformation characteristics vs. three curvature radii and velocity parameters were derived. A formula for predicting global deformation based on the depth of the depression and the depression boundary was proposed. The established expression can well describe the size effect and has a high prediction accuracy, and can provide reference for the design of impact load protection of large-sized curved thin shells in engineering. -
曲面薄壳结构具有优良的承载能力,被广泛应用于管道、土木工程、航空航天、海洋工程等领域,如建筑物的曲面穹顶结构、压力容器以及整流罩结构等[1]。冲击载荷是一种局部强动态载荷,容易导致曲面薄壳结构发生大变形,降低甚至失去承载能力,从而造成安全事故,因此,需要重点考虑冲击载荷作用下曲面薄壳防护结构的动态响应特性[2]。在众多不同曲率分布的曲面薄壳结构中,椭球壳能够有效地表征大跨度曲面薄壳结构局部的曲率分布,因此,开展薄壁椭球壳动态冲击响应研究具有一定的科学意义与工程应用价值。
球壳作为椭球壳的一种特殊情况,被广泛应用于工程设计中,除直接用于组成工程结构外,也被用于近似表征不同曲率分布的薄壳。针对球壳结构的动态响应,学者们已经开展了较广泛的研究:Updike等[3-4]、Kitching等[5]对半球壳在2块刚性平板间的轴向压缩进行了多次实验研究,探究了薄壁半球壳的大变形问题,得到了薄壁半球壳的变形模式;Gputa等[6]研究了半球壳及浅球壳在落锤冲击载荷作用下的坍塌行为及变形模式,并分析了壳厚度和曲率半径对壳体力学性能的影响;Wen[7]给出了球壳在钝头弹冲击下的塑性大变形计算方法,基于实验观测,将球壳的冲击问题等效为求解钝头弹冲击一定边界条件下的等效圆板问题;宁建国等[8]对圆柱弹丸撞击浅球壳的弹塑性变形进行了一系列实验和理论研究,针对变形较集中的边缘区域,提出了一种位移模式,并讨论了不同本构模型的影响;Li等[9]通过考虑径向位移对位移模式进行改进,分析了冲击球壳的射孔响应,此外,还讨论了理论模型中各参数对冲击壳变形和射孔响应的影响。当前关于球壳受冲击载荷作用方面的研究相对较完善,但由于球壳上曲率分布是均匀的,利用球壳能够在一定程度上近似任意曲率分布浅壳的局部特征,但无法反映曲率分布对大跨度薄曲壳冲击响应的影响。
不同于球壳,椭球壳表面曲率分布不均匀,能够有效地表征任意曲率曲面薄壳结构的局部几何特征,服务于复杂曲壳结构设计。同时,椭球壳本身也被广泛应用于工程设计,如压力容器的椭球封头等[10],因此,研究椭球壳的冲击响应具有更广泛的工程应用价值。目前,对椭球壳的研究集中在应力分析、屈曲等方面,Paliwal等[11]对置于弹性地基上的小椭圆度各向同性薄壳进行了均匀内压下的弯曲变形分析,采用渐近方法将一般旋转壳的方程简化为椭球壳,分析了均匀内压与边界载荷叠加下的应力分布;Patel等[12]通过内压实验研究了椭球形封头的弹塑性屈曲行为,并给出了一种预测屈曲压力的理论分析方法,发现屈曲压力随环面半径的增大而升高;Bushnell[13]认为,压力容器中椭球封头的屈曲为分叉屈曲,屈曲行为对几何非线性和材料非线性都很敏感,几何非线性导致屈曲压力升高,材料非线性导致屈曲压力降低;Chao等[14]通过解析分析确定了径向喷管与椭球壳交点处的临界应力,利用无量纲参数作为交点的应力集中系数,分析了壳结构各组分复杂的相互作用;Ross等[15]开展了玻璃钢半椭球壳在静水压力作用下的实验研究和理论分析,考虑了几何非线性和材料非线性预测屈曲压力,其理论结果与实验结果的一致性较好;Blachut等[16]采用数值模拟方法分析了椭球壳的缺陷形状、缺陷位置对其屈曲强度的影响;Smith等[17]分析了椭球壳弹塑性屈曲的理论临界条件,采用数值模拟与实验研究相结合的方法,重点分析了边界条件、材料性能以及壁厚对屈曲的影响。Liu等[18]基于等度量变换方法,对椭球壳在冲击载荷下的变形进行了理论分析,采用数值模拟方法求解了椭球壳在冲击作用下的位移模式;陈旭东等[19]通过哈密顿原理建立了中厚圆球壳及椭球壳的自由振动控制方程,采用数值模拟方法求解了不同边界条件的自振频率。上述研究与椭球壳的力学行为密切相关,但是针对椭球壳受局部冲击时的动态变形的研究较少,椭球壳受局部冲击载荷作用下的大变形在工程实际中具有重要的应用价值,对曲面薄壳结构的抗冲击设计具有重要的参考意义。在工程设计中,通过实验方法确定椭球壳的变形特征需要较高的时间成本和经济成本,通过数值模拟方法研究大尺寸椭球壳时需要平衡网格尺寸和计算精度,难以为工程设计提供及时有效的参考,而现有的冲击载荷下椭球壳动态响应模型给出的变形计算结果受参数选取的影响较大,且需要通过数值求解方法获得结果,并不能直接获得冲击载荷作用下椭球壳变形特征与曲壳几何、材料参数之间的定量关系,而定量关系往往能够在工程设计中提供较直观快速的指导。
基于此,本文中,首先,采用无量纲参数描述薄壁椭球壳在冲击载荷作用下的变形特征,利用二阶响应面拟合椭球壳的变形特征与曲率半径、冲击速度之间的定量关系,并讨论该模型对不同尺寸椭球薄壳的适用性。然后,针对薄壁椭球壳在局部冲击载荷下的动态大变形开展实验研究和数值模拟分析,通过实验获得椭球壳在不同速度圆柱形弹丸冲击下的终态变形、凹陷深度以及凹陷长短轴等变形特征,同时采用三维数字图像相关(digital image correlation,DIC)技术获得椭球壳从受冲击到最终达到稳态的全过程的变形情况;通过ABAQUS软件模拟对照实验的变形情况,分析一般椭球壳3个方向的曲率半径对中心点凹陷深度和凹陷边界的影响。最后,利用无量纲参数描述薄壁椭球壳的变形特征,考虑尺寸效应,提出一种包含3种曲率半径、壳体厚度和冲击速度的响应面模型,以期为冲击防护工程中自由弯曲壳体的设计提供参考。
1. 薄壁椭球壳的冲击实验与DIC分析
1.1 试件与实验方法
利用轻气枪加载平头圆柱弹丸开展薄壁椭球铝壳冲击加载实验,加载和测试装置如图1所示。测试系统由椭球壳夹具、激光测速仪、高速摄像机和轻气枪等组成。发射弹丸的速度范围为25.69~60.78 m/s。椭球壳的夹具如图2所示,采用多个螺栓将扁椭球壳紧固在2块厚度为1 cm的钢板上,钢板用螺栓固支在与轻气枪相连的工作台上,可以认为冲击扁椭球壳为固支边界条件。实验中采用的扁椭球壳是由完整椭球壳截取而来,完整椭球体的三轴(x、y、z方向)长度分别为300、200、100 mm,壳体厚度为1 mm,截取部分为z方向的顶部,表1给出了扁椭球壳的截面尺寸和5052Al壳的材料参数。表中:S1和S2为椭球壳截面半径,d为椭球壳截取深度,E、
ρ 和μ分别为壳体材料的杨氏模量、密度和泊松比。实验中使用的平头圆柱弹丸材料为45钢,半径为6 mm,长度为45 mm。实验中,采用Stereo-DIC技术记录椭球壳冲击过程,DIC技术能够重构椭球壳受冲击过程中的位移场变化云图,为观察椭球壳的变形过程提供了有效手段。在实验开始前,利用标定板标定左右相机的内外参数,用于变形过程的三维重构。
1.2 实验结果与分析
通过Stereo-DIC技术获取薄壁椭球壳在冲击载荷作用下的动态变形过程,图3给出了高速摄像机记录的椭球壳的变形过程。其中,v0为冲击速度。
中心凹陷深度(w)和凹陷边界是薄壁椭球壳受局部冲击过程的重要变形特征,实验结束后,拆下椭球壳测量中心凹陷深度和凹陷边界。将DIC方法得到的中心凹深与实验静态测量得到的中心凹深进行对比,其中DIC中心凹深受限于测量时间,选取振荡上下峰值的平均值,用于代替稳态凹深。
图4对比了静态测量椭球壳回收件的最大凹陷深度与DIC得到的最大凹陷深度,采用最小二乘法,通过线性拟合得到2条w-v0线,2种测量方法的平均差值为1.588 mm。实验结束后,试件被夹具约束,存在一定的残余变形,将试件从夹具卸载后进行静态测量,该部分变形被释放,而DIC结果直接从夹持的试件上获得,因此,两者之间存在一定的偏差。由于椭球壳自身几何特征类似于轴对称旋转壳,大部分残余变形被自身几何特征约束,因此,静态测量过程中释放的变形量极小,与DIC测试结果偏差不大。可见,DIC方法能够正确描述变形过程,证明了DIC测试方法的可靠性。静态测量变形后试件虽然可以得到中心点凹陷深度,但是难以得到椭球壳的全局凹陷变形,人工测量全局变形会引入很多不确定因素,因此,考虑采用DIC方法重构变形曲面,并分析动态变形过程。
通过对散斑点进行三维重构,去除离群点。由于冲击面圆柱弹体的遮挡,冲击点中心处局部很小的区域会出现图像缺失,但不影响对凹陷变形主要特征的识别,通过重构点云插值补全重构曲面,得到随时间变化的凹陷区域,对冲击速度为25.69 m/s的实验组进行分析。Stereo-DIC重构计算的位移场误差较小,图5(a)为初始未变形区域,各散斑均未发生位移;图5(b)~(f)为弹体撞击椭球壳不同时刻的离面位移云图。变形边界在不同时刻均为椭圆形边界,且随着冲击时间的增长,冲击波由中心向边界扩展,达到峰值位移后,椭球壳进行阻尼震荡,直至最终稳定。
观察位移云图可知,椭球壳凹陷边界区域呈现为一系列椭圆形包络线,其位移差值相对较小,将椭球壳冲击方向离面位移超过1 mm的区域定义为凹陷区域。在此基础上,对椭球壳凹陷变形过程中的椭圆形凹陷边界的长短半轴特征进行分析,对冲击后的椭球壳试件进行静态测量,得到椭球壳变形边界长半轴(L1)和短半轴(L2)随冲击速度的变化,如图6所示。
由图6可知,在弹丸未贯穿壳体时,壳体凹陷变形区域的椭圆边界长半轴和短半轴均随速度升高而增大,近似呈线性关系,且凹陷变形区域的椭圆长半轴与短半轴之比在1.4~1.5之间,接近未变形椭球壳初始长短轴比1.5。
2. 数值模拟分析
2.1 椭球壳冲击动态有限元模型
通过实验和DIC技术可以得到不同冲击速度下椭球壳在不同时间下的凹陷全局变形,但是受限于DIC子区域以及弹体对散斑的遮挡,对于变形较大的壳体,无法通过DIC方法完整获取长、短轴边界位置。因此,考虑采用数值模拟方法对实验工况进行模拟,分析长、短轴长度随时间的变化,并通过拓展工况研究曲率变化对薄壁椭球壳变形的影响。
利用ABAQUS进行实验工况的数值模拟,模型如图7(a)所示。曲壳几何尺寸与图2及表1所列参数一致,薄壁椭球壳边界条件设置为完全固支,弹体与壳体间设置通用接触,法向硬接触,切向无摩擦。子弹采用实体单元C3D8R,薄壁椭球壳采用壳单元S4R,S4R单元能够适用于薄壳单元及中厚度单元的模拟。材料参数如表1所示,椭球壳材料为5052Al,采用Johnson-Cook (J-C)本构模型,J-C本构模型参数如表2所示。其中:A为屈服应力,B为硬化模量,C为应变率强化系数,n为硬化指数,m为热软化指数。子弹材料为45钢,由于子弹在冲击过程中未发生显著变形,将子弹视为刚体。数值模拟分析前开展网格收敛性测试,结果如图7(b)所示。考虑到凹陷深度随网格尺寸变化,分别给出网格尺寸为1、2和4 mm情况下中心冲击点凹陷深度随冲击速度的变化,结果表明,弹体网格设置为2 mm时能兼顾计算效率和计算精度,为此壳体网格设置为全局2 mm网格。
2.2 数值模拟与实验结果的对比及分析
对实验工况进行验证,将实验组与模拟组进行全局凹陷变形对比。图8展示了5组实验与数值模拟全局变形形貌的对比,结果表明,数值模拟与实验得到的全局变形吻合较好。数值模拟全局变形展现出2种变形模式,分别为椭圆形凹陷、椭圆形凹陷且短轴两侧发生屈曲,静态测量试件的全局变形在这2种变形模式之外,还出现了另一种变形模式(v0=36.23 m/s),即发生椭圆形凹陷并于短轴一侧发生屈曲。考虑是因为实验中通过16颗螺栓紧固,预紧过程并不能保证完美的全固支边界条件,且椭球壳加工过程可能存在一定的内在缺陷,而数值模拟无法考虑到这样的因素。
图9给出了实验、DIC和数值模拟得到的凹陷深度和凹陷长短半轴的对比。从图9可以看出,对于凹陷深度,数值模拟结果与实验及DIC结果匹配较好,且数值模拟得到的凹陷长短半轴与实验得到的凹陷长短半轴匹配良好,椭圆度随冲击速度变化的趋势相同。因此,数值模拟结果能够较好地反映局部载荷下椭球薄壳的变形情况,可以采用数值模拟方法对椭球壳在冲击载荷下的变形规律进行进一步研究。
2.3 曲率半径对变形特征的影响
本文中,重点关注曲率半径变化对椭球壳变形特征如凹深、凹陷长短半轴的影响。固定厚度、深度、弹速以及弹体尺寸时,分别改变3个方向的曲率半径,研究其对凹深、凹陷长短轴的影响。其中R3为冲击方向(z方向)的曲率半径,R1和R2为垂直于冲击方向(x、y方向)的曲率半径,经过分析可知,R1、R2对中心凹深的影响完全相同,对各自方向凹陷半轴L1、L2的影响完全对称,即令R1,new=R2,old,R2,new=R1,old时中心凹深不变,凹陷长短半轴发生如下改变:L1,new=L2,old,L2,new=L1,old。因此,可以减少讨论对象,仅研究凹深、凹陷长短半轴与曲率半径R1、R3的关系。
固定厚度h=1.0 mm,截取深度d=40 mm,给出不同冲击速度(40、45和50 m/s)下中心凹深、凹陷长短半轴随曲率半径的变化,如图10所示。从图10可以看出,中心凹陷深度w受R1影响较小,在列出的几组R2、R3、v0组合中,中心凹陷深度w都随R1增大而先增大后减小,且变化范围很小;由于前述R1、R2在z方向冲击作用中地位相同,因此,中心凹陷深度w与R2的关系与中心凹陷深度w与R1的关系规律一致,不做讨论;中心凹陷深度w受R3的影响存在一定规律,在列出的几组R1、R2组合及不同冲击速度v0下,中心凹陷深度w均随R3增大而减小,且w-R3关系近似满足二次分布。
图11给出了不同冲击速度下凹陷边界L1、L2随曲率半径R1的变化。从图11可以看出,R1对L1影响很大,R1从200 mm增大至500 mm时,L1增大约60 mm,且L1与R1的关系可以用二阶多项式拟合;L2受R1影响较小,R1从200 mm增大至500 mm时,L2仅增大约5 mm。
图12给出了不同冲击速度下凹陷长短半轴随曲率半径R3的变化。从图12可以看出,在不同的R1、R2组合下,改变冲击速度v0,R3对L1、L2都存在较大影响,且在选择的2组R1、R2组合的多组不同速度冲击下,x、y方向的凹陷边界均随R3增大而减小,且L1-R3、L2-R3近似满足二次分布,从物理意义上理解,R3增大,意味着由z方向看,椭球壳逐渐收缩、尖锐,具备更强的抗冲击变形能力,因此,凹陷区域变小。
上述分析讨论了椭球壳变形特征受3个方向曲率半径单因素变化的影响,但是难以得到曲率半径之间的交互作用对变形特征的影响,需要通过量纲分析进行定量分析。
3. 量纲分析与响应面模型
3.1 量纲分析与参数敏感性分析
椭球壳在子弹冲击作用下的响应涉及多个物理量,可以通过量纲分析法得到各物理量之间的关系。薄壁椭球壳在冲击载荷作用下的最大挠度和凹陷长短半轴都是衡量壳体结构抗冲击性能的重要指标,表3给出了对中心凹陷深度w和凹陷长短半轴L1、L2存在影响的物理量参数及其量纲。表中:Y为壳体屈服强度,Ls为弹体长度,ms为弹体质量,Rs为弹体半径,Es为弹体弹性模量,μs为弹体泊松比。
表 3 物理量参数Table 3. Physical quantity parameters分析对象 参数 量纲 分析对象 参数 量纲 薄壁椭球壳 R1 L 圆柱弹丸 Ls L R2 L ms M R3 L Rs L ρ ML−3 Es ML−1T−2 E ML−1T−2 μs 1 Y ML−1T−2 v0 LT−1 μ 1 h L d L 基于表3中的15个物理量,中心凹陷深度w和凹陷长短半轴L1、L2与各物理量之间的关系可以表示为:
L1=f1(Ls,ms,Rs,Es,μs,v0,R1,R2,R3,ρe,E,Y,μ,h,d) (1) L2=f2(Ls,ms,Rs,Es,μs,v0,R1,R2,R3,ρe,E,Y,μ,h,d) (2) w=f3(Ls,ms,Rs,Es,μs,v0,R1,R2,R3,ρe,E,Y,μ,h,d) (3) 式(1)~(3)在LMT单位制里共有3个独立量纲,选取其中3个彼此独立的物理量作为基本量,本文中选取壳厚度h、椭球壳屈服强度Y和弹丸质量ms等3个具有独立量纲的物理量作为基本量,则根据π定理,式(1)~(3)可以改写成由12个无量纲量决定的因果关系:
L1h=f1(Lsh,Rsh,R1h,R2h,R3h,dh;μ,EY,μs,EsY,msv20Yh3,ρeh3ms) (4) L2h=f2(Lsh,Rsh,R1h,R2h,R3h,dh;μ,EY,μs,EsY,msv20Yh3,ρeh3ms) (5) wh=f3(Lsh,Rsh,R1h,R2h,R3h,dh;μ,EY,μs,EsY,msv20Yh3,ρeh3ms) (6) 式(4)~(6)为此类问题的通用计算模型,适用于求解不同壳体及弹体材料、不同几何参数、不同冲击速度下的最大挠度及凹陷长短半轴。
但是由于参数过多,难以给出显式的表达,从实验给出的曲壳动态响应过程来看,应力波传播过程对动态响应的影响不大,壳体材料泊松比和弹体材料泊松比对凹陷变形影响较小,壳体密度对凹陷变形影响较小,壳体材料和弹体材料的杨氏模量与变形的发展并没有紧密耦合在一起,因此,忽略以上几项参数,并保证部分尺寸参数的等比缩放,在Ls/h、Rs/h固定的前提下建立预测模型:
L1h=f1(R1h,R2h,R3h,dh,msv20Yh3) (7) L2h=f2(R1h,R2h,R3h,dh,msv20Yh3) (8) wh=f3(R1h,R2h,R3h,dh,msv20Yh3) (9) 选取固定厚度h=1.0 mm,R1、R2、R3均在200~320 mm范围内抽样,v0在30~45 m/s范围内抽样,d在50~80 mm范围内抽样,对以上因素进行正交实验设计,采用L16的5因素4水平进行参数敏感性分析。正交实验工况设计及数值模拟结果如表4所示。
表 4 正交实验工况设计及数值模拟结果Table 4. Orthogonal experimental design and simulation results工况 d/mm R1/mm R2/mm R3/mm v0/(m·s−1) w/mm L1/mm L2/mm 1 50 200 200 200 30 15.14 47.1 47.1 2 50 240 240 240 35 19.31 64.8 64.8 3 50 280 280 280 40 22.72 70.7 70.7 4 50 320 320 320 45 26.87 82.7 82.7 5 60 200 240 280 45 24.85 54.0 62.3 6 60 240 200 320 40 21.12 53.1 46.3 7 60 280 320 200 35 20.29 79.5 89.6 8 60 320 280 240 30 16.28 71.9 63.9 9 70 200 280 320 35 17.92 42.6 55.0 10 70 240 320 280 30 15.56 50.3 63.2 11 70 280 200 240 45 25.39 76.0 58.0 12 70 320 240 200 40 23.30 94.2 73.6 13 80 200 320 240 40 22.09 55.1 80.4 14 80 240 280 200 45 26.84 78.7 90.0 15 80 280 240 320 30 15.19 51.5 45.8 16 80 320 200 280 35 18.28 66.2 46.5 根据正交实验结果,各参数对凹深、凹陷长短轴的影响程度排序,极差分析如图13所示。从图13可以看出,对于中心凹深w,各参数的影响排序为:冲击速度v0为影响程度最大的因素,且远大于其他因素;R1、R2、R3的影响程度相当;深度d为最不敏感的参数。考虑到后续拟合需要,忽略深度d对凹深的影响;对于x方向凹陷边界,R1、v0、R3、R2、p的影响依次减小;对于y方向凹陷边界,R2、v0、R3、R1、d的影响依次减小。即对于凹陷边界,同方向的曲率半径为最大的影响因素,决定着抵抗能力的强弱,在凹陷边界的参数敏感性分析中,深度d同样是最不敏感参数,因此,出于减少参数的需求,可以忽略深度d对凹陷边界的影响。
椭球壳z方向的截取深度d在较大范围内对变形特征最大挠度、长短半轴长度影响较小,原因可能是截取深度d只会影响固支边界位置,不影响壳体形状。当固支边界区域的范围远大于变形区域时,可以忽略其作用,认为后续所建立的方程是限制在模型边界远大于凹陷边界这一条件内,即可选取d为定值。因此,后续工况设计中可以固定深度d为50 mm,即固定d/h=50。量纲分析式可以简化为:
L1h=f1(R1h,R2h,R3h,msv20Yh3) (10) L2h=f2(R1h,R2h,R3h,msv20Yh3) (11) wh=f3(R1h,R2h,R3h,msv20Yh3) (12) 3.2 曲率半径和冲击速度的响应面模型
2.3节中讨论了不同曲率半径单独变化时对中心凹深、凹陷长短轴的影响,本节中通过建立凹深、凹陷长短轴响应面模型,定量分析变形特征与模型参数的关系。
厚度h固定为1.0 mm,d固定为50 mm,采用正交抽样,R1、R2、R3在200~400 mm范围内抽样,v0在35~50 m/s范围内抽样,其中R1、R2、R3选取为3水平,v0为4水平,由于R1、R2相互替换后,凹深不变,凹陷长短半轴互换,所以R1、R2在z方向受局部冲击载荷作用时体现出相似的影响规律。选取R1≤R2的全工况下,数值模拟结果如表5所示。
表 5 响应面模型的计算工况Table 5. Calculation conditions of response surface modelR1/mm R2/mm R3/mm v0/(m·s−1) w/mm L1/mm L2/mm 200 200 200 35 18.45 52.9 52.9 200 300 200 35 18.95 55.5 77.0 200 400 200 35 19.15 57.4 99.4 300 300 200 35 20.37 85.0 85.0 300 400 200 35 20.64 87.0 113.0 200 200 200 40 21.90 58.3 58.3 200 300 200 40 22.37 60.6 84.0 200 400 200 40 22.57 63.5 107.0 300 300 200 40 24.09 92.7 92.7 300 400 200 40 24.82 95.7 123.5 200 200 200 45 25.43 63.0 63.0 200 300 200 45 25.87 65.4 90.6 200 400 200 45 25.91 68.0 115.6 300 300 200 45 27.90 99.6 99.6 300 400 200 45 28.82 102.5 132.9 200 200 200 50 29.05 68.2 68.2 200 300 200 50 29.38 70.0 97.1 200 400 200 50 29.26 72.0 122.5 300 300 200 50 31.63 106.0 106.0 300 400 200 50 32.98 109.6 141.8 200 200 300 35 17.67 42.5 42.5 200 300 300 35 18.00 44.4 60.0 200 400 300 35 18.08 45.6 76.4 300 300 300 35 19.43 66.9 66.9 300 400 300 35 19.89 67.5 86.4 400 400 300 35 20.98 94.4 94.4 200 200 300 40 20.96 46.7 46.7 200 300 300 40 21.29 48.5 66.3 200 400 300 40 21.26 50.0 83.6 300 300 300 40 23.10 73.6 73.6 300 400 300 40 23.53 76.0 97.7 400 400 300 40 24.77 102.8 102.8 200 200 300 45 24.38 51.0 51.0 200 300 300 45 24.64 52.83 72.2 200 400 300 45 24.51 54.0 89.4 300 300 300 45 26.70 79.5 79.5 300 400 300 45 27.19 81.8 104.6 400 400 300 45 28.56 110.8 110.8 200 200 300 50 27.83 54.7 54.7 200 300 300 50 28.15 56.8 78.0 200 400 300 50 27.80 57.3 97.1 300 300 300 50 30.42 85.3 85.3 300 400 300 50 30.86 87.6 111.2 400 400 300 50 32.63 118.6 118.6 200 200 400 35 17.11 36.3 36.3 200 300 400 35 17.36 37.5 50.7 200 400 400 35 17.14 38.5 62.8 200 200 400 40 20.37 40.1 40.1 200 300 400 40 20.60 41.7 56.5 200 400 400 40 20.31 42.3 70.0 200 200 400 45 23.63 43.5 43.5 200 300 400 45 23.92 45.0 61.3 200 400 400 45 23.51 45.7 76.2 200 200 400 50 27.04 46.8 46.8 200 300 400 50 27.37 48.8 66.4 200 400 400 50 26.81 48.8 81.9 300 300 400 35 18.70 56.7 56.7 300 400 400 35 19.14 58.8 74.6 300 300 400 40 22.22 62.5 62.5 300 400 400 40 22.72 64.8 81.9 300 300 400 45 25.87 68.2 68.2 300 400 400 45 26.30 70.2 88.8 300 300 400 50 29.52 73.7 73.7 300 400 400 50 29.84 75.3 94.6 2.3节中单独探究了不同方向曲率半径对凹深、凹陷长短轴的影响,凹深和曲率半径都可以近似使用一次分布、二次分布来描述,且凹陷长短轴和曲率半径也满足一次分布、二次分布形式,考虑到预测精度需要,对得到的数据进行二阶响应面拟合,结果如下:
wh=−4.95e−5(R1h)2−4.95e−5(R2h)2+2.19e−5(R3h)2−9.21e−6(msv20Yh3)2+9.5e−5R1hR2h−1.1e−5R1hR3h+8.58e−6×R1hmsv20Yh3−1.1e−5R2hR3h+8.58e−6R2hmsv20Yh3−8.48e−6R3hmsv20Yh3+0.01R1h+0.01R2h−0.012R3h+0.034msv20Yh3+3.99 (13) L1h=−6.3e−5(R1h)2−2.4e−4(R2h)2+4e−4(R3h)2−168.3(msv20Yh3)2+4e−4R1hR2h−5.3e−4R1hR3h+0.44×R1hmsv20Yh3−8.2e−5R2hR3h+0.064R2hmsv20Yh3−0.21R3hmsv20Yh3+0.17R1h+0.07R2h−0.14R3h+175msv20Yh3−16.6 (14) L2h=−2.4e−4(R1h)2−7e−5(R2h)2+4e−4(R3h)2−145(msv20Yh3)2+4e−4R1hR2h−8.2e−5R1hR3h+0.067R1hmsv20Yh3−5.3e−4R2hR3h+0.44R2hmsv20Yh3−0.21R3hmsv20Yh3+0.073R1h+0.17R2h−0.15R3h+159msv20Yh3−12.9 (15) 对所建立的二阶响应面模型进行误差分析,使用多重决定系数R2、相对平均绝对误差δRAAE、相对最大绝对误差δRMAE和均方根误差δRMS等4个指标同时对目标函数模型进行精确度评价,其中 R2、δRAAE和δRMS是对模型的全局误差进行评价,δRMAE对模型的局部精度进行评估。R2越接近1,模型的精度越高;δRAAE、δRMAE和δRMS越小,模型的误差越小。随机选取20组数据对模型进行误差分析,结果如表6所示。可以看出,所建立的关于椭球壳w/h、L1/h、L2/h的二阶响应面模型精确性良好。
表 6 响应面模型误差Table 6. Response surface model error变形特征 R2 δRAAE δRMAE δRMS w/h 0.995 0.015 0.032 0.019 L1/h 0.996 0.012 0.036 0.017 L2/h 0.996 0.018 0.050 0.025 为检验所建立的二阶响应面的尺寸效应及材料效应,建立h=2.0 mm和h=0.5 mm的数值模型,随机设计3组对应工况,前2组壳体材料为5052Al,厚度分别为0.5和2.0 mm,第3组壳体材料为45钢[21],厚度为0.5 mm。模型中无法考虑应变率效应,认为5052Al的屈服强度为121 MPa,45钢的屈服强度为507 MPa。表7中,ws、L1s和L2s为数值模拟得到的变形特征,wa、L1a和L2a为响应面预测变形特征。结果表明,对于5052Al材料椭球壳,响应面模型预测值与数值模拟计算值的相对误差较小,可以认为模型的尺寸效应良好,能够对椭球壳局部冲击凹陷变形特征进行预测;对于45钢材料椭球壳,响应面模型具备一定的预测能力,对中心凹陷深度的预测误差较大,对凹陷边界的预测误差相对较小,说明模型对于不同材料的椭球壳冲击响应具备一定的预测能力。误差较大的原因可能为,数值模拟中选用的材料本构模型为J-C本构,材料塑性阶段的应力不仅与屈服强度相关,还与应变及应变率效应相关。
表 7 数值模拟变形特征与响应面预测变形特征的对比Table 7. Comparison of simulated deformation characteristics with response surface predicted deformation features材料 h/mm R1/mm R2/mm R3/mm v0/(m·s−1) ws/mm wa/mm L1s/mm L1a/mm L2s/mm L2a/mm 5052Al 0.5 200 160 180 15 12.5 12.5 46.5 45.5 37.7 38.1 5052Al 2.0 450 500 480 80 27.1 29.2 90.8 96.8 97.2 104.9 45 steel 0.5 150 100 200 30 7.5 10.8 24.0 28.7 17.8 21.9 所建立的响应面模型是对一定范围内数据的拟合,对于v0处于35~50 m/s且R1、R2和R3处于200~400 mm范围内的工况具有良好的预测精度,但仅在上述工况范围内取样并不具有广泛的代表性,对所建立的响应面模型进行适用范围讨论。对于冲击速度v0、曲率半径R1、R2、R3进行一定范围的拓展验证。
以实际值上下各10%作为误差界限,不同冲击速度及曲率半径下的适用范围不同,由图14(a)可以得到,以较低速度冲击薄壁椭球壳时,x、y方向曲率半径R1、R2较大的情况下,预测值与真实值的误差很大,远远超出误差界限,原因是,当R3固定,R1、R2取较大值时,椭球壳截取部分近似接近于平板,当冲击速度较低时,不易形成塑性变形棱区,初始撞击造成的变形会发生大幅度的回弹,使得中心区域位移大幅减小,而本文中响应面模型基于的训练集的冲击速度及曲率半径范围均处于发生塑性变形棱区的情况,所以对接近低速冲击平板的情况预测不准确。对比图14(a)和(b)可以发现,其他条件不变,仅增大R3时,预测模型的预测范围更大,原因是,增大R3使得椭球壳截取部分更凸,在同一冲击速度下更易形成棱区变形。对比图14(a)、(c)、(e)或者(b)、(d)、(f),发现固定R3,仅提高冲击速度,当冲击速度由25 m/s升至55 m/s时,预测模型的预测范围逐渐增大,当冲击速度较高时,所建立的响应面模型预测范围可以进行拓展,曲率半径不局限于训练集的200~400 mm。
图15为预测模型凹陷边界在不同冲击速度及曲率半径下的适用范围,同样采用真实值的上下10%作为误差界限,预测凹陷边界的适用范围整体趋势与预测中心凹陷深度相同,原因同样是接近平板时塑性棱区不易形成,且同样随着冲击速度升高,预测范围增大。
根据得到的中心凹陷深度、凹陷边界描述全局变形区域,由于实验及数值模拟中得到的椭球壳变形区域边界为椭圆形,且边界椭圆度近似等于椭球壳初始椭圆度,假定椭球壳的凹陷变形区域始终为一个与原椭球壳长短轴之比相同的椭圆形,则椭球壳的离面位移投影到xy平面为一簇相似的椭圆环。给出描述全局离面位移w(x, y)与中心凹陷深度、凹陷边界之间的关系式:
w(x,y)=w[(x−x0)2L21+(y−y0)2L22−1]2 (16) 抽取工况进行验证,根据响应面模型计算中心凹陷深度、凹陷长短轴长,代入全局变形公式(式(16))对比数值模拟结果。
图16给出了响应面模型与数值模拟得到的1/4区域全局变形的对比。从图16可以看出,所提出的全局变形公式能够较好地描述薄壁椭球壳受局部冲击载荷作用下的凹陷变形规律,具有一定的物理意义,反映了椭球壳凹陷变形的扩展过程。
4. 结 论
开展了薄壁椭球壳受局部冲击载荷作用的实验和数值模拟分析,采用三维DIC技术记录实验中的动态变形过程,采用ABAQUS/Explicit进行模拟,数值模拟结果与实验结果吻合较好。通过数值模拟扩展工况,探究了曲率半径、冲击速度对变形特征的影响,通过量纲分析方法探讨了各因素对变形特征的共同作用,拟合得到了一种具有尺寸效应的响应面模型。得到以下主要结论。
(1) 椭球壳凹深受冲击方向曲率半径影响较大,受垂直于冲击方向的曲率半径影响较小;椭球壳凹陷长短半轴受冲击方向曲率半径和同轴方向曲率半径影响显著。
(2) 截取深度对变形特征影响较小,在关注的大尺寸壳体受局部冲击载荷问题中仅影响边界范围。
(3) 椭球壳在局部冲击载荷作用下以椭圆形凹陷为主。椭球壳全局变形分布可以用中心凹陷深度和凹陷变形边界来近似表达,从而将曲壳在局部冲击载荷作用下的全局变形分布问题转化为中心凹陷深度和凹陷变形边界。
(4) 所建立的响应面模型精度较高,模型具有尺寸效应和一定的延展性,可用于大尺寸椭球壳受冲击载荷作用的全局凹陷变形预测。
-
表 1 薄壁金属椭球壳几何与材料参数[20]
Table 1. Geometric and material parameters of thin-walled metal ellipsoidal shells[20]
S1/mm S2/mm d/mm E/GPa ρ/(kg∙m−3) μ 240 160 40 70 2 700 0.3 A/MPa B/MPa C n m 121 327 0.009 0.544 1.12 表 3 物理量参数
Table 3. Physical quantity parameters
分析对象 参数 量纲 分析对象 参数 量纲 薄壁椭球壳 R1 L 圆柱弹丸 Ls L R2 L ms M R3 L Rs L ρ ML−3 Es ML−1T−2 E ML−1T−2 μs 1 Y ML−1T−2 v0 LT−1 μ 1 h L d L 表 4 正交实验工况设计及数值模拟结果
Table 4. Orthogonal experimental design and simulation results
工况 d/mm R1/mm R2/mm R3/mm v0/(m·s−1) w/mm L1/mm L2/mm 1 50 200 200 200 30 15.14 47.1 47.1 2 50 240 240 240 35 19.31 64.8 64.8 3 50 280 280 280 40 22.72 70.7 70.7 4 50 320 320 320 45 26.87 82.7 82.7 5 60 200 240 280 45 24.85 54.0 62.3 6 60 240 200 320 40 21.12 53.1 46.3 7 60 280 320 200 35 20.29 79.5 89.6 8 60 320 280 240 30 16.28 71.9 63.9 9 70 200 280 320 35 17.92 42.6 55.0 10 70 240 320 280 30 15.56 50.3 63.2 11 70 280 200 240 45 25.39 76.0 58.0 12 70 320 240 200 40 23.30 94.2 73.6 13 80 200 320 240 40 22.09 55.1 80.4 14 80 240 280 200 45 26.84 78.7 90.0 15 80 280 240 320 30 15.19 51.5 45.8 16 80 320 200 280 35 18.28 66.2 46.5 表 5 响应面模型的计算工况
Table 5. Calculation conditions of response surface model
R1/mm R2/mm R3/mm v0/(m·s−1) w/mm L1/mm L2/mm 200 200 200 35 18.45 52.9 52.9 200 300 200 35 18.95 55.5 77.0 200 400 200 35 19.15 57.4 99.4 300 300 200 35 20.37 85.0 85.0 300 400 200 35 20.64 87.0 113.0 200 200 200 40 21.90 58.3 58.3 200 300 200 40 22.37 60.6 84.0 200 400 200 40 22.57 63.5 107.0 300 300 200 40 24.09 92.7 92.7 300 400 200 40 24.82 95.7 123.5 200 200 200 45 25.43 63.0 63.0 200 300 200 45 25.87 65.4 90.6 200 400 200 45 25.91 68.0 115.6 300 300 200 45 27.90 99.6 99.6 300 400 200 45 28.82 102.5 132.9 200 200 200 50 29.05 68.2 68.2 200 300 200 50 29.38 70.0 97.1 200 400 200 50 29.26 72.0 122.5 300 300 200 50 31.63 106.0 106.0 300 400 200 50 32.98 109.6 141.8 200 200 300 35 17.67 42.5 42.5 200 300 300 35 18.00 44.4 60.0 200 400 300 35 18.08 45.6 76.4 300 300 300 35 19.43 66.9 66.9 300 400 300 35 19.89 67.5 86.4 400 400 300 35 20.98 94.4 94.4 200 200 300 40 20.96 46.7 46.7 200 300 300 40 21.29 48.5 66.3 200 400 300 40 21.26 50.0 83.6 300 300 300 40 23.10 73.6 73.6 300 400 300 40 23.53 76.0 97.7 400 400 300 40 24.77 102.8 102.8 200 200 300 45 24.38 51.0 51.0 200 300 300 45 24.64 52.83 72.2 200 400 300 45 24.51 54.0 89.4 300 300 300 45 26.70 79.5 79.5 300 400 300 45 27.19 81.8 104.6 400 400 300 45 28.56 110.8 110.8 200 200 300 50 27.83 54.7 54.7 200 300 300 50 28.15 56.8 78.0 200 400 300 50 27.80 57.3 97.1 300 300 300 50 30.42 85.3 85.3 300 400 300 50 30.86 87.6 111.2 400 400 300 50 32.63 118.6 118.6 200 200 400 35 17.11 36.3 36.3 200 300 400 35 17.36 37.5 50.7 200 400 400 35 17.14 38.5 62.8 200 200 400 40 20.37 40.1 40.1 200 300 400 40 20.60 41.7 56.5 200 400 400 40 20.31 42.3 70.0 200 200 400 45 23.63 43.5 43.5 200 300 400 45 23.92 45.0 61.3 200 400 400 45 23.51 45.7 76.2 200 200 400 50 27.04 46.8 46.8 200 300 400 50 27.37 48.8 66.4 200 400 400 50 26.81 48.8 81.9 300 300 400 35 18.70 56.7 56.7 300 400 400 35 19.14 58.8 74.6 300 300 400 40 22.22 62.5 62.5 300 400 400 40 22.72 64.8 81.9 300 300 400 45 25.87 68.2 68.2 300 400 400 45 26.30 70.2 88.8 300 300 400 50 29.52 73.7 73.7 300 400 400 50 29.84 75.3 94.6 表 6 响应面模型误差
Table 6. Response surface model error
变形特征 R2 δRAAE δRMAE δRMS w/h 0.995 0.015 0.032 0.019 L1/h 0.996 0.012 0.036 0.017 L2/h 0.996 0.018 0.050 0.025 表 7 数值模拟变形特征与响应面预测变形特征的对比
Table 7. Comparison of simulated deformation characteristics with response surface predicted deformation features
材料 h/mm R1/mm R2/mm R3/mm v0/(m·s−1) ws/mm wa/mm L1s/mm L1a/mm L2s/mm L2a/mm 5052Al 0.5 200 160 180 15 12.5 12.5 46.5 45.5 37.7 38.1 5052Al 2.0 450 500 480 80 27.1 29.2 90.8 96.8 97.2 104.9 45 steel 0.5 150 100 200 30 7.5 10.8 24.0 28.7 17.8 21.9 -
[1] LIM H K, LEE J S. On the structural behavior of ship’s shell structures due to impact loading [J]. International Journal of Naval Architecture and Ocean Engineering, 2018, 10(1): 103–118. DOI: 10.1016/j.ijnaoe.2017.03.002. [2] MOHAMMAD Z, GUPTA P K, BAQI A, et al. Ballistic performance of monolithic and double layered thin-metallic hemispherical shells at normal and oblique impact [J]. Thin-Walled Structures, 2021, 159: 107257. DOI: 10.1016/j.tws.2020.107257. [3] UPDIKE D P, KALNINS A. Axisymmetric behavior of an elastic spherical shell compressed between rigid plates [J]. Journal of Applied Mechanics, 1970, 37(3): 635–640. DOI: 10.1115/1.3408592. [4] UPDIKE D P. On the large deformation of a rigid-plastic spherical shell compressed by a rigid plate [J]. Journal of Engineering for Industry, 1972, 94(3): 949–955. DOI: 10.1115/1.3428276. [5] KITCHING R, HOULSTON R, JOHNSON W. A theoretical and experimental study of hemispherical shells subjected to axial loads between flat plates [J]. International Journal of Mechanical Sciences, 1975, 17(11/12): 693–703. DOI: 10.1016/0020-7403(75)90072-7. [6] GUPTA N K, MOHAMED S N, VELMURUGAN R. Experimental and numerical investigations into collapse behaviour of thin spherical shells under drop hammer impact [J]. International Journal of Solids and Structures, 2007, 44(10): 3136–3155. DOI: 10.1016/j.ijsolstr.2006.09.014. [7] WEN H M. Large plastic deformation of spherical shells under impact by blunt-ended missiles [J]. International Journal of Pressure Vessels and Piping, 1997, 73(2): 147–152. DOI: 10.1016/S0308-0161(97)00043-4. [8] 宁建国, 杨桂通. 球形扁壳在冲击载荷作用下的超临界变形 [J]. 爆炸与冲击, 1992, 12(3): 206–212. DOI: 10.11883/1001-1455(1992)03-0206-7.NING J G, YANG G T. Supercritical deformations of shallow spherical shells under impact [J]. Explosion and Shock Waves, 1992, 12(3): 206–212. DOI: 10.11883/1001-1455(1992)03-0206-7. [9] LI J Q, REN H L, NING J G. Deformation and failure of thin spherical shells under dynamic impact loading: experiment and analytical model [J]. Thin-Walled Structures, 2021, 161: 107403. DOI: 10.1016/j.tws.2020.107403. [10] ZHENG J, LI K, LIU S, et al. Effect of shape imperfection on the buckling of large-scale thin-walled ellipsoidal head in steel nuclear containment [J]. Thin-Walled Structures, 2018, 124: 514–522. DOI: 10.1016/j.tws.2018.01.001. [11] PALIWAL D N, GUPTA R, JAIN A. Stress analysis of ellipsoidal shell on an elastic foundation [J]. International Journal of Pressure Vessels and Piping, 1993, 56(2): 229–242. DOI: 10.1016/0308-0161(93)90095-B. [12] PATEL P R, GILL S S. Experiments on the buckling under internal pressure of thin torispherical ends of cylindrical pressure vessels [J]. International Journal of Mechanical Sciences, 1978, 20(3): 159–175. DOI: 10.1016/0020-7403(78)90003-6. [13] BUSHNELL D. Nonsymmetric buckling of internally pressurized ellipsoidal and torispherical elastic-plastic pressure vessel heads [J]. Journal of Pressure Vessel Technology, 1977, 99(1): 54–63. DOI: 10.1115/1.3454520. [14] CHAO Y J, SUTTON M A. Stress analysis of ellipsoidal shell with radial nozzle [J]. International Journal of Pressure Vessels and Piping, 1985, 21(2): 89–108. DOI: 10.1016/0308-0161(85)90042-0. [15] ROSS C T F, HUAT B H, CHEI T B, et al. The buckling of GRP hemi-ellipsoidal dome shells under external hydrostatic pressure [J]. Ocean Engineering, 2003, 30(5): 691–705. DOI: 10.1016/s0029-8018(02)00039-2. [16] BLACHUT J, JAISWAL O R. On the choice of initial geometric imperfections in externally pressurized shells [J]. Journal of Pressure Vessel Technology, 1999, 121(1): 71–76. DOI: 10.1115/1.2883670. [17] SMITH P, BŁACHUT J. Buckling of externally pressurized prolate ellipsoidal domes [J]. Journal of Pressure Vessel Technology, 2008, 130(1): 011210. DOI: 10.1115/1.2834457. [18] LIU L, LI J Q. Dynamic deformation and perforation of ellipsoidal thin shell impacted by flat-nose projectile [J]. Materials, 2022, 15(12): 4124. DOI: 10.3390/ma15124124. [19] 陈旭东, 叶康生. 中厚椭球壳自由振动动力刚度法分析 [J]. 振动与冲击, 2016, 35(6): 85–90. DOI: 10.13465/j.cnki.jvs.2016.06.015.CHEN X D, YE K S. Free vibration analysis of moderately thick elliptical shells using the dynamic stiffness method [J]. Journal of Vibration and Shock, 2016, 35(6): 85–90. DOI: 10.13465/j.cnki.jvs.2016.06.015. [20] MA T B, SHEN Y, NING J G, et al. Analysis on dynamic shear fracture based on a novel damage evolution model [J]. International Journal of Impact Engineering, 2024, 183: 104810. DOI: 10.1016/j.ijimpeng.2023.104810. [21] 陈刚, 陈忠富, 徐伟芳, 等. 45钢的J-C损伤失效参量研究 [J]. 爆炸与冲击, 2007, 27(2): 131–135. DOI: 10.11883/1001-1455(2007)02-0131-05.CHEN G, CHEN Z F, XU W F, et al. Investigation on the J-C ductile fracture parameters of 45 steel [J]. Explosion and Shock Waves, 2007, 27(2): 131–135. DOI: 10.11883/1001-1455(2007)02-0131-05. -