姜婷婷，錢 坤，陳偉芳

(浙江大學(xué) 航空航天學(xué)院，杭州 310027)

高超聲速混合流動(dòng)的粒子模擬耦合算法*

姜婷婷，錢 坤，陳偉芳

(浙江大學(xué) 航空航天學(xué)院，杭州 310027)

通過(guò)采用基于三維非結(jié)構(gòu)網(wǎng)格的子網(wǎng)格技術(shù)，對(duì)適用于高超聲速混合流動(dòng)的自適應(yīng)時(shí)間步長(zhǎng)粒子模擬耦合算法(Improved Hybrid Particle Simulation Method, IHPSM)進(jìn)行了改進(jìn)，在保證算法計(jì)算效率的同時(shí)降低數(shù)值誤差。通過(guò)對(duì)三維雙曲鈍錐外形的數(shù)值仿真及與DSMC(Direct Simulation Monte Carlo)計(jì)算結(jié)果的對(duì)比分析，證明改進(jìn)的IHPSM算法具有較高的計(jì)算精度，且能較大幅度提高計(jì)算效率?；诟倪M(jìn)的IHPSM算法，文中針對(duì)雙曲鈍錐外形進(jìn)行了稀薄氣體效應(yīng)和飛行馬赫數(shù)對(duì)高超聲速流動(dòng)影響規(guī)律的研究。結(jié)果表明，氣體稀薄程度的增加會(huì)減緩雙曲鈍錐前端流場(chǎng)宏觀物理量的變化梯度，使流場(chǎng)激波結(jié)構(gòu)變厚；來(lái)流馬赫數(shù)的增大會(huì)使激波明顯增強(qiáng)，但對(duì)激波厚度與結(jié)構(gòu)的影響較小。

子網(wǎng)格技術(shù)；高超聲速混合流動(dòng)；粒子模擬耦合算法；IHPSM

0 引言

錢學(xué)森[1]根據(jù)Kn數(shù)的數(shù)值大小，將流動(dòng)按照氣體的稀薄程度劃分為連續(xù)流區(qū)、滑移流區(qū)、過(guò)渡流區(qū)與自由分子流區(qū)。在連續(xù)流區(qū)，通常通過(guò)數(shù)值求解N-S(Navier-Stokes)方程研究飛行器的氣動(dòng)特性；在氣體稀薄程度很大的區(qū)域，N-S方程失效，通常通過(guò)數(shù)值求解簡(jiǎn)化的Boltzmann方程或者DSMC(Direct Simulation Monte Carlo)方法來(lái)研究飛行器的氣動(dòng)特性。然而，高超聲速飛行器在近連續(xù)流區(qū)以及過(guò)渡流區(qū)飛行時(shí)，高速自由來(lái)流與飛行器控制面或翼面上的激波、邊界層相互作用，會(huì)出現(xiàn)連續(xù)流和稀薄流共存的現(xiàn)象[2-4]，即混合流動(dòng)現(xiàn)象。在混合流動(dòng)中，局部稀薄流的存在使得基于連續(xù)介質(zhì)模型假設(shè)的經(jīng)典連續(xù)性方程在理論上失效，而局部高密度區(qū)的存在，使得DSMC方法因其巨大的計(jì)算資源耗費(fèi)量亦失去應(yīng)用價(jià)值。近年來(lái)，N-S/DSMC耦合算法在求解這類流動(dòng)問(wèn)題中受到學(xué)術(shù)界廣泛關(guān)注[5-6]。但由于N-S方程與DSMC方法模擬本質(zhì)的巨大差異，N-S/DSMC耦合算法在分區(qū)并行計(jì)算、界面信息傳遞、復(fù)雜外形適應(yīng)性等方面仍面臨巨大挑戰(zhàn)。

