点源扩散燃烧时的形状函数研究

孟 杰，辛长范，马连泽，杨宇青，魏志芳，张树霞

(中北大学， 太原 030051)

1 引言

为了掌握枪炮膛内压力变化规律和弹丸速度变化规律[1]，了解燃烧时的形状函数变得尤为重要，目前在燃烧形状函数推导方面的相关研究很少，本研究提出了一种小尺寸、片状圆柱体的规则形状燃烧物，为了解其燃烧性能与形状函数，得到形状特征量，建立了燃烧模型，通过数值模拟计算，分析了燃烧过程，旨在为片状圆柱体的燃烧性能研究提供理论支持[2]。

2 形状函数

对于一定形状的易燃物来说，其燃烧具有一定的规律性。因此，已燃相对质量ψ、相对燃烧面σ和已燃相对厚度Z之间存在一定的函数关系。

将已燃相对厚度Z取为自变量，则：

ψ=f1(Z)

(1)

σ=f2(Z)

(2)

称为形状函数。

在以上2个形状函数关系式中，根据气体生成速率规律，只要知道其中的一个关系式，经过积分或者微分，便能得到另一个关系式[3]。

3 燃烧模型的建立

3.1 形状函数计算的基本假设

片状圆柱体结构如图1所示。

图1 片状圆柱体结构示意图Fig.1 Flake cylindrical structure diagram

为方便计算，做出如下假设[4]：

1) 圆柱体燃烧遵循几何燃烧定律；

2) 在t=0时刻，燃烧点在圆柱中心线上，以点源扩散的方式从燃烧点开始燃烧；

3) 燃烧物的结构为圆柱体，直径为D，高度为H[5]。

3.2 燃烧过程

某一圆柱体的直径为5 mm、高度为1 mm，由于圆柱体为轴对称图形，并且从燃烧点开始所燃烧过的体积和燃烧面的面积也为轴对称图形，为了便于计算，取圆柱体对称轴右侧部分进行分析，把片状圆柱体燃烧分为以下3个阶段[6]。

