CN112069745A - 一种固体推进剂废料切割处理的数值仿真方法及系统 - Google Patents

一种固体推进剂废料切割处理的数值仿真方法及系统 Download PDF

Info

Publication number
CN112069745A
CN112069745A CN202010944304.9A CN202010944304A CN112069745A CN 112069745 A CN112069745 A CN 112069745A CN 202010944304 A CN202010944304 A CN 202010944304A CN 112069745 A CN112069745 A CN 112069745A
Authority
CN
China
Prior art keywords
solid propellant
solid
particle
phase
establishing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010944304.9A
Other languages
English (en)
Other versions
CN112069745B (zh
Inventor
陈福振
严红
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202010944304.9A priority Critical patent/CN112069745B/zh
Publication of CN112069745A publication Critical patent/CN112069745A/zh
Application granted granted Critical
Publication of CN112069745B publication Critical patent/CN112069745B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Perforating, Stamping-Out Or Severing By Means Other Than Cutting (AREA)

Abstract

本发明涉及一种固体推进剂废料切割处理的数值仿真方法及系统。该方法包括：建立射流喷嘴几何模型和固体推进剂几何模型；建立前混合磨料水射流的液体‑固体颗粒两相流物理模型；建立固体推进剂的动力学控制方程；选取磨料、水、固体推进剂三种材料参数；采用SPH方法对所述前混合磨料水射流的液体‑固体颗粒两相流物理模型进行离散；采用SPH方法对所述固体推进剂的动力学控制方程进行离散；建立时间积分格式；根据所述射流喷嘴几何模型、固体推进剂几何模型、离散后的两相流物理模型、离散后的固体推进剂的动力学控制方程采用时间积分格式进行数值仿真计算，得到不同时刻的场变量。本发明能够迅速精确的确定场变量。

Description

一种固体推进剂废料切割处理的数值仿真方法及系统
技术领域
本发明涉及固体推进剂废料切割处理技术领域，特别是涉及一种固体推进剂废料切割处理的数值仿真方法及系统。
背景技术
在军工技术作业中，对固体火箭发动机的“处废”是一项难度极大、技术含量极高、安全系数较低的工作。如何能够有效的将发动机中的固体推进剂去除，同时又保障一定的安全系数，实现发动机壳体的有效回收利用，是科研和技术人员解决的难题之一。传统或者采用刀具切割的方法，需要人工近距离的作业，不仅推进剂的药物伤害身体，同时推进剂遭受机械作用后极易引发燃烧、爆炸等事故；或者采用高压水射流的方式进行切割，虽然降低了切割温度，但是该技术切割效率极低，易导致气溶胶集聚，增加次生灾害发生的几率。而前混合磨料水射流在此背景下逐渐发展起来，作为一种全新的现代切割技术，通过将高压水与磨料事先混合形成液固两相高能束流进行固体推进剂的切割，既有效降低了出口压力，又提高了切割过程的安全系数，是一种极具潜力的切割技术。
前混合磨料水射流切割固体推进剂废料过程中，磨料的配比如何调整、喷嘴的尺寸如何优化、出口的压力多少最为合适、固体推进剂的切割断面是什么样子的等等问题均需要深入研究。研究方法常用的主要有两种：一种是基于实验的直接观测的方法，该方法直接通过搭建实验台，采用高速摄影、粒子图像测速、激光纹影等方法，对磨料射流切割固体推进剂的全过程进行光学测量，获得切割过程中的细节以及切割的效果。该方法的主要缺点是：(1)需要专门的场地开展实验；(2)需要占用大量的人力、物力和财力；(3)实验的周期长，同时会经常出现失败情况，需要反复进行实验，进一步增加成本；(4)测试过程中很多不确定性因素无法控制，因此获得的实验结果有时与真实过程存在一定的差别。另一种是理论研究的方法，该方法通过分析磨料射流形成的机理、固体推进剂损伤破坏的机理、流动与传热的机理等建立理论预测模型，对改进优化后的结果进行预估。该方法存在的主要缺点是：(1)机理的分析依赖于实验的结果，实验无法捕捉到的现象，在进行机理分析的过程中就会存在盲点，影响理论模型的准确性；(2)理论模型中通常包含有很多人工参量，预测结果的精度与这些人工参量息息相关，影响结果的客观性；(3)通过理论模型只能获得最终的结果和结论，无法获得动力学过程中的细节，无法动态捕捉切割过程中的典型现象，预估的结果与实际经常存在较大的差别。
发明内容
本发明的目的是提供一种固体推进剂废料切割处理的数值仿真方法及系统，能够迅速精确的确定场变量。
为实现上述目的，本发明提供了如下方案：
一种固体推进剂废料切割处理的数值仿真方法，包括：
建立射流喷嘴几何模型和固体推进剂几何模型；
建立前混合磨料水射流的液体-固体颗粒两相流物理模型；
建立固体推进剂的动力学控制方程；
选取磨料、水、固体推进剂三种材料参数；
采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散；
采用SPH方法对所述固体推进剂的动力学控制方程进行离散；
建立时间积分格式；
根据所述射流喷嘴几何模型、固体推进剂几何模型、离散后的两相流物理模型、离散后的固体推进剂的动力学控制方程采用时间积分格式进行数值仿真计算，得到不同时刻的场变量。
可选地，所述建立前混合磨料水射流的液体-固体颗粒两相流物理模型，具体包括：
建立前混合磨料水射流的液体相控制方程组：
Figure BDA0002674711640000021
Figure BDA0002674711640000031
Figure BDA0002674711640000032
其中，下标l表示的是液体相的标志，与下面的固体颗粒相区分开；αl,ρl和 vl分别为液体的体积分数、密度和速度，
Figure BDA0002674711640000033
表示偏导数，
Figure BDA0002674711640000034
是指对时间的偏导数，Pl表示的是液体相的正压力，
Figure BDA0002674711640000035
表示的是梯度，
Figure BDA0002674711640000036
则表示的是散度，g为重力加速度，Rlp为液体相和固体颗粒相之间的相互作用力，液体的粘性作用对于颗粒的作用不可以忽略，τl为液体相的剪切应力张量；
建立前混合磨料水射流的固体颗粒相控制方程组：
Figure BDA0002674711640000037
Figure BDA0002674711640000038
Figure BDA0002674711640000039
Figure BDA00026747116400000310
其中，下标p表示的是固体颗粒相的标志，αp,ρp和vp分别为颗粒相体积分数、密度和速度，
Figure BDA00026747116400000311
为液体连续相压力梯度，
Figure BDA00026747116400000312
为固体颗粒相压力梯度，αpρpg为颗粒所受外部体积力，Rpl为相间相互作用力，hp为固体颗粒相的能量焓值及qp为固体颗粒相内部的热传导量，τp为颗粒相粘性应力张量。
可选地，所述建立固体推进剂的动力学控制方程，具体包括：
建立固体推进剂的动力学控制方程：
Figure BDA00026747116400000313
Figure BDA0002674711640000041
Figure BDA0002674711640000042
其中，下标s是固体的表征，ρs为固体推进剂的密度，vs α为固体推进剂的速度在α方向的分量，α取x、y以及z方向，x为位移，σs αβ为固体推进剂的全应力分量，
Figure BDA0002674711640000043
ps为固体推进剂所受到的正应力，
Figure BDA0002674711640000044
为固体推进剂所受到的剪应力，fs α为固体推进剂所受到的其他外部作用力，es为固体推进剂内能。
可选地，所述选取磨料、水、固体推进剂三种材料参数，具体包括：
设置水的材料参数为密度ρl＝998Kg/m3，粘度μl＝0.001003Kg/ms；
选取金刚砂作为磨料材料，金刚砂密度ρp＝3900Kg/m3，金刚砂粒径大小为80目，金刚砂与水的混合比1:10；
选取HTPB推进剂作为固体推进剂。
可选地，所述采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散，具体包括：
采用SPH方法对所述液体相控制方程组进行离散，得到离散后的液体相控制方程组：
Figure BDA0002674711640000045
Figure BDA0002674711640000046
Figure BDA0002674711640000047
其中，i,j分别是指的i粒子和j粒子，Wij为i粒子和j粒子之间的核函数的数值，W为核函数，h为光滑长度；
采用SPH方法对所述固体颗粒相控制方程组进行离散，得到离散后的固体颗粒相控制方程组：
Figure BDA0002674711640000051
Figure BDA0002674711640000052
Figure BDA0002674711640000053
Figure BDA0002674711640000054
其中，ρp,i为SPH粒子i的密度；ρp为颗粒的实际密度；vp,ij为速度矢量， vp,ij＝vp,i-vp,j，rp,ij为位移矢量，rp,ij＝rp,i-rp,j，Πij为粒子i和粒子j之间的人工粘性。
可选地，所述采用SPH方法对所述固体推进剂的动力学控制方程进行离散，具体包括：
采用SPH方法对所述固体推进剂的动力学控制方程进行离散，得到离散后的固体推进剂的动力学控制方程：
Figure BDA0002674711640000055
Figure BDA0002674711640000056
Figure BDA0002674711640000057
其中，
Figure BDA0002674711640000058
为应变率。
一种固体推进剂废料切割处理的数值仿真系统，包括：
几何模型建立模块，用于建立射流喷嘴几何模型和固体推进剂几何模型；
液体-固体颗粒两相流物理模型建立模块，用于建立前混合磨料水射流的液体-固体颗粒两相流物理模型；
动力学控制方程建立模块，用于建立固体推进剂的动力学控制方程；
材料参数确定模块，用于选取磨料、水、固体推进剂三种材料参数；
液体-固体颗粒两相流物理模型离散模块，用于采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散；
动力学控制方程离散模块，用于采用SPH方法对所述固体推进剂的动力学控制方程进行离散；
时间积分格式建立模块，用于建立时间积分格式；
不同时刻的场变量确定模块，用于根据所述射流喷嘴几何模型、固体推进剂几何模型、离散后的两相流物理模型、离散后的固体推进剂的动力学控制方程采用时间积分格式进行数值仿真计算，得到不同时刻的场变量。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明一方面在相比传统的实验研究和理论研究方面具有优势：仅需要电子计算机便可完成计算，无需实验开展所需要的实验台、光学测量装置、磨料材料、水、推进剂方坯等物品，因此大大减少人力、物力和财力的消耗，同时可以重复进行计算，不影响计算的结果，可以清晰的捕获切割过程中的每一个细节，是对开展实际实验的较好地补充；另外，本发明是从最本质的物理过程出发进行数值模拟，再现的是实际动力学过程中的各个细节，克服了理论预测将中间过程当成黑盒的不足，新技术不仅可以预测最终的切割性能，同时还可以深入揭示切割过程机理，改进理论预测模型，为高精度的理论预测提供支撑。另一方面在相比同类型的数值仿真方法上面具有优势：突破了传统界面追踪技术必须采用网格进行模拟的现状，网格模拟的过程中必须结合网格自适应才能捕获射流状态，而无网格粒子方法无需网格，更无需网格自适应，不仅对于液体水可以进行追踪，同时对于磨料颗粒来说也可以自然追踪；另外，新技术克服了传统拉格朗日网格法在计算固体变形破坏问题上的缺陷，不存在网格的扭曲和缠绕，可以自然地模拟固体材料的大变形、破坏和损伤。因此，本发明针对固体推进剂废料切割处理中同时存在的液体-颗粒两相流、磨料-水-固体推进剂多介质、固体材料大变形与损伤三个典型难题，均具有很好的解决，是一项非常具有实用价值的技术手段。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明固体推进剂废料切割处理的数值仿真方法流程图；
图2为本发明前混合磨料水射流装置喷嘴结构的粒子模型示意图；
图3为本发明固体推进剂的粒子模型示意图；
图4为本发明不同时刻磨料水射流混合物喷射过程示意图；
图5为本发明不同时刻磨料水射流混合物喷射过程示意图；
图6为本发明固体推进剂方坯冲击损坏过程与实验结果的对比图；
图7为本发明不同射流压力下切割时间变化曲线示意图；
图8为本发明不同射流压力下切割深度变化曲线示意图；
图9为本发明固体推进剂废料切割处理的数值仿真系统结构图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种固体推进剂废料切割处理的数值仿真方法及系统，能够迅速精确的确定场变量。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明固体推进剂废料切割处理的数值仿真方法流程图。如图1所示，一种固体推进剂废料切割处理的数值仿真方法包括：
步骤101：建立射流喷嘴几何模型和固体推进剂几何模型。
射流喷嘴几何模型和固体推进剂几何模型的建立均采用商业软件完成，通过商业软件建立三维几何模型，再导入网格划分软件，进行细致均匀的网格划分，最后网格文件导入到程序中按照一个网格转化成一个粒子的转化原则，全部转化为粒子，具体过程为：
1)首先建立实际前混合磨料水射流喷嘴几何模型和固体推进剂几何模型：采用达索公司旗下的SolidWorks子公司负责开发的三维软件SolidWorks软件完成。
2)在步骤1)的基础上，对于步骤1)建立的几何模型进行网格划分：对于模型网格生成，采用功能强大的CAE应用软件包—Hypermesh软件完成。
3)将步骤2)划分后形成的网格文件导入到自编程序中进行网格到粒子的转化：按照一个网格对应一个粒子的原则，采用任意六面体体积计算方法(现有技术)计算六面体网格的体积即为SPH粒子的体积，采用任意六面体质心计算方法(现有技术)计算六面体网格的质心即为SPH粒子的质心也就是SPH 粒子的初始位置，这样就获得了结构的粒子模型。
SPH单个粒子的体积决定了SPH单个粒子的质量，在下面步骤105中会涉及，是算法计算的核心数据，SPH粒子的质心，也就是SPH粒子的初始位置，直接决定了物质的初始位置。图2为本发明前混合磨料水射流装置喷嘴结构的粒子模型示意图，图3为本发明固体推进剂的粒子模型示意图。其中喷嘴长度6cm，孔径1mm，固体推进剂方坯尺寸为：20cm*10cm*10cm。
步骤102：建立前混合磨料水射流的液体-固体颗粒两相流物理模型。
前混合磨料水射流中的水为液体相，磨料为固体颗粒相，根据两相混合性质，基于双流体模型(现有技术)，建立描述前混合磨料水射流两相流动特性的物理模型，公式如下：
液体相控制方程组：
液体相为连续相，物理模型共包括质量守恒方程、动量守恒方程、能量守恒方程三个，如下
Figure BDA0002674711640000091
Figure BDA0002674711640000092
Figure BDA0002674711640000093
其中，下标l表示的是液体相的标志，与下面的固体颗粒相区分开；αl,ρl和 vl分别为液体的体积分数、密度和速度，
Figure BDA0002674711640000094
表示偏导数，
Figure BDA0002674711640000095
是指对时间的偏导数，Pl表示的是液体相的正压力，
Figure BDA0002674711640000096
表示的是梯度，
Figure BDA0002674711640000097
则表示的是散度，g为重力加速度，Rlp为液体相和固体颗粒相之间的相互作用力。液体的粘性作用对于颗粒的作用不可以忽略，τl为液体相的剪切应力张量，用下式计算
Figure BDA0002674711640000098
式中，
Figure BDA0002674711640000099
为液体的剪切速率张量，
Figure BDA00026747116400000910
为液体的广义粘性，K为稠度系数，n为流动指数。对于牛顿流体，n＝1，
Figure BDA00026747116400000911
对于幂律型流体，n＜1，
Figure BDA00026747116400000912
随剪切速率
Figure BDA00026747116400000913
的增大而减小。本发明中水为牛顿流体，则n＝1，
Figure BDA00026747116400000914
将公式(4)及剪切速率表达式
Figure BDA00026747116400000915
代入剪切力引起的加速度公式
Figure BDA00026747116400000916
中可得，Fl (v)为由液体剪切力引起的液体相的加速度值
Figure BDA00026747116400000917
公式(3)中的hl为液体的能量焓值及ql为液体相内部的热传导量，公式为
Figure BDA00026747116400000918
Tref代表参考温度值，T为现在时刻液体的温度值，cp,i为液体的比热容，kl为液体相内部之间的热传导系数，ε为液体相和固体颗粒相之间的热传导系数，Tp和Tl分别为固体颗粒相和液体相的温度。
固体颗粒相控制方程组：
固体颗粒相为离散相，物理模型共包括质量守恒方程、动量守恒方程、能量守恒方程以及拟温度守恒方程四个，如下
Figure BDA0002674711640000101
Figure BDA0002674711640000102
Figure BDA0002674711640000103
Figure BDA0002674711640000104
其中，下标p表示的是固体颗粒相的标志，αp,ρp和vp分别为颗粒相体积分数、密度和速度，
Figure BDA0002674711640000105
为液体连续相压力梯度，
Figure BDA0002674711640000106
为固体颗粒相压力梯度，αpρpg为颗粒所受外部体积力，Rpl为相间相互作用力。hp为固体颗粒相的能量焓值及qp为固体颗粒相内部的热传导量。τp为颗粒相粘性应力张量
Figure BDA0002674711640000107
其中，μp和λp为颗粒相的剪切粘度和体粘度，I为单位张量。
引入对于颗粒相压力Pp和粘性应力τp的描述，从而实现对颗粒相控制方程(8)(9)的封闭。采用颗粒动理学理论获知，颗粒相压力Pp和粘性应力τp与颗粒速度脉动的最大值相关，而颗粒的速度脉动由颗粒拟温度描述，颗粒拟温度守恒方程如式(10)，其中(-ppI+τp):
Figure BDA0002674711640000108
为由颗粒相应力产生的能量，
Figure BDA0002674711640000109
为能量耗散项，kp为能量耗散系数，Ncθp为颗粒间碰撞产生的能量耗散项，具体公式为
Figure BDA00026747116400001010
Figure BDA0002674711640000111
Pp＝αpρp[1+2(1+e)αpg0]θp (14)
Figure BDA0002674711640000112
式中，dp为颗粒的直径，epp为颗粒之间的碰撞恢复系数，ξp为由颗粒碰撞产生的颗粒相有效容积粘度，为中间变量。g0为颗粒的径向恢复系数
Figure BDA0002674711640000113
αp,max为可压缩条件下，颗粒相可达到的最大体积分数值。
公式(10)中，φlp为连续相与颗粒相间的能量交换，一般取
φlp＝-3βlpθp (17)
其中βlp为气体相与颗粒相间的曳力系数，由公式(19)和(23)决定。
液体相和固体颗粒相之间的相互作用力模型：
作用于单颗粒上的曳力可由动量交换系数β和两相间滑移速度vl-vp表示
Rlp＝βlp(vl-vp) (18)
这里的Rlp等同于公式2和公式8中的R。大量研究表明，颗粒相的体积分数对于决定颗粒群运动的曳力来说具有重要的影响。本发明采用Gidaspow 等提出的公式，即对于密相的计算采用Ergun方程以及对于稀相的计算采用 Wen-Yu方程：
Figure BDA0002674711640000121
曳力系数CD为
Figure BDA0002674711640000122
相对雷诺数Rep定义为
Figure BDA0002674711640000123
为消除两个方程间的不连续性，引入松弛因子对过渡区域中的动量交换系数进行光滑处理
Figure BDA0002674711640000124
因此，动量交换系数β可以表示为
Figure BDA0002674711640000125
由此，可得作用于单位质量颗粒上的曳力Rl′p为
Figure BDA0002674711640000126
步骤103：建立固体推进剂的动力学控制方程。
固体推进剂为固体材料结构，其变形与破坏过程需要采用具有材料强度的动力学控制方程进行描述，具体如下
Figure BDA0002674711640000127
Figure BDA0002674711640000128
Figure BDA0002674711640000131
其中，下标s是固体的表征，ρs为固体推进剂的密度，vs α为固体推进剂的速度在α方向的分量，α可以取x、y以及z方向。x为位移。σs αβ为固体推进剂的全应力分量，
Figure BDA0002674711640000132
ps为固体推进剂所受到的正应力，
Figure BDA0002674711640000133
为固体推进剂所受到的剪应力，fs α为固体推进剂所受到的其他外部作用力，es为固体推进剂内能。
公式
Figure BDA00026747116400001316
中的正压力ps采用Mie-Gruneisen状态方程计算
Figure BDA0002674711640000134
其中
Figure BDA0002674711640000135
Figure BDA0002674711640000136
Figure BDA0002674711640000137
式中各系数为Γ＝1.99，CS＝3940，SS＝1.489。
公式
Figure BDA0002674711640000138
中的剪应力
Figure BDA0002674711640000139
的计算采用以下方式，假定应力为应变和应变率的函数，对于各向异性剪切应力，若假设小位移，则应力率与应变率成正比，比例系数为剪切模量
Figure BDA00026747116400001310
G为剪切模量，
Figure BDA00026747116400001311
为应力率偏张量，εαβ为应变率张量，为
Figure BDA00026747116400001312
Figure BDA00026747116400001313
为εαβ的剪切变形部分。
计算获得
Figure BDA00026747116400001314
进而获得
Figure BDA00026747116400001315
之后，计算冯米塞斯等效应力
Figure BDA0002674711640000141
然后再比较冯米塞斯等效应力与固体推进剂的本构模型计算得到的屈服强度σeq大小，若J≤σeq则剪切应力
Figure BDA0002674711640000142
保持不变，若J＞σeq则剪切应力
Figure BDA0002674711640000143
按照比例退回到屈服面上
Figure BDA0002674711640000144
接下就阐述屈服强度σeq的计算：固体推进剂的本构模型采用含损伤的 Johnson-Cook本构模型描述，有效刻画固体推进剂材料的屈服应力及损伤演化，该模型中将材料的屈服强度表示为损伤变量、等效应变、等效应变率和温度的函数：
Figure BDA0002674711640000145
其中A,B,C,n,m是材料常数，D为损伤变量，D＝0表示材料没有损伤，D＝1 表示材料完全失效。r是累积损伤塑性应变，r＝(1-D)ζ，
Figure BDA0002674711640000146
ζ是累积塑性应变，
Figure BDA0002674711640000147
是自定义参考应变率。温度T*＝(T-T0)/(Tm-T0)，T0是室温，Tm是材料熔点。
损伤变量D是累积塑性应变ζ的函数，当D＝1时发生损伤破坏：
D＝∑Δζ/ζf (35)
ζf是断裂塑性应变，与材料的应力三轴度、应变率和温度相关。本构模型中的剪切损伤演化模型将ζf描述如下：
Figure BDA0002674711640000148
其中D1-D5为材料常数，σ*＝σm/σeq为应力三轴度，σm＝(σx+σy+σz)/3为平均正应力。
步骤104：选取磨料、水、固体推进剂三种材料参数，具体包括：
1)设置水的材料参数为密度ρl＝998Kg/m3，粘度μl＝0.001003Kg/ms。
2)选取金刚砂作为磨料材料，金刚砂密度ρp＝3900Kg/m3，金刚砂粒径大小为80目，金刚砂与水的混合比1:10。
3)选取HTPB推进剂作为固体推进剂。
选取了这些初始参数之后，才能采用步骤105和106进行计算。具体的，选择这些参数之后，直接代入到步骤105和106之中的公式中，开始计算，获得密度、速度等数值的变化量，再采用时间积分更新下一时刻的数值。具体的参数如表1所示：
表1参数表
Figure BDA0002674711640000151
步骤105：采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散。
所述采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散，具体包括：
采用SPH方法对所述液体相控制方程组(1)、(2)、(3)进行离散，得到离散后的液体相控制方程组：
Figure BDA0002674711640000152
Figure BDA0002674711640000153
Figure BDA0002674711640000154
其中，i,j分别是指的i粒子和j粒子，其他的符合在前面都已经解释，Wij为i粒子和j粒子之间的核函数的数值，W为核函数，h为光滑长度。
采用SPH方法对所述固体颗粒相控制方程组(7)(8)(9)(10)进行离散，得到离散后的固体颗粒相控制方程组：
Figure BDA0002674711640000161
Figure BDA0002674711640000162
Figure BDA0002674711640000163
Figure BDA0002674711640000164
ρp,i为SPH粒子i的密度(即颗粒相有效密度)；ρp为颗粒的实际密度；速度矢量vp,ij＝vp,i-vp,j，位移矢量rp,ij＝rp,i-rp,j。
Πij为粒子i和粒子j之间的人工粘性，用以下公式计算
Figure BDA0002674711640000165
式中，
Figure BDA0002674711640000166
ε＝0.01用于防止粒子相互靠近时产生的数值发散，常数α在模拟激波时一般设定为1，在本发明中模拟流体和固体结构冲击破坏等问题时，α最小为0.02，可以保证计算稳定。
由于SPH在对颗粒相进行计算时，每个SPH粒子表征一系列具有一定粒径分布的颗粒，因此作用于SPH粒子上单位质量曳力
Figure BDA0002674711640000167
及对流换热量
Figure BDA0002674711640000168
为
Figure BDA0002674711640000169
Figure BDA00026747116400001610
其中
Figure BDA0002674711640000171
Figure BDA0002674711640000172
为作用于颗粒k上的曳力。
Figure BDA0002674711640000173
为作用于颗粒k上的对流换热量。N为 SPH粒子所表征的颗粒的数量。
步骤106：采用SPH方法对所述固体推进剂的动力学控制方程进行离散，具体包括：
采用SPH方法对所述固体推进剂的动力学控制方程(25)(26)(27)进行离散，得到离散后的固体推进剂的动力学控制方程：
Figure BDA0002674711640000174
Figure BDA0002674711640000175
Figure BDA0002674711640000176
其中，
Figure BDA0002674711640000177
为应变率。对于应变率
Figure BDA0002674711640000178
的计算，采用以下SPH离散式
Figure BDA0002674711640000179
公式中的ps和ss计算分别采用如上公式(28)(32)计算。
步骤107：建立时间积分格式。
对于建立的磨料水射流的液体-固体颗粒两相流物理模型以及固体推进剂的考虑材料强度的动力学模型，如公式(37)-(43)和(47)-(49)，左端项中均包含有时间的导数项，求得时间导数项，才可以对场变量ρ、v、e、h、θp以及位移x进行时间积分，从而获得不同时刻的场变量，采用的时间积分格式如下：
Figure BDA00026747116400001710
xi(t+δt)＝xi(t)+vi(t+δt/2)δt (52)
式中
Figure BDA0002674711640000181
表示粒子i的密度ρ，速度v，内能e，焓值h、拟温度值θp；xi为粒子i处的位置坐标。
Figure BDA0002674711640000182
是通过方程组(37)-(43)和(47)-(49)计算获得的dvi(t)和dρi(t)；δt为采用公式(53)-(55)计算获得的时间步长。
对SPH离散方程采用显式时间积分求解，采用蛙跳积分方法，对时间具有二阶精度，并且存储量低，计算效率高。对于蛙跳积分，时间步长必须满足稳定性条件，这里应用柯朗－弗里德里奇－列维(Courant-Friedrich-Lewy，简称CFL)条件对时间步长进行估计，具体表达式为：
Figure BDA0002674711640000183
Figure BDA0002674711640000184
Figure BDA0002674711640000185
其中，f为作用于单位质量上的外力，μ为流体的动力粘度，最终取式(53)～(55)中最小值作为SPH计算的时间步长。αΠβΠ为无量纲参量，可调参数值，φij为粒子i和粒子j之间的粘性。
步骤108：根据所述射流喷嘴几何模型、固体推进剂几何模型、离散后的两相流物理模型、离散后的固体推进剂的动力学控制方程采用时间积分格式进行数值仿真计算，得到不同时刻的场变量。
在以上步骤步骤105-107的基础上，自编程序实现以上算法，然后进行编译，采用OpenMP并行计算的方式进行计算，计算前混合磨料水射流经由射流喷嘴几何模型喷出后形成射流，进一步与固体推进剂几何模型相互作用，沿固体推进剂方坯中心线方向以恒定速度切割，获得不同时间节点上的磨料水混合物的速度v、位移x、密度ρ以及固体推进剂方坯的密度ρ、速度v、损伤变量 D、位移x等场变量数据。
将上述场变量数据结果进行后处理，采用商业软件Tecplot完成，将步骤 108计算获得的数据导入到Tecplot软件中，生成各个时刻磨料水射流和固体推进剂物质在空间中的分布状况图片，给出速度矢量、密度、损伤变量等数据在图中的显示情况，同时选取模型中特定的粒子为研究对象，提取该粒子的相关参量随时间的变化数据，通过Tecplot软件生成该粒子的相关变量的时间历程曲线。
图4为本发明不同时刻磨料水射流混合物喷射过程示意图。图5为本发明不同时刻磨料水射流混合物喷射过程示意图。图6为本发明固体推进剂方坯冲击损坏过程与实验结果的对比图，可以看出两者在射流侵彻形态、推进剂破坏形态上都吻合较好，验证了本发明在计算固体推进剂废料切割处理问题上的有效性。
图7为本发明不同射流压力下切割时间变化曲线示意图，随着射流压力的增加，切割时间逐渐降低，因为在相同的切割深度的情况下，压力越大，射流速度越大，切割相同深度时间越短，所以增加了横向移动的时间，减小了切割的时间。图8为本发明不同射流压力下切割深度变化曲线示意图，随着射流压力的增大，在相同的横向移动速度下，切割深度会增大，因为射流压力增加，则射流速度增大，射流的动能增加，切割深度增大。
通过对不同时刻固体推进剂被切割处理的细节以及不同参数影响下切割特性的变化曲线进行分析，深入揭示磨料水射流形成的物理机理、推进剂切割破坏机理、磨料射流切割过程中的飞溅运动机理、推进剂断面形成机理等，获得切割规律，进一步改进和完善理论预测模型；另一方面，可直接在改进优化磨料射流配比参数、喷嘴几何参数、喷注压力等基础上，进行模拟仿真，获得这些改进优化措施的有效性，最终使得磨料水射流切割性能达到最优。
传统对于磨料射流切割固体推进剂这种特种材料性能预测大多采用实验和理论预测的方式，实验方法费时费力，同时受实验测量技术的制约很多切割细节无法捕获；理论预测仅能获得最终切割特性与初始参数之间的关系，无法有效揭示磨料射流切割过程机理；而本发明中采用数值模拟可以较好解决实验和理论计算存在的不足，但传统数值模拟大多采用基于网格的数值模拟技术，对于本发明所针对的问题中不仅存在磨料与水这种混合多相问题，同时存在固体结构受多相流切割而发生大变形破坏问题，这种涉及复杂多介质和大变形问题，传统网格数值模拟技术无能为力，本发明独辟蹊径引入无网格粒子仿真技术很好的解决该问题。
本发明提出了一种计算液体-固体颗粒两相流动问题的新思路，新思路基于无网格粒子仿真方法，考虑液相的界面演化，考虑两相之间的相互耦合，考虑固体材料的大变形损伤破坏等因素，计算过程中无需进行网格重画，计算量得到了控制，同时该技术属于拉格朗日方法，对于多相流和多介质的界面追踪比较自然，无需引入特别的界面追踪技术，为解决含液体-颗粒多相流和材料损伤破坏问题提供了新的途径。
图9为本发明固体推进剂废料切割处理的数值仿真系统结构图。如图9 所示，一种固体推进剂废料切割处理的数值仿真系统包括：
几何模型建立模块201，用于建立射流喷嘴几何模型和固体推进剂几何模型。
液体-固体颗粒两相流物理模型建立模块202，用于建立前混合磨料水射流的液体-固体颗粒两相流物理模型。
动力学控制方程建立模块202，用于建立固体推进剂的动力学控制方程。
材料参数确定模块203，用于选取磨料、水、固体推进剂三种材料参数。
液体-固体颗粒两相流物理模型离散模块204，用于采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散。
动力学控制方程离散模块205，用于采用SPH方法对所述固体推进剂的动力学控制方程进行离散。
时间积分格式建立模块206，用于建立时间积分格式。
不同时刻的场变量确定模块207，用于根据所述射流喷嘴几何模型、固体推进剂几何模型、离散后的两相流物理模型、离散后的固体推进剂的动力学控制方程采用时间积分格式进行数值仿真计算，得到不同时刻的场变量。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (7)