不同于N-S/DSMC耦合算法的思想，Macrossan[7]提出了一種粒子模擬耦合算法(Hybrid Particle Simulation Method, HPSM)。該方法采用DSMC方法計(jì)算流場(chǎng)的局部稀薄流區(qū)，采用基于DSMC方法發(fā)展而來(lái)的平衡粒子模擬方法(Equilibrium Particle Simulation Method, EPSM)[8]計(jì)算流場(chǎng)的局部連續(xù)流區(qū)。與DSMC方法類似，EPSM對(duì)仿真分子進(jìn)行數(shù)值模擬。但不同的是，EPSM省去了復(fù)雜的碰撞過(guò)程，而是在遵守能量守恒和動(dòng)量守恒的基礎(chǔ)上，通過(guò)從Maxwell平衡分布中隨機(jī)抽樣來(lái)得到仿真分子碰撞之后的狀態(tài)，因此提高了計(jì)算效率?；贛acrossan的思想，陳偉芳[9-10]、朱榮麗[11]等研究者進(jìn)一步在計(jì)算模型、計(jì)算效率等方面發(fā)展了粒子模擬耦合算法，并通過(guò)相關(guān)高超聲速典型算例驗(yàn)證了算法的有效性。由于在整個(gè)計(jì)算區(qū)域內(nèi)均采用粒子模擬方法進(jìn)行數(shù)值求解，HPSM具有物理模型清晰、區(qū)域邊界條件自然、界面信息傳遞簡(jiǎn)單等優(yōu)點(diǎn)。然而，盡管如此，目前國(guó)內(nèi)外的研究者針對(duì)HPSM所做的研究仍不充分，若要實(shí)現(xiàn)其廣泛的應(yīng)用，仍需解決諸多難題。首先，傳統(tǒng)的連續(xù)介質(zhì)假設(shè)失效判據(jù)在HPSM中并不適用，因此亟需建立一套合理的劃分DSMC計(jì)算區(qū)域和EPSM計(jì)算區(qū)域的流場(chǎng)區(qū)域劃分判據(jù)。其次，雖然HPSM已經(jīng)在一定程度上較DSMC方法提高了計(jì)算效率，但距離工程應(yīng)用的需求仍有較大差距。此外，粒子模擬方法(無(wú)論是DSMC方法還是EPSM)對(duì)計(jì)算網(wǎng)格的尺寸較為敏感[12-13]，當(dāng)計(jì)算網(wǎng)格尺寸遠(yuǎn)大于當(dāng)?shù)仄骄肿幼杂沙虝r(shí)(這種情況在高密度的局部連續(xù)流區(qū)不可避免)，會(huì)導(dǎo)致該方法產(chǎn)生較大的數(shù)值誤差。在前期研究工作[14]中，本文作者針對(duì)前兩個(gè)問(wèn)題，構(gòu)建了適用于HPSM算法的平衡態(tài)失效判據(jù)(Equilibrium Breakdown Parameter)和自適應(yīng)時(shí)間步長(zhǎng)方法，從而建立了改進(jìn)的HPSM算法(Improved Hybrid Particle Simulation Method, IHPSM)。驗(yàn)證算例表明，IHPSM在保證計(jì)算精度的同時(shí)較大程度的提高了計(jì)算效率。

本文針對(duì)第3個(gè)問(wèn)題，在IHPSM算法的基礎(chǔ)上引入三維非結(jié)構(gòu)子網(wǎng)格技術(shù)，以減小算法對(duì)網(wǎng)格尺寸的依賴性。同時(shí)，本文通過(guò)對(duì)超聲速雙曲鈍錐外形的數(shù)值計(jì)算來(lái)驗(yàn)證引入子網(wǎng)格技術(shù)的IHPSM的可行性。此外，本文采用IHPSM算法開(kāi)展了稀薄氣體效應(yīng)和飛行速度對(duì)高超聲速雙曲鈍錐外形流場(chǎng)結(jié)構(gòu)的影響規(guī)律研究。

1 自適應(yīng)時(shí)間步長(zhǎng)IHPSM算法

IHPSM[14]通過(guò)平衡態(tài)失效判據(jù)劃分計(jì)算區(qū)域，然后在滿足連續(xù)介質(zhì)模型假設(shè)且達(dá)到局部平衡態(tài)分布的流動(dòng)區(qū)域使用EPSM進(jìn)行計(jì)算，其他區(qū)域使用DSMC方法進(jìn)行計(jì)算[15]，從而實(shí)現(xiàn)對(duì)整個(gè)流場(chǎng)的數(shù)值求解。

1.1 平衡態(tài)失效判據(jù)

