孟 杰，辛長范，馬連澤，楊宇青，魏志芳，張樹霞

(中北大學(xué)， 太原 030051)

1 引言

為了掌握槍炮膛內(nèi)壓力變化規(guī)律和彈丸速度變化規(guī)律[1]，了解燃燒時的形狀函數(shù)變得尤為重要，目前在燃燒形狀函數(shù)推導(dǎo)方面的相關(guān)研究很少，本研究提出了一種小尺寸、片狀圓柱體的規(guī)則形狀燃燒物，為了解其燃燒性能與形狀函數(shù)，得到形狀特征量，建立了燃燒模型，通過數(shù)值模擬計算，分析了燃燒過程，旨在為片狀圓柱體的燃燒性能研究提供理論支持[2]。

2 形狀函數(shù)

對于一定形狀的易燃物來說，其燃燒具有一定的規(guī)律性。因此，已燃相對質(zhì)量ψ、相對燃燒面σ和已燃相對厚度Z之間存在一定的函數(shù)關(guān)系。

將已燃相對厚度Z取為自變量，則：

ψ=f1(Z)

(1)

σ=f2(Z)

(2)

稱為形狀函數(shù)。

在以上2個形狀函數(shù)關(guān)系式中，根據(jù)氣體生成速率規(guī)律，只要知道其中的一個關(guān)系式，經(jīng)過積分或者微分，便能得到另一個關(guān)系式[3]。

3 燃燒模型的建立

3.1 形狀函數(shù)計算的基本假設(shè)

片狀圓柱體結(jié)構(gòu)如圖1所示。

圖1 片狀圓柱體結(jié)構(gòu)示意圖Fig.1 Flake cylindrical structure diagram

為方便計算，做出如下假設(shè)[4]：

1) 圓柱體燃燒遵循幾何燃燒定律；

2) 在t=0時刻，燃燒點在圓柱中心線上，以點源擴(kuò)散的方式從燃燒點開始燃燒；

3) 燃燒物的結(jié)構(gòu)為圓柱體，直徑為D，高度為H[5]。

3.2 燃燒過程

某一圓柱體的直徑為5 mm、高度為1 mm，由于圓柱體為軸對稱圖形，并且從燃燒點開始所燃燒過的體積和燃燒面的面積也為軸對稱圖形，為了便于計算，取圓柱體對稱軸右側(cè)部分進(jìn)行分析，把片狀圓柱體燃燒分為以下3個階段[6]。