1.一种固体推进剂废料切割处理的数值仿真方法，其特征在于，包括：
建立射流喷嘴几何模型和固体推进剂几何模型；
建立前混合磨料水射流的液体-固体颗粒两相流物理模型；
建立固体推进剂的动力学控制方程；
选取磨料、水、固体推进剂三种材料参数；
采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散；
采用SPH方法对所述固体推进剂的动力学控制方程进行离散；
建立时间积分格式；
根据所述射流喷嘴几何模型、固体推进剂几何模型、离散后的两相流物理模型、离散后的固体推进剂的动力学控制方程采用时间积分格式进行数值仿真计算，得到不同时刻的场变量。
2.根据权利要求1所述的固体推进剂废料切割处理的数值仿真方法，其特征在于，所述建立前混合磨料水射流的液体-固体颗粒两相流物理模型，具体包括：
建立前混合磨料水射流的液体相控制方程组：
Figure FDA0002674711630000011
Figure FDA0002674711630000012
Figure FDA0002674711630000013
其中，下标l表示的是液体相的标志，与下面的固体颗粒相区分开；αl、ρl和vl分别为液体的体积分数、密度和速度，
Figure FDA0002674711630000014
表示偏导数，
Figure FDA0002674711630000015
是指对时间的偏导数，Pl表示的是液体相的正压力，
Figure FDA0002674711630000016
表示的是梯度，
Figure FDA0002674711630000017
则表示的是散度，g为重力加速度，Rlp为液体相和固体颗粒相之间的相互作用力，液体的粘性作用对于颗粒的作用不可以忽略，τl为液体相的剪切应力张量；
建立前混合磨料水射流的固体颗粒相控制方程组：
Figure FDA0002674711630000021
Figure FDA0002674711630000022
Figure FDA0002674711630000023
Figure FDA0002674711630000024
其中，下标p表示的是固体颗粒相的标志，αp,ρp和vp分别为颗粒相体积分数、密度和速度，
Figure FDA0002674711630000025
为液体连续相压力梯度，
Figure FDA0002674711630000026
为固体颗粒相压力梯度，αpρpg为颗粒所受外部体积力，Rpl为相间相互作用力，hp为固体颗粒相的能量焓值及qp为固体颗粒相内部的热传导量，τp为颗粒相粘性应力张量。
3.根据权利要求1所述的固体推进剂废料切割处理的数值仿真方法，其特征在于，所述建立固体推进剂的动力学控制方程，具体包括：
建立固体推进剂的动力学控制方程：
Figure FDA0002674711630000027
Figure FDA0002674711630000028
Figure FDA0002674711630000029
其中，下标s是固体的表征，ρs为固体推进剂的密度，vs α为固体推进剂的速度在α方向的分量，α取x、y以及z方向，x为位移，σs αβ为固体推进剂的全应力分量，
Figure FDA00026747116300000210
ps为固体推进剂所受到的正应力，
Figure FDA00026747116300000211
为固体推进剂所受到的剪应力，fs α为固体推进剂所受到的其他外部作用力，es为固体推进剂内能。
4.根据权利要求1所述的固体推进剂废料切割处理的数值仿真方法，其特征在于，所述选取磨料、水、固体推进剂三种材料参数，具体包括：
设置水的材料参数为密度ρl＝998Kg/m3，粘度μl＝0.001003Kg/ms；
选取金刚砂作为磨料材料，金刚砂密度ρp＝3900Kg/m3，金刚砂粒径大小为80目，金刚砂与水的混合比1:10；
选取HTPB推进剂作为固体推进剂。
5.根据权利要求2所述的固体推进剂废料切割处理的数值仿真方法，其特征在于，所述采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散，具体包括：
采用SPH方法对所述液体相控制方程组进行离散，得到离散后的液体相控制方程组：
Figure FDA0002674711630000031
Figure FDA0002674711630000032
Figure FDA0002674711630000033
其中，i,j分别是指的i粒子和j粒子，Wij为i粒子和j粒子之间的核函数的数值，W为核函数，h为光滑长度；
采用SPH方法对所述固体颗粒相控制方程组进行离散，得到离散后的固体颗粒相控制方程组：
Figure FDA0002674711630000034
Figure FDA0002674711630000035
Figure FDA0002674711630000041
Figure FDA0002674711630000042
其中，ρp,i为SPH粒子i的密度；ρp为颗粒的实际密度；vp,ij为速度矢量，vp,ij＝vp,i-vp,j，rp,ij为位移矢量，rp,ij＝rp,i-rp,j，Πij为粒子i和粒子j之间的人工粘性。
6.根据权利要求3所述的固体推进剂废料切割处理的数值仿真方法，其特征在于，所述采用SPH方法对所述固体推进剂的动力学控制方程进行离散，具体包括：
采用SPH方法对所述固体推进剂的动力学控制方程进行离散，得到离散后的固体推进剂的动力学控制方程：
Figure FDA0002674711630000043
Figure FDA0002674711630000044
Figure FDA0002674711630000045
其中，
Figure FDA0002674711630000046
为应变率。
7.一种固体推进剂废料切割处理的数值仿真系统，其特征在于，包括：
几何模型建立模块，用于建立射流喷嘴几何模型和固体推进剂几何模型；
液体-固体颗粒两相流物理模型建立模块，用于建立前混合磨料水射流的液体-固体颗粒两相流物理模型；
动力学控制方程建立模块，用于建立固体推进剂的动力学控制方程；
材料参数确定模块，用于选取磨料、水、固体推进剂三种材料参数；
液体-固体颗粒两相流物理模型离散模块，用于采用SPH方法对所述前混合磨料水射流的液体-固体颗粒两相流物理模型进行离散；
动力学控制方程离散模块，用于采用SPH方法对所述固体推进剂的动力学控制方程进行离散；
时间积分格式建立模块，用于建立时间积分格式；
不同时刻的场变量确定模块，用于根据所述射流喷嘴几何模型、固体推进剂几何模型、离散后的两相流物理模型、离散后的固体推进剂的动力学控制方程采用时间积分格式进行数值仿真计算，得到不同时刻的场变量。
CN202010944304.9A 2020-09-10 2020-09-10 一种固体推进剂废料切割处理的数值仿真方法及系统 Active CN112069745B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010944304.9A CN112069745B (zh) 2020-09-10 2020-09-10 一种固体推进剂废料切割处理的数值仿真方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010944304.9A CN112069745B (zh) 2020-09-10 2020-09-10 一种固体推进剂废料切割处理的数值仿真方法及系统