在N-S/DSMC耦合算法中，通常采用連續(xù)介質(zhì)模型假設(shè)失效判據(jù)來(lái)劃分N-S方程與DSMC計(jì)算區(qū)域。而EPSM適用于氣體狀態(tài)達(dá)到近平衡態(tài)分布時(shí)的情況，它相當(dāng)于通過(guò)數(shù)值模擬仿真分子的運(yùn)動(dòng)與碰撞過(guò)程來(lái)求解Euler方程。因此，EPSM的適用條件更為嚴(yán)苛。在IHPSM中，首先通過(guò)連續(xù)介質(zhì)假設(shè)失效判據(jù)來(lái)劃分局部連續(xù)流區(qū)和局部稀薄流區(qū)，然后在局部連續(xù)流區(qū)通過(guò)“平衡態(tài)失效參數(shù)[14]”判斷網(wǎng)格內(nèi)分子狀態(tài)是否處于局部平衡態(tài)分布，即是否滿足EPSM的使用條件。

在此，采用的連續(xù)介質(zhì)假設(shè)失效判據(jù)為當(dāng)?shù)乜伺瓟?shù)[2]，其定義為

(1)

式中λ為平均分子自由程；ρ為氣體密度；T為平動(dòng)溫度；u為速度；a為聲速。

當(dāng)網(wǎng)格內(nèi)KnGLL,max>0.05時(shí)，為稀薄流動(dòng)區(qū)域，采用DSMC方法來(lái)模擬該區(qū)域的碰撞過(guò)程；當(dāng)KnGLL,max<0.05時(shí)，為連續(xù)流動(dòng)區(qū)域，此時(shí)需要繼續(xù)判斷該網(wǎng)格是否達(dá)到平衡分布狀態(tài)?；贐ird[16-17]的研究，選取網(wǎng)格中的分子碰撞數(shù)作為判斷網(wǎng)格內(nèi)的流動(dòng)是否達(dá)到平衡態(tài)的準(zhǔn)則。分子碰撞數(shù)f的定義：

(2)

式中 Δt為計(jì)算時(shí)間步長(zhǎng)；γ為分子碰撞頻率。

對(duì)于VHS模型來(lái)說(shuō)，分子碰撞頻率γ的表達(dá)式：

式中dref為氣體分子在參考溫度為0 ℃時(shí)的參考直徑；n為氣體分子數(shù)密度；kB為Boltzmann常數(shù)；ms為氣體分子質(zhì)量；Tref為參考溫度；ω為粘性指數(shù)。

在一個(gè)時(shí)間步長(zhǎng)內(nèi)，若網(wǎng)格中的平均分子碰撞數(shù)f≥10，則網(wǎng)格內(nèi)的分子運(yùn)動(dòng)狀態(tài)將達(dá)到平衡態(tài)，即可采用EPSM進(jìn)行計(jì)算。

1.2 自適應(yīng)時(shí)間步長(zhǎng)方法

為了提高計(jì)算效率，在IHPSM中，每個(gè)網(wǎng)格的時(shí)間步長(zhǎng)不再保持統(tǒng)一，而是在計(jì)算過(guò)程中不斷調(diào)整，以適應(yīng)流場(chǎng)物理量的變化。

在自適應(yīng)時(shí)間步長(zhǎng)方法[14]中，定義3種與時(shí)間步長(zhǎng)相關(guān)的變量：基時(shí)間步長(zhǎng)Δt、時(shí)間區(qū)Ti(Ti≥1，且為整數(shù)，i為網(wǎng)格單元編號(hào))以及網(wǎng)格時(shí)間步長(zhǎng)Δti。基于DSMC方法對(duì)時(shí)間步長(zhǎng)的要求，選取初始基時(shí)間步長(zhǎng)[11]Δt：

(3)

式中i為網(wǎng)格編號(hào)；Δxi、Δyi、Δzi為網(wǎng)格i在3個(gè)方向的參考長(zhǎng)度；u∞、v∞、w∞為來(lái)流氣體運(yùn)動(dòng)速度；γ∞為來(lái)流氣體碰撞頻率。

Δt將作為整個(gè)流場(chǎng)推進(jìn)過(guò)程中的最大時(shí)間步長(zhǎng)，且不隨著推進(jìn)過(guò)程發(fā)生改變。在高超聲速流場(chǎng)中，大多數(shù)網(wǎng)格的分子碰撞頻率隨著時(shí)間推進(jìn)過(guò)程而增加。因此，這個(gè)基時(shí)間步長(zhǎng)對(duì)于大多數(shù)網(wǎng)格來(lái)講不會(huì)出現(xiàn)因時(shí)間步長(zhǎng)過(guò)小而降低計(jì)算效率的情況。

