本文關(guān)鍵詞：推進(jìn)技術(shù)雜志，由筆耕文化傳播整理發(fā)布。

《推進(jìn)技術(shù)雜志》2015年第十二期

摘要：

通過數(shù)值仿真方法針對(duì)液體火箭發(fā)動(dòng)機(jī)內(nèi)的氣相化學(xué)動(dòng)力學(xué)與振蕩聲場(chǎng)的熱聲耦合機(jī)理進(jìn)行了研究。采用任意拉格朗日算法解耦流動(dòng)與化學(xué)源項(xiàng)間的剛性。采用的多步總包反應(yīng)機(jī)理考慮了底層的準(zhǔn)穩(wěn)態(tài)組分脈動(dòng)。通過入口流量邊界的流量脈動(dòng)向燃燒室中引入縱向聲波，并建立了冷流聲學(xué)相似場(chǎng)模型以分析熱聲耦合效應(yīng)的強(qiáng)度。研究發(fā)現(xiàn)：在線性小振幅聲場(chǎng)中，氣相化學(xué)動(dòng)力學(xué)控制的釋熱系統(tǒng)與聲學(xué)振蕩無明顯耦合激勵(lì)；在非線性有限幅值聲場(chǎng)中，燃燒室壓力與釋熱波動(dòng)出現(xiàn)“突躍”并表現(xiàn)為陡峭前沿波，氣相化學(xué)動(dòng)力學(xué)控制的釋熱系統(tǒng)與聲學(xué)振蕩發(fā)生耦合激勵(lì)，反應(yīng)流較其冷流聲學(xué)相似場(chǎng)的壓力振蕩振幅增強(qiáng)約300%。最后分析了耦合激勵(lì)發(fā)生的可能原因，提出了氣相化學(xué)動(dòng)力學(xué)體系的“釋熱分岔”假說。

關(guān)鍵詞：

熱聲耦合；化學(xué)動(dòng)力學(xué)；任意拉格朗日歐拉算法；冷流聲學(xué)相似場(chǎng)；數(shù)值仿真

1引言

液體火箭發(fā)動(dòng)機(jī)高頻燃燒不穩(wěn)定是一種熱聲耦合現(xiàn)象［1］，可由Rayleigh準(zhǔn)則進(jìn)行描述，這一準(zhǔn)則早在19世紀(jì)就被發(fā)現(xiàn)且至今仍在燃燒不穩(wěn)定研究領(lǐng)域享有基礎(chǔ)性地位。Rayleigh準(zhǔn)則本身是簡(jiǎn)單的，將其置于液體火箭發(fā)動(dòng)機(jī)具體的環(huán)境中卻立刻變得捉摸不清。液體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)是一個(gè)復(fù)雜的物理化學(xué)系統(tǒng)，其包含霧化、蒸發(fā)、湍流、燃燒等多個(gè)子過程，其中哪些子過程可以作為燃燒不穩(wěn)定的控制機(jī)制至今未有定論。為了研究釋熱與聲場(chǎng)的耦合機(jī)理，有必要對(duì)燃燒當(dāng)中的各子過程與聲場(chǎng)的相互作用機(jī)制進(jìn)行單獨(dú)研究［2～4］。長(zhǎng)期以來，特征時(shí)間較長(zhǎng)的噴霧和蒸發(fā)過程一直是燃燒不穩(wěn)定領(lǐng)域的重點(diǎn)研究對(duì)象［5，，6］，兩者在釋熱鏈中是熱釋放速率的決定性因素且與聲學(xué)過程的固有頻率相當(dāng)。除此之外，近年來相關(guān)文獻(xiàn)通過在仿真中采用LES或DES湍流模型捕捉瞬態(tài)湍流結(jié)構(gòu)來研究湍流對(duì)燃燒不穩(wěn)定的影響，其中比較有代表性的是普度大學(xué)Smith等的研究［7］，他們認(rèn)為大尺度的渦撞擊破碎引發(fā)的周期性釋熱脈動(dòng)可能成為燃燒不穩(wěn)定的驅(qū)動(dòng)機(jī)制。理論上說，液體火箭發(fā)動(dòng)機(jī)內(nèi)流場(chǎng)流動(dòng)燃燒的諸多子過程均是驅(qū)動(dòng)燃燒不穩(wěn)定的潛在因素，然而熱釋放的最終途徑化學(xué)動(dòng)力學(xué)過程由于特征時(shí)間極短而被長(zhǎng)期忽視，特別是針對(duì)氣相化學(xué)動(dòng)力學(xué)與聲場(chǎng)耦合機(jī)理的仿真研究相對(duì)匱乏。在十分有限的研究文獻(xiàn)當(dāng)中，均未對(duì)化學(xué)動(dòng)力學(xué)釋熱過程進(jìn)行充分的孤立，因而難以得到較為明確的結(jié)論，仍存在一系列重要的問題亟待被解答。鑒于此，本文構(gòu)建了可進(jìn)行考慮了底層基元反應(yīng)的多步總包化學(xué)動(dòng)力學(xué)計(jì)算的任意拉格朗日算法平臺(tái)，并在此基礎(chǔ)上研究了氣相化學(xué)反應(yīng)與聲學(xué)振蕩的耦合機(jī)理，提出了冷流聲學(xué)相似場(chǎng)的分析模型，探討了化學(xué)動(dòng)力學(xué)過程作為高頻燃燒不穩(wěn)定的控制因素的可能性。

