考虑细砂冲击压缩特性和吸热相变的数学模型及建模方法

摘要

本发明适用于爆炸与冲击数值模拟技术领域，公开了一种考虑细砂冲击压缩特性和吸热相变的数学模型及建模方法，该方法包括以下步骤：基于连续介质假设，建立双流体模型控制方程；根据爆炸类型建立气相状态方程封闭模型；采用凝聚介质状态方程描述细砂的冲击压缩特性，根据材料冲击状态方程实验数据确定状态方程参数，建立固相状态方程封闭模型；根据多相流流动状态，建立相间阻力作用模型；根据流场对流换热特性和颗粒介质熔化、汽化特性，建立相间热交换率和相间传质模型。本发明提出了描述冲击波与充砂结构相互作用物理过程的封闭数学模型，为爆炸加载下充砂结构消波吸能效应数值模拟技术研究提供了理论模型。

说明书

技术领域

本发明属于爆炸与冲击数值模拟技术领域，涉及一种考虑细砂冲击压缩特性和吸热相变的双流体数学模型及其建模方法。

背景技术

充砂结构作为一种经济型爆炸防护结构，用于降低爆炸冲击载荷对于人员和各类目标的损伤。通过对充砂结构设计形式进行合理优化，可以最大程度地提升结构的防护性能、降低结构主体颗粒群的高速运动冲击效应。在爆炸防护技术研究领域，数值研究是解决充砂结构科学优化设计问题的主要技术手段之一。

充砂结构由细砂组成，在冲击波加载下，细砂将呈现拟流体的冲击压缩特性，二相流场发展将经历堆积颗粒流动、稠密气固二相流动、稀疏气固二相流动等多种流动型态变化，而且高速高温流体与细砂之间将产生强烈的质量、动量和能量的耦合作用。在爆炸与冲击数值模拟技术领域，适用于解决冲击波与充砂结构相互作用问题的数值模拟方法为多介质欧拉方法。目前，在多介质欧拉坐标体系下所发展的、适于求解该类问题的多介质数学模型，大致可以分为两类：一类假设流体-细砂主体间具有清晰物质界面，数学模型考虑物质界面处动量、能量的耦合作用，但无法考虑界面处因传热引起的温升相变效应，如：主流商业软件LS-DYNA、AUTODYN所集成的多介质欧拉求解器；另一类假设介质间相互渗透、体积共占，数学模型全面考虑相间的阻力作用、传热及传质过程，但并不考虑介质的冲击压缩特性，如：气固多相流领域所发展的双流体模型。为了解决现有多介质欧拉方法数学模型无法全面考虑细砂主体的冲击压缩特性和吸热相变效应的问题，需要建立一种适用于在欧拉坐标体系下进行数值求解的更为完善的数学模型。

发明内容

本发明的目的是提出一种考虑细砂冲击压缩特性和吸热相变的双流体数学模型及其建模方法，能够合理而全面地描述冲击波加载下二相流场内介质间的相互作用物理过程，弥补了现有冲击波与充砂结构相互作用多介质欧拉数值计算模型的不足，对于冲击波与充砂结构相互作用多介质欧拉数值计算模拟技术的发展具有重要的促进作用，对于开展充砂结构爆炸防护性能数值研究具有重要的工程应用价值。

本发明的技术方案是提供一种考虑细砂冲击压缩特性和吸热相变的数学模型，其特殊之处在于：包括双流体模型控制方程、气相状态方程封闭模型、固相状态方程封闭模型、相间阻力作用模型、相间热交换率模型及相间传质模型。

进一步地，双流体模型控制方程如下：

气相质量守恒方程：

式中，ε

气相动量守恒方程：

式中，

气相能量守恒方程：

式中，E

固相质量守恒方程：

式中，ε

固相x、y、z方向动量守恒方程：

式中，p

固相能量守恒方程：

式中，E

固相粒子数守恒方程：

气相和固相体积分数之间所满足的关系：

ε

固相粒子数密度和固相体积分数满足如下关系：

式中，d

进一步地，所述固相状态方程封闭模型如下式：

固相单位质量总能量E

固相分压p

