An improved finite particle method for discontinuous interface problems
-
摘要: 有限粒子法(finite particle method,FPM)作为SPH(smoothed particle hydrodynamics)方法的重要改进,有效提高了边界区域粒子的近似精度,但是当FPM处理多物理场时,在不连续界面附近的计算精度会大大降低,并且FPM必须满足的矩阵非奇异性也提高了对界面处理的要求。本文中基于DSPH(discontinuous SPH)方法,提出了一种考虑界面不连续的改进FPM—DSFPM(discontinuous special FPM)法,旨在改善FPM在界面不连续处的计算精度,从而进一步提高其计算效率和稳定性。首先,分析了DSFPM的核近似精度。其次,根据不同的工程问题,给出DSFPM处理小变形和大变形问题的算法流程。利用DSFPM、DSPH和FPM等3种方法对弹性铝块小变形碰撞冲击算例进行了模拟,通过对比分析铝块的速度和应力以及计算时间验证了DSFPM算法在非连续界面处计算精度和计算效率的优势。最后,通过结合DSFPM和DFPM(discontinuous FPM)实现了对于大变形问题的模拟。Abstract: The finite particle method(FPM) is an important improvement for the smoothed particle hydrodynamics(SPH) method, which effectively improves the calculation accuracy of boundary particles. However, when the discontinuous physical field is solved by the FPM, the accuracy in the vicinity of the discontinuous interface is greatly reduced, and the non-singularity of the matrix must be satisfied in the FPM, which requires an elaborate handling of the interface. Based on the discontinuous SPH(DSPH) method, this paper proposed an improved FPM-discontinuous special FPM(DSFPM), which considers the discontinuous interface, aiming to improve the computational accuracy at the interface and further improve the efficiency and stability of the FPM. In this paper, the estimation accuracy of the DSFPM was analyzed firstly, and then the algorithm flow diagram of the DSFPM to deal with the small deformation and large deformation problems was demonstrated. Next, the DSFPM, DSPH and FPM were used to simulate the small deformation problem-elastic aluminum blocks impact. By comparing the velocity and stress of the aluminum blocks and computational time, we verified the accuracy and computational efficiency of the DSFPM. Finally, the simulation of the large deformation problem was realized by a combining method with the DSFPM and DFPM.
-
Key words:
- FPM /
- interface /
- DSFPM /
- accuracy /
- computational cost
-
爆轰产物组成和爆轰参数是决定炸药爆炸性能重要参数, 也是设计和改进炸药的重要依据。从20世纪40年代开始, 很多国外学者在爆轰产物状态方程研究领域作了大量的工作, 提出了很多的数学物理模型[1], 对爆轰产物组成和爆轰参数进行研究, 如BKW[2]、LJD、JCZ3及WCA模型, 并发展了有关的计算方法及计算程序, 如BKW、APPEGE、RUBY、LAMINEUR、TIGER、CHEQ[3]程序都能较成功地解决这一领域的一些问题, 其中BKW、APPEGE、RUBY、LAMINEUR、TIGER等Fortran程序中应用了BKW状态方程, 而CHEQ程序中主要应用了液体状态方程的微扰理论, 并且还考虑了碳的石墨-金刚石-液碳三相组成。
从20世纪60年代至今, 随着计算机科学技术的发展, 中国学者在爆轰产物状态方程研究领域也作了大量的研究。X.Wu[4]利用BKW状态方程模型, 从热力学着手, 确定各种可能混合物的状态方程, 计算混合物的自由能, 找出具有最小自由能的组成, 最终求出爆轰CJ点爆速、爆压及爆轰温度, 但计算爆轰温度和实验结果差别比较大, 而爆速的误差也达到500 m/s。李德华等[5]利用WAC状态方程作为爆轰气相产物的物态方程, 考虑化学平衡条件和化学计算条件, 基于Gibbs自由能最小原理, 编程对化学平衡状态的炸药进行数值模拟计算, 对几种单质炸药的爆轰参数作预言, 得到的爆轰CJ点的爆速与实验值得误差小于3%、爆压与实验值的误差小于4%、爆轰温度与实验值的误差小于5%。赵艳红等[6]采用van der Waals等效单组分流体模型和Ross硬球微扰理论软球修正模型计算爆轰气相产物的状态方程; 用石墨相、金刚石相、类石墨液相和类金刚石液相4种相态描述凝聚成分, 由Gibbs自由能最小确定不同状态下的凝聚产物相态, 对爆轰产物混合系统采用Gibbs自由能最小原理, 通过化学平衡方程组求解炸药爆轰产物系统的平衡组分, 计算结果与BKW和LJD的结果相近。
虽然针对爆轰产物组成的研究和编程很多, 但很少有人公开自己的程序, 且用基于最小自由能原理建立的非线性方程组研究爆轰产物和爆轰参数时, 需要计算出各产物的最小自由能, 求解过程复杂。本文中严格按照数学推导过程, 自行设计更易于求解的理论计算模型求解多元非线性方程, 选用BKW状态方程, 依据碳的相图和文献[8]选定碳的4种单质相态, 将一般爆轰压力(0≤p≤600 GPa)和温度(300 K≤T≤15 000 K)下游离态的碳当作金刚石处理, 采用自编程序对爆轰产物组成和爆轰参数进行计算。
1. 数值分析求解最小自由能组分
1.1 BKW状态方程
BKW状态方程最初是由Becker提出的方程, 后来经过Kistiakowsky和Wilson等的修正, 最后确定为以下形式:
pVmRT=σ(X)=1+XeβX (1) 式中:p为压力; Vm为爆轰气体的摩尔体积; R为理想气体常数; T为爆轰气体的热力学温度; X=κ∑xiki/[Vm(T+θ)α], κ、α和θ为经验常数, ki为第i种物质的余容, xi为第i种物质的摩尔分数。BKW状态方程4个常数α、β、κ和θ的值通常如表 1所示, 本文中对欧洲民用爆炸物品系列标准[9]取值进行修正, 修正前α、β、θ和κ对应的取值分别是0.5、0.298、6 620和10.50。通过自编程序计算发现, 若κ取9.50, 计算的爆轰参数与实验值更接近, 本文中采用欧洲标准修正后的取值作为BKW状态方程参数, 如表 1所示, 然后计算爆轰参数关系。
表 1 BKW状态方程参数取值Table 1. The parameter values of BKW equition of state炸药 α β θ κ RDX 0.5 0.16 400 10.91 工业炸药 0.5 0.16 2 890 10.91 本文中 0.5 0.298 6 620 9.50 1.2 求解爆轰产物组分的推导过程
以往的研究人员根据化学平衡方程组建立的非线性方程组, 根据化学平衡条件推导出爆炸产物组成与自由能关系。在假定压力和温度下由状态方程或者高压状态下分子作用模型确定假定爆轰产物组分总自用能[10-11], 根据求得假定组分自由能求出一组可能的爆轰产物组成, 由热力学关系确定假定爆轰组分的压力和温度, 再利用求出的可能爆轰产物组成求出爆轰组分总自由能, 如此利用程序反复迭代, 直到求得爆轰产物组分总自由能、压力和温度在迭代程序中达到一定精度, 最终得到爆轰产物组成、压力和温度。
本文中利用化学平衡条件和数学转化推导出非线性方程组, 根据质量守恒和物质自由能公式可得:
N0k=Nk (2) m∑i=1aiknig+n∑j=1ajknjs=Nkk=1,2,⋯,l (3) G=m∑i=1ginig+n∑j=1gjnjs (4) 式中:Nk0为每千克炸药配方含有第k种元素的总物质的量, Nk为每千克爆炸产物含有第k种元素的总物质的量, aik为每摩尔第i种气相组分含有k元素的物质的量, nig为每千克爆炸产物含有第i种气相组分的物质的量, ajk为每摩尔第j种凝聚相组分含有k元素的物质的量, njs为每千克爆炸产物含有第j种凝聚相组分的物质的量, G为爆炸产物总自由能, gi为每摩尔i气相组分的吉布斯自由能, gj为每摩尔j凝聚相组分的吉布斯自由能。
利用拉格朗日乘数法把式(2)~(4)的条件极值转化为非条件极值, 根据拉乌尔定律和牛顿迭代法处理得到以下公式:
x(ν+1)=x(ν)−F′(x(ν))−1F(x(ν)) (5) F′(x)=[∂f1(x)∂n1g⋯∂f1(x)∂nmg∂f1(x)∂n1s⋯∂f1(x)∂nns∂f1(x)∂λ1⋯∂f1(x)∂λl∂f2(x)∂n1g⋯∂f2(x)∂nmg∂f2(x)∂n1s⋯∂f2(x)∂nns∂f2(x)∂λ1⋯∂f2(x)∂λl⋮⋮⋮⋮⋮⋮∂fm+n+l(x)∂n1g⋯∂fm+n+l(x)∂nmg∂fm+n+l(x)∂n1s⋯∂fm+n+l(x)∂nns∂fm+n+l(x)∂λ1⋯∂fm+n+l(x)∂λl] (6) F=[∂G(n,λk)∂n1g,⋯,∂G(n,λk)∂nig,⋯,∂G(n,λk)∂nmg,∂G(n,λk)∂n1s,⋯,∂G(n,λk)∂njs,⋯,∂G(n,λk)∂nns,∂G(n,λk)∂λ1,⋯,∂G(n,λk)∂λk,⋯,∂G(n,λk)∂λl]T (7) 式中:x=(n1g, …, nig, …, nmg, n1s, …, njs, …, nns, λ1, …, λk, …, λl)T, λk为拉格朗日因子; fi(n, λk)=∂G(n, λk)/∂nig, i=1, …, m; fj+m(n, λk)=∂G(n, λk)/∂njs, j=1, …, n; fk+m+n(n, λk)=∂G(n, λk)/∂λk, k=1, …, l; x(ν)代表第ν次迭代的爆轰组分物质量;
G(n,λk)=m∑i=1(g0iRT+lnp+lnnigng)nig+n∑j=1g0jRTnjs+l∑k=1λk(Nk−m∑i=1aiknig−n∑j=1ajknjs) (8) gi0和gj0分别为第i种气体和第j种凝聚相组分的吉布斯自由能。
利用数学的方法设计了计算爆轰产物组成的方法, 利用数学方法计算爆炸产物组成, 不需要考虑单个组分的自由能。前期学者研究爆轰产物组成考虑单个组分的自由能, 由于分子的作用势和自由能转换计算复杂, 在一定程度上会给计算带来误差。
2. 爆轰参数的确定
由能量守恒法则知爆炸前后系统能量守恒可得[9]:
E0=E=Eg+Es (9) 式中:E0为炸药内能, E为爆炸产物内能, Eg为爆炸产物气体内能, Es为爆炸产物固体内能。
{Eg=m∑i=1nigE0i+EimpEs=n∑j=1njsE0j (10) 式中: Ei0为第i种气体的内能, Eimp为BKW状态方程对应的第i种气体的非理想条件项, Ej0为第j种固体的内能。根据BKW状态方程, 可以推导出Eimp表达式[9]:
Eimp=ngRTαTT+θ(σ−1) (11) 炸药爆轰气相产物分子体系必须满足Hugniot关系[5-7]:
E1−E0=(pH+p0)(V0−VH)/2 (12) 式中:E1为爆轰产物内能, pH爆轰产物CJ条件的爆轰压力, p0为初始压力, VH为爆轰产物CJ条件下的体积, V0为初始体积。爆速表达式为:
D=V0√(pH−p0)/(V0−VH) (13) 由于CJ点的爆速是最小的, 利用抛物线最小法由Hugoniot曲线上的3个点便可求得爆轰CJ点的爆速DCJ。
3. PETN的计算结果
为了验证本文计算方法的准确性, 采用具有实验数据的密度为1.77 g/cm3的PETN进行验证。本文中依据纯物质热化学数据手册[12], 根据(1)~(13)式, 自编程序求得PETN在CJ爆轰点的爆轰参数。密度为1.77 g/cm3、初始摩尔内能为-399.572 kJ/mol的PETN炸药的计算结果见表 2和表 3。本文中得到的压力、爆速和温度的误差分别为0.01%、0.51%和0.21%。利用本文方法还能确定密度为1.77 g/cm3的PETN在CJ爆轰点的产物密度为2.43 g/cm3, 与以往学者认为的爆轰产物密度范围2~3 g/cm3相符[6]。
表 2 每摩尔PETN炸药CJ点爆轰产物组成Table 2. Compositions of detonation products per mole of PETN explosive at CJ point方法 n(CO2)/mol n(CO)/mol n(N2)/mol n(O2)/mol n(H2)/mol BKW[1, 2] 3.950 9 0.096 2 1.999 2 1.359 1×10-5 0.568 9×10-4 LJD[1, 2] 3.779 0 0.446 9 1.969 5 0 0.480 4×10-1 WCA[13] 3.934 9 0.154 6 1.991 7 4.494 3×10-5 3.989 5×10-4 本文 3.999 5 0.001 1 1.999 9 1.190 4×10-5 2.077 7×10-4 方法 n(H2O)/mol n(NH3)/mol n(NO)/mol n(CH4)/mol n(C)/mol BKW[1, 2] 3.998 6 1.327 5×10-4 2.244 1×10-4 0.537 3×10-5 0.951 6 LJD[1, 2] 3.929 1 0 0.619 5×10-1 1.197 9×10-2 0.759 8 WCA[13] 3.975 2 1.629 0×10-2 3.141 7×10-4 0 0.910 5 本文 3.999 7 7.354 8×10-8 4.530 7×10-6 5.482 0×10-6 0.999 1 4. 结论
(1) 本文中基于最小自由能原理, 自行设计理论计算模型, 求解爆轰产物和爆轰参数, 选取BKW方程作为产物状态方程, 状态方程参数采用欧洲标准取值, 并对欧洲标准参数取值进行修正。利用数值分析的方法, 推导求解爆轰产物的方程组, 在计算爆轰产物化学平衡组成时不需要单独计算每个产物组分的自由能, 具有优化计算减小误差的优点, 且容易编程实现。
(2) 根据自行推导的计算方法, 自编程序实现了爆轰参数计算, 用密度为1.77 g/cm3的PETN验证计算方法的准确性, 所得CJ点对应爆压、爆温和爆速与实验值的误差均小于1%, 并得到密度为1.77 g/cm3的PETN在CJ爆轰点的产物密度为2.43 g/cm3。
-
-
[1] MONAGHAN J J. Simulating free surface flows with SPH[J]. Journal of Computational Physics, 1994, 110:399-406. DOI: 10.1006/jcph.1994.1034. [2] 杨秀峰, 刘谋斌.SPH方法Delaunay三角刨分与自由液面重构[J].计算力学学报, 2016, 33(4):594-598. DOI: 10.7511/jslx201604027.YANG Xiufeng, LIU Moubin. Delaunay triangulation and free surface extraction for SPH method[J]. Chinese Journal of Computational Mechanics, 2016, 33(4):594-598. DOI: 10.7511/jslx201604027. [3] 龙厅, 胡德安, 韩旭.FE-ISPH与FE-WCSPH模拟流固耦合问题的比较研究[C]//中国计算力学大会.贵阳, 2014: 547. http://cpfd.cnki.com.cn/Article/CPFDTOTAL-AGLU201408004093.htm [4] 刘谋斌, 宗智, 常建忠.光滑粒子动力学方法的发展与应用[J].力学进展, 2011, 41(2):219-236. DOI: 10.6052/1000-0992-2011-2-lxjzJ2010-078.LIU Moubin, ZONG Zhi, CHANG Jianzhong. Developments and applications of smoothed particle hydrodynamics[J]. Advances in Mechanics, 2011, 41(2):219-236. DOI: 10.6052/1000-0992-2011-2-lxjzJ2010-078. [5] 傅学金, 强洪夫, 杨月诚.固体介质中SPH方法的拉伸不稳定性问题研究进展[J].力学进展, 2007, 37(3):375-388. DOI: 10.3321/j.issn:1000-0992.2007.03.005.FU Xuejin, QIANG Hongfu, YANG Yuecheng. Advances in the tensile instability of smoothed particle hydrodynamics applied to solid dynamics[J]. Advances in Mechanics, 2007, 37(3):375-388. DOI: 10.3321/j.issn:1000-0992.2007.03.005. [6] LIU W K, JUN S, LI S, et al. Reproducing kernel particle methods for structure dynamics[J]. International Journal for Numerical Methods in Engineering, 1995, 38(10):1655-1679. DOI: 10.1002/nme.1620381005. [7] CHEN J K, BERAUN J E. A generalized smoothed particle hydrodynamics method for nonlinear dynamic problem[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 190:225-239. DOI: 10.1016/S0045-7825(99)00422-3. [8] 章杰, 苏少卿, 郑宇, 等.改进SPH方法在陶瓷材料层裂数值模拟中的应用[J].爆炸与冲击, 2013, 33(4):401-407. DOI: 10.3969/j.issn.1001-1455.2013.04.011.ZHANG Jie, SU Shaoqing, ZHENG Yu, et al. Application of modified SPH method to numerical simulation of ceramic spallation[J]. Explosion and Shock Waves, 2013, 33(4):401-407. DOI: 10.3969/j.issn.1001-1455.2013.04.011. [9] LIU M B, LIU G R. Restoring particle consistency in smoothed particle hydrodynamics[J]. Applied Numerical Mathematics, 2006, 56(1):19-36. DOI: 10.1016/j.apnum.2005.02.012. [10] 郑兴, 段文洋.K2_SPH方法及其对二维非线性水波的模拟[J].计算物理, 2011, 28(5):659-666. DOI: 10.3969/j.issn.1001-246X.2011.05.004.ZHENG Xing, DUAN Wenyang. K2_SPH Method and application for 2D nonlinear water wave simulation[J]. Chinese Journal of Computational Physics, 2011, 28(5):659-666. DOI: 10.3969/j.issn.1001-246X.2011.05.004. [11] 刘谋斌, 杨秀峰, 邵家儒.高精度SPH方法及其在海洋工程中的应用[C]//颗粒材料计算力学会议论文集.兰州, 2014: 39-41. http://cpfd.cnki.com.cn/Article/CPFDTOTAL-AGLU201408001007.htm [12] ADAMI S, HU X Y, ADAMS N A. A generalized wall boundary condition for smoothed particle hydrodynamics[J]. Journal of Computational Physics, 2012, 231(21):7057-7075. DOI: 10.1016/j.jcp.2012.05.005. [13] LIU M B, SHAO J R, CHANG J Z. On the treatment of solid boundary in smoothed particle hydrodynamics[J]. Science China:Technological Sciences, 2012, 55(1):244-254. DOI: 10.1007/s11431-011-4663-y. [14] LIU M B, LIU G R, LAM K Y. A one-dimensional meshfree particle formulation for simulating shock waves[J]. Shock Wave, 2003, 13:201-211. DOI: 10.1007/s00193-003-0207-0. [15] XU F, ZHAO Y, YAN R, et al. Multi-dimensional discontinuous SPH method and its application to metal penetration analysis[J]. International Journal for Numerical Methods in Engineering, 2013, 93:1125-1146. DOI: 10.1002/nme.4414. [16] 闫蕊, 徐绯, 张岳青.DSPH方法的有效性验证及应用[J].爆炸与冲击, 2013, 33(2):133-139. DOI: 10.3969/j.issn.1001-1455.2013.02.004.YAN Rui, XU Fei, ZHANG Yueqing. Validation of DSPH method and its application to physical problems[J]. Explosion and Shock Waves, 2013, 33(2):133-139. DOI: 10.3969/j.issn.1001-1455.2013.02.004. [17] 宋俊豪, 张超英, 梁朝湘, 等.RDSPH:一种适用于一维非连续条件的新SPH方法[J].广西师范大学学报(自然科学版), 2009, 27(3):9-13. DOI: 10.3969/j.issn.1001-6600.2009.03.003.SONG Junhao, ZHANG Chaoying, LIANG Chaoxiang, et al. A new one-dimensional smoothed particle hydrodynamics method in simulating discontinuous problem[J], Journal of Guangxi Normal University (Natural Science Edition), 2009, 27(3):9-13. DOI: 10.3969/j.issn.1001-6600.2009.03.003. [18] YANG Yang, XU Fei, ZHANG Meng, et al. An effective improved algorithm for finite particle method[J]. International Journal of Computational Methods, 2016, 13(4):1641009. DOI: 10.1142/S0219876216410097. [19] MONAGHAN J J. Smoothed particle hydrodynamic[J]. Annual Review of Astronomy and Astrophysics, 1992, 30(1):543-574. DOI: 10.1146/annurev.aa.30.090192.002551. [20] MONAGHAN J J, KAJTAR J. SPH particle boundary forces for arbitrary boundaries[J]. Computer Physics Communications, 2009, 180(10):1811-1820. DOI: 10.1016/j.cpc.2009.05.008. 期刊类型引用(6)
1. 伍俊英,姚雨乐,郑富德,刘嘉锡,李钧剑,王健宇,陈朗. 多脉冲飞秒激光切割加工炸药装药过程的热响应. 兵工学报. 2024(10): 3462-3473 . 百度学术
2. 吴俊浩,郭子如,杜宝强,刘伟,张金元. 炸药爆炸产物组成的确定及其对爆热的影响. 工程爆破. 2022(03): 91-96 . 百度学术
3. 李根,卢芳云,李翔宇,李志斌. 基于气固两相反应流的温压炸药能量释放规律数值模拟及实验验证. 火炸药学报. 2021(02): 195-204 . 百度学术
4. 刘飞,杨超志,夏明,贾鑫,汪剑辉. 钢筋混凝土板爆炸动态响应研究进展. 防护工程. 2020(05): 1-9 . 百度学术
5. 杜明燃,张金龙,梁琳娜,黄亮亮,韩体飞,刘锋. 炸药爆轰热力学数据非线性拟合方法研究. 爆破器材. 2019(03): 38-43 . 百度学术
6. 何伟平,黄菊,陈厚和,刘晓静,王德堂. 基于BKW状态方程的爆轰产物及参数的改进算法. 火炸药学报. 2017(03): 53-59 . 百度学术
其他类型引用(5)
-