2數(shù)值方法與模型

2.1數(shù)值算法想要深入研究氣相化學(xué)反應(yīng)釋熱過程對(duì)燃燒不穩(wěn)定的作用機(jī)制，就必須要引入多步化學(xué)反應(yīng)機(jī)理模型。然而，多步化學(xué)反應(yīng)的時(shí)間尺度遠(yuǎn)小于流場(chǎng)的流動(dòng)時(shí)間尺度，導(dǎo)致流動(dòng)控制方程中的組分、能量輸運(yùn)方程出現(xiàn)了剛性化學(xué)動(dòng)力學(xué)源項(xiàng)。這種流動(dòng)與化學(xué)源項(xiàng)之間的剛性耦合可以通過運(yùn)算分離算法消除［8］。本文采用任意拉格朗日歐拉（ALE）運(yùn)算分離算法作為CFD核心算法，此算法將一個(gè)時(shí)間步長(zhǎng)分為A，B，C三個(gè)階段進(jìn)行計(jì)算。前兩個(gè)階段為拉格朗日計(jì)算，不計(jì)算對(duì)流項(xiàng)，其中A階段只進(jìn)行化學(xué)動(dòng)力學(xué)源項(xiàng)的計(jì)算，B階段采用類似SIMPLE算法的半隱式迭代法求解壓力-速度耦合場(chǎng)；最后C階段為歐拉計(jì)算，求解對(duì)流項(xiàng)［9］。燃燒室內(nèi)流場(chǎng)采用歐拉坐標(biāo)系下的完全N-S方程進(jìn)行描述，N-S控制方程的通用形式參見文獻(xiàn)［10］。方程的空間離散為二階精度，時(shí)間離散采用一階向前差分。

2.2物理模型為了保證對(duì)化學(xué)動(dòng)力學(xué)釋熱過程的孤立分析，湍流模型采用標(biāo)準(zhǔn)k-ε湍流模型。大渦模擬、分離渦模擬等瞬態(tài)湍流模型可以捕捉到諸如漩渦脫落這種強(qiáng)瞬態(tài)湍流結(jié)構(gòu)，而這類瞬態(tài)結(jié)構(gòu)可能和聲學(xué)振蕩發(fā)生周期性耦合誘發(fā)不穩(wěn)定從而影響仿真結(jié)果的分析［7］，因而選取了雷諾時(shí)均的湍流模型�；瘜W(xué)動(dòng)力學(xué)模型采用通用有限速率化學(xué)反應(yīng)模型。湍流與燃燒的相互作用作如下處理：流動(dòng)為燃燒輸運(yùn)組分和參數(shù)，燃燒為流動(dòng)提供組分源項(xiàng)和能量源項(xiàng)�；瘜W(xué)反應(yīng)機(jī)理采用正十二烷多步總包機(jī)理，此機(jī)理包含了34組分、31步的總包機(jī)理及底層的52組分、217步的骨架機(jī)理，由正十二烷的詳細(xì)反應(yīng)機(jī)理［11］通過直接關(guān)系圖法（DRG）結(jié)合計(jì)算奇異值攝動(dòng)法（CSP）的機(jī)理簡(jiǎn)化方法獲得，可以考慮底層基元反應(yīng)的準(zhǔn)穩(wěn)態(tài)（QSS）組分的濃度脈動(dòng)。在求解QSS組分濃度時(shí)采用LQSSA方法，解耦QSS組分間的非線性耦合，并通過拓?fù)浞治龅慕Y(jié)果使用高斯直接消除法求解QSS濃度線性方程組［12］。采用剛性積分求解器VODE解化學(xué)動(dòng)力學(xué)方程組（1），并采用基于粒度細(xì)化的動(dòng)態(tài)負(fù)載平衡技術(shù)搭建MPI并行環(huán)境以提升化學(xué)動(dòng)力學(xué)計(jì)算效率［13］。

