張群峰， 閆盼盼， 黎軍

(1.北京交通大學(xué) 土木工程學(xué)院， 北京 100044; 2.中國航空工業(yè)集團(tuán)公司 沈陽飛機(jī)設(shè)計(jì)研究所， 遼寧 沈陽 110035)

內(nèi)埋式彈艙與彈體相互影響的精細(xì)模擬

張群峰1， 閆盼盼1， 黎軍2

(1.北京交通大學(xué) 土木工程學(xué)院， 北京 100044; 2.中國航空工業(yè)集團(tuán)公司 沈陽飛機(jī)設(shè)計(jì)研究所， 遼寧 沈陽 110035)

為了研究內(nèi)埋式武器艙彈體投放過程中下落彈體與艙體之間強(qiáng)烈的相互作用，利用基于Menter SSTk-ω湍流模型的分離渦模擬方法，結(jié)合六自由度剛體動(dòng)力學(xué)方程和重疊網(wǎng)格技術(shù)對(duì)某一簡(jiǎn)化開式彈艙和彈體模型的三維流場(chǎng)進(jìn)行了非定常計(jì)算，利用非平穩(wěn)信號(hào)處理方法平滑偽Wigner-Ville分布分析了艙內(nèi)測(cè)點(diǎn)壓力的時(shí)頻特性。研究結(jié)果表明：在彈體出艙過程中，剪切層被破壞，導(dǎo)致艙體內(nèi)流動(dòng)狀態(tài)發(fā)生較大改變，不再呈現(xiàn)開式艙體自持振蕩特性，艙體內(nèi)噪聲的各階單調(diào)聲模態(tài)不明顯，受彈體的影響，艙體內(nèi)強(qiáng)度較高的渦結(jié)構(gòu)主要集中于艙體后緣，使艙體后緣噪聲水平升高；彈體出艙后剪切層迅速重建，自持振蕩回路再次形成，艙體流動(dòng)狀態(tài)恢復(fù)為典型的純空艙流動(dòng)狀態(tài)，但受彈體頭部加速氣流的影響，剪切層不穩(wěn)定性增強(qiáng)，導(dǎo)致艙體內(nèi)部噪聲水平升高。彈體在出艙過程中由于剪切層影響受到抬頭力矩的作用，使得彈體具有較大的迎角。彈體出艙后，在迎角的影響下，彈體開始受豎直向上的力和低頭力矩的作用，豎直向上的力會(huì)阻礙彈體的下落。同時(shí)，彈體下落過程受艙體強(qiáng)非定常流場(chǎng)的影響，彈體受力也存在強(qiáng)烈的波動(dòng)。研究結(jié)果為內(nèi)埋武器艙優(yōu)化設(shè)計(jì)和制定武器安全投放控制規(guī)律提供參考。

流體力學(xué)； 內(nèi)埋艙； 武器分離； 重疊網(wǎng)格； SST分離渦模擬； 平滑偽Wigner-Ville分布

0 引言

武器裝載于內(nèi)埋艙有助于戰(zhàn)斗機(jī)降低阻力和減小雷達(dá)探測(cè)面積，其逐漸替代外掛式而被廣泛采用。然而當(dāng)艙門開啟時(shí)，開式流動(dòng)內(nèi)埋式武器艙內(nèi)部存在極其復(fù)雜的流動(dòng)物理現(xiàn)象，彈艙內(nèi)會(huì)形成大幅度的壓力脈動(dòng)，造成劇烈震蕩并產(chǎn)生刺耳的噪聲，對(duì)控制武器的電子系統(tǒng)及相關(guān)結(jié)構(gòu)造成疲勞破壞[1]。

為了研究內(nèi)埋式彈艙的復(fù)雜流場(chǎng)，研究人員通常把彈艙簡(jiǎn)化為矩形空腔來研究流場(chǎng)的穩(wěn)態(tài)和動(dòng)態(tài)特性，并已經(jīng)進(jìn)行了大量的試驗(yàn)及數(shù)值模擬研究，研究表明：腔體的長深比[2-3]、腔體形狀[4]、來流馬赫數(shù)[5]、來流邊界層厚度[6-7]及有無控制措施[8]等因素均會(huì)對(duì)腔體內(nèi)噪聲環(huán)境產(chǎn)生影響。基于簡(jiǎn)化純空腔流動(dòng)的研究結(jié)果，文獻(xiàn)[9-10]分別對(duì)艙體帶彈和彈體固定在艙體不同位置時(shí)的流場(chǎng)進(jìn)行了研究，給出了彈體固定時(shí)艙體流場(chǎng)特性。值得提出的是，飛機(jī)在進(jìn)行武器投放時(shí)，彈體下落過程中，彈艙剪切層存在破壞和重建過程。這種強(qiáng)烈的相互作用一方面使得空艙具有更強(qiáng)的非定常流動(dòng)特性；另一方面，在彈體穿越剪切層的過程中，受艙體內(nèi)強(qiáng)非定常流場(chǎng)及剪切層的影響，彈體受力出現(xiàn)較大幅度的波動(dòng)，彈體的下落姿態(tài)也發(fā)生較大的變化，這一階段很大程度上決定了整個(gè)發(fā)射過程的成敗。為了保護(hù)艙內(nèi)的機(jī)載設(shè)備和投放安全，同時(shí)保證彈體分離品質(zhì)，對(duì)彈體投放過程內(nèi)埋式艙體與彈體間相互影響進(jìn)行深入研究顯得尤為重要。

現(xiàn)階段采用風(fēng)洞試驗(yàn)技術(shù)研究內(nèi)埋投放主要有可控軌跡試驗(yàn)(CTS)和風(fēng)洞自由飛試驗(yàn)兩種，其中CTS基于準(zhǔn)定常假設(shè)，對(duì)于強(qiáng)耦合條件下的分離問題預(yù)測(cè)，無法獲得強(qiáng)非定常非線性效應(yīng)。自由飛技術(shù)可以克服上述問題，但是難以精確控制投放分離的初始條件，同時(shí)無法獲得瞬態(tài)流場(chǎng)的精細(xì)化流動(dòng)結(jié)構(gòu)。由此可見，由于受試驗(yàn)設(shè)備和試驗(yàn)安全的限制，采用試驗(yàn)手段研究彈體投放過程中彈體和彈艙的強(qiáng)耦合流動(dòng)特性具有一定的局限性。而數(shù)值模擬方法則有其優(yōu)越性，尤其是一些新型湍流模擬方法的出現(xiàn)，使得計(jì)算精度有較大幅度的提高。本文將通過求解三維非定常Navier-Stokes方程，采用基于SSTk-ω模型的SST-IDDES方法，并利用重疊網(wǎng)格技術(shù)[11-12]和求解剛體動(dòng)力學(xué)方程實(shí)現(xiàn)彈體的運(yùn)動(dòng)，研究艙門打開后到彈體發(fā)射完成整個(gè)過程中，彈體與彈艙之間強(qiáng)烈的相互作用。