確定了基時(shí)間步長(zhǎng)以后，網(wǎng)格內(nèi)當(dāng)?shù)貢r(shí)間步長(zhǎng)隨計(jì)算過(guò)程的自動(dòng)調(diào)整步驟可描述如下：

(1)將初始時(shí)間區(qū)設(shè)為Ti=1；

(2)每隔一個(gè)基時(shí)間步長(zhǎng)Δt，計(jì)算網(wǎng)格 的分子碰撞數(shù)fi，同時(shí)將時(shí)間區(qū)更新為Ti=int(fi)+1；

(3)令

(4)

1.3 子網(wǎng)格計(jì)算技術(shù)

1.4 IHPSM算法的主要計(jì)算步驟

IHPSM算法的主要計(jì)算步驟描述如下：

(1)按照式(3)求出初始基時(shí)間步長(zhǎng)Δt，并賦予各網(wǎng)格初始時(shí)間區(qū)Ti=1，則各網(wǎng)格的初始時(shí)間步長(zhǎng)均為Δt。

(2)在一個(gè)時(shí)間步長(zhǎng)內(nèi)，求出每個(gè)網(wǎng)格中的當(dāng)?shù)乜伺瓟?shù)KnGLL，max和碰撞數(shù)fi，同時(shí)定義碰撞數(shù)累計(jì)值CNi=CNi+fi。

(3)當(dāng)KnGLL，max≥0.05時(shí)：若fi≤1，則令網(wǎng)格屬性IDi=1；若fi>1，則令I(lǐng)Di=-1。

當(dāng)KnGLL，max<0.05時(shí)：若fi≤1，則令I(lǐng)Di=1；若1

為了更加清晰地描述該步驟，圖1給出了較詳細(xì)的計(jì)算流程。

(4)每個(gè)基時(shí)間步長(zhǎng)之后，對(duì)仿真分子進(jìn)行重新排序，同時(shí)統(tǒng)計(jì)仿真分子的物理信息，以便得到宏觀物理信息。

(5)重復(fù)步驟(2)～(4)，直到計(jì)算結(jié)束。

2 數(shù)值模擬結(jié)果及分析

2.1 高馬赫數(shù)雙曲鈍錐外形驗(yàn)證算例

為了驗(yàn)證自適應(yīng)時(shí)間步長(zhǎng)IHPSM計(jì)算程序，對(duì)高超聲速過(guò)渡流條件下的三維雙曲鈍錐縮比模型進(jìn)行數(shù)值計(jì)算，并與DSMC計(jì)算結(jié)果進(jìn)行對(duì)比。其中，DSMC計(jì)算結(jié)果由課題組內(nèi)部已經(jīng)經(jīng)過(guò)驗(yàn)證的DSMC計(jì)算軟件[18]計(jì)算得到。三維雙曲鈍錐外形雙曲面半角為42.5°，頭部半徑為0.013 62 m。來(lái)流介質(zhì)為空氣，計(jì)算條件：

圖3為IHPSM方法計(jì)算得到的流場(chǎng)馬赫數(shù)云圖、壓力云圖和溫度云圖。從流場(chǎng)馬赫數(shù)云圖可看出，氣流經(jīng)過(guò)激波后開(kāi)始減速，從超聲速流動(dòng)減為亞聲速流動(dòng)。流場(chǎng)壓力云圖表明，激波后壓力增大，在駐點(diǎn)處壓力達(dá)到最大值。激波后的溫度變化則是先增大而后在駐點(diǎn)處又稍微降低，這主要是由等溫壁條件所致。

圖4給出了IHPSM方法、DSMC方法駐點(diǎn)線溫度、密度的分布曲線，通過(guò)對(duì)比可看出，IHPSM方法獲得了與DSMC計(jì)算方法基本一致的計(jì)算結(jié)果，表明IHPSM方法可較精確地描述高超聲速混合流動(dòng)的流場(chǎng)結(jié)構(gòu)。