第1阶段：燃烧面未到达圆柱体下底面阶段(0

图2 第1阶段燃烧平面示意图Fig.2 First stage combustion diagram

第2阶段：燃烧面到达圆柱体下底面但未到达圆柱体侧面阶段(H/Rmax

图3 第2阶段燃烧示意图Fig.3 Second stage combustion diagram

第3阶段：燃烧面到达圆柱体侧面但圆柱体未燃尽(D/2Rmax

图4 第三阶段燃烧平面示意图Fig.4 Third stage combustion diagram

4 形状函数推导

从图2、图3、图4可知，燃烧部分为一定形状的旋转体，所以从旋转体入手推导形状函数。

4.1 实际公式推导

实际公式的推导采用微元法。

4.1.1 旋转体体积微元

假设某一小矩形块长为dx、宽为y，距离旋转轴Y的距离为x，如图5(a)所示。让该小矩形块绕Y轴旋转一周，得到的是一个柱壳，如图5(b)所示。接下来要做的是求一个柱壳的体积。柱壳的半径为x、高为y、厚度为dx。将它展开平铺成一个矩形状金属片，因为壁厚很小，可以将它考虑成一个长方体。理想柱壳的展开图，其厚度为dx、高度为柱壳的高度y、其长度等于柱壳的底圆周长2πx，所以柱壳的体积与2πxydx很接近[7]。

图5 旋转体体积微元示意图Fig.5 The volume element diagram of a rotating body

4.1.2 旋转体表面积微元

假设s是距离Y轴为x的圆弧段上的小段圆弧，如图6(a)中蓝色曲线部分，将该圆弧绕Y轴旋转一周，得到一个曲边环，如图6(b)所示。

为了计算方便，用直边代替曲边，如图6(c)蓝色线段所示。设线段在Y轴上的投影为L，蓝色线段与其投影线的夹角为θ，则线段的长度为L/cosθ。用该直线代替弧段旋转时，得到的是一个直边环。此时，虽然得到的是直边环，但环的边并不平行于Y轴，所以得到的环实际上是圆锥表面的一部分。这个表面积是可以计算的，但是比较麻烦。因此为了简化，可以进一步近似，近似成为一个边长都相等的环，即这个环是圆柱形的[8]。最终的结果，得到了一个半径为x、宽度为L/cosθ的圆柱形环，如图6(d)所示，其表面积为2πxL/cosθ。

图6 旋转体表面积微元示意图Fig.6 The surface area element diagram of a rotating body

实际公式下的形状函数ψ-Z与σ-Z关系由上述2个部分推导得到。

4.2 理论公式推导

理论公式的推导采用球的常用表面积、体积公式，用来与实际公式形成对比。

如果孩子通过语言等表达不满，仍然无法阻止对方的欺负行为，那就要马上离开，跑去一个安全的区域，比如小朋友很多的地方，或者家长聊天的区域，或者老师身边。即便一时找不到自己熟悉的小朋友或者成人也没关系。越是恶劣的、有主观意图的欺负甚至校园霸凌，都越是隐藏在那些阴暗的死角不能曝光，因此跑到其他人能注意到的地带，就相对安全很多。

4.2.1 旋转体体积公式

当燃烧到达第3阶段时，燃烧二维状态图如图7所示。其中A点为燃烧面与圆柱下表面交线上的一点，B点为燃烧面与圆柱右侧面交线上的一点，角θ1为燃烧点和B点的连线与上表面的夹角，角θ为燃烧点和A点的连线与上表面的夹角。

图7 第3阶段时的燃烧二维状态图Fig.7 The combustion state diagram of the third stage

假设此时燃烧面的半径为R，则θ1角处的圆面半径为Rcosθ1，当角度增加dθ时，其所增加的高度h为Rdθcosθ1。因此体积微元为[9]：

dV=π(Rcosθ1)2Rcosθ1dθ

(3)

角θ1与角θ所夹得这部分体积即为：

(4)

积分得：

(5)

根据式(5)，便可以得到燃烧时每个阶段的体积公式[10]，即：

(6)

4.2.2 旋转体表面积公式

dS=2πRcosθ1Rdθ=2πR2cosθ1dθ

(7)

角θ1与角θ所夹得这部分表面积即为：

(8)

燃烧时每个阶段的表面积公式为：

(9)

理论公式下的形状函数ψ-Z与σ-Z关系由式(6)和式(9)所得。

5 数值模拟

片状圆柱体的直径为5 mm、高度为1 mm，定义沿着径向为X轴，沿着高度方向为Y轴，在二维平面上划分网格(网格单位长度可以根据燃烧速率划分)，沿着X轴方向网格以0.01 mm为一个单位长度，沿着Y轴方向网格以0.01 mm为一个单位长度[12]。

假设物体燃烧速率为一个单位时间燃烧一个网格单位长度，t时刻时燃烧面为图8中的曲线部分，将已燃烧完的每个小网格标记为0，未燃烧或未燃烧完的小网格标记为1。

那么，实际公式下旋转体的体积就转换为由图8中标记为0的每一个相同于图5中小矩形块绕Y轴旋转一周的体积相加所得，取dx为网格的长、y为网格的宽，则：

(10)

式(10)中：i为标记为0的网格数；n为标记为0的网格总数。将图8中圆弧按每一行划分成多个圆弧段，实际公式下的旋转体表面积就转换为由相同于图6中每一个圆弧段绕Y轴旋转一周的表面积相加所得，取L为网格的宽，则：

(11)

式(11)中：j为燃烧面所经过的网格行数；m为燃烧面所燃烧过的总行数。

图8所示的燃烧状态标记图则需要在Matlab中进行一系列编程得到。

图8 t时刻燃烧状态标记图Fig.8 Burning status diagram at time t

5.1 形状特征量

在Matlab中编写程序，以燃烧面半径每增加一单位时的厚度除以最大燃烧面半径，即已燃相对厚度Z作为横轴；分别以所燃烧过的体积除以片状圆柱的总体积，即已燃相对质量ψ和以燃烧处表面积除以片状圆柱的总表面积，即相对燃烧面σ作为纵轴画曲线，并拟合得出其对应函数关系式[13]。

图9反映了形状函数ψ-Z关系与拟合关系，其中黑色线为ψ-Z关系，红色线为拟合关系，图10反映了σ-Z关系。图10中绿点为第1阶段过渡到第2阶段的特殊标注点，蓝点为第2阶段过渡到第3阶段的特殊标注点。

图9 形状函数ψ-Z关系与拟合关系曲线Fig.9 The shape function ψ-Z relation and fitting relation diagram

图10 形状函数σ-Z关系曲线Fig.10 The shape function σ-Z relation diagram

图9拟合所得出的已燃相对质量ψ与已燃相对厚度Z之间的函数关系为：

ψ=0.220 1Z3+0.981 3Z2-0.087 7Z

(10)

同时对应

ψ=χZ(1+λZ+μZ2)

(11)

便可求得3个形状特征量χ、λ、μ，即：

(12)

由于2个形状函数式存在积分或者微分关系，只需求出一个即可。

5.2 实际与理论的对比分析

将实际形状函数推导与理论形状函数推导所得关系进行画图比较[14]。其中，图11中黑色线为实际推导形状函数ψ-Z关系，红色线为理论推导形状函数ψ-Z关系；图12中黑色线为实际推导形状函数σ-Z关系，红色线为理论推导形状函数σ-Z关系[15]。

图11 形状函数ψ-Z实际与理论曲线Fig.11 Comparison diagram of actual and theoretical shape function ψ-Z

图12 形状函数σ-Z实际与理论曲线Fig.12 Comparison diagram of actual and theoretical shape function σ-Z

由图11、图12可看出，实际推导所得曲线与理论推导所得曲线几乎重合，由此本文的实际公式是正确的。

6 结论

1) 实际燃烧形状函数与理论分析吻合较好，为了得到更准确的形状函数及其形状特征量，可以在划分网格时将网格划分的更小更细；

2) 可以考虑改变程序，通过改变程序中步长，即相当于改变燃烧速率，得到相应的形状函数，所以此研究可以对计算多种燃烧速率下的形状函数提供参考；

3) 由于实际推导下的体积与表面积求解方法为微元法，因此可以对求解多种规则形状物体以及不规则形状物体的形状函数的数值模拟与计算提供参考。