1 數(shù)值計(jì)算方法

內(nèi)埋式艙體和彈體的內(nèi)外流場(chǎng)控制方程包括連續(xù)方程、動(dòng)量方程、能量方程、氣體狀態(tài)方程、湍流模型方程和幾何守恒方程(SCL)[13-14]。限于篇幅，這里僅介紹本研究中所用到的SSTk-ω模型和SST DES模型、控制方程的離散方法、重疊網(wǎng)格法及計(jì)算方法驗(yàn)證。

1.1 Menter SST k-ω模型

剪應(yīng)力修正的SSTk-ω模型混合了k-ω模型和k-ε模型的優(yōu)勢(shì)，在近壁面處使用k-ω模型，而在邊界層外部區(qū)域采用k-ε模型。SSTk-ω模型包含了修正的湍流黏性公式，考慮了湍流剪切應(yīng)力的效應(yīng)，能更精確地模擬逆壓梯度引起的分離和分離區(qū)大小，SSTk-ω模型公式[15-19]為

(1)

(2)

式中：ρ為密度；k為湍動(dòng)能；t為時(shí)間；U為速度向量；Ug為網(wǎng)格移動(dòng)速度，可以通過SCL方程進(jìn)行求解得到；μ為分子黏性系數(shù)；μt為湍流黏性系數(shù)

(3)

(4)

Pk為湍動(dòng)能生成項(xiàng)；ω為湍流耗散比，β=0.075；α=α1F1+α2(1-F1)，α1=5/9,α2=0.44；σk1=0.85；σω1=0.85；S為應(yīng)變率張量；

(5)

y為網(wǎng)格格心到最近壁面的距離，ν為運(yùn)動(dòng)黏度，

(6)

Pk=min(μtS2,10Cμρkω)；

(7)

Cμ=0.09；σω2=0.856.

1.2 SST DES模型

Spalart將大渦模擬(LES)方法和一方程 (S-A)湍流模型相結(jié)合，用S-A湍流模型和LES方法分別模擬以耗散和大渦輸運(yùn)為主要流動(dòng)特征的區(qū)域，并通過當(dāng)?shù)鼐W(wǎng)格尺度與壁面距離的比較和判斷，實(shí)現(xiàn)這兩種模擬方法的自動(dòng)轉(zhuǎn)換，構(gòu)建了分離渦模擬(DES)方法，被習(xí)慣性稱為DES97[20-21]。由于網(wǎng)格分布不當(dāng)時(shí)，會(huì)產(chǎn)生因雷諾應(yīng)力不匹配而形成的?；瘧?yīng)力損耗，甚至出現(xiàn)非物理的分離。為了解決這一問題，2006年，Spalart等通過對(duì)邊界層參數(shù)進(jìn)行控制，延遲了RANS湍流模型的作用，提出了Delayed-DES(DDES)模型[22]。

2008年Shur等[23]通過引入壁面?；鬁u模擬(WMLES)成功克服了對(duì)數(shù)律不匹配問題，并重新定義亞格子長度尺度為

ΔIDDES=min{max[cwdw,cwhmax,hwn],hmax},

(8)

提出了Improved-DDES。式中：hwn是垂直壁面方向的網(wǎng)格步長；dw為到壁面距離；cw為經(jīng)驗(yàn)常數(shù)，取0.15；hmax為hwn的最大值。Strelets[24]提出以SST模型代替S-A模型的SST-DES模型，使得新的混合模型兼有了SST模型的諸多優(yōu)點(diǎn)。2012年Gritskevich[25]對(duì)SST DDES和IDDES模型的常數(shù)進(jìn)行了進(jìn)一步的精確校準(zhǔn)，下述文中的SST DES均指的是SST-IDDES模型。

1.3 平滑偽Wigner-Ville分布時(shí)頻分析法

無論是彈體投放中彈體出艙穿越剪切層還是彈體出艙后頭部加速區(qū)與剪切層相互干擾，都會(huì)使得彈艙內(nèi)部壓力脈動(dòng)發(fā)生較大變化而改變艙體聲學(xué)特性。由于彈體下落的影響，艙體內(nèi)部脈動(dòng)壓力數(shù)據(jù)屬于非平穩(wěn)信號(hào)，因此不能采用傳統(tǒng)傅里葉變換來處理彈體下落過程中艙內(nèi)的壓力脈動(dòng)，此時(shí)必須采用時(shí)頻分析技術(shù)進(jìn)行艙內(nèi)監(jiān)測(cè)點(diǎn)壓力脈動(dòng)的頻譜特性分析。

對(duì)于非平穩(wěn)信號(hào)，現(xiàn)在已有許多時(shí)頻分析技術(shù)，如短時(shí)傅里葉變換(STFT)、小波變換(WT)、Wigner-Ville分布(WVD)及平滑偽Wigner-Ville分布(SPWVD)。STFT方法存在時(shí)間分辨率和頻率分辨率的矛盾，導(dǎo)致時(shí)頻聚集性很差。WT時(shí)頻聚集性比STFT優(yōu)異，但是計(jì)算復(fù)雜耗時(shí)很長，不利于較大數(shù)據(jù)量的計(jì)算。WVD屬于二次型時(shí)頻描述法，由于二次型時(shí)頻方法固有特性，嚴(yán)重的交叉干擾對(duì)其應(yīng)用范圍有較大限制。SPWVD方法繼承了WVD的優(yōu)點(diǎn)同時(shí)對(duì)交叉干擾進(jìn)行了平滑處理，因此在計(jì)算量和精度上都是適當(dāng)?shù)?。本文即選用SPWVD方法對(duì)彈體下落過程中艙體內(nèi)部壓力非平穩(wěn)信號(hào)進(jìn)行時(shí)頻分析，以得到艙體聲學(xué)特性變化。信號(hào)x(t)的SPWVD定義[26]為