2.3仿真對(duì)象選用文獻(xiàn)［14］的縮比實(shí)驗(yàn)平臺(tái)進(jìn)行仿真研究，選取其中編號(hào)為15D1的工況作為仿真工況，在實(shí)驗(yàn)當(dāng)中，此工況表現(xiàn)出了典型的縱向不穩(wěn)定燃燒。縮比燃燒室結(jié)構(gòu)如圖1所示。燃燒室頭部為單噴嘴氣液噴霧燃燒，中心為氧化劑入口，氧化劑由比重為42%的氧氣與58%的水氣構(gòu)成，外側(cè)環(huán)縫為煤油的旋流入口。實(shí)驗(yàn)當(dāng)中的重要參數(shù)及結(jié)果如表1所示。在仿真時(shí)并未考慮液相過程，采用氣氣噴注的方式組織燃燒。氣氣噴注不同于氣液噴注，會(huì)導(dǎo)致仿真流場(chǎng)的溫度分布與實(shí)驗(yàn)情形差異明顯，這會(huì)影響到燃燒室內(nèi)釋熱與聲學(xué)的時(shí)空交互作用從而可能無法獲得實(shí)驗(yàn)當(dāng)中的不穩(wěn)定現(xiàn)象。然而，本文的目的并非復(fù)現(xiàn)實(shí)驗(yàn)現(xiàn)象，而是僅將此實(shí)驗(yàn)工況作為一種基礎(chǔ)工況，用氣氣燃燒取代氣液燃燒從而排除了液相的噴霧、蒸發(fā)等可能影響系統(tǒng)穩(wěn)定性的因素，加之k-ε湍流模型的使用，最終孤立了化學(xué)動(dòng)力學(xué)過程，使之成為燃燒釋熱的控制過程與決定系統(tǒng)穩(wěn)定性的最重要因素。圖2是本文所用的計(jì)算網(wǎng)格，網(wǎng)格數(shù)約為11萬，經(jīng)網(wǎng)格獨(dú)立性驗(yàn)證，此計(jì)算網(wǎng)格滿足仿真要求。監(jiān)測(cè)點(diǎn)選在燃燒劇烈的上游區(qū)域。入口邊界采用質(zhì)量流量邊界，出口邊界外推。初始時(shí)刻高溫氮?dú)獬涑庹麄�(gè)計(jì)算區(qū)域以達(dá)到點(diǎn)火的目的。

3無聲場(chǎng)激勵(lì)下的仿真結(jié)果及分析

仿真計(jì)算過程中，為了減小分離誤差的影響，仿真時(shí)間步長(zhǎng)選取為1.0×10-7s。圖3是監(jiān)測(cè)點(diǎn)的壓力變化曲線，仿真初始階段監(jiān)測(cè)點(diǎn)壓力曲線振蕩上升，在3ms附近達(dá)到峰值后開始振蕩下降，同時(shí)振蕩幅值逐漸衰減。在11ms后曲線趨于平穩(wěn)，最終的穩(wěn)態(tài)壓力約為2.37MPa，與實(shí)驗(yàn)的穩(wěn)態(tài)數(shù)據(jù)相差較�。ㄒ姳�1），但仿真工況表現(xiàn)出了線性穩(wěn)定性。系統(tǒng)表現(xiàn)出穩(wěn)定性的原因在于系統(tǒng)當(dāng)中的聲能激勵(lì)因素（熱聲耦合效應(yīng)）未能克服聲能阻滯因素（邊界及粘性對(duì)聲能的耗散），后者關(guān)系到燃燒室的聲學(xué)特性，本文的重點(diǎn)在于考察前者，即單純的化學(xué)動(dòng)力學(xué)過程控制的釋熱能否與聲波發(fā)生耦合互激從而成為燃燒不穩(wěn)定的潛在誘發(fā)機(jī)制。為了分析化學(xué)動(dòng)力學(xué)過程作為燃燒不穩(wěn)定的激勵(lì)因素的可能性，本文通過向燃燒室中引入聲場(chǎng)來考察氣相化學(xué)反應(yīng)與聲學(xué)振蕩的耦合機(jī)理。