p

式中，p

采用凝聚介质状态方程描述细砂的冲击压缩响应特性，p

式中：

进一步地，相间阻力作用模型如下：

采用牛顿阻力公式描述相间阻力作用，相间阻力F

式中，C

C

颗粒雷诺数Re

式中，μ为流体动力粘性系数。

进一步地，相间热交换率模型如下：

相间热量交换率Q

式中，努赛尔数(Nusselt)Nu用于表征对流换热与导热之比，λ是流体的导热系数，T

进一步地，相间传质模型以相间质量交换率Γ

相间质量交换率Γ

式中：c

进一步地，固相温度T

式中：T

进一步地，气相状态方程封闭模型，具体为：

气相单位质量总能量E

气相热力学温度与单位质量比内能满足如下关系：

T

气相分压p

γ

式中：γ

本发明还提供一种上述考虑细砂冲击压缩特性和吸热相变的数学模型的建模方法，其特殊之处在于，基于连续介质假设，按照下述步骤建立以下模型：

步骤1、建立双流体模型控制方程；

步骤2、根据爆炸类型建立气相状态方程封闭模型；

采用凝聚介质状态方程描述细砂的冲击压缩特性，根据材料冲击状态方程实验数据确定状态方程参数，建立固相状态方程封闭模型；

根据多相流流动状态，建立相间阻力作用模型；

根据流场对流换热特性和颗粒介质熔化、汽化特性，建立相间热交换率模型；

根据流场对流换热特性和颗粒介质熔化、汽化特性，建立相间传质模型。

进一步地，假设组成充砂结构的主体细砂为一种拟流体，采用双流体模型描述冲击波与充砂结构相互作用物理过程，在欧拉-欧拉坐标系下描述气相和固相。文献专著中所推导双流体模型理论控制方程极为复杂，数学模型封闭难、求解难，不适于工程应用，因此，本发明在建模中仅关注重点影响因素和机制，具体为：

(1)假设充砂结构主体细砂由单一球形惰性颗粒介质组成；

(2)假设气固二相流动体系绝热，与外界不存在表面传热、体加热和外力做功；

(3)不考虑气固二相流动湍流效应；

(4)考虑气体和颗粒群之间的热传导、辐射传热过程；

(5)考虑吸热引起惰性颗粒的熔化、汽化过程；

(6)假设相间传质过程中颗粒粒子数守恒；

(7)考虑气体和颗粒群之间的相间阻力作用；

(8)考虑细砂的冲击压缩过程；

(9)不考虑细砂的剪切效应。

进一步地，步骤1建立的双流体模型控制方程如下：

气相质量守恒方程：

气相动量守恒方程：

气相能量守恒方程：

固相质量守恒方程：

固相x、y、z方向动量守恒方程：

固相能量守恒方程：

固相粒子数守恒方程：

气相和固相体积分数之间所满足的关系：

ε

固相粒子数密度和固相体积分数满足如下关系：

进一步地，步骤2中采用凝聚介质状态方程描述细砂的冲击压缩特性，根据细砂冲击状态方程实验，建立固相状态方程封闭模型，具体如下：

固相单位质量总能量E

固相分压p

p

式中：p

采用凝聚介质状态方程描述充砂结构的冲击压缩响应特性，p

上式右端第一项主要反映“冷状态”的性质，由介质体积变形引起；第二项主要反映“热状态”的性质。

根据材料冲击状态方程实验数据确定状态方程参数c

进一步地，步骤2中根据多相流流动状态，建立相间阻力作用模型，当流体-颗粒之间保持相对运动时，将受到来自流体的广义相间阻力和侧向力，前者包括：阻力、附加质量力、Basset力，后者包括：升力、Magnus力和Staffman力等。建模中仅考虑发挥主要作用的流体阻力，并采用牛顿阻力公式描述相间阻力作用，给出步骤4中相间阻力F

采用牛顿阻力公式描述相间阻力作用，相间阻力F

颗粒群阻力系数C

颗粒雷诺数Re