SPWVDx(t,f)=

(9)

式中：g(u)h(τ)為實(shí)對(duì)稱窗函數(shù)，其長度可以獨(dú)立進(jìn)行選擇。

1.4 重疊網(wǎng)格和離散格式

重疊網(wǎng)格法中存在兩套或多套相互交疊的網(wǎng)格，分別稱為背景網(wǎng)格區(qū)域和重疊網(wǎng)格區(qū)域。在計(jì)算中通過將背景網(wǎng)格中被重疊區(qū)遮擋的網(wǎng)格隱藏，通常稱之為“挖洞”，形成一套計(jì)算網(wǎng)格。如果重疊網(wǎng)格區(qū)域發(fā)生移動(dòng)，則按一定規(guī)律重新進(jìn)行“挖洞”，來得到新的計(jì)算網(wǎng)格，即可實(shí)現(xiàn)重疊網(wǎng)格區(qū)域的運(yùn)動(dòng)。背景網(wǎng)格與重疊網(wǎng)格相交的邊界區(qū)域通過一定的方法進(jìn)行插值來實(shí)現(xiàn)變量值的傳遞，并保證數(shù)值解在交界處光滑過渡。

控制方程對(duì)流項(xiàng)的空間離散采用具有2階精度的Roe格式，擴(kuò)散項(xiàng)采用中心差分格式。選取修正的Venkatakrishnan[27]限制器保證2階精度插值且具有TVD性質(zhì)，同時(shí)又具有較小的數(shù)值耗散。非定常計(jì)算的時(shí)間離散涉及控制方程中每一個(gè)與Δt有關(guān)的項(xiàng)。物理量φ隨時(shí)間變化描述為

F(φ)包含所有的空間離散項(xiàng)。計(jì)算中引入雙重時(shí)間步法，即在控制方程中引入虛擬時(shí)間項(xiàng)，利用物理時(shí)間步求解真實(shí)解，而每一物理時(shí)間步通過虛擬時(shí)間迭代達(dá)到收斂。

1.5 數(shù)值方法驗(yàn)證

1.5.1 對(duì)空腔噪聲計(jì)算方法的驗(yàn)證

開式流動(dòng)腔體會(huì)在腔體唇口處形成剪切層，剪切層由于自身不穩(wěn)定性或外界擾動(dòng)形成脫落渦，脫落渦向后傳播并與腔后壁碰撞，產(chǎn)生強(qiáng)烈的噪聲形成聲波，聲波通過空腔內(nèi)部向上游傳遞，與腔體前壁相撞后激發(fā)了剪切層內(nèi)新的渦生成、發(fā)展與脫落，再次與后壁相撞，形成二次聲波，在滿足一定的空氣動(dòng)力條件和相位條件時(shí)，腔內(nèi)流動(dòng)形成自激振蕩。自激振蕩產(chǎn)生后會(huì)形成聲波反饋回路，誘發(fā)腔內(nèi)產(chǎn)生強(qiáng)烈噪聲，出現(xiàn)多個(gè)聲壓級(jí)峰值的激振頻率。

為了研究空腔氣動(dòng)聲學(xué)特性，歐洲QinetiQ組織對(duì)M219空腔進(jìn)行了多輪風(fēng)洞試驗(yàn)得到了可信的試驗(yàn)數(shù)據(jù)。M219腔體為典型開式空腔[28]，尺寸：長L為0.508 m，寬W和深D均為0.101 6 m，長深比L/D=5∶1. 模型幾何尺寸如圖1所示。

圖1 M219模型幾何形狀及監(jiān)測(cè)點(diǎn)分布Fig.1 Sketch of M219 cavity anddistribution of monitoring points

模型的計(jì)算網(wǎng)格數(shù)為1 300萬，腔體內(nèi)網(wǎng)格尺寸為1 mm，第一層邊界層網(wǎng)格取2×10-3mm,以保證y+值在1左右。來流馬赫數(shù)Ma=0.85，來流靜溫T=266.53 K，來流靜壓p=63 000 Pa. 與試驗(yàn)相同，在腔體底面等間距的布置了10個(gè)壓力監(jiān)測(cè)點(diǎn)。計(jì)算時(shí)間步選取1×10-5s，總的計(jì)算步為20 000步，總采樣時(shí)長為0.2 s. 對(duì)空腔流動(dòng)，可用Rossiter[29]、Heller等[30-31]改進(jìn)的Rossiter半經(jīng)驗(yàn)公式來估算各階振蕩頻率：

(10)

式中：m=1,2,3，…；r、k為常數(shù)分別取0.25、0.57；Ma∞為遠(yuǎn)方來流馬赫數(shù)；U∞為遠(yuǎn)方來流速度。由(10)式求得的各階頻率為148 Hz、363 Hz、578 Hz.

對(duì)得到的壓力時(shí)域分布曲線進(jìn)行傅里葉變換，得到各測(cè)點(diǎn)的聲壓級(jí)頻譜圖，如圖2所示。從圖2中可以看出求得的聲壓級(jí)頻譜曲線與試驗(yàn)測(cè)得的結(jié)果及公式計(jì)算得到的結(jié)果三者吻合得很好，證明本文所采用的數(shù)值方法和網(wǎng)格分布可以較為準(zhǔn)確地求解空腔噪聲問題。

圖2 腔底聲壓級(jí)分布Fig.2 SPL along the center line of cavity bottom

