0 引言
为有效模拟复杂的动态响应行为,物质点法(material point method, MPM)作为一种无网格的数值模拟工具,展现了其处理大变形、复杂接触及材料破坏等问题的显著优势[6-7]。与传统有限元法(finite element method, FEM)相比,MPM能够避免网格畸变,特别是在高速冲击导致的极端大变形问题中,MPM展现了较高的计算精度和稳定性[8]。这使得MPM特别适用于模拟弹体与目标结构的复杂相互作用。尽管MPM在处理这些极端条件下的动态响应问题中表现优异,但是材料参数、冲击速度、角度等输入条件的不确定性依然会导致系统响应的显著随机性和敏感性[9]。如何有效评估结构对这些不确定性输入参数变化的敏感性,成为进一步提升预测精度和设计优化的关键。
文中将MPM和sPCE结合,提出一种新颖的敏感性分析方法,用于定量研究弹体高速冲击下结构的动态响应[17]。通过MPM建立高速冲击的数值模型,模拟弹体在不同冲击速度和角度下的冲击过程,分析结构在此过程中产生的非线性变形、损伤及破坏行为。然后通过sPCE进行敏感性分析,评估系统对关键参数变化的响应特性,定量分析材料参数、冲击速度和角度等因素对结构动态响应的影响程度。
1 自适应敏感性分析方法
结合物质点法(MPM)和稀疏多项式混沌展开(sPCE)的自适应敏感性分析方法的流程如图1所示。
首先通过MPM模拟结构在动态载荷下的响应,建立物质点模型。接着利用蒙特卡洛采样生成输入变量样本,并通过MPM计算获得一系列输出样本。基于这些样本构建稀疏多项式混沌展开(sPCE)代理模型,最终实现大规模响应计算和敏感性分析。所有算法均使用Matlab编写,在CPU为Intel(R) Core(TM) i7-7700、内存为128 GB的PC机上运行。6个具体步骤为:1)物质点法建模;2)确定性计算;3)正交多项式混沌展开;4)展开系数求解;5)概率响应计算;6)敏感性分析。
在步骤2)中选择蒙特卡洛采样随机生成大量样本,根据需要生成服从均匀分布或正态分布的样本,能够广泛覆盖输入变量的分布,适合处理高维复杂问题和输入变量的不确定性。生成N个输入变量样本ξ后,在每一个采样点上使用原MPM模型执行确定性计算,从而获得结构的动态响应f(ξ)。
在步骤3)中根据各输入变量概率分布函数,选择对应的正交多项式作为基函数构造PCE[18]:
y=f(ξ)= cαψα(ξ)
式中:α=[α1 α2…αd]表示满足 =α1+α2+…+αd的多阶下标;ψα(ξ)为一个多元多项式基函数;cα为对应于ψα(ξ)的未知展开系数。一般来说,随机变量服从正态分布选择Hermite正交多项式,服从均匀分布选择Legendre正交多项式。Hermite多项式在区间(-∞,∞)上表示为{Hh(ξ) ,而H0(ξ)=1,H1(ξ)=ξ,…,Hh+1(ξ)=ξHh(ξ)-hHh-1(ξ)。
Legendre多项式在区间[-1,1]上表示为{Ll(ξ) ,而L0(ξ)=1,L1(ξ)=ξ,…,Ll+1(ξ)=(2l+1)/(l+1)ξLl(ξ)-l/(l+1)Ll-1(ξ)。
其中,1≤h≤o,1≤l≤o,o为多项式展开阶数。
$\left.C_{\mathrm{V}}=\sqrt{\frac{\sum_{j=1}^{N} \mathcal{L}\left(\boldsymbol{\xi}^{(j)}\right)}{N}} \right\rvert\, \frac{\sum_{j=1}^{N} f\left(\boldsymbol{\xi}^{(j)}\right)}{N}$
$R^{2}=1-\frac{\sum_{j=1}^{N} \mathcal{L}\left(\boldsymbol{\xi}^{(j)}\right)}{\sum_{j=1}^{N}\left(f\left(\boldsymbol{\xi}^{(j)}\right)-\bar{f}(\boldsymbol{\xi})\right)^{2}}$
式中:$\mathcal{L}(\boldsymbol{\xi})=\left(f^{o}(\boldsymbol{\xi})-f(\boldsymbol{\xi})\right)^{2}$,fo(ξ)为代理模型计算的响应值,f(ξ)为原始模型的响应值; (ξ)为响应值f(ξ)的均值。若CV<5%,则表明模型误差已经在可接受范围内,达到了较好的精度平衡,这一标准是基于工程领域的经验,既保证了模型的逼近精度,又避免了过度计算。为进一步验证模型的拟合质量,还采用决定系数R2进行拟合验证,R2越接近1,表明模型对真实数据的拟合程度越高,从而确保了模型的收敛性和鲁棒性。
在步骤5)中,通过sPCE直接生成输出样本,由输出响应样本统计得到其概率分布和概率统计信息。
在步骤6)中,由步骤4)中得到的sPCE展开系数可直接得到展开式,敏感度可通过求目标函数f(ξ)的偏导数得到:
= β=1,2,…,d
式中0≤ξβ≤1。
2 弹体碰撞分析
每个时间步长的应变能定义为:
U= upVp
式中:up为粒子p的应变能密度,up=1/2σp,ijεp,ij,σp,ij为应力张量分量,εp,ij为应变张量分量。up的表达式为:
up=
式中:μ为剪切模量;κ为科洛索夫常数。
图3展示了t取0.15~0.4 s时应力分布(σxy)的演化过程,表明了材料在外载条件下的破坏机制。从t取0.15~0.2 s,材料中心部位开始出现应力集中,红色区域表明拉应力达到最大值。随着时间推进至t取0.25~0.4 s,拉应力进一步扩展,材料逐渐破裂,局部应力集中现象显著,最终形成不规则的裂纹分布。颜色条显示应力的大小范围,从-3 000~1 300 Pa,说明在破坏过程中存在较大的应力梯度。这一过程揭示了材料内部应力集中导致的动态破坏行为。
通过MPM建立的物理模型通常是理想状态模型,因此需要考虑实际情况的不确定性因素,包括材料特性、外界载荷等。此算例研究了弹体材料属性的不确定性对模型响应的影响,将材料的杨氏模量E和密度ρ作为不确定性因素,在动态碰撞过程中,取应变能的极值作为响应值。表1给出了两个不确定性分析方案,方案1中E的不确定范围设定在500~1 500 MPa,方案2中ρ的不确定范围设定在5~15 g/cm3。
表1 模型参数的不确定范围Table 1 Uncertainty ranges for model parameters |
| Case | E/MPa | ρ/(g/cm3) | ν | κ | G/MPa |
|---|---|---|---|---|---|
| C1 | [500, 1 500] | 10 | 0.3 | 2.8 | 384.62 |
| C2 | 1 000 | [5, 15] | 0.3 | 2.8 | 384.62 |
为了验证所提方法的有效性,对以下4种算法进行对比测试:1)MCs,以基于简单随机采样的MCs抽取样本;2)GPR,以高斯过程回归(Gaussian process regression)作为代理模型;3)DT,以决策树(decision tree)作为代理模型;4)sPCE,即通过文中的自适应敏感性分析方法,以稀疏多项式混沌展开作为代理模型。通过MCs在不确定范围内分别进行随机取样,杨氏模量E服从正态分布ξ~N(1 000,2 502),密度服从正态分布ξ~N(10,2.52),采样次数分别为500和100。500个数据通过MPM计算,结果作为衡量其他三种算法计算结果准确性的对比标准,100个数据通过GPR、DT和sPCE计算,在这两种方案中,sPCE所采用的正交多项式均选择了Hermite多项式,并且通过文中的自适应敏感性分析方法在阶数o=3处进行了截断。
表2 4种算法的概率统计量和评价指标对比(方案1)Table 2 Comparison of probability statistics and evaluation metrics for 4 algorithms (C1) |
| Algorithm | μ/J | σ/J | CV/% | R2 |
|---|---|---|---|---|
| MCs | 57.948 | 30.659 2 | — | — |
| GPR | 57.948 | 31.426 0 | 5.373 2 | 0.954 5 |
| DT | 57.948 | 33.327 9 | 7.157 7 | 0.908 3 |
| sPCE | 57.948 | 30.656 0 | 1.796 7 | 0.998 8 |
表3 4种算法的概率统计量和评价指标对比(方案2)Table 3 Comparison of probability statistics and evaluation metrics for 4 algorithms (C2) |
| Algorithm | μ/J | σ/J | CV/% | R2 |
|---|---|---|---|---|
| MCs | 53.932 | 13.241 7 | — | — |
| GPR | 53.932 | 13.978 2 | 6.871 2 | 0.923 6 |
| DT | 53.932 | 14.647 5 | 7.346 5 | 0.902 4 |
| sPCE | 53.932 | 13.241 1 | 3.778 2 | 0.984 5 |
除了不同算法的对比分析外,研究还进一步比较了不同阶数下sPCE与MPM的计算结果。如图4所示,3阶sPCE与MPM的结果具有良好的拟合度,表明基于自适应敏感性分析方法构建的sPCE能够更精确地刻画应变能随杨氏模量和密度变化的非线性特征。
图5展示了两种方案的敏感性结果。从方案1中可以观察到,敏感度随杨氏模量E的变化呈现非线性趋势。在500~900 MPa之间,敏感度逐渐增大并在900 MPa附近达到最大值,此后敏感度开始下降,在接近1 400 MPa时趋于0。这表明材料在该范围内对杨氏模量的变化表现出最大敏感性,而在极端值处(500 MPa和1 500 MPa)材料的响应趋于平稳。此外,嵌入图展示了应变能随杨氏模量的变化趋势,表明应变能在杨氏模量较高时增加较快,显示出非线性增长。
方案2展示了敏感度与密度的非线性关系,密度ρ为5 g/cm3时敏感度达到峰值1.3左右,在密度ρ为8 g/cm3时接近零,随后开始下降至零以下,呈现负相关。该趋势表明,材料对密度变化的敏感性在特定区间内显著增强。嵌入图则显示了应变能随密度的变化,随着密度的增大,应变能呈现出先增后减的趋势,最大值出现在8 g/cm3左右,表明此时材料的能量积累达到最大。
3 弹体侵彻分析
表4 AISI 52-100铬钢材料模型参数Table 4 AISI 52-100 chromium steel material model parameters |
| E/MPa | ρ/(g/cm3) | ν | G/MPa |
|---|---|---|---|
| 20 000 | 0.007 85 | 0.3 | 7 692.31 |
表5 6061-T6铝材料模型参数Table 5 6061-T6 aluminum material model parameters |
| E/MPa | ρ/(g/cm3) | ν | κ | G/MPa |
|---|---|---|---|---|
| 782 000 | 0.002 7 | 0.3 | 7 218.43 | 30 076.92 |
图7展示了不同入射角度(θ=0°,15°,30°,45°)下,在不同时刻(t=5 μs,15 μs,35 μs,50 μs)材料在冲击作用下的von Mises应力分布演变情况。随着时间的推移,从t=5 μs到t=50 μs,可以观察到材料内部的应力逐渐从冲击点向外扩散,表现出由应力集中引发的破坏逐步扩大。对于垂直冲击(θ=0°),应力分布相对对称,破坏主要集中在冲击点正下方,随着时间的推进,破坏区域逐渐扩展,但破坏形态相对规整。然而,当入射角度增大(θ=15°,30°,45°)时,冲击过程中的应力分布呈现出显著的非对称性。应力集中区域随着角度增大而向一侧偏移,尤其在较大的角度(θ=45°)下,应力集中和破坏主要集中在冲击点一侧,并伴随着材料的局部剥离和撕裂。此外,入射角度越大,材料的破坏区域扩展越广泛,损伤的非对称性越明显,破坏的剧烈程度也越高。这表明入射角度对材料在冲击载荷下的应力分布和破坏模式有显著影响,较大的角度导致更严重的局部破坏与应力集中现象。
图8展示了在不同冲击速度(v=660 m/s,910 m/s,1160 m/s,1410 m/s)以及不同入射角度(θ=0°,15°,30°,45°)下,材料应变能随时间的演变过程。总体来看,随着冲击速度的增加,应变能显著升高,且不同角度下应变能的增长趋势有所差异。在低速冲击(660 m/s)时,θ=0°的垂直冲击导致的应变能较其他角度略高,但随时间推进,所有角度的应变能逐渐趋于一致。
然而,随着冲击速度的增加,在更高速的冲击下(910 m/s及更高),θ=45°的斜向冲击显著增加了应变能,特别是在应变能峰值出现时,表现出明显高于其他角度的特征。这表明较大入射角和高速冲击会引发更显著的应变能累积和材料破坏。随着冲击时间的推移,各角度的应变能先快速上升,达到峰值后逐渐回落,呈现出应力集中导致的破坏扩展和能量释放的非线性演化规律。
在此基础上,研究冲击载荷的不确定性对模型响应的影响,将弹体的冲击速度v和入射角度θ作为不确定性因素,在动态碰撞过程中,取应变能的极值作为响应值。入射角度不确定范围设定在0°~45°,冲击速度的不确定范围设定在500~1 500 m/s。
通过MCs在不确定范围内分别进行随机取样,入射角度服从均匀分布ξ~U(0, 45),冲击速度服从均匀分布ξ~U(500, 1500),采样次数分别为1 000和200。1 000个数据通过MPM计算,结果作为衡量其他3种算法计算结果准确性的对比标准;200个数据通过MPM计算响应值作为GPR、DT和sPCE的训练样本,结果如图9所示。sPCE所采用的正交多项式选择了Legendre多项式,并且通过文中的自适应敏感性分析方法在阶数为5处进行了截断。
表6展示了4种算法在计算时间成本(N=1000)和评价指标上的对比结果。
表6 4种算法的时间成本和评价指标对比Table 6 Comparison of time costs and evaluation metrics for 4 algorithms |
| Algorithm | Tpre/s | Ttot/s | CV/% | R2 |
|---|---|---|---|---|
| MCs | — | 172 350 | — | — |
| GPR | 3.45 | 34 473.45 | 7.26 | 0.969 |
| DT | 4.79 | 34 474.79 | 6.39 | 0.971 |
| sPCE | 3.27 | 34 473.27 | 1.74 | 0.997 |
MCs的总计算时间成本Ttot主要由生成1 000个输入样本及MPM的确定性计算两部分组成,而其余3种代理模型的总计算时间成本Ttot则包括生成200个输入样本、MPM确定性计算(200组样本共耗时34 470 s)和预测1 000组输出响应的时间Tpre。与MCs方法相比,这3种代理模型显著降低了计算时间,显著提升了计算效率,其中sPCE的表现最为优异。在精度方面,sPCE以最低的变异系数(1.74%)和最高的决定系数(0.997)展现出更强的稳定性和拟合能力,优于GPR和DT。综合来看,sPCE在时间成本和精度方面均具有显著优势,是最为高效可靠的选择。
图10展示了使用5阶sPCE的应变能的连续表面分布。相较于MPM结果,sPCE能够很好地预测整体应变能随冲击速度和角度变化的宏观趋势。特别是在冲击速度较高(v>10 km/s)和角度较大的区域(θ>30°),应变能呈现显著的增加。
图11展示了角度敏感度随冲击速度和入射角度的变化,表现出非线性的波动趋势。在中等速度(v≈10 km/s)时,角度敏感度达到最大值,尤其是在θ接近20°左右,系统对角度变化最为敏感。在较低或较高速度条件下,角度敏感度降低,表明系统对角度变化的敏感性随速度的变化而有所减弱。总体来看,角度敏感度在速度和角度的中间值处达到峰值,而在两端趋于平缓。
图12展示了速度敏感度随冲击速度和角度的变化。与角度敏感度相比,速度敏感度表现出更为平缓的趋势,尤其是在速度较低时(v<10 km/s),敏感度接近于零,表明在低速情况下,系统对速度变化不敏感。然而,随着速度增加到中等值,速度敏感度略有上升,尤其在较大的入射角度下(θ接近40°),表现出轻微的负敏感度。这表明在高速且高角度冲击下,系统对速度变化的敏感性略有减弱。综合来看,角度敏感度和速度敏感度均展现出随冲击速度和角度变化的复杂非线性特征,但角度敏感度的变化幅度和波动性更为显著。
图13展示了入射角度θ和速度v的二阶敏感度随冲击速度和入射角度的变化趋势,表明系统对这两个参数之间交互效应的敏感性。图中的二阶敏感度呈现出显著的波动性,尤其是在中等速度(v≈10 km/s)和较大角度(θ>20°)时,敏感度达到较高的正值,显示出速度和角度之间强烈的耦合作用。在低速和低角度区域,二阶敏感度值较小且较为平稳,表明此时速度和角度的交互效应对系统响应的影响较弱。相较于一阶敏感度,二阶敏感度在中高速度和较大角度条件下展现出更剧烈的变化,表明在这些条件下,速度和角度的共同变化对系统的影响更为复杂,二者之间的耦合效应对材料的整体行为起到重要作用。
4 结论
通过结合物质点法和稀疏多项式混沌展开,提出了一种自适应敏感性分析方法,以深入探讨弹体在高速冲击下的结构动态响应,通过算例验证了方法的有效性。研究结果表明,杨氏模量E和密度ρ对弹体的动态响应表现出显著的非线性敏感性。在特定范围内,杨氏模量E的变化对材料的应变能影响较大。入射角度θ和冲击速度v的敏感性分析表明,系统对入射角度θ的变化在中等速度下最为敏感,尤其在接近20°时,敏感度达到峰值。相较之下,速度敏感度在低速条件下接近零,随着速度的增加而略有上升,尤其在较大入射角度下表现出负敏感度。二阶敏感度分析揭示了入射角度θ和冲击速度v之间的强耦合作用,尤其在中等速度和较大角度条件下,敏感度的波动性显著。该方法的有效性和高精度为理解弹体在高速冲击下的动态响应提供了坚实的理论基础,尤其是在需要优化设计以提高结构抗冲击性能的领域,为未来的工程应用和优化设计提供了指导。