进一步地，本发明建模中综合考虑了相间热传导和辐射传热传热机制。当高温气流流经颗粒表面，如果相间存在温度差，那么在温度梯度驱动下相间将产生热传导过程，单位时间内流体通过颗粒表面传递给颗粒的热量正比于温度差和颗粒表面积，建模中采用牛顿热传导方程描述相间传热率；除此之外，流体-颗粒表面还可以以热辐射方式进行相间传热，建模中采用Stefan-Boltzmann定律描述颗粒向周围流体的辐射热流强度。最终给出步骤5中相间热量交换率Q

相间热量交换率Q

进一步地，高温高速气流冲刷作用下，细砂吸热产生温升，当颗粒温度达到熔化温度直至汽化温度时，颗粒自表面产生熔化、汽化相变过程，将与气相产生传质过程。在强迫对流环境中，采用考虑强迫对流效应的液滴蒸发公式描述相间传质过程，给出步骤6中相间质量交换率Γ

相间质量交换率Γ

进一步地，步骤5和步骤6中，在高温气流作用下，考虑固相颗粒介质的熔化、汽化相变过程时，固相温度T

进一步地，步骤2根据爆炸类型，建立气相状态方程封闭模型。当爆炸源区质量能量密度不同时，需要采用特定的产物状态方程描述源区膨胀过程中热力学参数的变化，具体为：

气相单位质量总能量E

气相热力学温度与单位质量比内能满足如下关系：

T

气相分压p

γ

本发明的有益效果是：

本发明提出一种考虑细砂冲击压缩特性和吸热相变的双流体数学模型，在统一欧拉-欧拉型计算模型框架下，能够综合考虑冲击波加载下细砂的冲击压缩过程、辐射与对流机制引起的吸热相变过程，弥补了现有冲击波与充砂结构作用多介质欧拉数值计算模型无法合理描述充砂结构冲击压缩过程的不足，对于冲击波与充砂结构相互作用多介质欧拉数值模拟技术的发展具有重要的促进作用，对于开展充砂结构爆炸防护性能数值研究具有重要的工程应用价值。

附图说明

图1是实施例的方法流程图。

图2是发明实例图，(a)为管道内爆炸源示意图，(b)为管道内充砂结构示意图。

图3是t＝0.2ms时刻，气固两相流场参数云图分布，(a)为固相体积分数云图分布，(b)为气相、固相压力云图分布，(c)为气相、固相速度模量云图分布，(e)为气相、固相温度云图分布。

图4是t＝5.5ms、6.0ms、6.5ms、7.0ms、7.5ms时刻，固相体积分数、固相压力分布云图，其中(a)是固相体积分数分布云图，(b)是固相压力p

图5是流场能量、质量守恒性检验，其中(a)是能量守恒检验，(b)是质量守恒检验，(c)是两相动能占总能量比例变化曲线，(d)是两相内能占总能量比例变化曲线。

图6是压力时间历程数值模拟结果与实验测试结果对比，其中(a)是管道壁面0.621m处压力时间历程数值模拟结果与实验测试结果对比，(b)是管道壁面0.691m处压力时间历程数值模拟结果与实验测试结果对比。

具体实施方式

为了使本发明的目的、数学模型的建立方法及优点更加清楚明白，下面结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。

针对现有多介质欧拉模拟方法中数学模型未同时考虑介质冲击响应特性和传热机制的问题，本发明提供了一种考虑细砂冲击压缩特性和吸热相变的双流体数学模型，下面结合附图对本发明作详细的描述。

如图1所示，本实施例以二维轴对称数学模型为例，进行说明，具体建模过程如下：

S101：基于连续介质假设，建立二维轴对称双流体模型控制方程，形式如下；

气相质量守恒方程：

气相径向R、轴向Z的动量守恒方程：

气相能量守恒方程：

固相质量守恒方程：

固相径向R、轴向Z的动量守恒方程：

固相能量守恒方程：

固相粒子数守恒方程：

气相和固相体积分数之间所满足的关系：

ε

上述各式中：下表g、p分别代表气相和固相；ε为相k(＝g，p)的体积分数，代表单位体积中相k所占的体积分数；ρ为相密度，代表相k的质量与相k所占体积之比；u、v为相k的径向R、轴向Z的速度分量；n