利用SPWVD時(shí)頻分析技術(shù)對(duì)驗(yàn)證算例中計(jì)算結(jié)果進(jìn)行時(shí)頻分析。取后緣監(jiān)測(cè)點(diǎn)10的壓力時(shí)域曲線進(jìn)行處理，得到其能量密度時(shí)頻分布如圖3所示。圖3中顏色條即表示能量密度SPWVDx(t,f)，其單位是Pa2/(s·Hz)，較高的能量密度意味著較強(qiáng)的噪聲水平。由圖3可知能量密度主要集中分布在3個(gè)頻率附近，分別對(duì)應(yīng)著M219腔體的1～3階模態(tài)，其中1階模態(tài)能量密度較弱，2階、3階模態(tài)能量密度較強(qiáng)，說明腔體內(nèi)噪聲以2階、3階模態(tài)為主導(dǎo)，1階模態(tài)噪聲水平較弱，這與驗(yàn)證算例得到的聲壓級(jí)頻譜圖相一致。同時(shí)可以看出各階模態(tài)的能量密度隨著時(shí)間增加并不保持為常數(shù)：在0～0.04 s時(shí)間內(nèi)能量密度主要分布在1階、2階模態(tài)，而0.04 s之后，1階、2階模態(tài)能量密度減弱，轉(zhuǎn)為集中在第3階模態(tài)。在0.14 s，2階模態(tài)能量密度又有所增加、3階模態(tài)能量密度有所減弱。這一現(xiàn)象即文獻(xiàn)[32]中所提到的腔體噪聲主導(dǎo)模態(tài)隨時(shí)間的轉(zhuǎn)換。通過以上分析可以看出，腔體內(nèi)的壓力脈動(dòng)信號(hào)嚴(yán)格意義上是非平穩(wěn)信號(hào)，因此通常采用的傅里葉變換是一種近似處理，而SPWVD方法可以精確地分析腔體內(nèi)氣動(dòng)聲學(xué)時(shí)頻特性。

圖3 M219腔體底部能量密度時(shí)頻分布圖Fig.3 Time-frequency distribution of energy density on the bottom of M219 cavity

1.5.2 彈體投放問題的驗(yàn)證

機(jī)翼/掛架/帶舵外掛物(WPFS)模型是美國發(fā)起的第一次分離投放計(jì)算流體力學(xué)計(jì)算驗(yàn)證時(shí)使用的試驗(yàn)?zāi)Ｐ?，CTS試驗(yàn)測(cè)試由阿諾德工程發(fā)展中心(AEDC)完成，具體的幾何尺寸及結(jié)果數(shù)據(jù)詳見參考文獻(xiàn)[33]。本文以該模型試驗(yàn)來驗(yàn)證采用本文所選用數(shù)值方法并應(yīng)用重疊網(wǎng)格技術(shù)和剛體動(dòng)力學(xué)方程模擬彈體投放過程的準(zhǔn)確性。

WPFS模型幾何外形如圖4所示，模型網(wǎng)格劃分為兩個(gè)區(qū)域，包含機(jī)翼及掛架的背景網(wǎng)格區(qū)域和包含彈體的重疊網(wǎng)格區(qū)域。背景網(wǎng)格和重疊網(wǎng)格區(qū)域交界處網(wǎng)格尺寸之比保持在1.0～1.2，以保證插值的精度。背景網(wǎng)格區(qū)域網(wǎng)格總數(shù)為1 150萬，重疊網(wǎng)格區(qū)域網(wǎng)格總數(shù)為320萬。壁面第一層網(wǎng)格尺度為2×10-3mm. 計(jì)算時(shí)間步取1×10-3s.

圖4 WPFS模型幾何外形Fig.4 Wing-pylon-store geometry

彈體發(fā)射時(shí)受到兩個(gè)彈射力的作用，大小分別為10 679.4 N和42 717.5 N，彈射力作用距離為0.1 m，作用時(shí)間約為0.05 s，詳細(xì)參數(shù)見表1. 總的計(jì)算時(shí)間為0.3 s.

表1 彈體參數(shù)及發(fā)射條件Tab.1 Missile parameters and launching conditions

圖5為3個(gè)方向彈體質(zhì)心位移隨時(shí)間變化曲線。從圖5中可以看出，計(jì)算得到的彈體位移變化與試驗(yàn)符合較好，由于彈射力作用，彈體下落方向位移要遠(yuǎn)大于另外兩個(gè)方向的位移。圖6為彈體滾轉(zhuǎn)、偏航、俯仰姿態(tài)角隨時(shí)間變化曲線，同樣由于彈體彈射力作用，在下落初期，彈體俯仰角度增長最快，但撤去彈射力后，在重力和氣動(dòng)力共同作用下俯仰角又開始減小，其余兩個(gè)方向姿態(tài)角保持持續(xù)增加的趨勢(shì)?？梢钥闯霰疚牟扇〉臄?shù)值方法可以準(zhǔn)確地模擬彈體下落運(yùn)動(dòng)。

圖5 彈體質(zhì)心位移隨時(shí)間變化Fig.5 Displacement of missile Cg versus time

圖6 彈體姿態(tài)角隨時(shí)間變化Fig.6 Missile attitude angle versus time

2 計(jì)算模型和網(wǎng)格劃分

2.1 計(jì)算模型

為了更接近真實(shí)地模擬內(nèi)埋艙彈體投放過程，本文選取的艙體計(jì)算模型為等比放大的M219腔體模型。艙體尺寸為4.57 m×0.914 m× 0.914 m，彈艙長深比保持5∶1. 彈體模型即為WPFS中彈體，其長度為3.38 m，彈徑0.508 m，質(zhì)心位于距彈頭1.42 m處，彈體質(zhì)量為907.2 kg，x、y、z軸3個(gè)方向的轉(zhuǎn)動(dòng)慣量分別為27 kg·m2、488 kg·m2、488 kg·m2.

選擇的坐標(biāo)系為x軸為逆航向，z軸向上，y軸向右，則按此坐標(biāo)系定義，彈射力為負(fù)，彈體抬頭為正、左偏航為正。

2.2 計(jì)算條件

來流馬赫數(shù)Ma=0.85，靜壓p=6.3×104Pa，靜溫T=266.5 K，基于每米的雷諾數(shù)Re=1.35×107. 來流迎角為α=0°，入口條件設(shè)置為遠(yuǎn)場(chǎng)自由來流條件，艙體壁面均采用無滑移壁面條件。彈體投放采用彈射投放方式，彈體所受的彈射力作用規(guī)律與WPFS中彈射力作用規(guī)律一致，詳見表1.

2.3 艙體網(wǎng)格尺寸選擇

為了能準(zhǔn)確地模擬大尺寸艙體內(nèi)噪聲問題，需要適當(dāng)?shù)臄?shù)值方法和恰當(dāng)?shù)木W(wǎng)格密度。對(duì)于本文所采用的數(shù)值方法已經(jīng)得到了驗(yàn)證，但是對(duì)于大尺寸艙體內(nèi)網(wǎng)格尺度的選擇，還需進(jìn)一步確定。本文將通過試算的方式來確定恰當(dāng)?shù)木W(wǎng)格尺度。在1.5.1節(jié)已經(jīng)看到，通過Rossiter公式可以準(zhǔn)確地預(yù)測(cè)艙體的各階模態(tài)頻率，因此本文將以Rossiter公式求得的頻率作為依據(jù)來檢驗(yàn)網(wǎng)格尺寸是否滿足要求，根據(jù)Rossiter公式求得各階模態(tài)的頻率分別為16 Hz、40 Hz、64 Hz. 本文分別選取了艙體內(nèi)4種網(wǎng)格尺寸，分別為16 mm、12 mm、8 mm和6 mm. 壁面第一層網(wǎng)格尺度保持與驗(yàn)證算例相一致。最終4種情況總網(wǎng)格數(shù)目分別為350萬、700萬、2 400萬和4 200萬。