4聲場(chǎng)的激勵(lì)方式

本文在穩(wěn)態(tài)流場(chǎng)的基礎(chǔ)上采用入口流量邊界引入脈動(dòng)量的方法在燃燒室內(nèi)激勵(lì)聲波［7］，這種方法比較容易在邊界設(shè)定擾動(dòng)條件，且誘導(dǎo)效果明顯，并對(duì)算法的穩(wěn)定性和收斂性影響較小。為激勵(lì)縱向聲波，在氧化劑與燃料入口加入同頻率、同相位的擾動(dòng)。

5釋熱系統(tǒng)在有限幅值擾動(dòng)激勵(lì)下的響應(yīng)

燃燒室在仿真工況下的一階縱向振型約為1300Hz，此聲學(xué)振型在實(shí)驗(yàn)當(dāng)中表現(xiàn)出了最不穩(wěn)定的性質(zhì)，以此振頻作為擾動(dòng)信號(hào)的頻率；流量擾動(dòng)法的最大幅值可達(dá)到穩(wěn)態(tài)流量的100%，以此作為擾動(dòng)信號(hào)的振幅，考察此時(shí)燃燒室的壓力振蕩與釋熱響應(yīng)。圖4是燃燒室的壓力與釋熱振蕩曲線。壓力曲線在第四周期發(fā)生“突躍”，在“突躍”后振蕩逐漸增大直至以較大幅值穩(wěn)定振蕩，在“突躍”后，壓力與釋熱每一周期的波動(dòng)均存在陡峭前沿，這種陡峭前沿波更接近于真實(shí)燃燒不穩(wěn)定發(fā)生時(shí)所觀測(cè)到的波形［1］。釋熱曲線呈現(xiàn)出了高頻小幅的振蕩，這種振蕩應(yīng)為基元反應(yīng)中的瞬態(tài)組分脈動(dòng)所引發(fā)。此外，由Rayleigh準(zhǔn)則可知：熱聲耦合發(fā)生正向激勵(lì)的必要條件是釋熱與壓力的振蕩保持相位一致性。由于多步化學(xué)動(dòng)力學(xué)系統(tǒng)在基元反應(yīng)中存在著鏈傳遞和壓力敏感反應(yīng)等機(jī)制，因而在聲場(chǎng)波動(dòng)影響下的釋熱波動(dòng)相較于聲場(chǎng)波動(dòng)發(fā)生了小幅相位偏移，但基本保持了相位一致，這種同相性意味著單純由氣相化學(xué)動(dòng)力學(xué)控制的釋熱過程與聲學(xué)過程存在著能夠發(fā)生耦合激勵(lì)的基礎(chǔ)。為進(jìn)一步明確化學(xué)動(dòng)力學(xué)過程是否增強(qiáng)了聲學(xué)振蕩，需要一個(gè)對(duì)比參照的系統(tǒng)，這個(gè)系統(tǒng)應(yīng)該提供聲波的原始振蕩強(qiáng)度，并且應(yīng)保持除化學(xué)動(dòng)力學(xué)過程以外的最大程度的近似，以較好地剝離化學(xué)動(dòng)力學(xué)因素的影響從而方便對(duì)其單獨(dú)研究。為此，本文建立了冷流聲學(xué)相似場(chǎng)的對(duì)比模型。對(duì)應(yīng)于反應(yīng)流流場(chǎng)，在相同質(zhì)量流量邊界的基礎(chǔ)上，考慮化學(xué)動(dòng)力學(xué)的組分源項(xiàng)但不加入化學(xué)動(dòng)力學(xué)能量源項(xiàng)從而建立冷流流場(chǎng)，兩者最大的區(qū)別在于是否存在釋熱。為了使冷流的流場(chǎng)平均聲速和平均壓力與反應(yīng)流相當(dāng)，提升入口邊界的初始溫度。如此一來，兩者的固有振型一致、質(zhì)量流量和平均室壓相同，唯一的區(qū)別在于反應(yīng)流通過化學(xué)動(dòng)力學(xué)釋熱獲得溫度而冷流相似場(chǎng)在入口處直接給定。這樣，就可以量化釋熱過程對(duì)聲學(xué)過程的激勵(lì)效應(yīng)的大小。圖5對(duì)比了相同振幅擾動(dòng)下的反應(yīng)流壓力響應(yīng)和其冷流聲學(xué)相似場(chǎng)的壓力響應(yīng)。冷流聲學(xué)相似場(chǎng)壓力振蕩的最終幅值與反應(yīng)流突躍前的振蕩幅值基本一致，因此可以斷定化學(xué)動(dòng)力學(xué)的釋熱過程確實(shí)增強(qiáng)了壓力振蕩。比較圖4中的壓力曲線和圖5（a）可以發(fā)現(xiàn)，在反應(yīng)流的振蕩達(dá)到穩(wěn)定后，較其冷流聲學(xué)相似場(chǎng)的振幅增強(qiáng)約300%。多步化學(xué)動(dòng)力學(xué)過程與聲學(xué)過程耦合發(fā)生了相互激勵(lì)，而這種相互激勵(lì)的產(chǎn)生可能與“突躍”和“陡峭前沿”有著密切的聯(lián)系，也有可能與瞬態(tài)組分脈動(dòng)所引發(fā)的高頻釋熱脈動(dòng)有關(guān)。