固相粒子数密度和固相体积分数满足如下关系：

S102：根据爆炸类型，建立气相状态方程封闭模型。

图2为发明实例示意图，图中：爆炸源区为一种化学装药，在管道中心爆炸，管道两侧各设置13道砂墙。给定以下形式状态方程：

气相单位质量总能量E

气相热力学温度与单位质量比内能满足如下关系：

气相分压p

式(39)～(41)中：γ

S103：采用凝聚介质状态方程描述细砂的冲击压缩特性，根据细砂材料状态方程实验数据，确定固相状态方程封闭模型。

固相单位质量总能量与单位质量比内能满足如下关系：

固相分压与气相分压之间满足如下关系：

p

式中：p

采用凝聚介质状态方程描述细砂的冲击压缩响应特性，p

式(19)中：

图2所示充砂结构由石英细砂组成，参阅文献Laine等(2001)、Brown等(2007)给出的石英干砂材料状态方程实验测试数据，确定状态方程参数c

c

石英细砂其他材料参数取值为：

S104：根据多相流流动状态，建立相间阻力作用模型，具体为：

采用牛顿阻力公式描述相间阻力作用，相间阻力F

其中：

C

颗粒雷诺数计算公式如下：

式中：d

S105：所述根据流场对流换热特性和颗粒介质熔化、汽化特性，建立相间热交换率模型。

相间热量交换率Q

式中：右端第一项表示通过热传导机制引起的相间热量交换，右端第二项表示通过辐射传热机制引起的热量交换；努赛尔特数(Nusselt)Nu用于表征对流换热与导热之比；λ是流体的导热系数；θ

采用著名的Ranz-Marshall公式计算Nu，形式如下：

Nu＝2.0+0.6Re

式中：普朗特(Prandtl)数Pr表征定压比热与定容比热之比。

S106：所述根据流场对流换热特性和颗粒介质熔化、汽化特性，建立相间质量交换率模型。

相间质量交换率Γ

式中：q

S105和S106：在高温气流作用下，考虑固相颗粒介质的熔化、汽化相变过程时，固相温度T

式中：T

石英砂其他热力学参数见表1：

表1石英砂热力学参数

图3所示为t＝0.2ms时刻气固两相流场参数云图分布，如图所示：利用所研制二维轴对称气固二相数值模拟程序，可以获取任意时刻流场中砂墙体积分数、气固二相压力、二相速度、二相温度分布等参数的细致数值模拟结果。爆炸冲击波驱动下，高速气流曳力作用下，离散颗粒介质最终将在容器末端产生堆积行为。图4中(a)所示为t＝5.5ms～7.5ms时刻，容器末端固相体积分数分布云图，如图所示：固相体积分数经历了先增后减的变化趋势：反映了气相介质作用下离散砂粒运动至容器末端被压缩、被稀疏的运动过程。图4中(b)所示为t＝5.5ms～7.5ms时刻，离散砂粒压缩稀疏过程中固相压力分布云图，如图所示：随固相堆积状态变化，固相介质压缩正应力也将随之产生或增强或衰减的变化。图4说明，当前固相压力模型可以描述离散颗粒在压缩过程中，内部压缩正应力的变化过程。

图5中(a)所示为流场总能量、固相总能量、气相总能量随时间变化曲线。图5中(b)所示为流场总质量、固相总质量、气相总质量随时间变化曲线。如图所示：流场总能量、总质量保持良好守恒性。图5中(c)为气固两相总动能占总能量比例随时间变化曲线，图5中(d)为两相内能占总能量比例随时间变化曲线，根据两相内能、动能变化曲线可以对充砂结构的吸能效应做出预估。

图6中(a)～(b)所示为管道壁面0.621m和0.691m处冲击波压力时间历程曲线数值模拟结果与实验测试结果对比曲线。如图所示：数值模拟结果所得压力时间历程曲线整体分布趋势与实验测试结果基本一致，压力达到准静态数值模拟结果与实测结果基本一致，而且数值模拟结果可以捕捉到波系传播演化过程中的详细细节。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。