計(jì)算時(shí)間步長為5×10-5s，總計(jì)算步為20 000步，總物理時(shí)間為1 s. 計(jì)算中按照相同的方式在艙體底面均勻布置10個(gè)壓力監(jiān)測(cè)點(diǎn)，來流條件保持與驗(yàn)證算例相同。取監(jiān)測(cè)點(diǎn)10將得到的結(jié)果進(jìn)行傅里葉變換，得到聲壓級(jí)頻譜如圖7所示。

圖7 不同網(wǎng)格密度下艙底聲壓級(jí)分布Fig.7 SPLs along the center line of cavity bottom at different mesh density

通過圖7中的對(duì)比可以看出，當(dāng)網(wǎng)格尺度從16 mm，減小到6 mm的過程中，計(jì)算得到的聲壓級(jí)均與理論公式求得的聲壓級(jí)基本相符。然而對(duì)于低頻區(qū)域的寬頻噪聲，16 mm和12 mm網(wǎng)格得到的寬頻噪聲偏大，8 mm與6 mm網(wǎng)格的寬頻噪聲結(jié)果較小且相差不大。綜合考慮計(jì)算精度、計(jì)算資源和計(jì)算時(shí)間，選取 8 mm的網(wǎng)格作為基準(zhǔn)網(wǎng)格尺度。本文后續(xù)艙體網(wǎng)格劃分及計(jì)算時(shí)間步長選取均采用本小節(jié)選取的尺度。

彈體網(wǎng)格劃分時(shí)保證交界處重疊網(wǎng)格尺寸與背景網(wǎng)格尺寸之比小于1.2. 彈體壁面第一層網(wǎng)格厚度同樣保持為2×10-3mm，總網(wǎng)格數(shù)目為450萬。

3 結(jié)果分析

3.1 彈體投放對(duì)艙體聲學(xué)特性影響

彈體下落算例中，在建立起流場(chǎng)之后，又對(duì)艙體帶彈狀態(tài)進(jìn)行了0.3 s的非定常計(jì)算，之后彈體開始下落，下落時(shí)長為0.45 s，總計(jì)算時(shí)間為0.75 s. 取艙體內(nèi)部前緣監(jiān)測(cè)點(diǎn)2和后緣監(jiān)測(cè)點(diǎn)8的壓力數(shù)據(jù)進(jìn)行SPWVD時(shí)頻分析，得到能量密度時(shí)頻分布圖，如圖8所示。

圖8 艙體底部監(jiān)測(cè)點(diǎn)能量密度時(shí)頻分布圖Fig.8 Energy density distributions on the bottom of cavity

通過能量密度時(shí)頻分布圖可以看出，在彈體下落之前的0.3 s內(nèi)，能量密度集中在各階模態(tài)對(duì)應(yīng)的3個(gè)頻率附近，由于艙體帶彈的影響，頻率值與Rossiter公式求得的3個(gè)模態(tài)的頻率16 Hz、40 Hz、64 Hz有小幅偏差。

在0.3～0.6 s這一階段，整個(gè)艙體內(nèi)氣動(dòng)聲學(xué)特性完全改變，從前后緣兩個(gè)測(cè)點(diǎn)都可以看出，能量密度不再分布于各階模態(tài)對(duì)應(yīng)的頻率上，而是集中于低頻區(qū)域，后緣監(jiān)測(cè)點(diǎn)的能量密度幅值較彈體下落之前存在著大幅度的增加，前緣點(diǎn)能量密度幅值雖然有所增加但幅度較小。這說明由于彈體的下落，艙體內(nèi)噪聲水平有所增高，后緣點(diǎn)噪聲增加幅度更大，噪聲環(huán)境更惡劣。通過圖9的截面馬赫數(shù)分布云圖可以看出，這段時(shí)間為彈體穿越艙體唇口剪切層出艙過程，在此階段艙體剪切層被破壞，開式空艙自持振蕩的回路被阻斷。因此艙體內(nèi)氣動(dòng)聲學(xué)特性發(fā)生顯著改變，不再呈現(xiàn)典型的開式艙體氣動(dòng)特性。同時(shí)由于彈體對(duì)剪切層的干擾以及彈體下落后艙體容積改變產(chǎn)生的容積效應(yīng)，使得艙體內(nèi)壓力脈動(dòng)更加劇烈，噪聲強(qiáng)度增大。圖10為艙體內(nèi)部Q值等值面分布圖。從圖10中可以看出，由于彈體下落的影響，艙體內(nèi)渦結(jié)構(gòu)主要分布在艙體后緣區(qū)域，與之相比艙體前部流動(dòng)較為穩(wěn)定，因此與艙體后緣相比前緣噪聲水平較低。彈體下落后，艙體后緣噪聲水平升高，增強(qiáng)的噪聲前傳使得前緣噪聲水平高于彈體下落前噪聲水平。

圖10 Q準(zhǔn)則等值面分布圖Fig.10 Iso-surface contribution of Q criterion

在0.60 s之后，能量密度又再次集中在3個(gè)模態(tài)對(duì)應(yīng)的頻率附近。但是艙體后緣點(diǎn)各階頻率上的能量密度較彈體下落前都有所增大。當(dāng)彈體頭部穿過剪切層之后，彈體完成出艙過程，艙體剪切層迅速重建，自持振蕩回路再次形成，艙體流動(dòng)狀態(tài)恢復(fù)為典型的純空艙流動(dòng)狀態(tài)。從圖9(c)的馬赫數(shù)分布云圖中可以看出，由于彈體前緣頭部的加速氣流與剪切層相互作用，使得剪切層不穩(wěn)定性增強(qiáng)，同時(shí)會(huì)使得剪切層內(nèi)渦傳播速度增加，脫落渦具有更大的動(dòng)能與后緣撞擊，導(dǎo)致艙體后緣壓力脈動(dòng)更加劇烈，因此艙體后部噪聲水平升高。

