SPH-HLLC coupled method for one-dimentional elastic-perfectly plastic model
-
摘要: 通过弹塑性波分析求得HLLC(Harten-Lax-van Leer-contact)近似黎曼解,提出了SPH(smoothed particle hydrodynamics)与一维理想弹塑性体模型下近似的HLLC黎曼求解器耦合的一种构造简单的算法。在SPH计算中,支持域内每个粒子对都存在一个黎曼间断问题,它的黎曼解被代入控制方程中计算。其中一维理想弹塑性体的HLLC近似黎曼解的思想是:先假设整体处于弹性状态计算黎曼解,然后对计算结果进行塑性条件修正,最后用修正后的物理变量计算HLLC近似黎曼解。将提出的SPH-HLLC耦合算法与传统SPH算法在一维算例下的计算结果进行对比,结果表明,该算法能有效模拟一维理想弹塑性体材料的碰撞,并能有效抑制在不同材料之间的压强和偏应力震荡,这是传统SPH方法很难做到的。Abstract: A 1D SPH (smoothed particle hydrodynamics) and approximate HLLC (Harten-Lax-van Leer-contact) Riemann solver coupled method for elastic-perfectly plastic model is proposed through elastic and plastic wave analysis. In SPH simulations, each particle pair in the supporting domain generates a Riemann problem, whose solutions are substituted into governing equations. The philosophy of HLLC approximate Riemann solver is to divide the procedure into three steps: assume the whole state in elastic deformation and compute Riemann problem, and then reconstruct flux under von Mises yielding conditions and compute the final HLLC Riemann solution with reconstructed fluxes. We compare the new SPH-HLLC method with the traditional SPH method in several numerical tests, which show that this method can effectively simulate collision and reflected rarefaction waves between the materials, and it can profoundly suppress oscillations of pressure and deviatoric stress at contact interface between different materials, which the traditional SPH method finds difficult to realize. Moreover, the new SPH-HLLC scheme shows better energy performance than the traditional SPH method in 2D test case where initial kinetic energy is successfully transformed into internal energy with new SPH-HLLC scheme while total energy significantly decreases with time using the traditional SPH method.
-
测试冲击波信号不失真的条件是测试系统幅频特性的平直段是否覆盖被测信号的频谱, 如果传感器的幅频特性平直段不能覆盖冲击波信号频谱, 则被测信号的波形就会发生畸变, 会造成较大的动态测试误差[1]。因此分析冲击波测试传感器的动态特性是测试冲击波的前提, 是评价和分析测试精度的基础。当前, 多数冲击波测试传感器只是按照理论设计给出传感器的固有频率, 假设其为二阶系统, 根据固有频率作为动态响应的依据。对传感器进行动态校准的也是当作二阶系统, 从动态校准曲线直接求取阻尼和固有频率等参数, 进而得到幅、相频特性。事实上, 多数冲击波传感器的数学模型并不是典型的二阶系统, 按照二阶系统得到的传感器动态特性就不准确, 采用传统的低通滤波方法处理测试数据就会带来较大的误差。以压阻式压力传感器为例, 选用模型辨识方法来获取冲击波测试传感器的传递函数和动态特性指标, 并判断是否满足冲击波的测试要求。当传感器动态特性不满足测试要求时, 提出采用动态特性补偿方法进行数据处理来提高测试精度。
对于传感器的动态特性补偿, 国外相应的方法有单参数滤波法和反卷积补偿法[2-3]。在国内, 有频率域修正法、数值微分法、叠加积分法、零极点相消原理补偿法和基于数学反演算法的多参数模型反滤波方法[1, 4-5]。单参数反卷积方法、反卷积补偿法和多参数模型反滤波方法不能采用递归实现, 运算量较大。频率域修正法、数值微分法和叠加积分法均没有考虑误差对处理结果的稳定性问题, 对实测信号是难于适用的, 因此本文中选用零极点相消原理的补偿法。
1. 冲击波传感器的指标分析方法
1.1 传感器动态标定
冲击波信号的频率高达100 kHz, 传感器在使用中会受到具有微秒级的冲击波压力上升前沿的激励, 要求用于标定的激励信号频谱分量高于冲击波信号。激波管可产生上升沿小于1 μs、平台保持时间大于5 ms的阶跃压力信号, 是理想的激励信号源[6], 选用激波管作为系统的动态标定装置, 如图 1所示。图 2为美国ENDEVCO公司生产的某支8510C-50型压阻式压力传感器的动态标定曲线。
1.2 传感器模型辨识
选用K.Steiglitz等[7]提出的一种迭代方法(GLS(SM)), 即特殊白化滤波器广义最小二乘迭代法, 它具有收敛速度快、辨识精度高等特点。选取6阶进行辨识, 建立压力传感器的动态数学模型, 得出其离散传递函数形式为:
G(z)=b0+b1Z−1+b2Z−2+b3Z−3+b4Z−4+b5Z−5+b6Z−61+a1Z−1+a2Z−2+a3Z−3+a4Z−4+a5Z−5+a6Z−6 传递函数各系数值分别为: b0=0.219 545, b1=-0.704 395, b2=0.943 348, b3=-0.997 469, b4=1.116 571, b5=-0.807 486, b6=0.230 211, a1=-4.896 553, a2=10.618 130, a3=-13.357 093, a4=10.368 354, a5=-4.653 092, a6=0.920 576。
对所建模型输入阶跃信号进行仿真计算, 得到仿真输出。仿真输出与原始输出对比如图 3所示, 从图中可以看出, 二者吻合较好。
模型在实际应用之前, 需要进行模型检验。本文中通过检验系统模型输出与原始输出残差的白噪性来判定模型的好坏。思路为求取并判断残差, 如果残差是(或近似是)零均值序列, 则认为所获取的模型是较好的模型, 否则可能是不好的。上述模型的残差均值为-0.000 14, 其值较小, 说明残差具有白噪性, 辨识的数学模型是实用的。
1.3 传感器指标分析
1.3.1 频率域性能指标
在传递函数中, 令z=ejω, 可得到传感器的频率特性和相频特性, 并得出工作频带指标。传感器的幅频特性见图 4, 从图中可以看出, 主谐振频率为332 kHz, 在低频段还有小的谐振峰, 受其影响, 系统幅频特性显著降低, 幅值误差为±2%时的工作频带为10 kHz; 幅值误差为±5%时的工作频带为16 kHz; 幅值误差为±10%时的工作频带为21.5 kHz; 幅值误差为±3 dB时的工作频带为34.5 kHz, 不能满足冲击波信号有效带宽100 kHz(±3 dB)的要求[8], 如果直接采用低通滤波器的数据处理方法会带来较大误差。
1.3.2 时间域性能指标
从图 3中仿真响应曲线可以求出传感器的时间域指标:超调量为57.0%, 上升时间为3.5 μs, 峰值时间为11.0 μs, 稳态响应时间为382.0 μs(±5%), 可见传感器的稳态响应时间较长。
从指标分析可知, 需要对传感器动态特性进行补偿, 以使其工作频带能够满足冲击波测试要求。
2. 冲击波传感器的动态特性补偿方法
测试系统的动态特性与其传递函数的极点位置密切相关, 对传递函数的零极点进行分析, 考察它们对动态性能的影响, 然后采用零极点相消补偿原理方法, 即串接一个补偿环节, 重新调节加入的极点位置, 而将原来不符合要求的极点消去, 使系统动态特性得以改善, 这就是零极点相消补偿法补偿测试系统动态特性的基本思路[4]。
按零极点相消补偿原理方法设计步骤[1, 4], 设计了补偿滤波器补偿GC(z), 补偿后等效系统形式为:
H(z)=G(z)GC(z)=b0z4+b1z3+b2z2+b3z+b4z4+a1z3+a2z2+a3z+a4 式中各系数分别为:b0=0.001 925, b1=0.005 086, b2=0.006 639, b3=0.005 723, b4=0.002 244, a1=-2.754 480, a2=2.946 422, a3=-1.447 633, a4=0.005 723。
该传感器补偿前后幅频特性对比如图 4所示。从图中可知, 补偿后系统工作频带由34.5 kHz(±3 dB)提高到100 kHz(±3 dB)。传感器补偿前后阶跃响应曲线对比如图 5所示, 从图中可以看出, 补偿后时间域指标为:超调量为6.1%, 上升时间为3.5 μs, 峰值时间为7.0 μs, 稳态响应时间为8.0 μs(±5%)。可见, 传感器动态特性补偿后频域和时域指标得到显著提高。
使用上面设计的补偿滤波器对校准曲线进行修正, 得到修正后的校准曲线, 如图 6所示。从图中可以看出, 修正后的校准曲线由原来过冲64.5%变为4.3%。直接采用截止频率为100 kHz的二阶巴特沃斯低通滤波器对传感器的动态校准曲线进行处理, 处理结果如图 7所示。从图中可以看出, 处理后传感器过冲从64.5%变为21.3%, 可见, 直接采用低通滤波数据处理方法误差较大。
3. 补偿滤波器在实验测试数据处理中的应用
测试系统建立时, 需对每支冲击波测试传感器进行动态标定, 并设计相应的补偿滤波器。图 8为某型战斗部静爆威力某点位实测冲击波超压曲线, 采用上述方法对数据进行修正处理, 结果如图 8所示。从图中可以看出, 原始压力峰值为80.3 kPa, 修正后值为63.3 kPa, 更加接近真实值。
4. 结论
(1) 用于冲击波测试的压力传感器不能简单地看作二阶系统, 应当进行动态标定, 然后通过模型辨识方法来建立数学模型, 分析动态特性。
(2) 对测试系统和实测数据的成果应用表明, 动态特性补偿方法可以拓宽测试系统的工作频带, 提高冲击波测试数据的准确度。
-
表 1 弹塑性碰撞测试
Table 1. Parameters related Al impact
算例 ρ/(kg·m−3) u/(m·s−1) p/MPa S/MPa 坐标区间/m t/ms 铝碰撞 左侧 2700 200 0 −200 −1~0 0.1 右侧 2700 −200 0 0 0~1 表 2 Wilkins算例的相关参数
Table 2. Parameters related to Wilkins test
算例 ρ/(kg·m−3) u/(m·s−1) Y0/MPa μ/GPa a0/(m·s−1) ρ0/(kg·m−3) s Γ0 坐标区间/mm t/μs Wilkins 左侧 2785 800 300 27.6 5328 2785 1.338 2 0 ~ 5 5 右侧 2785 0 300 27.6 5328 2785 1.338 2 5 ~ 50 表 3 铜铝材料撞击算例的相关参数
Table 3. Parameters related to Cu/Al collision
算例 ρ/(kg·m−3) u/(m·s−1) Y0/MPa μ/GPa a0/(m·s−1) ρ0/(kg·m−3) s Γ0 坐标区间/mm t/μs 铜铝材料撞击 左侧 8930 60 90 45.0 3940 8930 1.490 2 0 ~ 25 2 右侧 2785 0 300 27.6 5328 2785 1.338 2 25~50 表 4 二维厚壁圆筒铍壳体碰撞算例的相关参数
Table 4. Parameters related to collapse of a thick-walled cylindrical beryllium shell
算例 ρ0/(kg·m−3) |u|/(m·s−1) Y0/MPa μ/GPa a0/(m·s−1) s Γ t/μs 铍壳体 1845 417.1 330 151.9 12870 1.124 2 130 表 5 计算前后总能比(结束时刻能量/初始时刻能量)
Table 5. Total energy ratio (final total energy/initial total energy)
方法 总能比/% l0 = 0.2 mm l0 = 0.125 mm l0 = 0.08 mm 传统 90.40 92.75 94.69 SPH-HLLC 100.59 100.72 100.64 -
[1] MONAGHAN J J. Smoothed particle hydrodynamics [J]. Annual Review of Astronomy and Astrophysics, 1992, 30: 543–574. DOI: 10.1146/annurev.aa.30.090192.002551. [2] VILA J P. On particle weighted methods and smooth particle hydrodynamics [J]. Mathematical Models and Methods in Applied Sciences, 1999, 9(2): 161–209. DOI: 10.1142/S0218202599000117. [3] PARSHIKOV A N, MEDIN S A. Smoothed particle hydrodynamics usinginterparticle contact algorithms [J]. Journal of Computational Physics, 2002, 180(1): 358–382. DOI: 10.1006/jcph.2002.7099. [4] LIBERSKY L D, RANDLES P W. Shocks and discontinuities in particle methods [J]. AIP Conference Proceedings, 2006, 845(1): 1089–1092. DOI: 10.1063/1.2263512. [5] MEHRA V, CHATURVEDI S. High velocity impact of metal sphere on thin metallic plates: a comparative smooth particle hydrodynamics study [J]. Journal of Computational Physics, 2006, 212(1): 318–337. DOI: 10.1016/j.jcp.2005.06.020. [6] LIN X, BALLMANN J. A Riemann solver and a second-order Godunov method for elastic-plastic wave propagation in solids [J]. International Journal of Impact Engineering, 1993, 13(3): 463–478. DOI: 10.1016/0734-743X(93)90118-Q. [7] 姚成宝, 付梅艳, 韩峰, 等. 基于多介质Riemann问题的流体-固体耦合数值方法及其在爆炸与冲击问题中的应用 [J]. 兵工学报, 2021, 42(2): 340–355. DOI: 10.3969/j.issn.1000-1093.2021.02.012.YAO C B, FU M Y, HAN F, et al. A numerical scheme for fluid-solid interactions based on multi-medium Riemann problem and its application in explosion andimpact problems [J]. Acta Armamentarii, 2021, 42(2): 340–355. DOI: 10.3969/j.issn.1000-1093.2021.02.012. [8] CHENG J B. Harten-Lax-van Leer-contact (HLLC) approximation Riemann solver with elastic waves for one-dimensional elastic-plastic problems [J]. Applied Mathematics and Mechanics, 2016, 37(11): 1517–1538. DOI: 10.1007/s10483-016-2104-9. [9] GAO S, LIU T G. 1D exact elastic-perfectly plastic solid Riemann solver and its multi-material application [J]. Advances in Applied Mathematics and Mechanics, 2017, 9(3): 621–650. DOI: 10.4208/aamm.2015.m1340. [10] GAO S, LIU T G, YAO C B. A complete list of exact solutions for one-dimensional elastic-perfectly plastic solid Riemann problem without vacuum [J]. Communications in Nonlinear Science and Numerical Simulation, 2018, 63: 205–227. DOI: 10.1016/j.cnsns.2018.02.030. [11] LI X, ZHAI J Y, SHEN Z J. An HLLC-type approximate Riemann solver for two-dimensional elastic-perfectly plastic model [J]. Journal of Computational Physics, 2022, 448: 110675. DOI: 10.1016/j.jcp.2021.110675. [12] LIU M B, LIU G R. Smoothed particle hydrodynamics (SPH): an overview and recent developments [J]. Archives of Computational Methods in Engineering, 2010, 17(1): 25–76. DOI: 10.1007/s11831-010-9040-7. [13] HUI W H, KUDRIAKOV S. On wall overheating and other computational difficulties of shock-capturing methods [J]. Computational Fluid Dynamics Journal, 2001, 10(2): 192–209. [14] TORO E F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction [M]. New York: Springer, 1997. [15] CHEN Q, LI L, QI J, et al. A cell-centeredLagrangian scheme with an elastic-perfectly plastic solid Riemann solver for wave propagations in solids [J]. Advances in Applied Mathematics and Mechanics, 2022, 14(3): 703–724. DOI: 10.4208/aamm.OA-2020-0344. [16] WILKINS M L. Calculation of elastic-plastic flow: NSA-18-002406 [R]. Livermore: Lawrence Radiation Laboratory, 1963. [17] LIU L, CHENG J B, LIU Z. A multi-material HLLC Riemann solver with both elastic and plastic waves for 1D elastic-plastic flows [J]. Computers & Fluids, 2019, 192: 104265. DOI: 10.1016/j.compfluid.2019.104265. [18] HOWELL B P, BALL G J. A free-Lagrange augmented Godunov method for the simulation of elastic-plastic solids [J]. Journal of Computational Physics, 2002, 175(1): 128–167. DOI: 10.1006/jcph.2001.6931. [19] MAIRE P H, ABGRALL R, BREIL J, et al. A nominally second-order cell-centered Lagrangian scheme for simulating elastic-plastic flows on two-dimensional unstructured grids [J]. Journal of Computational Physics, 2013, 235: 626–665. DOI: 10.1016/j.jcp.2012.10.017. [20] QUINLAN N J, BASA M, LASTIWKA M. Truncation error in mesh-free particle methods [J]. International Journal for Numerical Methods in Engineering, 2006, 66(13): 2064–2085. DOI: 10.1002/nme.1617. [21] ZHANG Z L, LIU M B. Smoothed particle hydrodynamics with kernel gradient correction for modeling high velocity impact in two- and three-dimensional spaces [J]. Engineering Analysis with Boundary Elements, 2017, 83: 141–157. DOI: 10.1016/j.enganabound.2017.07.015. -