第1階段：燃燒面未到達(dá)圓柱體下底面階段(0

圖2 第1階段燃燒平面示意圖Fig.2 First stage combustion diagram

第2階段：燃燒面到達(dá)圓柱體下底面但未到達(dá)圓柱體側(cè)面階段(H/Rmax

圖3 第2階段燃燒示意圖Fig.3 Second stage combustion diagram

第3階段：燃燒面到達(dá)圓柱體側(cè)面但圓柱體未燃盡(D/2Rmax

圖4 第三階段燃燒平面示意圖Fig.4 Third stage combustion diagram

4 形狀函數(shù)推導(dǎo)

從圖2、圖3、圖4可知，燃燒部分為一定形狀的旋轉(zhuǎn)體，所以從旋轉(zhuǎn)體入手推導(dǎo)形狀函數(shù)。

4.1 實際公式推導(dǎo)

實際公式的推導(dǎo)采用微元法。

4.1.1 旋轉(zhuǎn)體體積微元

假設(shè)某一小矩形塊長為dx、寬為y，距離旋轉(zhuǎn)軸Y的距離為x，如圖5(a)所示。讓該小矩形塊繞Y軸旋轉(zhuǎn)一周，得到的是一個柱殼，如圖5(b)所示。接下來要做的是求一個柱殼的體積。柱殼的半徑為x、高為y、厚度為dx。將它展開平鋪成一個矩形狀金屬片，因為壁厚很小，可以將它考慮成一個長方體。理想柱殼的展開圖，其厚度為dx、高度為柱殼的高度y、其長度等于柱殼的底圓周長2πx，所以柱殼的體積與2πxydx很接近[7]。

圖5 旋轉(zhuǎn)體體積微元示意圖Fig.5 The volume element diagram of a rotating body

4.1.2 旋轉(zhuǎn)體表面積微元

假設(shè)s是距離Y軸為x的圓弧段上的小段圓弧，如圖6(a)中藍(lán)色曲線部分，將該圓弧繞Y軸旋轉(zhuǎn)一周，得到一個曲邊環(huán)，如圖6(b)所示。

為了計算方便，用直邊代替曲邊，如圖6(c)藍(lán)色線段所示。設(shè)線段在Y軸上的投影為L，藍(lán)色線段與其投影線的夾角為θ，則線段的長度為L/cosθ。用該直線代替弧段旋轉(zhuǎn)時，得到的是一個直邊環(huán)。此時，雖然得到的是直邊環(huán)，但環(huán)的邊并不平行于Y軸，所以得到的環(huán)實際上是圓錐表面的一部分。這個表面積是可以計算的，但是比較麻煩。因此為了簡化，可以進(jìn)一步近似，近似成為一個邊長都相等的環(huán)，即這個環(huán)是圓柱形的[8]。最終的結(jié)果，得到了一個半徑為x、寬度為L/cosθ的圓柱形環(huán)，如圖6(d)所示，其表面積為2πxL/cosθ。

圖6 旋轉(zhuǎn)體表面積微元示意圖Fig.6 The surface area element diagram of a rotating body

實際公式下的形狀函數(shù)ψ-Z與σ-Z關(guān)系由上述2個部分推導(dǎo)得到。

4.2 理論公式推導(dǎo)

理論公式的推導(dǎo)采用球的常用表面積、體積公式，用來與實際公式形成對比。

如果孩子通過語言等表達(dá)不滿，仍然無法阻止對方的欺負(fù)行為，那就要馬上離開，跑去一個安全的區(qū)域，比如小朋友很多的地方，或者家長聊天的區(qū)域，或者老師身邊。即便一時找不到自己熟悉的小朋友或者成人也沒關(guān)系。越是惡劣的、有主觀意圖的欺負(fù)甚至校園霸凌，都越是隱藏在那些陰暗的死角不能曝光，因此跑到其他人能注意到的地帶，就相對安全很多。

4.2.1 旋轉(zhuǎn)體體積公式

當(dāng)燃燒到達(dá)第3階段時，燃燒二維狀態(tài)圖如圖7所示。其中A點為燃燒面與圓柱下表面交線上的一點，B點為燃燒面與圓柱右側(cè)面交線上的一點，角θ1為燃燒點和B點的連線與上表面的夾角，角θ為燃燒點和A點的連線與上表面的夾角。

圖7 第3階段時的燃燒二維狀態(tài)圖Fig.7 The combustion state diagram of the third stage

假設(shè)此時燃燒面的半徑為R，則θ1角處的圓面半徑為Rcosθ1，當(dāng)角度增加dθ時，其所增加的高度h為Rdθcosθ1。因此體積微元為[9]：

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

(3)

角θ1與角θ所夾得這部分體積即為：

(4)

積分得：

(5)

根據(jù)式(5)，便可以得到燃燒時每個階段的體積公式[10]，即：

(6)

4.2.2 旋轉(zhuǎn)體表面積公式

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

(7)

角θ1與角θ所夾得這部分表面積即為：

(8)

燃燒時每個階段的表面積公式為：

(9)

理論公式下的形狀函數(shù)ψ-Z與σ-Z關(guān)系由式(6)和式(9)所得。

5 數(shù)值模擬

片狀圓柱體的直徑為5 mm、高度為1 mm，定義沿著徑向為X軸，沿著高度方向為Y軸，在二維平面上劃分網(wǎng)格(網(wǎng)格單位長度可以根據(jù)燃燒速率劃分)，沿著X軸方向網(wǎng)格以0.01 mm為一個單位長度，沿著Y軸方向網(wǎng)格以0.01 mm為一個單位長度[12]。

假設(shè)物體燃燒速率為一個單位時間燃燒一個網(wǎng)格單位長度，t時刻時燃燒面為圖8中的曲線部分，將已燃燒完的每個小網(wǎng)格標(biāo)記為0，未燃燒或未燃燒完的小網(wǎng)格標(biāo)記為1。

那么，實際公式下旋轉(zhuǎn)體的體積就轉(zhuǎn)換為由圖8中標(biāo)記為0的每一個相同于圖5中小矩形塊繞Y軸旋轉(zhuǎn)一周的體積相加所得，取dx為網(wǎng)格的長、y為網(wǎng)格的寬，則：

(10)

式(10)中：i為標(biāo)記為0的網(wǎng)格數(shù)；n為標(biāo)記為0的網(wǎng)格總數(shù)。將圖8中圓弧按每一行劃分成多個圓弧段，實際公式下的旋轉(zhuǎn)體表面積就轉(zhuǎn)換為由相同于圖6中每一個圓弧段繞Y軸旋轉(zhuǎn)一周的表面積相加所得，取L為網(wǎng)格的寬，則：

(11)

式(11)中：j為燃燒面所經(jīng)過的網(wǎng)格行數(shù)；m為燃燒面所燃燒過的總行數(shù)。

圖8所示的燃燒狀態(tài)標(biāo)記圖則需要在Matlab中進(jìn)行一系列編程得到。

圖8 t時刻燃燒狀態(tài)標(biāo)記圖Fig.8 Burning status diagram at time t

5.1 形狀特征量

在Matlab中編寫程序，以燃燒面半徑每增加一單位時的厚度除以最大燃燒面半徑，即已燃相對厚度Z作為橫軸；分別以所燃燒過的體積除以片狀圓柱的總體積，即已燃相對質(zhì)量ψ和以燃燒處表面積除以片狀圓柱的總表面積，即相對燃燒面σ作為縱軸畫曲線，并擬合得出其對應(yīng)函數(shù)關(guān)系式[13]。

圖9反映了形狀函數(shù)ψ-Z關(guān)系與擬合關(guān)系，其中黑色線為ψ-Z關(guān)系，紅色線為擬合關(guān)系，圖10反映了σ-Z關(guān)系。圖10中綠點為第1階段過渡到第2階段的特殊標(biāo)注點，藍(lán)點為第2階段過渡到第3階段的特殊標(biāo)注點。

圖9 形狀函數(shù)ψ-Z關(guān)系與擬合關(guān)系曲線Fig.9 The shape function ψ-Z relation and fitting relation diagram

圖10 形狀函數(shù)σ-Z關(guān)系曲線Fig.10 The shape function σ-Z relation diagram

圖9擬合所得出的已燃相對質(zhì)量ψ與已燃相對厚度Z之間的函數(shù)關(guān)系為：

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

(10)

同時對應(yīng)

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

(11)

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

(12)

由于2個形狀函數(shù)式存在積分或者微分關(guān)系，只需求出一個即可。

5.2 實際與理論的對比分析

將實際形狀函數(shù)推導(dǎo)與理論形狀函數(shù)推導(dǎo)所得關(guān)系進(jìn)行畫圖比較[14]。其中，圖11中黑色線為實際推導(dǎo)形狀函數(shù)ψ-Z關(guān)系，紅色線為理論推導(dǎo)形狀函數(shù)ψ-Z關(guān)系；圖12中黑色線為實際推導(dǎo)形狀函數(shù)σ-Z關(guān)系，紅色線為理論推導(dǎo)形狀函數(shù)σ-Z關(guān)系[15]。

圖11 形狀函數(shù)ψ-Z實際與理論曲線Fig.11 Comparison diagram of actual and theoretical shape function ψ-Z

圖12 形狀函數(shù)σ-Z實際與理論曲線Fig.12 Comparison diagram of actual and theoretical shape function σ-Z

由圖11、圖12可看出，實際推導(dǎo)所得曲線與理論推導(dǎo)所得曲線幾乎重合，由此本文的實際公式是正確的。

6 結(jié)論

1) 實際燃燒形狀函數(shù)與理論分析吻合較好，為了得到更準(zhǔn)確的形狀函數(shù)及其形狀特征量，可以在劃分網(wǎng)格時將網(wǎng)格劃分的更小更細(xì)；

2) 可以考慮改變程序，通過改變程序中步長，即相當(dāng)于改變?nèi)紵俾剩玫较鄳?yīng)的形狀函數(shù)，所以此研究可以對計算多種燃燒速率下的形狀函數(shù)提供參考；

3) 由于實際推導(dǎo)下的體積與表面積求解方法為微元法，因此可以對求解多種規(guī)則形狀物體以及不規(guī)則形狀物體的形狀函數(shù)的數(shù)值模擬與計算提供參考。