3.2 艙體對(duì)彈體投放的影響

內(nèi)埋彈體投放過程中存在彈體穿越剪切層過程，此時(shí)彈體受到的氣動(dòng)力呈現(xiàn)出強(qiáng)烈的非定常特性，對(duì)彈體的姿態(tài)產(chǎn)生較大的影響。

圖11、圖12分別為彈體豎直方向受力及速度隨時(shí)間變化曲線，圖13、圖14、圖15為彈體下落過程俯仰力矩、俯仰角速度及俯仰角隨時(shí)間變化曲線。從0.30 s～0.35 s彈體受彈射力作用，此時(shí)彈體所受合力較大，具有很大的加速度，短時(shí)間內(nèi)速度由0增加到-3.5 m/s，由于彈射力的作用，使得彈體在這一階段內(nèi)受抬頭力矩的作用，抬頭角速度迅速增加，但由于彈射力作用時(shí)間短，彈體迎角約1.6°. 0.35 s時(shí)彈射力消失，彈體在重力和氣動(dòng)力共同作用下下落出艙。在0.35～0.42 s時(shí)，彈體受力較平穩(wěn)，向下合力維持在10 kN左右，彈體下落的加速度減小，速度增幅下降但繼續(xù)保持穩(wěn)定增長。圖16為0.38 s艙體中心截面馬赫數(shù)分布云圖及彈體表面壓力分布云圖，圖16 (a)中截面為馬赫數(shù)分布，彈體表面為壓力分布，圖16 (b)為彈體下表面壓力分布云圖，圖16 (c)為彈體上表面壓力分布云圖。圖17為彈體上下表面中心線壓力系數(shù)分布曲線。由圖17可知：彈體此時(shí)具有一定的迎角，彈體頭部所在區(qū)域流動(dòng)速度較低，壓力較高；彈體尾部區(qū)域更靠近剪切層，所以流動(dòng)速度較高，壓力較低。因此彈體此階段受到較小抬頭力矩的作用，彈體俯仰角速度小幅度增加，彈體迎角持續(xù)增大，到0.42 s時(shí)，彈體迎角達(dá)到7°左右。

圖11 彈體豎直方向受力隨時(shí)間變化Fig.11 Vertical force of missile versus time

圖13 彈體俯仰力矩隨時(shí)間變化Fig.13 Pitch moment of missile versus time

圖14 彈體俯仰角速度隨時(shí)間變化Fig.14 Pitch angular velocity of missile versus time

圖15 彈體俯仰角隨時(shí)間變化Fig.15 Pitch angle of missile versus time

圖16 t=0.38 s時(shí)艙體截面馬赫數(shù)及彈體表面壓力分布云圖Fig.16 Mach number contours of cavity section and pressure distribution of missile for t=0.38 s

圖17 t=0.38 s時(shí)彈體上下表面對(duì)稱線壓力系數(shù)分布Fig.17 Pressure coefficients of symmetric lines on upper and lower surfaces for t=0.38 s

圖18 t=0.49 s時(shí)截面馬赫數(shù)及彈體表面壓力分布云圖Fig.18 Mach number contours of cavity section and pressure distribution of missile for t=0.49 s

圖19 t=0.49 s時(shí)彈體上下表面對(duì)稱線壓力系數(shù)分布Fig.19 Pressure coefficient distribution of symmetric lines on upper and lower surfaces for t=0.49 s

從0.42 s起，彈體開始穿越剪切層，圖18 (a)為0.49 s時(shí)艙體中心截面馬赫數(shù)分布云圖及彈體表面壓力分布云圖。圖18(b)、圖18(c)分別為彈體下表面和上表面壓力分布云圖。圖19為彈體上下表面中心線壓力系數(shù)分布曲線。從圖18 (a)中可知，彈體頭部位置剪切層受到下落彈體的阻礙，撞擊在彈體下表面，一部分動(dòng)能轉(zhuǎn)變?yōu)閴耗?。從圖18 (b)中可以看出，彈體頭部下表面受到剪切層沖擊的區(qū)域局部壓力顯著增高。從圖18(b)、圖18(c)及圖19中可以看出彈體上表面壓力明顯低于下表面。彈體頭部下表面的局部高壓區(qū)使得彈體俯仰力矩保持為抬頭力矩，上下表面較大的壓力差導(dǎo)致氣動(dòng)力合力方向?yàn)樨Q直向上，且隨著彈體出艙過程不斷增大。因此在彈體穿越剪切層階段，彈體俯仰角速度及抬頭角度持續(xù)增加，且彈體下落受到較大的阻礙，下落速度增速減緩。

隨著彈體進(jìn)一步下落，彈體頭部穿過剪切層，作用在彈體上的抬頭力矩減??；同時(shí)彈體受抬頭角速度的影響，彈體迎角持續(xù)增加，0.52 s時(shí)，彈體迎角增加到14°，彈體后部下表面壓力繼續(xù)增加，尤其是彈翼下表面壓力及高壓區(qū)面積也明顯增大。受此影響，最終使彈體合力矩變?yōu)榈皖^力矩，且隨著彈體的迎角增加、低頭力矩持續(xù)增大。當(dāng)俯仰力矩變?yōu)榈皖^力矩后，彈體抬頭角速度開始減小，但迎角持續(xù)增加。到0.53 s時(shí)向下的合力減小為0，隨后變?yōu)橄蛏系暮狭?，彈體下落速度不增反降。

直到0.67 s之后，彈體俯仰角速度減小到0并進(jìn)一步變?yōu)榈皖^角速度，彈體迎角達(dá)到最大值后開始減小，彈體的豎直方向合力開始下降，迎角的減小也使得彈體受到的低頭力矩逐漸降低。

4 結(jié)論

1) 彈體投放過程中內(nèi)埋艙體與彈體存在強(qiáng)烈的相互作用。彈體下落出艙穿越剪切層會(huì)使艙體唇口剪切層被破壞，導(dǎo)致開式艙體流動(dòng)狀態(tài)發(fā)生改變。艙體自持振蕩回路被阻斷，艙內(nèi)噪聲不再呈現(xiàn)明顯的各階單調(diào)聲模態(tài)。受彈體的影響，艙體內(nèi)渦結(jié)構(gòu)主要集中于艙體后緣，且強(qiáng)度較高使艙體后緣噪聲水平升高。