IHPSM方法和DSMC方法的測(cè)試計(jì)算機(jī)配置均為2 GHz的AMD處理器，內(nèi)存24 G。2種方法均采用單核進(jìn)行計(jì)算。流場(chǎng)達(dá)到穩(wěn)定狀態(tài)時(shí)，DSMC方法所耗費(fèi)的計(jì)算時(shí)間約為720 h，而本文改進(jìn)的IHPSM方法花費(fèi)的計(jì)算時(shí)間僅為310 h。因此，在相同的仿真分子運(yùn)動(dòng)時(shí)間內(nèi)，IHPSM方法在取得較高預(yù)測(cè)精度的前提下，顯著地提高了計(jì)算效率。

2.2 雙曲鈍錐外形流場(chǎng)結(jié)構(gòu)的變化規(guī)律研究

本節(jié)采用IHPSM算法開(kāi)展氣體稀薄效應(yīng)與飛行速度對(duì)三維雙曲鈍錐外形流場(chǎng)結(jié)構(gòu)的影響規(guī)律的研究。

首先，在來(lái)流馬赫數(shù)為20、物面溫度為500 K的條件下，取飛行高度分別為75、80、85、90 km的計(jì)算狀態(tài)進(jìn)行數(shù)值仿真，并對(duì)流場(chǎng)結(jié)果進(jìn)行對(duì)比分析。不同飛行高度下計(jì)算得到的駐點(diǎn)線密度和溫度分布曲線如圖5所示。

從圖5中各物理量變化趨勢(shì)可看出，隨著飛行高度的增大，氣體越來(lái)越稀薄，鈍錐前端流場(chǎng)中密度、溫度的變化梯度逐漸減小，流場(chǎng)激波結(jié)構(gòu)越來(lái)越厚。這是因?yàn)?，隨著來(lái)流氣體稀薄程度的增大，分子的碰撞頻率減小，粘性作用也因此而減小，從而導(dǎo)致流場(chǎng)物理量變化率減緩，宏觀表現(xiàn)便是激波的厚度增加。

此外，在飛行高度80 km、物面溫度為500 K的條件下，通過(guò)改變來(lái)流速度來(lái)研究來(lái)流馬赫數(shù)變化對(duì)流動(dòng)的影響。來(lái)流馬赫數(shù)分別取6、10、15、20。不同馬赫數(shù)狀態(tài)下計(jì)算得到的駐點(diǎn)線密度和溫度分布曲線如圖6所示。從圖6中可看出，隨著來(lái)流馬赫數(shù)增大，流體的動(dòng)能在頭部駐點(diǎn)區(qū)轉(zhuǎn)化為氣體內(nèi)能，激波明顯增強(qiáng)，導(dǎo)致激波后駐點(diǎn)區(qū)密度越來(lái)越大，而溫度急劇升高。但從密度與溫度分布來(lái)看，來(lái)流馬赫數(shù)對(duì)激波厚度的影響不大。

3 結(jié)論

(1)在IHPSM算法中加入三維非結(jié)構(gòu)子網(wǎng)格技術(shù)，降低了算法由于對(duì)計(jì)算網(wǎng)格的依賴性而產(chǎn)生的數(shù)值誤差。

(2)通過(guò)對(duì)三維雙曲鈍錐外形的數(shù)值仿真以及與DSMC計(jì)算結(jié)果的比對(duì)分析，表明改進(jìn)的IHPSM方法能同時(shí)滿足高超聲速混合流動(dòng)問(wèn)題對(duì)計(jì)算精度和計(jì)算效率的需求。

(3)針對(duì)雙曲鈍錐外形,進(jìn)行了稀薄氣體效應(yīng)和來(lái)流馬赫數(shù)對(duì)高超聲速流動(dòng)的影響規(guī)律研究。結(jié)果表明,隨著來(lái)流氣體稀薄程度的增加，鈍錐前端流場(chǎng)宏觀物理量的變化梯度會(huì)減緩，流場(chǎng)激波結(jié)構(gòu)越來(lái)越厚；隨著來(lái)流速度的增大，激波會(huì)明顯增強(qiáng)，但激波的厚度與結(jié)構(gòu)的變化不大。

[1] Tsien H S.Superaerodynamics,mechanics of rarefied gases[J].Journal of the Aeronautical Sciences,1946,13(12):653-664.

[2] Boyd I D,Chen G,Candler G V.Predicting failure of the continuum fluid equations in transitional hypersonic flows[J].Physics of Fluids,1995,7(1):210-219.