Publications (2)

Publication Number Publication Date
CN112069745A true CN112069745A (zh) 2020-12-11
CN112069745B CN112069745B (zh) 2022-02-08

Family

ID=73663306

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010944304.9A Active CN112069745B (zh) 2020-09-10 2020-09-10 一种固体推进剂废料切割处理的数值仿真方法及系统

Country Status (1)

Country Link
CN (1) CN112069745B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114254572A (zh) * 2021-12-16 2022-03-29 西北工业大学太仓长三角研究院 考虑污染物沉积的航发压气机流场性能预测方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090320443A1 (en) * 2008-05-09 2009-12-31 Geisler Robert L Propulsion system, opposing grains rocket engine, and method for controlling the burn rate of solid propellant grains
CN103294855A (zh) * 2013-05-15 2013-09-11 西安近代化学研究所 固体推进剂羽流特性虚拟试验及羽流数据结构网格化方法
CN104143027A (zh) * 2014-08-01 2014-11-12 北京理工大学 一种基于sph算法的流体热运动仿真系统
CN110283354A (zh) * 2019-05-29 2019-09-27 武汉理工大学 一种针对报废端羟基聚丁二烯复合固体推进剂的降解溶剂和应用
CN111475937A (zh) * 2020-04-03 2020-07-31 中国地质科学院地质力学研究所 一种流固二相流流化滑坡的模拟仿真方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090320443A1 (en) * 2008-05-09 2009-12-31 Geisler Robert L Propulsion system, opposing grains rocket engine, and method for controlling the burn rate of solid propellant grains
CN103294855A (zh) * 2013-05-15 2013-09-11 西安近代化学研究所 固体推进剂羽流特性虚拟试验及羽流数据结构网格化方法
CN104143027A (zh) * 2014-08-01 2014-11-12 北京理工大学 一种基于sph算法的流体热运动仿真系统
CN110283354A (zh) * 2019-05-29 2019-09-27 武汉理工大学 一种针对报废端羟基聚丁二烯复合固体推进剂的降解溶剂和应用
CN111475937A (zh) * 2020-04-03 2020-07-31 中国地质科学院地质力学研究所 一种流固二相流流化滑坡的模拟仿真方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
XIAO, LQ等: "Simulation of solid propellant microstructures by combining the collective rearrangement method with the discrete element method", 《AIP ADVANCES》 *
常武军: "复合固体推进剂细观损伤及其数值仿真研究", 《万方学位论文》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114254572A (zh) * 2021-12-16 2022-03-29 西北工业大学太仓长三角研究院 考虑污染物沉积的航发压气机流场性能预测方法及系统
CN114254572B (zh) * 2021-12-16 2024-01-02 西北工业大学太仓长三角研究院 考虑污染物沉积的航发压气机流场性能预测方法及系统