2)彈體出艙后，艙體唇口處剪切層迅速重建，艙體恢復(fù)典型空腔流動(dòng)特性。但是受彈體頭部加速氣流的影響，剪切層不穩(wěn)定性增強(qiáng)，脫落渦強(qiáng)度增高，與艙體后壁撞擊產(chǎn)生的噪聲更大，使得艙體內(nèi)部聲壓級(jí)水平總體上升。

3)彈體下落過程受艙體強(qiáng)非定常流場(chǎng)的影響，彈體受力及力矩均出現(xiàn)較大幅的波動(dòng)。彈體穿越剪切層時(shí)頭部受剪切層撞擊產(chǎn)生局部高壓區(qū)，使彈體產(chǎn)生抬頭力矩導(dǎo)致彈體姿態(tài)角發(fā)生改變出現(xiàn)抬頭。彈體姿態(tài)角的抬頭又使彈體受力發(fā)生改變，氣動(dòng)力阻礙彈體下落，影響飛機(jī)與彈體安全分離。

References)

[1] 馮必鳴, 聶萬勝, 車學(xué)科. 超聲速條件下內(nèi)埋式武器分離特性的數(shù)值分析[J]. 飛機(jī)設(shè)計(jì), 2009, 29(4): 1-5. FENG Bi-ming, NIE Wan-sheng, CHE Xue-ke. Simulation of the store separation from a cavity at supersonic speed[J]. Aircraft Design, 2009, 29(4)：1-5. (in Chinese)

[2] 楊黨國, 范召林, 李建強(qiáng), 等. 超聲速空腔流激振蕩與聲學(xué)特性研究[J]. 航空動(dòng)力學(xué)報(bào), 2010, 25(7): 1567-1572. YANG Dang-guo, FAN Zhao-lin, LI Jian-qiang. Cavity flow oscillations and aeroacoustic characteristics at supersonic speeds [J]. Journal of Aerospace Power, 2010, 25(7): 1567-1572. (in Chinese)

[3] 謝露, 艾俊強(qiáng), 李權(quán), 等. 長深比對(duì)空腔流動(dòng)與聲學(xué)特性的影響分析[J]. 航空工程進(jìn)展, 2014, 5(1): 18-24. XIE Lu, AI Jun-qiang, LI Quan, et al. Effect of length-to-depth ratio on flow and aero acoustic characteristics of cavity [J]. Advances in Aeronautical Science and Engineering, 2014, 5(1): 18-24. (in Chinese)

[4] 徐路, 桑為民, 雷熙薇. 三維內(nèi)埋式彈艙流動(dòng)特性及形狀影響數(shù)值分析[J]. 應(yīng)用力學(xué)學(xué)報(bào), 2011, 28(1): 85-89. XU Lu, SANG Wei-min, LEI Xi-wei. Numerical analysis of flow characteristics and shape effect of three dimensional internal weapons bay [J]. Chinese Journal of Applied Mechanics, 2011, 28(1): 85-89. (in Chinese)

[5] 王一丁, 陳濱琦, 郭亮, 等. 空腔噪聲非線性數(shù)值模擬[J]. 國防科技大學(xué)學(xué)報(bào), 2015, 37(4): 151-157. WANG Yi-ding, CHEN Bin-qi, GUO Liang, et al. Nonlinear numerical simulation of cavity noise [J]. Journal of National University of Defense Technology, 2015, 37(4): 151-157. (in Chinese)

[6] 張群峰,閆盼盼,黎軍.邊界層厚度對(duì)腔體氣動(dòng)聲學(xué)特性影響數(shù)值模擬[J].航空動(dòng)力學(xué)報(bào), 2016, 31(3): 717-725. ZHANG Qun-feng, YAN Pan-pan, LI Jun. Numerical simulation on influence of boundary-layer thickness to cavity aeroacoustic characteristics [J]. Journal of Aerospace Power, 2016, 31(3):717-725. (in Chinese)