[3] Gnoffo P A.CFD validation studies for hypersonic flow prediction [R].AIAA 2001-1025.

[4] Hash D B,Hassan H A.Assessment of schemes for coupling monte carlo and navier-stokes solution methods[J].Journal of thermophysics and Heat Transfer,1996,10(2):242-249.

[5] Schwartzentruber T E,Boyd I D.A modular particle-continuum numerical method for hypersonic nonequilibrium gas flows[J].Journal of Computational Physics,2007,225:1159-1174.

[6] Lian Y Y,Wu J S,Cheng G,et al.Development of a parallel hybrid method for the DSMC and NS solver[R].AIAA 2005-0435.

[7] Macrossan M N.A particle-only hybrid method for near-continuum flows[C]//American Institute of Physics,2001,388-395.

[8] Pullin D I.Direct simulation methods for compressible ideal gas flow[J].Journal of Computational Physics,1980,34:231-244.

[9] 陳偉芳,吳明巧,任兵.DSMC/EPSM混合算法研究[J].計(jì)算力學(xué)學(xué)報(bào),2003,20(3):274-278.

[10] 吳明巧,陳偉芳,任兵.超聲速平頭圓柱繞流DSMC/EPSM混合算法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2002,20(2):206- 210.

[11] 朱榮麗,曹義華.高超聲速飛行器DSMC/EPSM自適應(yīng)當(dāng)?shù)貢r(shí)間步長(zhǎng)混合算法 [J].空氣動(dòng)力學(xué)學(xué)報(bào),2008,26(3):394-399.

[12] Macrossan M N.The equilibrium flux method for the calculation of flows with non-equilibrium chemical reactions[J].Journal of Computational Physics,1989,80(1):204-231.

[13] Breuer K S,Piekos E S,Gonzales D A.DSMC simulations of continuum flows[R].AIAA 1995-2088.

[14] Jiang T T,Xia C C,Chen W F.An improved hybrid particle scheme for hypersonic rarefied-continuum flow[J].Vacuum,2016,124:76-84.

[15] Laux M.Local time stepping with automatic adaptation for the DSMC method[R].AIAA 1998-2670.

[16] Bird G A.Near continuum impact of an under-expanded jet[C]//Proc.AIAA Computational Dynamics Conference,1973:103.

[17] Bird G A.Molecular gas dynamics and the direct simulation of gas flow[M].Clarendon:Oxford,1994.

[18] 吳雄,陳偉芳,石于中.過(guò)渡區(qū)球頭氣動(dòng)力特性DSMC仿真研究[J].航天返回與遙感,2002,23(3):1-5.

(編輯：呂耀輝)

Study on hybrid particle simulation method for hypersonic mixed flows

JIANG Ting-ting,QIAN Kun,CHEN Wei-fang

(School of Aeronautics and Astronautics,Zhejiang University,Hangzhou 310027,China)

A hybrid all-particle scheme(Improved Hybrid Particle Simulation Method,IHPSM),which is intended for hypersonic mixed rarefied-continuum flows,is extended for improved efficiency and wider applicability by employing the three-dimensional unstructured sub-cell technique.To assess the validity of the improved hybrid particle simulation scheme,a typical test case for three dimensional hyperbolic blunt cone is simulated.The results of IHPSM method for flow field are in good agreement with the standard DSMC results.Notably,the hybrid procedures allow a large reduction in overall simulation expense relative to DSMC.In addition,the relationship between the hypersonic flow field and the rarefied gas or Mach number is studied by simulating a series of test cases of the hyperbolic blunt cone.As the inflow gas becomes rarer,the gradient of the macroscopic properties in front of the hyperbolic blunt cone becomes smaller,and the shock wave is thicker.The increase of Mach number increases the values of the macroscopic properties in the shock wave,but has little influence on the shock wave structure.

sub-cell technique;hypersonic mixed flows;hybrid all-particle scheme;IHPSM

2016-03-09；

2016-03-21。

國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(2014CB340201)；國(guó)家自然科學(xué)基金(51575487)。

姜婷婷(1989—)，女，博士生，流體力學(xué)專業(yè)。E-mail：jiangyan8140@163.com

V411

A

1006-2793(2017)03-0283-06

10.7673/j.issn.1006-2793.2017.03.003