6釋熱系統(tǒng)在小振幅擾動(dòng)激勵(lì)下的響應(yīng)

取擾動(dòng)幅值為0.5%，此小振幅下的聲學(xué)過程基本是線性的［15］。圖6對(duì)比了線性聲場(chǎng)中的反應(yīng)流場(chǎng)與其冷流聲學(xué)相似場(chǎng)的壓力波動(dòng)及頻譜。反應(yīng)流與其冷流聲學(xué)相似場(chǎng)對(duì)線性擾動(dòng)的壓力波動(dòng)響應(yīng)基本一致，化學(xué)動(dòng)力學(xué)釋熱過程未能與聲學(xué)過程發(fā)生明顯的耦合激勵(lì)。在線性擾動(dòng)下，燃燒室的壓力波形更加趨近于標(biāo)準(zhǔn)的正弦波形，不再具有“突躍”及“陡峭前沿峰”等非線性現(xiàn)象。圖7對(duì)比了線性與非線性擾動(dòng)下的系統(tǒng)釋熱響應(yīng)。線性擾動(dòng)下與非線性擾動(dòng)下的釋熱曲線均存在著高頻小幅振蕩，然而線性擾動(dòng)下并未發(fā)生明顯的熱聲耦合激勵(lì)，因此，高頻小幅的釋熱振蕩并不是聲波放大的誘因。由于這種高頻小幅的釋熱振蕩能量太低且頻率太高導(dǎo)致了壓力波無法對(duì)其進(jìn)行充分的響應(yīng)，因而與高頻脈動(dòng)的釋熱曲線相對(duì)的壓力曲線十分光滑。

7分析與討論

在有限幅值擾動(dòng)激勵(lì)的非線性聲場(chǎng)中，壓力波動(dòng)“突躍”后發(fā)展成為具有陡峭前沿的壓力振蕩，化學(xué)反應(yīng)釋熱與聲波發(fā)生了耦合激勵(lì)；而在小幅值擾動(dòng)激勵(lì)的線性聲場(chǎng)中，化學(xué)反應(yīng)釋熱與聲波并未發(fā)生明顯互激。在兩種情形下，兩者的釋熱過程均是以小幅高頻脈動(dòng)進(jìn)行的，前者存在著后者所沒有的非線性現(xiàn)象，前者發(fā)生熱聲激勵(lì)耦合而后者沒有，因此可以基本認(rèn)為：化學(xué)動(dòng)力學(xué)釋熱與聲學(xué)過程的相互激勵(lì)作用來自于“突躍”和“陡峭前沿”，瞬態(tài)組分的高頻脈動(dòng)對(duì)聲學(xué)過程并無影響。在非線性聲學(xué)振蕩下的系統(tǒng)的釋熱響應(yīng)之所以會(huì)發(fā)生“突躍”，可能與氣相化學(xué)動(dòng)力學(xué)控制的釋熱系統(tǒng)自身蘊(yùn)含的分岔特性的相關(guān)。化學(xué)動(dòng)力學(xué)系統(tǒng)存在著“點(diǎn)火延遲”、“可燃極限”等概念，這些概念當(dāng)中蘊(yùn)含著發(fā)生釋熱分岔的可能性。在化學(xué)動(dòng)力學(xué)釋熱過程中，存在起鏈、鏈傳遞等子過程，因此會(huì)出現(xiàn)以溫度迅速躍升為標(biāo)志的點(diǎn)火點(diǎn)。圖8給出了在混合比為1.0，初溫1200K條件下，環(huán)境壓力分別為2MPa和4MPa時(shí)，簡(jiǎn)化機(jī)理在0-D勻質(zhì)反應(yīng)器中的溫度對(duì)時(shí)間的變化曲線。由圖可知，在不同壓力下，溫度曲線會(huì)出現(xiàn)平移，壓力越高，曲線向左平移幅度越大。以拉格朗日法的觀點(diǎn)來看，圖中的時(shí)間軸近似代表了燃燒室中的空間位置。如此一來，在豎直方向劃一直線交不同環(huán)境壓力下的曲線于兩個(gè)交點(diǎn)，則這兩點(diǎn)代表了同一空間位置在壓力波動(dòng)下的溫度響應(yīng)，而這種溫度響應(yīng)直接反映了釋熱響應(yīng)。若兩個(gè)交點(diǎn)分別處于各自曲線的未燃區(qū)和已燃區(qū)，則壓力的波動(dòng)就造成了燃燒室中某處的燃燒行為發(fā)生了分岔。圖中兩曲線所圍成的近似矩形的區(qū)域?yàn)榇藟毫φ袷幭碌尼専犴憫?yīng)敏感區(qū)域，壓力越高曲線向右平移越大，也就意味著敏感區(qū)域會(huì)隨著壓力波動(dòng)的增大而擴(kuò)大；同時(shí)，在同一空間區(qū)域的釋熱波動(dòng)也會(huì)增大。在非線性聲場(chǎng)下出現(xiàn)了的“突躍”，可能就源于上游某區(qū)域由于壓力的大幅振蕩所導(dǎo)致的局部釋熱分岔，這種分岔反過來激勵(lì)了壓力波的振蕩，從而形成正反饋，出現(xiàn)了熱聲耦合激勵(lì)的現(xiàn)象。之所以線性聲場(chǎng)下熱聲效應(yīng)不明顯，其原因可能在于小波動(dòng)線性條件下壓力接近的兩條曲線的敏感區(qū)域基本不存在從而無法發(fā)生劇烈的釋熱分岔。

8結(jié)論

（1）氣相化學(xué)動(dòng)力學(xué)控制的釋熱系統(tǒng)與線性小振幅聲波無法發(fā)生明顯的耦合激勵(lì)。（2）聲學(xué)過程對(duì)化學(xué)動(dòng)力學(xué)的高頻小幅釋熱脈動(dòng)無響應(yīng)。（3）氣相化學(xué)動(dòng)力學(xué)控制的釋熱系統(tǒng)與非線性有限振幅聲波可以發(fā)生明顯的耦合激勵(lì)，發(fā)生耦合激勵(lì)后的壓力振蕩較冷流聲學(xué)相似場(chǎng)振幅增大了約300%。這一激勵(lì)產(chǎn)生的原因可能是化學(xué)動(dòng)力學(xué)釋熱系統(tǒng)的釋熱分岔效應(yīng)所致。（4）氣相化學(xué)動(dòng)力學(xué)釋熱過程無法成為線性不穩(wěn)定的控制因素，然而在燃燒室發(fā)生非線性振蕩時(shí)（自激或外激），它可能會(huì)促進(jìn)熱聲耦合的激勵(lì)效應(yīng)。

作者：薛誠(chéng)尤 聶萬勝 何博 鄭剛 單位：中國(guó)人民解放軍裝備學(xué)院 研究生院 中國(guó)人民解放軍裝備學(xué)院 航天裝備系

  本文關(guān)鍵詞：推進(jìn)技術(shù)雜志，由筆耕文化傳播整理發(fā)布。