[7] 楊黨國,羅新福,李建強(qiáng).來流邊界層厚度對(duì)開式空腔氣動(dòng)聲學(xué)特性的影響分析[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2011, 29(4): 486-490. YANG Dang-guo, LUO Xin-fu, LI Jian-qiang. Analysis of aeroacoustic characteristics in open cavities influenced by boundary-layer thickness[J]. Acta Aerodynamica Sinica, 2011, 29(4): 486-490. (in Chinese)

[8] 黎軍, 李天, 張群峰. 開式流動(dòng)腔體的流動(dòng)機(jī)理與控制[J]. 實(shí)驗(yàn)流體力學(xué),2008,22(1):80-83. LI Jun, LI Tian, ZHANG Qun-feng. The mechanism and control of open cavity flow[J].Journal of Experiments in Fluid Mechanics, 2008, 22(1): 80-83. (in Chinese)

[9] 黃長強(qiáng), 國海峰, 唐上欽, 等. 超聲速帶彈武器艙氣動(dòng)特性數(shù)值研究[J]. 兵工學(xué)報(bào),2013, 34(8): 975-980. HUANG Chang-qiang, GUO Hai-feng, TANG Shang-qin, et al. Numerical research on aerodynamic characteristics of cavity with store at supersonic speeds[J].Acta Armamentarii, 2013, 34(8): 975-980. (in Chinese)

[10] 司海青, 王同光. 數(shù)值模擬有外掛物的空腔流動(dòng)[J].空氣動(dòng)力學(xué)學(xué)報(bào),2007, 25(3): 404-409. SI Hai-qing, WANG Tong-guang. Numerical simulations of the cavity with a store [J]. Acta Aerodynamica Sinica,2007, 25(3):404-409. (in Chinese)

[11] 朱自強(qiáng). 應(yīng)用計(jì)算流體力學(xué)[M]. 北京：北京航空航天大學(xué)出版社, 1998: 173-174. ZHU Zi-qiang. The application of computational fluid dynamics [M]. Beijing: Beijing University of Aeronautics and Astronautics Press, 1998: 173-174. (in Chinese)

[12] 閻超. 計(jì)算流體力學(xué)方法及應(yīng)用[M].北京：北京航空航天大學(xué)出版社，2006: 197-217. YAN Chao. The computational fluid dynamics method and its application [M]. Beijing, Beijing University of Aeronautics and Astronautics Press, 2006: 197-217. (in Chinese)

[13] Tamura Y, Fujii K. Conservation law for moving and transformed grids[C]∥Proceedings of 6th National Symposium on Computational Fluid Dynamics. Tokyo, Japan: Japan Society of Computational Fluid Dynamics, 1993: 519-522.

[14] 劉偉. 細(xì)長機(jī)翼搖滾機(jī)理的非線性動(dòng)力學(xué)分析及數(shù)值模擬方法研究[D]. 長沙: 國防科學(xué)技術(shù)大學(xué), 2004. LIU Wei. Nonlinear dynamics analysis for mechanism of slenderwing rock and study of numerical simulation method [D]. Changsha： National University of Defense Technology, 2004. (in Chinese)

[15] Durbin P A, Pettersson B A. Statistical theory and modeling for turbulent flows[M]. 2nd ed. New York: A John Wiley Sons Ltd, 2011.

[16] Wilcox D C. Turbulence modeling for CFD[M]. La Canada, CA: DCW Industries, 1998: 142-148.

[17] Menter F R. Two-equation eddy-viscosity turbulence modeling for engineering applications [J]. AIAA Journal, 1994, 32(8): 1598-1605.

[18] Menter F R, Kuntz M. A daptation of eddy viscosity turbulence models to unsteady separated flow behind vehicles[M]∥McCallen R, Browand F, Ross J. The aerodynamics of heavy vehicles: trucks, buses and trains. Berlin: Springer Berlin Heidelberg, 2004:339-352.

[19] Menter F R, Kuntz M, Langtry R. Ten years of industrial experience with the SST turbulence model [C]∥ Proceedings of the Fourth International Symposium on Turbulence, Heat and Mass Transfer. Danbury, CT: Begell House, 2003: 625-632.

[20] Spalart P R. Detached eddy simulation [J]. Annual Review of Fluid Mechanics, 2009, 41(1):203-229.

[21] Spalart P R, Jou W H, Strelets M, et al. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[C]∥1st AFOSR international Conference on DNS/LES. Columbus, OH: AFOSR, 1997.

[22] Spalart P R, Deck S, Shur M L, et al. A new version of detached-eddy simulation, resistant to ambiguous grid densities[J]. Theoretical and Computational Fluid Dynamics, 2006, 20(3): 181-195.

[23] Shur M L, Spalart P R, Strelets M K, et al. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities [J]. International Journal of Heat and Fluid Flow, 2008, 29(6): 1638-1649.

[24] Strelets M. Detached eddy simulation of massively separated flows [C]∥39th AIAA Aerospace Sciences Meeting and Exhibit. Reno, NV US: AIAA, 2001.

[25] Gritskevich M. Development of DDES and IDDES formulations for thek-ωshear stress transport model [J]. Flow, Turbulence and Combustion, 2012, 88(3): 431-449.

[26] Auger F, Flandrin P. Improving the readability of time-frequency and time-scale representations by the reassignment method [J]. IEEE Transactions on Signal Processing, 1995, 43(5): 1068-1089.

[27] Venkatakrishnan V. On the accuracy of limiters and convergence to steady state solutions[C]∥ 31st Aerospace Sciences Meeting and Exhibit. Reno, NV, US: AIAA, 1994:152-164.

[28] Henshaw M J C. M219 cavity case: verification and validation data for computational unsteady aerodynamics, RTO-TR-26[R]. London: QinetiQ, 2002.

[29] Rossiter J E. Wind tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds [R]. Farnborough: Aeronautical Research Council, 1964.

[30] Heller H H, Holmes D G, Covert E E. Flow induced pressure oscillations in shallow cavities [J]. Journal of Sound and Vibration, 1971, 18(4):545-546.

[31] Heller H H, Bliss D B. The physical mechanism of flow-induced pressure fluctuations in cavities and concepts for their suppression [C]∥2nd AIAA Aeroacoustics Conference. Hampton,VA,US: AIAA, 1975.

[32] Zhuang N. Experimental investigation of supersonic cavity flow and their control [D]. Tallahassee, Florida, US: Florida State University, 2007.

[33] Heim E R. CFD wing/pylon/finned store mutual interference wind tunnel experiment, AEDC-TSR-91-P4 [R]. Tennessee: Arnold Engineering Development Center, 1991.

Elaborate Simulation of Interaction Effect between Internal Weapon Bay and Missile

ZHANG Qun-feng1，YAN Pan-pan1，LI Jun2

(1.School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China;2.Shenyang Aircraft Design and Research Institute,Aviation Industry Corporation of China，Shenyang 110035, Liaoning, China)

To study the strong interaction effect between internal weapon bay and missile during the weapon release, the unsteady flow around a simplified weapon bay and missile model is simulated using SST k-ω improved delayed detached eddy simulation (IDDES) method, six-degrees-of-freedom rigid body dynamics equations and overset mesh method. The fluctuation of pressure in the weapon bay is analyzed using the smooth pseudo Winger-Vile distribution (SPWVD) method, and the time-frequency characteristics of pressure are obtained. The research results show that, when a missile is delivered from the weapon bay, the shear layer is destroyed to lead to the change in the flow structure, the self-sustained oscillation disappears in the weapon bay, and there is no obvious cavity tone. Stronger vortexes concentrate in the trailing edge of cavity, and thus the sound pressure level (SPL) in the trailing edge of the weapon bay increases. After the missile is delivered from the weapon bay, the shear layer is rapidly reestablished and the self-sustained oscillation shows up again. The impact of the accelerated flow near the head of the missile leads to enhance the instability in shear layer. Accordingly, the sound pressure level in the weapon bay is enhanced. The missile flies at a large angle of attack since it is affected by upward force moment during passing through the shear layer. After the missile is delivered from the weapon bay, it begins to be affected by vertical upward force and pitch down moment. The vertical upward force can hinder the missile from falling. Meanwhile, the force and moment acting on the missile fluctuate strongly under the effect of the unsteady flow in the weapon bay.

fluid mechanics; internal weapon bay; weapon release; overset mesh; SST DES model; smooth pseudo Winger-Vile distribution

2016-05-03

國家自然科學(xué)基金項(xiàng)目(11172283)

張群峰(1972—)，男，講師，博士。E-mail: zhangqunfeng@263.net

V211.3

A

1000-1093(2016)12-2366-11

10.3969/j.issn.1000-1093.2016.12.024