Also Published As

Publication number Publication date
CN112069745B (zh) 2022-02-08

Similar Documents

Publication Publication Date Title
Sun et al. Numerical simulation of two-phase flows in complex geometries by using the volume-of-fluid/immersed-boundary method
Kleb et al. Sketch-to-solution: An exploration of viscous CFD with automatic grids
Premnath et al. Lattice Boltzmann model for axisymmetric multiphase flows
Ma et al. Euler–Lagrange simulations of bubble cloud dynamics near a wall
Slone et al. A finite volume unstructured mesh approach to dynamic fluid–structure interaction: an assessment of the challenge of predicting the onset of flutter
Zheleznyakova et al. Molecular dynamics-based unstructured grid generation method for aerodynamic applications
Janssen et al. Modeling of wave breaking andwave-structure interactions by coupling of fully nonlinear potential flow and lattice-Boltzmann models
Jiang et al. The SPH method for simulating a viscoelastic drop impact and spreading on an inclined plate
Chen et al. GPU-accelerated smoothed particle hydrodynamics modeling of granular flow
Pederson et al. The Sedov blast wave as a radial piston verification test
US20230195978A1 (en) Performance prediction method and system for whole atomization process of aeroengine fuel
Candler et al. Current status and future prospects for the numerical simulation of hypersonic flows
CN112069745B (zh) 一种固体推进剂废料切割处理的数值仿真方法及系统
Fei et al. A multi-scale model for coupling strands with shear-dependent liquid
Zhang et al. Continuous contact force model with an arbitrary damping term exponent: model and discussion
Bašić et al. Simulation of water entry and exit of a circular cylinder using the ISPH method
Larsen Impact loads on circular cylinders
Horton et al. Benchmarking of computational fluid methodologies in resolving shear-driven flow fields
JP2005316614A (ja) 最適化方法及び最適化プログラム
Lygidakis et al. Numerical analysis of flow over the NASA Common Research Model using the academic Computational Fluid Dynamics code Galatea
Xu et al. SPH and ALE formulations for sloshing tank analysis
Hall et al. A computational model for the flow of resin in self-healing composites
Mierke et al. GPU-accelerated large-eddy simulation of ship-ice interactions
Shi et al. SPH method with space-based variable smoothing length and its applications to free surface flow
Ge et al. Particle methods for multiscale simulation of complex flows

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant