• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    高超聲速飛行器電離反應(yīng)流場(chǎng)特性研究

    2020-07-06 08:25:26衛(wèi)麗云王學(xué)德
    彈道學(xué)報(bào) 2020年2期
    關(guān)鍵詞:物面電離飛行器

    衛(wèi)麗云,王學(xué)德

    (南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 南京 210094)

    隨著人類對(duì)高空稀薄流領(lǐng)域探索的不斷深入,高超聲速飛行器的性能和用途也得到了極大的改善,尤其在軍事方面,高超聲速飛行器得到了全面發(fā)展,比如,高超聲速飛行器能有效地進(jìn)行高空預(yù)警,偵查和高空打擊工作,極大地拓寬了作戰(zhàn)空間。為此,高超聲速飛行器水平已逐漸成為各國(guó)軍事實(shí)力的表現(xiàn)。但隨著飛行條件的不斷升高,稀薄氣體電離成為高超聲速飛行器不得不考慮的問題,由其導(dǎo)致的“黑障現(xiàn)象”[1]很大程度上影響了再入飛行器與地面的通信交流,且電離反應(yīng)的吸放熱也會(huì)影響流場(chǎng)中飛行器的氣動(dòng)特性。而對(duì)于電離反應(yīng)流場(chǎng)的研究,必須在非電離反應(yīng)流場(chǎng)的基礎(chǔ)上引入新的組分,即非電離流場(chǎng)五組分所對(duì)應(yīng)的離子和游離電子。離子和電子的加入使流場(chǎng)中的能量分配產(chǎn)生新的變化,且因?yàn)榉磻?yīng)機(jī)理的不同,電離反應(yīng)的實(shí)現(xiàn)方法也一直是研究的熱點(diǎn)。

    直接模擬蒙特卡羅方法(DSMC)是一種基于統(tǒng)計(jì)學(xué)的數(shù)值模擬方法,該方法自提出以來被廣泛應(yīng)用于稀薄流問題的研究[2-3]。將電離反應(yīng)引入DSMC一直以來被視為最可能的研究電離反應(yīng)的途徑。早期的研究大多將電離反應(yīng)幾率與Arrhenius 速率系數(shù)函數(shù)相關(guān)聯(lián),而Arrhenius速率函數(shù)依賴于單溫度速率常數(shù)的宏觀信息,且建立在化學(xué)反應(yīng)平衡的假設(shè)上,這一思想在文獻(xiàn)[4-5]中都有應(yīng)用。但是電離反應(yīng)溫度可達(dá)上萬開爾文,目前并沒有如此高溫下完備的單溫度速率數(shù)據(jù),因此,在這一溫度范圍內(nèi),宏觀單溫度速率常數(shù)并不適用。盡管文獻(xiàn)[6]在此基礎(chǔ)上引入振動(dòng)能的影響,使這種方法更接近真實(shí)效應(yīng),但對(duì)電離反應(yīng)發(fā)生的幾率的處理依舊依賴于氣體的宏觀性質(zhì)。BIRD[7]提出了一種基于量子動(dòng)力學(xué)的模型—Q-K模型,該方法擺脫了對(duì)宏觀速率方程和氣體平衡假設(shè)的依賴,將所有的反應(yīng)與微觀粒子的振動(dòng)能相關(guān)聯(lián)。但由于電離反應(yīng)的特殊性,主導(dǎo)電離反應(yīng)發(fā)生的往往不是粒子的振動(dòng)能的激發(fā),而是相對(duì)應(yīng)的電子能的激發(fā),因此,僅考慮粒子振動(dòng)能的Q-K模型并不完全適用于電離反應(yīng)。LIECHTY等[8]基于粒子微觀動(dòng)力學(xué)理論,將電子能級(jí)引入電離反應(yīng)模型,認(rèn)為電子能在流場(chǎng)中的計(jì)算應(yīng)該包括3個(gè)部分:電子能作為電離反應(yīng)碰撞能的一部分,必須從能量分布中取值;電子能在粒子碰撞的過程中也需和其他能量一樣,有合理的分配過程;電子能的產(chǎn)生和消耗必須有合理的轉(zhuǎn)換速率。

    由于研究采用模型的局限性或者研究的側(cè)重點(diǎn)不同,學(xué)者們?cè)诔绦蛑锌紤]的電離反應(yīng)也不盡相同。LIECHTY等[9]應(yīng)用分子動(dòng)力學(xué)原理將聯(lián)合電離反應(yīng)和置換電離反應(yīng)引入DSMC程序中,重點(diǎn)研究了將這一方法用于計(jì)算極高溫度下的速率常數(shù)的可能性。SHEVYRIN等[10]根據(jù)流場(chǎng)電中性假設(shè),不將電子作為粒子進(jìn)行模擬,直接把網(wǎng)格離子數(shù)密度之和作為電子數(shù)密度。這一方法雖然簡(jiǎn)化了程序,但直接忽略了電子與原子碰撞發(fā)生的電子激發(fā)電離反應(yīng)。FANG等[11]采用不同的權(quán)重因子將中性粒子與帶電粒子進(jìn)行區(qū)分,重點(diǎn)考慮了含有電子參與的電離反應(yīng),忽略了包含電荷轉(zhuǎn)移的置換電離反應(yīng)。

    本文從分子動(dòng)力學(xué)理論出發(fā),將電離反應(yīng)分類進(jìn)一步細(xì)化,考慮更多的電離反應(yīng),并根據(jù)電離反應(yīng)發(fā)生機(jī)理的不同采用不同的反應(yīng)發(fā)生判據(jù)。首先通過對(duì)RAMC-Ⅱ外形的飛行器的飛行模擬,驗(yàn)證了本文所采用二維軸對(duì)稱電離程序的有效性。隨后又應(yīng)用這一模型,研究了電離反應(yīng)流場(chǎng)與非電離反應(yīng)流場(chǎng)的不同以及不同種類電離反應(yīng)對(duì)稀薄流場(chǎng)飛行器特性的影響程度。

    1 電離反應(yīng)及實(shí)現(xiàn)方法

    1.1 電子激發(fā)能的計(jì)算

    在電離反應(yīng)中,伴隨著電子的產(chǎn)生或轉(zhuǎn)移,粒子的電子態(tài)受到相應(yīng)的激發(fā),會(huì)產(chǎn)生一部分能量,定義為電子激發(fā)能。電子激發(fā)態(tài)從相應(yīng)的碰撞實(shí)效溫度(由碰撞中分子的相對(duì)平動(dòng)能和電子能之和求得)的平衡態(tài)取樣得到,而一個(gè)粒子的電子激發(fā)能由下式計(jì)算得到:

    Ee=∑(Nj/N)ej

    (1)

    式中:ej為電子能級(jí)j的能量,Nj為能級(jí)j中的粒子數(shù),N為總粒子數(shù),Nj/N可以根據(jù)平衡分布得出,如下式所示:

    (2)

    式中:k為玻爾茲曼常數(shù);T為碰撞實(shí)效溫度,由碰撞對(duì)粒子的碰撞數(shù)量平衡分布取樣而得;gj為能級(jí)的簡(jiǎn)并度,具體的gj和ej的值可參考文獻(xiàn)[4]。

    1.2 電離反應(yīng)的DSMC實(shí)現(xiàn)

    根據(jù)電離反應(yīng)發(fā)生機(jī)理的不同,電離反應(yīng)大概可分為3類:聯(lián)合電離反應(yīng)、電子激發(fā)電離反應(yīng)、置換電離反應(yīng)。其中,置換電離反應(yīng)因?yàn)榘l(fā)生置換的粒子的種類不同,起主導(dǎo)作用的能量形式會(huì)有相應(yīng)的變化,可進(jìn)一步分為電子交換電離反應(yīng)和含離子的置換反應(yīng)兩類。

    聯(lián)合電離反應(yīng)的過程可以看做兩步反應(yīng),如式(3)所示,其中e代表電子。第一步反應(yīng)和復(fù)合反應(yīng)相似,生成的AB*可以看作絡(luò)合物,極不穩(wěn)定。第二步反應(yīng)中,反應(yīng)初始的碰撞能應(yīng)該包括離解能,再根據(jù)L-B振動(dòng)松弛理論進(jìn)行能量分配后,若此時(shí)AB*的平動(dòng)能和電子激發(fā)能之和大于其電離活化能,則聯(lián)合電離反應(yīng)發(fā)生。本文所考慮的聯(lián)合電離反應(yīng)如式(4)~式(6)所示。

    A+B→AB*→AB++e

    (3)

    (4)

    (5)

    N+O→NO++e

    (6)

    電子激發(fā)電離是單原子分子和電子碰撞生成離子的反應(yīng),其一般反應(yīng)式可以寫為式(7)的形式。其中電子沒有電子激發(fā)能,碰撞對(duì)的碰撞總能記為碰撞對(duì)的相對(duì)平動(dòng)能和單原子分子的電子激發(fā)能之和。該類電離反應(yīng)是否發(fā)生主要取決于碰撞能與電離活化能的相對(duì)大小,當(dāng)該碰撞對(duì)的碰撞能大于單原子分子的電離活化能時(shí),電子激發(fā)電離反應(yīng)發(fā)生。本文所考慮的電子激發(fā)電離反應(yīng)如式(8)~式(9)所示。

    A+e→A++e+e

    (7)

    N+e→N++e+e

    (8)

    O+e→O++e+e

    (9)

    置換電離反應(yīng)是由一個(gè)離子和一個(gè)中性粒子碰撞而發(fā)生的反應(yīng)。根據(jù)反應(yīng)機(jī)理又可細(xì)分為2類,這里稱兩粒子間碰撞得失電子的反應(yīng)為電子交換電離反應(yīng),而兩粒子間碰撞得失原子的反應(yīng)為含離子的置換反應(yīng)。

    電子交換電離反應(yīng)的過程可以看作中性粒子受電子態(tài)激發(fā)失去電子,之后由于電荷間的相互作用,離子得到電子,其一般式可寫為式(10)。對(duì)于這類化學(xué)反應(yīng)中的吸熱反應(yīng),采用L-B振動(dòng)松弛理論進(jìn)行能量分配后,若此時(shí)中性粒子A的電子能級(jí)仍大于反應(yīng)的活化電子能級(jí),則反應(yīng)發(fā)生。而對(duì)于這類反應(yīng)中的放熱反應(yīng),判據(jù)與吸熱反應(yīng)相同,只是在進(jìn)行能量分配時(shí),除了本身的碰撞總能外,還需考慮反應(yīng)能。本文所考慮的電子交換電離反應(yīng)如式(11)~式(14)所示。

    A+B+→A++B

    (10)

    (11)

    (12)

    (13)

    N+NO+N++NO

    (14)

    對(duì)于含離子的置換反應(yīng),其一般式可寫為式(15)。這類化學(xué)反應(yīng)與中性粒子的置換反應(yīng)機(jī)理類似,當(dāng)碰撞對(duì)的碰撞總能達(dá)到一定值,雙原子粒子的化學(xué)鍵發(fā)生斷裂,其中一個(gè)原子會(huì)與碰撞對(duì)中另一反應(yīng)原子結(jié)合,發(fā)生置換電離反應(yīng)。所以,對(duì)于這一類化學(xué)反應(yīng),主要的影響因素還是雙原子粒子的振動(dòng)能級(jí)。對(duì)于吸熱反應(yīng),對(duì)碰撞能采用L-B振動(dòng)松弛理論進(jìn)行能量分配,若此時(shí)雙原子粒子的振動(dòng)能級(jí)仍高于反應(yīng)的活化振動(dòng)能級(jí),則反應(yīng)發(fā)生。而對(duì)于放熱反應(yīng),在進(jìn)行能量分配時(shí)需考慮碰撞能和反應(yīng)能的總和,之后若反應(yīng)生成的中性粒子的振動(dòng)能級(jí)高于活化能級(jí),則反應(yīng)發(fā)生。本文所考慮的含離子的置換反應(yīng)如式(16)~式(19)所示。

    A+BC+→AB(AB+)+C+(C)

    (15)

    N+NO+N2+O+

    (16)

    (17)

    (18)

    NO+O+O2+N+

    (19)

    2 計(jì)算方法

    2.1 DSMC簡(jiǎn)介

    稀薄流領(lǐng)域的非連續(xù)性導(dǎo)致N-S方程失效,基于N-S方程的數(shù)值模擬方法也就不再適用于該領(lǐng)域問題的研究,BIRD提出的基于Boltzmann方程的DSMC方法成為目前研究稀薄流領(lǐng)域問題最常用的方法。該方法通過對(duì)Boltzmann方程所描述物理過程的模擬,將粒子的運(yùn)動(dòng)和碰撞解耦,解決了Boltzmann方程碰撞項(xiàng)難解的問題,用一個(gè)模擬粒子代表一定數(shù)量的真實(shí)粒子,各模擬粒子之間以及模擬粒子與物面之間不斷進(jìn)行著能量交換,最后通過采樣統(tǒng)計(jì)得到計(jì)算域內(nèi)計(jì)算網(wǎng)格的宏觀特性分布。

    2.2 軸對(duì)稱流場(chǎng)的DSMC實(shí)現(xiàn)方法

    二維的DSMC程序由于忽略了z方向的粒子運(yùn)動(dòng),計(jì)算結(jié)果會(huì)與實(shí)際情況有較大的偏差,但是三維的DSMC程序又會(huì)極大地增加計(jì)算難度,降低計(jì)算效率。因此,本文引入二維軸對(duì)稱程序[12],將z方向的運(yùn)動(dòng)并入y方向,提高二維程序計(jì)算精度的同時(shí)沒有增加很大的計(jì)算量。

    將z方向的運(yùn)動(dòng)疊加到y(tǒng)方向,實(shí)際上就是將Oyz面內(nèi)的運(yùn)動(dòng)作極坐標(biāo)轉(zhuǎn)化,程序中的粒子位置的y坐標(biāo)取Oyz面的極坐標(biāo)值,相應(yīng)的Δt時(shí)間內(nèi)y坐標(biāo)的變動(dòng)也需考慮y和z2個(gè)方向上的運(yùn)動(dòng)。

    此外,作軸對(duì)稱處理后沿網(wǎng)格體積徑向方向變化較大,而DSMC程序中的碰撞取樣是按網(wǎng)格體積均分取樣的,這就會(huì)導(dǎo)致網(wǎng)格取樣不均的問題,當(dāng)靠近軸線的網(wǎng)格達(dá)到計(jì)算要求的取樣數(shù)時(shí),外圍的網(wǎng)格取樣數(shù)會(huì)過多,極大地提高了計(jì)算難度。在此引入半徑權(quán)重系數(shù)WF,設(shè)定一個(gè)參考半徑RF,RF的取值要比網(wǎng)格平均尺寸小得多。WF的定義為

    WF=r/RF

    (20)

    式中:r為粒子所處半徑,即粒子的y坐標(biāo);WF的值最小取1。網(wǎng)格內(nèi)的每個(gè)模擬粒子所代表的真實(shí)粒子的數(shù)量由原來的F變?yōu)閃FF,這樣大大減小了外圍網(wǎng)格與靠近軸線網(wǎng)格之間粒子分布數(shù)量的差值,取樣數(shù)量也得到平均。

    2.3 電子運(yùn)動(dòng)在DSMC中的實(shí)現(xiàn)方法

    電子運(yùn)動(dòng)的處理也是電離反應(yīng)的一大難點(diǎn)。電子的質(zhì)量為9.11×10-31kg,相較于流場(chǎng)中的其他重離子,電子要輕大約5個(gè)數(shù)量級(jí),而粒子的熱運(yùn)動(dòng)速度與粒子質(zhì)量的平方根成反比[13],這就表示電子的運(yùn)動(dòng)速度要比其他重粒子快大約2個(gè)數(shù)量級(jí)。在適合捕捉其他粒子運(yùn)動(dòng)的時(shí)間步長(zhǎng)內(nèi),無法準(zhǔn)確跟蹤電子運(yùn)動(dòng),但若以電子的運(yùn)動(dòng)速度取時(shí)間步長(zhǎng),又極大地提高了DSMC計(jì)算的難度。

    目前處理電子的運(yùn)動(dòng)主要有3種思路:①人為增大電子質(zhì)量3個(gè)數(shù)量級(jí)[4];②捆綁法[14];③平均速度法[15]。這3種方法都可以滿足流場(chǎng)基本電中性的假設(shè),且在國(guó)內(nèi)外均有應(yīng)用。

    本文采用以上的第一種方法,人為增大電子質(zhì)量,將電子質(zhì)量認(rèn)為是9.11×10-28kg。為了保證碰撞過程的質(zhì)量守恒,必須同時(shí)改變離子的質(zhì)量。因?yàn)榧词闺娮拥馁|(zhì)量增大3個(gè)數(shù)量級(jí),還是比離子的質(zhì)量小大約2個(gè)數(shù)量級(jí),占離子質(zhì)量的比重很小,若將離子質(zhì)量減少相應(yīng)比重,對(duì)其熱運(yùn)動(dòng)速度的影響很小,可忽略不計(jì)。表1為離子質(zhì)量改變的詳細(xì)信息。

    3 驗(yàn)證算例

    3.1 程序有效性驗(yàn)證

    本文的驗(yàn)證模型取RAMC-Ⅱ試驗(yàn)中采用的飛行器外形,具體的尺寸如圖1所示。沿機(jī)身方向分別布置了4個(gè)微型反射器,微型反射器用于測(cè)量所在位置垂直物面方向的最大電子密度,具體的位置也在圖1中標(biāo)出。計(jì)算條件取自文獻(xiàn)[10],如表2所示。

    圖1 RAMC-Ⅱ試驗(yàn)中飛行器外形及探測(cè)點(diǎn)分布

    表2 計(jì)算條件

    本文的網(wǎng)格采用非結(jié)構(gòu)貼體網(wǎng)格,分子碰撞采用VHS(變徑硬球)模型,碰撞能量分配采用Larsen-Borgnakke模型,分子與物面相互作用采用完全漫反射模型,物面溫度設(shè)置為恒溫1 000 K。由于飛行器的攻角為0°,只計(jì)算一半流場(chǎng)。

    圖2為不同高度(H)下得到的飛行器周圍密度云圖。飛行器頭部有很明顯的弓形激波,激波后的氣體密度急劇增大,且隨著飛行高度的增加,流場(chǎng)的密度逐漸減小。

    圖2 不同高度流場(chǎng)密度云圖

    圖3為不同高度下飛行器周圍溫度云圖。激波后的氣體溫度迅速升高,但沿貼近壁面方向溫度逐漸降低。隨著飛行高度的升高,激波的厚度有所增加。以上密度和溫度云圖的變化趨勢(shì)與文獻(xiàn)[16]一致(二者的計(jì)算條件相差不大),從定性角度分析,所采用的電離程序是符合變化規(guī)律的。

    圖3 不同高度流場(chǎng)溫度云圖

    圖4為不同探測(cè)點(diǎn)(已在圖1中標(biāo)出)垂直于物面方向上的最大電子數(shù)密度隨飛行高度變化的對(duì)比。圖中,H為飛行器運(yùn)行的海拔高度,ne為垂直物面方向上的最大電子數(shù)密度,對(duì)ne取對(duì)數(shù),即lgne。由于DSMC中統(tǒng)計(jì)是以網(wǎng)格單元為單位,而在粒子運(yùn)動(dòng)過程中,各網(wǎng)格中的粒子數(shù)會(huì)不斷發(fā)生變化,取樣數(shù)也會(huì)發(fā)生相應(yīng)的變化,造成統(tǒng)計(jì)過程中存在一定的漲落誤差,一般認(rèn)為,電子密度計(jì)算誤差允許在3倍以內(nèi)[17]。從圖中可以看出,無論哪一探測(cè)點(diǎn),本文所采用的程序得到的計(jì)算結(jié)果基本與文獻(xiàn)[10]的結(jié)果以及實(shí)驗(yàn)數(shù)據(jù)處于同一數(shù)量級(jí),就相差倍數(shù)來說,與文獻(xiàn)[10]之間的差距更小一些,具有良好的一致性。這就從定量的角度說明本文所采用電離程序的正確性。

    圖4 不同探測(cè)點(diǎn)的最大電子數(shù)密度

    3.2 網(wǎng)格無關(guān)性驗(yàn)證

    在DSMC中,網(wǎng)格主要起2種作用:①對(duì)宏觀流場(chǎng)性質(zhì)進(jìn)行空間離散,便于統(tǒng)計(jì)流場(chǎng)信息;②便于分子碰撞對(duì)的選擇。網(wǎng)格尺寸過大,會(huì)造成流場(chǎng)宏觀特性的估計(jì)誤差增大;網(wǎng)格尺寸過小又會(huì)極大地增加程序的計(jì)算量,降低工作效率。因此,合理的網(wǎng)格設(shè)計(jì)對(duì)DSMC計(jì)算結(jié)果的可信度有重要影響。對(duì)于網(wǎng)格尺寸的選擇,學(xué)者們采用的標(biāo)準(zhǔn)不一,但原則上最大不超過來流分子的平均自由程。本文引入一個(gè)無量綱參數(shù)Z:

    (21)

    式中:s為網(wǎng)格單元的平均尺寸,λ為來流分子平均自由程,Kn為來流的克努森數(shù),L為流場(chǎng)特征長(zhǎng)度。

    以90 km飛行高度條件為例,采用非結(jié)構(gòu)貼體網(wǎng)格,取Z=0.33,Z=0.50,Z=0.80,Z=1.00分別計(jì)算了飛行器附近流場(chǎng)的特性。Z值越大,網(wǎng)格單元平均尺寸與分子自由程越相近,網(wǎng)格分布也就越稀疏。圖5為不同網(wǎng)格密度下物面附近流場(chǎng)特性對(duì)比。圖中,Cp為物面壓力系數(shù),CF為物面阻力系數(shù),q為物面熱流密度。為了方便比較,將橫坐標(biāo)x作無量綱化處理,此處的流場(chǎng)特征長(zhǎng)度L取飛行器的模擬長(zhǎng)度。

    從圖5中可以看出,4種不同尺寸的網(wǎng)格計(jì)算得到的結(jié)果吻合較高,特別是從物面壓力系數(shù)分布和熱流密度分布來看,這4種尺度的網(wǎng)格得到的結(jié)果幾乎一致。但從摩阻系數(shù)分布來看,不同網(wǎng)格尺寸得到的摩阻系數(shù)峰值稍有不同,其中Z=0.50和Z=0.80兩種網(wǎng)格得到的計(jì)算結(jié)果始終保持較高的吻合度,同時(shí)為了減小計(jì)算量,本文以下的計(jì)算中,統(tǒng)一采用Z=0.80的網(wǎng)格,即網(wǎng)格的平均尺寸取為來流分子平均自由程的4/5。

    4 電離反應(yīng)對(duì)流場(chǎng)的影響

    4.1 電離流場(chǎng)特性

    為了方便看出電離反應(yīng)的影響,將5組分化學(xué)反應(yīng)程序與11組分含電離的化學(xué)反應(yīng)程序得到的計(jì)算結(jié)果進(jìn)行對(duì)比,飛行條件取80 km海拔處的來流條件,如表3所示。

    表3 來流條件

    圖6為80 km飛行高度下不同組分流場(chǎng)的溫度對(duì)比。

    圖6 不同組分下流場(chǎng)溫度云圖

    可以看出,11組分流場(chǎng)的整體溫度要比5組分低,其中5組分溫度場(chǎng)的最大值為16 989 K,而11組分為14 380 K,相差2 500 K左右。從激波位置來看,11組分流場(chǎng)中的激波更靠近壁面,更薄一些。這是因?yàn)樵诒疚牡难芯恐?11組分與5組分相比,多考慮了21個(gè)電離反應(yīng),而電離反應(yīng)的活化能與非電離反應(yīng)相比較高,要使化學(xué)反應(yīng)充分發(fā)生且達(dá)到平衡所需消耗的熱量也就較大,也就會(huì)使11組分流場(chǎng)整體溫度有所下降?;诹鲌?chǎng)中氣體的熱脹冷縮效應(yīng),在相同的壁面溫度下,溫度較低的流場(chǎng)的壓縮性更強(qiáng)一些,以致11組分流場(chǎng)中的氣體在壁面附近受到更大的壓縮,形成的激波厚度也就相對(duì)較薄。

    圖7為80 km飛行高度下不同組分流場(chǎng)中的物面氣動(dòng)特性對(duì)比,x/R為流場(chǎng)橫坐標(biāo)與飛行器頭部半徑的比值。由圖可知,物面壓力系數(shù)的變化趨勢(shì)不隨流場(chǎng)組分的變化而變化,都是在駐點(diǎn)區(qū)域達(dá)到最大值,后沿x軸方向急劇下降,直到飛行器肩部位置后趨于平穩(wěn)。而物面摩阻系數(shù)在駐點(diǎn)附近最小,之后急劇升高,到x/R約等于0.2附近達(dá)到峰值,之后的變化趨勢(shì)同壓力系數(shù)相似,快速下降到肩部附近后趨于平穩(wěn),且不同組分流場(chǎng)相差不大。由此可以得出,電離反應(yīng)對(duì)飛行器表面的氣動(dòng)特性影響不大。壁面壓力系數(shù)和摩阻系數(shù)分別取決于物面入射分子與反射分子的法向動(dòng)量差和切向動(dòng)量差,在相同來流條件下,雖然中性粒子在物面附近有可能電離產(chǎn)生對(duì)應(yīng)的離子和電子,使飛行器周圍的粒子數(shù)密度有所增大,但電子質(zhì)量較小,即使本文將電子質(zhì)量人為增大1 000倍,仍與重粒子相差2個(gè)數(shù)量級(jí),因此電子對(duì)物面附近入射分子與反射分子的動(dòng)量影響不大,因而會(huì)出現(xiàn)圖7所示的情況。

    圖7 80 km飛行高度下不同組分流場(chǎng)物面氣動(dòng)特性變化

    圖8為80 km飛行高度下不同組分流場(chǎng)中物面熱流密度對(duì)比。與氣動(dòng)特性曲線一樣,不同組分流場(chǎng)的沿物面熱流密度變化趨勢(shì)一致,但5組分流場(chǎng)中的熱流密度峰值明顯高于11組分流場(chǎng)中的熱流密度峰值,這是因?yàn)殡婋x反應(yīng)大多是吸熱反應(yīng),會(huì)吸收大量熱量,從而使飛行器表面與周圍氣體的熱量交換減少。而這一差別主要體現(xiàn)在駐點(diǎn)附近,這是因?yàn)殡婋x反應(yīng)大多集中在駐點(diǎn)附近。電離反應(yīng)的發(fā)生主要是因?yàn)楦邷叵码娮幽軕B(tài)的激發(fā),從圖6可以看出,駐點(diǎn)附近溫度是流場(chǎng)中的高溫集中區(qū)域,因而發(fā)生的電離反應(yīng)數(shù)量也就較多。

    圖8 80 km飛行高度下不同組分流場(chǎng)物面熱流密度變化

    4.2 不同電離反應(yīng)對(duì)流場(chǎng)的影響

    為了研究本文所考慮的3種電離反應(yīng)的主次關(guān)系,依次將3種電離反應(yīng)從程序中抹去,計(jì)算了這3種情況下的流場(chǎng)特性。飛行高度取為80 km,其他計(jì)算條件與表3相同。

    本文來流速度為7 650 m/s,對(duì)應(yīng)所含的能量大約為9 eV,根據(jù)文獻(xiàn)[13],此時(shí)空氣中產(chǎn)生的主導(dǎo)電離反應(yīng)主要為聯(lián)合電離反應(yīng)。從本文所考慮的所有電離反應(yīng)式來看,其他有帶電粒子參與的電離反應(yīng)的發(fā)生都是以聯(lián)合電離反應(yīng)產(chǎn)生為前提的。因而如果不考慮聯(lián)合電離反應(yīng),則其他電離反應(yīng)也不發(fā)生,計(jì)算得到的結(jié)果同5組分流場(chǎng)的計(jì)算結(jié)果相近,在以下的討論中不作贅述。為方便敘述,將考慮所有電離反應(yīng)的計(jì)算稱為工況1,不計(jì)電子激發(fā)電離反應(yīng)的計(jì)算稱為工況2,不計(jì)置換電離反應(yīng)的計(jì)算稱為工況3。由前面的討論可知,電離反應(yīng)對(duì)飛行器表面氣動(dòng)特性的影響不是很大,所以后續(xù)討論的重點(diǎn)為不同電離反應(yīng)對(duì)電子密度及飛行器表面氣動(dòng)熱的影響。

    圖9為3種不同工況下計(jì)算所得的飛行器頭部電子密度的對(duì)比圖。由圖9可知,工況1和工況3的電子密度分布很相近,工況2的電子層相對(duì)較薄,密度也就小,不過3種工況的電子密度基本還在一個(gè)數(shù)量級(jí)上。從化學(xué)平衡的角度出發(fā),若將電子激發(fā)電離反應(yīng)從程序中抹去,即只考慮置換電離反應(yīng)和聯(lián)合電離反應(yīng),根據(jù)比較結(jié)果可以看出,置換電離反應(yīng)對(duì)生成游離電子的聯(lián)合電離反應(yīng)有抑制作用,從而使流場(chǎng)中的電子密度有所降低,這是因?yàn)殡m然置換電離反應(yīng)會(huì)消耗一部分聯(lián)合電離反應(yīng)所生成的離子,看似會(huì)促進(jìn)聯(lián)合電離反應(yīng)的正向發(fā)生,但由于置換電離反應(yīng)較多,生成的離子種類也多,有的離子生成物與聯(lián)合電離反應(yīng)的離子生成物相同,一定程度上依然會(huì)抑制聯(lián)合電離反應(yīng)的正向發(fā)生,從計(jì)算結(jié)果來看,置換電離反應(yīng)的抑制作用比促進(jìn)作用大。若將置換電離反應(yīng)從程序中抹去,即只考慮電子激發(fā)電離反應(yīng)和聯(lián)合電離反應(yīng),可以看到流場(chǎng)的電子密度整體偏大。這是因?yàn)殡娮蛹ぐl(fā)電離反應(yīng)消耗聯(lián)合電離反應(yīng)生成的游離電子的同時(shí)不會(huì)產(chǎn)生與其相同的離子產(chǎn)物,對(duì)聯(lián)合電離的發(fā)生只有促進(jìn)作用;另一方面,除去置換反應(yīng)也使聯(lián)合電離反應(yīng)受到的抑制作用有所減弱。

    圖9 電子數(shù)密度對(duì)比(單位:m-3)

    圖10為不同工況下沿x方向探測(cè)點(diǎn)的垂直方向上的最大電子數(shù)密度比較。從圖中可以看出,電子密度沿x軸方向逐漸減小,這與沿x軸溫度逐漸降低,反應(yīng)數(shù)減少有關(guān)。從3條曲線兩兩之間的差值來看,工況2下的電子數(shù)密度始終最小,工況3下的電子數(shù)密度始終最大,與密度云圖所展示的差別一致,且工況2與工況1之間的差距比工況3與工況1之間的差距大。由此可以得出結(jié)論,本文所考慮的置換電離反應(yīng)對(duì)聯(lián)合電離反應(yīng)的抑制作用比電子激發(fā)電離反應(yīng)對(duì)聯(lián)合電離反應(yīng)的促進(jìn)作用更明顯。

    圖10 不同探測(cè)點(diǎn)垂直方向上的最大電子數(shù)密度

    圖11為3種工況下飛行器頭部熱流密度的分布情況。無論抹去哪種電離反應(yīng)都會(huì)使物面的熱流密度有所升高,這種差距沿x軸方向逐漸減小,到飛行器肩部以后基本重合。因?yàn)殡婋x反應(yīng)大多是吸熱反應(yīng),抹去其中一種電離反應(yīng)都會(huì)使電離反應(yīng)發(fā)生次數(shù)減少,從而減小了氣動(dòng)熱在電離反應(yīng)上的耗散,但由于在80 km處空氣較為稀薄,總的電離反應(yīng)數(shù)本來就小,所以會(huì)出現(xiàn)圖11中工況2和工況3的熱流密度雖有提升但差距不大的情況。

    圖11 不同工況飛行器頭部熱流密度的分布

    從定量的角度分析,3種工況下物面最大熱流密度保持在同一數(shù)量級(jí),其中,工況2的物面最大熱流密度為825 kW/m2,工況3的物面最大熱流密度為870 kW/m2,以工況1為參照進(jìn)行比較,差值百分比分別為4.43%和10.13%,也就是說,抹去任一種電離反應(yīng)帶來的熱流密度變動(dòng)不超過11%,其中抹去電子激發(fā)電離反應(yīng)對(duì)熱流密度的的影響尤其小。表明電子激發(fā)電離反應(yīng)消耗的熱量較少,這是因?yàn)楸疚乃紤]的電子激發(fā)電離反應(yīng)都為單原子分子與電子的反應(yīng),單原子分子本身的電離活化能就比雙原子分子的低,因而反應(yīng)所需吸收的熱量也就較小。

    5 結(jié)束語(yǔ)

    本文從電離反應(yīng)發(fā)生機(jī)理出發(fā),對(duì)不同電離反應(yīng)采用不同反應(yīng)發(fā)生判據(jù),采用人為增大電子質(zhì)量1 000倍的方法處理電子運(yùn)動(dòng),計(jì)算了高超聲速稀薄流域的二維軸對(duì)稱流場(chǎng),研究了電離反應(yīng)對(duì)高超聲速飛行器流場(chǎng)特性的影響以及不同種類電離反應(yīng)對(duì)流場(chǎng)的影響,結(jié)論如下:

    ①通過模擬RAMC-Ⅱ試驗(yàn)中飛行器的高空運(yùn)行情況,驗(yàn)證了本文所采用電離模塊的正確性,表明從分子動(dòng)力學(xué)角度出發(fā)考慮電離反應(yīng)是可行的。

    ②相同飛行高度下,與只考慮非電離反應(yīng)的情況對(duì)比,考慮電離反應(yīng)會(huì)使流場(chǎng)溫度整體下降,激波厚度變薄,飛行器表面的熱流密度也會(huì)有所降低,但電離反應(yīng)對(duì)高超聲速下飛行器表面氣動(dòng)特性沒有明顯的影響。

    ③聯(lián)合電離反應(yīng)是所有電離反應(yīng)的基礎(chǔ),是生成游離電子的主要反應(yīng),本文所考慮的置換電離反應(yīng)對(duì)稀薄流流場(chǎng)游離電子的生成有抑制作用,而電子激發(fā)電離反應(yīng)對(duì)游離電子的生成有促進(jìn)作用,且置換電離反應(yīng)的抑制作用比電子激發(fā)電離反應(yīng)的促進(jìn)作用更明顯。置換反應(yīng)和電子激發(fā)電離反應(yīng)均會(huì)使流場(chǎng)的熱流密度峰值減小,但影響不超過11%,其中電子激發(fā)電離反應(yīng)的影響尤其小。

    猜你喜歡
    物面電離飛行器
    電離與離子反應(yīng)高考探源
    高超聲速飛行器
    激波/湍流邊界層干擾壓力脈動(dòng)特性數(shù)值研究1)
    水的電離平衡問題解析
    復(fù)雜飛行器的容錯(cuò)控制
    電子制作(2018年2期)2018-04-18 07:13:25
    如何復(fù)習(xí)“水的電離”
    讓吸盤掛鉤更牢固
    新型單面陣自由曲面光學(xué)測(cè)量方法成像特性仿真
    神秘的飛行器
    彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究
    亚洲aⅴ乱码一区二区在线播放| 夫妻性生交免费视频一级片| 又爽又黄a免费视频| 亚洲av综合色区一区| 男人爽女人下面视频在线观看| 亚洲欧美日韩卡通动漫| 美女福利国产在线 | 22中文网久久字幕| 人人妻人人澡人人爽人人夜夜| 一个人免费看片子| 久久久国产一区二区| 国产精品秋霞免费鲁丝片| 街头女战士在线观看网站| 少妇人妻 视频| 久久人人爽人人爽人人片va| 亚洲欧美日韩卡通动漫| 免费黄色在线免费观看| 新久久久久国产一级毛片| 久久毛片免费看一区二区三区| 一个人看视频在线观看www免费| 亚洲欧美一区二区三区黑人 | 2021少妇久久久久久久久久久| 亚洲成人手机| 少妇猛男粗大的猛烈进出视频| 女人十人毛片免费观看3o分钟| 搡老乐熟女国产| 亚州av有码| 亚洲精品乱久久久久久| 久久久精品94久久精品| 搡老乐熟女国产| 看免费成人av毛片| 99热这里只有是精品50| 国产精品蜜桃在线观看| 99久久精品一区二区三区| 网址你懂的国产日韩在线| av在线观看视频网站免费| 最近最新中文字幕免费大全7| 在线观看国产h片| 男人和女人高潮做爰伦理| 久久精品久久久久久久性| 成年av动漫网址| 久久婷婷青草| 免费看光身美女| 国产黄片美女视频| 欧美区成人在线视频| 在线播放无遮挡| 狠狠精品人妻久久久久久综合| 男人爽女人下面视频在线观看| 欧美zozozo另类| 亚洲色图综合在线观看| 91精品国产国语对白视频| 久久 成人 亚洲| 老司机影院成人| 日韩一区二区三区影片| 欧美 日韩 精品 国产| 国产欧美日韩一区二区三区在线 | 欧美人与善性xxx| 久久久久久久久久人人人人人人| 日本与韩国留学比较| 国产精品久久久久久精品古装| 国产免费视频播放在线视频| videossex国产| 偷拍熟女少妇极品色| 久久精品国产亚洲av涩爱| 久久毛片免费看一区二区三区| 97超碰精品成人国产| 欧美日韩一区二区视频在线观看视频在线| 亚洲精品日韩av片在线观看| 国产精品爽爽va在线观看网站| 中文字幕精品免费在线观看视频 | 国产片特级美女逼逼视频| 国产深夜福利视频在线观看| a级一级毛片免费在线观看| 伦精品一区二区三区| 尾随美女入室| 99精国产麻豆久久婷婷| av天堂中文字幕网| 91午夜精品亚洲一区二区三区| 国产成人免费观看mmmm| 一级毛片久久久久久久久女| 99久久综合免费| 99热国产这里只有精品6| 热re99久久精品国产66热6| 国产老妇伦熟女老妇高清| 91精品国产九色| 老女人水多毛片| 丰满少妇做爰视频| 久久久a久久爽久久v久久| 如何舔出高潮| 伊人久久精品亚洲午夜| 日韩制服骚丝袜av| 亚洲成人中文字幕在线播放| 欧美精品一区二区免费开放| 精品一区二区免费观看| 少妇人妻 视频| 99热国产这里只有精品6| 国产女主播在线喷水免费视频网站| 我的老师免费观看完整版| 另类亚洲欧美激情| 色视频www国产| 乱系列少妇在线播放| 国产一区二区在线观看日韩| 国产一区二区在线观看日韩| h日本视频在线播放| 国产免费一区二区三区四区乱码| 麻豆国产97在线/欧美| av国产免费在线观看| 日韩欧美 国产精品| 不卡视频在线观看欧美| 18+在线观看网站| 如何舔出高潮| 人妻系列 视频| 热99国产精品久久久久久7| 黄色怎么调成土黄色| 久久久久久伊人网av| 精品久久久久久久久亚洲| 亚洲av中文av极速乱| 制服丝袜香蕉在线| 国产精品久久久久久久电影| 九草在线视频观看| 亚洲在久久综合| av视频免费观看在线观看| 欧美日韩综合久久久久久| 晚上一个人看的免费电影| 女人十人毛片免费观看3o分钟| 99热网站在线观看| 哪个播放器可以免费观看大片| 91午夜精品亚洲一区二区三区| 成人毛片a级毛片在线播放| 久久久久久久久久人人人人人人| 亚洲欧美日韩东京热| 午夜免费观看性视频| 亚洲国产毛片av蜜桃av| 晚上一个人看的免费电影| 亚洲欧美日韩另类电影网站 | 欧美日韩精品成人综合77777| 亚洲av欧美aⅴ国产| 黑丝袜美女国产一区| 看免费成人av毛片| 国产精品久久久久久久电影| 欧美精品人与动牲交sv欧美| 黄片无遮挡物在线观看| 男人舔奶头视频| 亚洲av国产av综合av卡| 亚洲一区二区三区欧美精品| 国产成人精品久久久久久| 免费观看在线日韩| 久久久久久久久久人人人人人人| h日本视频在线播放| 一个人看的www免费观看视频| 日韩,欧美,国产一区二区三区| 人妻制服诱惑在线中文字幕| 久久这里有精品视频免费| 亚洲人成网站在线观看播放| 黄色欧美视频在线观看| 久久精品熟女亚洲av麻豆精品| 夜夜爽夜夜爽视频| 青春草视频在线免费观看| 久久久久久久久久久丰满| 99热这里只有精品一区| av线在线观看网站| videos熟女内射| 精品亚洲成a人片在线观看 | 纵有疾风起免费观看全集完整版| 日本免费在线观看一区| 99视频精品全部免费 在线| 最近中文字幕2019免费版| 大香蕉久久网| 国产在线免费精品| 亚洲av二区三区四区| 免费看光身美女| 亚洲精品一二三| 男女无遮挡免费网站观看| 久久精品国产鲁丝片午夜精品| 国产无遮挡羞羞视频在线观看| 美女高潮的动态| 亚洲人成网站在线播| 老熟女久久久| 国产精品麻豆人妻色哟哟久久| 丝袜喷水一区| 国产午夜精品一二区理论片| 成人一区二区视频在线观看| a级一级毛片免费在线观看| 一级毛片aaaaaa免费看小| 人人妻人人爽人人添夜夜欢视频 | 亚洲国产精品成人久久小说| 在线观看av片永久免费下载| 尾随美女入室| 欧美日韩视频高清一区二区三区二| 国产白丝娇喘喷水9色精品| 男人狂女人下面高潮的视频| 91精品国产九色| 久久这里有精品视频免费| 国产精品一区二区性色av| 十分钟在线观看高清视频www | 丰满乱子伦码专区| 在线 av 中文字幕| 黄色欧美视频在线观看| 熟女电影av网| 久久久成人免费电影| 久久 成人 亚洲| 一级黄片播放器| 黄色欧美视频在线观看| 天堂中文最新版在线下载| 九九久久精品国产亚洲av麻豆| 久久久久人妻精品一区果冻| 亚洲四区av| 国产大屁股一区二区在线视频| 国产在线免费精品| 免费高清在线观看视频在线观看| 大码成人一级视频| 一区二区三区乱码不卡18| 在线亚洲精品国产二区图片欧美 | 日本欧美视频一区| 亚洲色图av天堂| 日韩中字成人| 视频中文字幕在线观看| 久久精品夜色国产| 在线观看免费视频网站a站| 熟女电影av网| 免费不卡的大黄色大毛片视频在线观看| 啦啦啦视频在线资源免费观看| 内地一区二区视频在线| 伊人久久精品亚洲午夜| 亚洲欧美精品专区久久| 夫妻午夜视频| 亚洲四区av| 亚洲av成人精品一区久久| 高清黄色对白视频在线免费看 | 七月丁香在线播放| 最近手机中文字幕大全| 狠狠精品人妻久久久久久综合| 99re6热这里在线精品视频| 最黄视频免费看| 欧美日韩视频高清一区二区三区二| 亚洲精品中文字幕在线视频 | 成人美女网站在线观看视频| 日本免费在线观看一区| 18禁裸乳无遮挡动漫免费视频| 婷婷色综合www| 国产精品伦人一区二区| 午夜福利在线观看免费完整高清在| 久久精品国产亚洲av天美| 内射极品少妇av片p| 国产 一区精品| 国产视频首页在线观看| 久久久色成人| 久热这里只有精品99| 久久久亚洲精品成人影院| 伦理电影免费视频| 国产精品一二三区在线看| 色5月婷婷丁香| 亚洲天堂av无毛| 天美传媒精品一区二区| 久久久久国产精品人妻一区二区| 在线观看一区二区三区激情| 99精国产麻豆久久婷婷| 亚州av有码| 欧美xxxx黑人xx丫x性爽| 国产爽快片一区二区三区| 日韩国内少妇激情av| 亚洲熟女精品中文字幕| 啦啦啦中文免费视频观看日本| 美女视频免费永久观看网站| 亚洲无线观看免费| 亚洲成人中文字幕在线播放| 美女内射精品一级片tv| 又黄又爽又刺激的免费视频.| 身体一侧抽搐| 午夜激情久久久久久久| 在现免费观看毛片| 久久精品熟女亚洲av麻豆精品| 免费人妻精品一区二区三区视频| 欧美少妇被猛烈插入视频| 男人爽女人下面视频在线观看| 国产免费又黄又爽又色| 日韩av不卡免费在线播放| 少妇裸体淫交视频免费看高清| 免费黄网站久久成人精品| 波野结衣二区三区在线| 国产精品国产三级国产av玫瑰| 大香蕉97超碰在线| 精品视频人人做人人爽| 日韩人妻高清精品专区| 国产视频首页在线观看| 亚洲成人一二三区av| www.av在线官网国产| 美女国产视频在线观看| 国内精品宾馆在线| h视频一区二区三区| 日韩一本色道免费dvd| 青青草视频在线视频观看| 在线亚洲精品国产二区图片欧美 | 少妇熟女欧美另类| 亚洲性久久影院| 日韩不卡一区二区三区视频在线| 欧美xxxx黑人xx丫x性爽| 精品久久久噜噜| 亚洲av欧美aⅴ国产| 在线观看美女被高潮喷水网站| 18禁动态无遮挡网站| 国产在线免费精品| 国产成人免费观看mmmm| 天堂8中文在线网| 亚洲欧美一区二区三区国产| 少妇精品久久久久久久| 久久久久久久久久成人| 免费观看a级毛片全部| 国模一区二区三区四区视频| 男女无遮挡免费网站观看| 日产精品乱码卡一卡2卡三| 综合色丁香网| 高清毛片免费看| 午夜福利在线观看免费完整高清在| 有码 亚洲区| 18禁动态无遮挡网站| 免费av不卡在线播放| 一本—道久久a久久精品蜜桃钙片| 亚洲欧美中文字幕日韩二区| 午夜免费男女啪啪视频观看| av国产免费在线观看| 插阴视频在线观看视频| 赤兔流量卡办理| 成人国产麻豆网| 亚洲av免费高清在线观看| 国产伦精品一区二区三区四那| 亚洲欧美日韩无卡精品| 少妇猛男粗大的猛烈进出视频| av在线app专区| 久久久亚洲精品成人影院| 一二三四中文在线观看免费高清| tube8黄色片| 久久99蜜桃精品久久| 毛片一级片免费看久久久久| 国产日韩欧美在线精品| 日本av免费视频播放| 日本一二三区视频观看| 91精品一卡2卡3卡4卡| 视频区图区小说| 九九爱精品视频在线观看| 精品久久久噜噜| 久久热精品热| 午夜福利视频精品| 久久精品熟女亚洲av麻豆精品| 国产 一区精品| 一区二区三区免费毛片| 欧美变态另类bdsm刘玥| 国产精品无大码| 一本久久精品| 欧美激情极品国产一区二区三区 | 深爱激情五月婷婷| 精品亚洲乱码少妇综合久久| 男女免费视频国产| 嫩草影院新地址| 尤物成人国产欧美一区二区三区| 纵有疾风起免费观看全集完整版| 18禁在线播放成人免费| 老熟女久久久| 国产亚洲5aaaaa淫片| 欧美日韩国产mv在线观看视频 | 黑人猛操日本美女一级片| 熟妇人妻不卡中文字幕| 午夜福利网站1000一区二区三区| 欧美一区二区亚洲| 亚洲精品国产av蜜桃| 国产女主播在线喷水免费视频网站| 久久ye,这里只有精品| 中文字幕制服av| 中文天堂在线官网| 偷拍熟女少妇极品色| av黄色大香蕉| 国产亚洲5aaaaa淫片| 一级毛片我不卡| 欧美日本视频| 亚洲天堂av无毛| 国产黄色视频一区二区在线观看| 午夜免费男女啪啪视频观看| 午夜福利在线在线| 一区二区三区免费毛片| 尾随美女入室| 久久久久久久久大av| 天堂俺去俺来也www色官网| 国产精品一区二区在线观看99| 日韩av免费高清视频| 尤物成人国产欧美一区二区三区| 久久久午夜欧美精品| videossex国产| 日韩欧美一区视频在线观看 | 欧美日韩精品成人综合77777| 老司机影院成人| 不卡视频在线观看欧美| 午夜激情久久久久久久| 国产69精品久久久久777片| 插逼视频在线观看| 欧美性感艳星| 又爽又黄a免费视频| 只有这里有精品99| 毛片女人毛片| 人妻一区二区av| 国产亚洲欧美精品永久| 中文字幕制服av| 日韩av在线免费看完整版不卡| 天天躁日日操中文字幕| 亚洲欧美成人精品一区二区| 男女下面进入的视频免费午夜| 精品一区二区三区视频在线| 五月开心婷婷网| 91久久精品国产一区二区三区| 国产久久久一区二区三区| 22中文网久久字幕| 一区二区三区乱码不卡18| 伦精品一区二区三区| 国产男人的电影天堂91| 免费看日本二区| 成人毛片a级毛片在线播放| 亚洲aⅴ乱码一区二区在线播放| 小蜜桃在线观看免费完整版高清| 国产精品av视频在线免费观看| av女优亚洲男人天堂| 大香蕉久久网| 欧美国产精品一级二级三级 | 特大巨黑吊av在线直播| 交换朋友夫妻互换小说| 午夜免费鲁丝| 中文字幕久久专区| 日韩不卡一区二区三区视频在线| 人妻系列 视频| 你懂的网址亚洲精品在线观看| 久久99热这里只频精品6学生| 99热这里只有是精品50| 婷婷色综合大香蕉| 国产中年淑女户外野战色| 亚洲成人av在线免费| 最近中文字幕高清免费大全6| 国产成人午夜福利电影在线观看| 99久久中文字幕三级久久日本| 成人二区视频| 中文字幕av成人在线电影| 日本一二三区视频观看| 少妇高潮的动态图| 超碰av人人做人人爽久久| 午夜免费观看性视频| av免费在线看不卡| 嘟嘟电影网在线观看| 亚洲精品自拍成人| 中文乱码字字幕精品一区二区三区| 久久国产精品大桥未久av | 六月丁香七月| 少妇熟女欧美另类| 在线精品无人区一区二区三 | 一级毛片我不卡| 欧美激情极品国产一区二区三区 | 国产精品一二三区在线看| 日韩av在线免费看完整版不卡| 性色avwww在线观看| 日本黄色日本黄色录像| 国产精品爽爽va在线观看网站| 精品视频人人做人人爽| 日本vs欧美在线观看视频 | 欧美高清性xxxxhd video| 汤姆久久久久久久影院中文字幕| 国产老妇伦熟女老妇高清| 插逼视频在线观看| 久久韩国三级中文字幕| 久久毛片免费看一区二区三区| 自拍欧美九色日韩亚洲蝌蚪91 | 99久久人妻综合| 欧美日韩精品成人综合77777| 久久人人爽人人片av| 永久免费av网站大全| 国产精品久久久久成人av| 国产成人freesex在线| 亚洲精品国产成人久久av| 在线观看一区二区三区激情| 久久久成人免费电影| 狠狠精品人妻久久久久久综合| 亚洲美女黄色视频免费看| 韩国av在线不卡| av在线蜜桃| 大香蕉97超碰在线| 欧美日韩一区二区视频在线观看视频在线| 午夜免费鲁丝| 18+在线观看网站| 有码 亚洲区| 亚洲丝袜综合中文字幕| 各种免费的搞黄视频| 欧美 日韩 精品 国产| 日日啪夜夜撸| 亚洲精品第二区| 日韩人妻高清精品专区| 日韩亚洲欧美综合| 国产精品久久久久久久电影| 18禁在线无遮挡免费观看视频| 亚洲中文av在线| 国产成人免费无遮挡视频| 色综合色国产| 国精品久久久久久国模美| 联通29元200g的流量卡| 精品一区二区免费观看| 特大巨黑吊av在线直播| 久久精品久久久久久噜噜老黄| 99热全是精品| av在线观看视频网站免费| 秋霞在线观看毛片| 2022亚洲国产成人精品| 一级黄片播放器| 国产亚洲5aaaaa淫片| 97精品久久久久久久久久精品| 午夜福利高清视频| 狂野欧美激情性bbbbbb| 欧美三级亚洲精品| 日本欧美国产在线视频| 久久久久久久久久人人人人人人| 天美传媒精品一区二区| 亚洲在久久综合| 久久99热这里只频精品6学生| 久久精品国产亚洲av涩爱| 天美传媒精品一区二区| 久久久久久久久久成人| 人人妻人人添人人爽欧美一区卜 | 成人黄色视频免费在线看| 建设人人有责人人尽责人人享有的 | 欧美日本视频| 伦理电影免费视频| 18+在线观看网站| 又粗又硬又长又爽又黄的视频| 午夜免费观看性视频| 国产亚洲午夜精品一区二区久久| 青春草视频在线免费观看| 久久久国产一区二区| 在线观看免费视频网站a站| 国产精品免费大片| 久热这里只有精品99| 精品久久国产蜜桃| av视频免费观看在线观看| av一本久久久久| 亚洲精品一二三| 天堂8中文在线网| 日韩国内少妇激情av| 啦啦啦视频在线资源免费观看| 精品酒店卫生间| 好男人视频免费观看在线| 国产男女超爽视频在线观看| 最近手机中文字幕大全| 全区人妻精品视频| 大香蕉97超碰在线| 两个人的视频大全免费| 熟女人妻精品中文字幕| 亚洲精品自拍成人| 日韩一区二区三区影片| 日本av免费视频播放| 欧美极品一区二区三区四区| 日韩制服骚丝袜av| 丝袜喷水一区| 国产精品欧美亚洲77777| 一区二区三区四区激情视频| 欧美激情国产日韩精品一区| 97热精品久久久久久| 亚洲av欧美aⅴ国产| 热99国产精品久久久久久7| 日韩av在线免费看完整版不卡| 乱系列少妇在线播放| 婷婷色麻豆天堂久久| 91精品国产国语对白视频| 免费观看av网站的网址| 中文乱码字字幕精品一区二区三区| 亚洲精品日韩在线中文字幕| 免费在线观看成人毛片| av免费观看日本| 乱码一卡2卡4卡精品| 亚洲精品国产成人久久av| 亚洲av中文字字幕乱码综合| 色视频在线一区二区三区| 国产精品99久久久久久久久| 国产免费一区二区三区四区乱码| 在线天堂最新版资源| 欧美日韩国产mv在线观看视频 | 777米奇影视久久| 免费看不卡的av| 人人妻人人添人人爽欧美一区卜 | 欧美丝袜亚洲另类| 久久久亚洲精品成人影院| 最近2019中文字幕mv第一页| 一级毛片黄色毛片免费观看视频| 国内精品宾馆在线| av网站免费在线观看视频| 天堂俺去俺来也www色官网| 国产精品无大码| 在线亚洲精品国产二区图片欧美 | 小蜜桃在线观看免费完整版高清| 91久久精品国产一区二区三区| 精品少妇黑人巨大在线播放| 女人十人毛片免费观看3o分钟| 波野结衣二区三区在线| 在线 av 中文字幕| 日产精品乱码卡一卡2卡三| 国产精品蜜桃在线观看| 亚洲国产av新网站| 免费少妇av软件| 亚洲精品国产色婷婷电影| 国产在线视频一区二区| 大片免费播放器 马上看| 身体一侧抽搐| 日韩一区二区三区影片| 男女下面进入的视频免费午夜| 久久久a久久爽久久v久久| 免费久久久久久久精品成人欧美视频 | 99久久综合免费| 人人妻人人澡人人爽人人夜夜| 热99国产精品久久久久久7| 国产黄色免费在线视频| 欧美日韩亚洲高清精品| 一本—道久久a久久精品蜜桃钙片| 在线观看免费视频网站a站| 色视频在线一区二区三区|