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

    三維等離子體MHD氣動熱環(huán)境數(shù)值模擬

    2017-11-20 03:12:33丁明松江濤董維中高鐵鎖劉慶宗
    航空學(xué)報(bào) 2017年8期
    關(guān)鍵詞:磁流體駐點(diǎn)超聲速

    丁明松 , 江濤, 董維中, 高鐵鎖, 劉慶宗

    中國空氣動力研究與發(fā)展中心 計(jì)算空氣動力研究所, 綿陽 621000

    三維等離子體MHD氣動熱環(huán)境數(shù)值模擬

    丁明松*, 江濤, 董維中, 高鐵鎖, 劉慶宗

    中國空氣動力研究與發(fā)展中心 計(jì)算空氣動力研究所, 綿陽 621000

    電磁流動控制技術(shù)是一個多學(xué)科交叉融合的重要研究方向,在高超聲速飛行器氣動特性優(yōu)化、氣動熱環(huán)境減緩、邊界層轉(zhuǎn)捩和等離子體分布等流動控制方面顯示出廣闊的應(yīng)用前景??紤]高超聲速飛行器繞流流場中發(fā)生的離解、復(fù)合、電離和置換等化學(xué)反應(yīng),氣體分子振動能激發(fā)以及化學(xué)非平衡效應(yīng),耦合電磁場作用并基于低磁雷諾數(shù)假設(shè),通過數(shù)值模擬求解三維非平衡Navier-Stokes流場控制方程和Maxwell電磁場控制方程,建立磁場與三維等離子體流場耦合數(shù)值模擬方法及程序,采用典型算例進(jìn)行考核。在此基礎(chǔ)上,開展不同條件下磁場對再入三維等離子體流場以及氣動熱環(huán)境影響分析。研究表明:建立的高超聲速飛行器的等離子體流場與磁場耦合計(jì)算方法及程序,其數(shù)值模擬結(jié)果與文獻(xiàn)符合,外加磁場使飛行器頭部弓形激波外推,磁場強(qiáng)度越強(qiáng),激波面外推距離越大;不同磁場強(qiáng)度環(huán)境下,流場中溫度峰值大小略有變化,變化幅度較??;磁場對絕大部分區(qū)域的熱流有減緩作用,作用的大小與飛行高度、馬赫數(shù)以及磁場的配置緊密相關(guān);當(dāng)前的計(jì)算條件下,飛行的高度越高,磁場的作用越明顯。

    MHD; 等離子體; 化學(xué)非平衡; 數(shù)值模擬; 氣動熱環(huán)境

    高超聲速飛行器再入大氣層的過程中,速度可高達(dá)5 km/s以上甚至達(dá)到第一、第二宇宙速度,此時與空氣發(fā)生強(qiáng)烈的相互作用,飛行器激波層內(nèi)溫度高達(dá)8 000 K以上[1],高溫流場中發(fā)生復(fù)雜的物理化學(xué)現(xiàn)象,如離解、復(fù)合、電離和交換等,氣體分子振動能量模態(tài)激活,形成具有弱導(dǎo)電性等離子體鞘套?;诖帕黧w動力學(xué)(MHD),利用機(jī)載磁場發(fā)生裝置向弱導(dǎo)電性的等離子體鞘套引入適當(dāng)?shù)膭恿亢湍芰繉Ω叱曀亠w行器流場進(jìn)行控制,這種電磁流動控制技術(shù)可以有效改進(jìn)和提高飛行器的氣動特性、通信性能、隱身與突防能力,在高超聲速飛行器氣動力控制、氣動熱環(huán)境減緩、邊界層轉(zhuǎn)捩控制和等離子體電子密度分布調(diào)整等方面顯示出廣闊的應(yīng)用前景。

    由于電磁場變化的特征時間和流動的特征時間存在數(shù)量級差別,同時高超聲速飛行器外流場等離子體導(dǎo)電性一般較弱,因此,數(shù)值模擬高超聲速飛行器外流場時常采用“低磁雷諾數(shù)近似”:在電導(dǎo)率較小時,感應(yīng)磁場相對于外加磁場很弱,基本上可以忽略,此時可以假設(shè)外磁場未受流動干擾,舍去電磁交叉項(xiàng),使方程形式得到簡化。采用低磁雷諾數(shù)MHD方程組,避免了數(shù)值模擬時的奇性和剛性問題,減小了非必要的繁復(fù)計(jì)算,因而計(jì)算效率相對較高。

    高溫空氣電導(dǎo)率是高超聲速電磁流動控制數(shù)值模擬最重要的參數(shù)之一。影響混合氣體電導(dǎo)率的因素很多,要得到較為準(zhǔn)確的氣體電導(dǎo)率,依賴于高超聲速流場的準(zhǔn)確模擬,包括混合氣體電離程度、溫度、密度、組分和熱力學(xué)狀態(tài)等參數(shù)的精確模擬。因此,在高超聲速飛行器磁流體數(shù)值模擬時,必須考慮高溫流場中發(fā)生的復(fù)雜物理化學(xué)現(xiàn)象,包括高溫氣體熱化學(xué)反應(yīng)、氣體分子能量模態(tài)激活及其非平衡效應(yīng),即高溫氣體非平衡效應(yīng),得到詳細(xì)的流場參數(shù)。

    近20年來,隨著人工電離技術(shù)、超導(dǎo)磁體技術(shù)以及計(jì)算機(jī)性能的突飛猛進(jìn),高超聲速飛行器磁流體控制技術(shù)迎來了新的研究熱潮。2005年,Maccormack[2]考慮低磁雷諾數(shù)近似,對半球圓筒的磁流體高超聲速流動進(jìn)行了數(shù)值模擬;2010年,Lee等[3]采用低磁雷諾數(shù)假設(shè)和高溫氣體凍結(jié)與平衡假設(shè),數(shù)值模擬了飛行高度60 km、馬赫數(shù)10的半球圓柱流場,分析了磁場強(qiáng)度對流場氣動熱環(huán)境的影響;2011年,Bisek和Poggie[4]考慮低磁雷諾數(shù)假設(shè),利用燒蝕引入“種子”(假設(shè)它不參加化學(xué)反應(yīng),只提高氣體電導(dǎo)率),數(shù)值模擬了帶尾舵鈍頭橢圓錐體馬赫數(shù)Ma=8和10的高超聲速流場,分析了外磁場對鈍頭、尾舵前緣熱流的影響;2012年,Chernyshev等[5]采用兩溫度模型,將等離子體分為電子、原子和離子3種成分,考慮高溫氣體非平衡效應(yīng)和低磁雷諾數(shù)近似,分析了磁場對尖錐柱體斜激波的干擾;2013年,F(xiàn)ujino和Ishikawa[6]考慮高溫氣體效應(yīng),數(shù)值模擬了沿再入軌道的二維大鈍頭低磁雷諾數(shù)磁流體流動;2015年,Takahashi等[7]數(shù)值模擬了火星大氣二維軸對稱球錐低磁雷諾MHD流動,發(fā)現(xiàn)通過電磁流動控制技術(shù)可以有效減速,降低熱流;2015年,Masuda等[8]考慮高溫空氣化學(xué)反應(yīng),對三維球錐外形低磁雷諾數(shù)MHD進(jìn)行了數(shù)值模擬,發(fā)現(xiàn)傾斜磁場或者非零迎角情況下,磁場控制會引入側(cè)向力。

    國內(nèi)的MHD研究早期多集中于冶金和宇宙物理學(xué)領(lǐng)域,高超聲速飛行器電磁流體控制研究起步較晚。2005年,赫新等[9]采用無波動、無自由參數(shù)的耗散(Non-oscillatory,containing No free parameters and Dissipative, NND)格式,數(shù)值模擬了理想磁流體一維MHD激波管流動和二維噴管MHD流動;2007年,潘勇[10]采用非結(jié)構(gòu)網(wǎng)格對二維理想MHD繞流逆風(fēng)格式解法進(jìn)行了研究;2013年,孫曉輝[11]采用Euler方程,考慮化學(xué)反應(yīng),開展了新型沖壓推進(jìn)系統(tǒng)的波系結(jié)構(gòu)及其電磁流動控制的研究;2014年,張洪量[12]考慮低磁雷諾數(shù)假設(shè)和電子束誘導(dǎo)電離,發(fā)現(xiàn)MHD控制能有效提高飛行器的飛行性能。

    可以看出,在高超聲速飛行器外場的磁流動控制研究方面,國外很重要的一個發(fā)展趨勢是:在低磁雷諾數(shù)假設(shè)下,考慮高溫氣體熱化學(xué)非平衡效應(yīng),研究不同磁場配置條件(如磁場類別、強(qiáng)度、方向等)對各種外形高超聲速飛行器流動的影響規(guī)律,達(dá)到飛行器氣動特性控制的目的。但這些研究大多較為零散,系統(tǒng)地分析磁場強(qiáng)度、飛行高度以及飛行速度等因素對高超聲速飛行器外場磁流體控制影響規(guī)律的研究較少。而國內(nèi)高超聲速飛行器外場的磁流體控制研究,數(shù)值模擬多以二維或簡化氣體模型(如理想氣體、平衡氣體等)為主,考慮高溫氣體非平衡效應(yīng)的三維磁流體控制研究并不多見。

    筆者所在團(tuán)隊(duì)對高超聲速飛行器非平衡流場進(jìn)行了較為廣泛的研究[13-17]。本文在此基礎(chǔ)上,耦合考慮高超聲速飛行器高溫氣體非平衡效應(yīng)與電磁場效應(yīng),建立三維低磁雷諾數(shù)化學(xué)非平衡MHD數(shù)值模擬方法和程序,開展高超聲速等離子體流場與電磁場耦合效應(yīng)研究,較為系統(tǒng)地分析不同磁場強(qiáng)度、飛行高度以及飛行速度條件下磁場對再入三維等離子體流場以及氣動熱環(huán)境的影響。

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

    在飛行器實(shí)際的飛行外場中,流體介質(zhì)一般符合低電導(dǎo)率特征(電導(dǎo)率最大值σmax一般為102S/m 量級),通常滿足低磁雷諾數(shù)(Rem?1)假設(shè)[18]。此時,感應(yīng)磁場相對于外磁場基本可以忽略,控制方程與常規(guī)流體力學(xué)控制方程類似。右端出現(xiàn)電磁源項(xiàng)WMHD,無量綱形式為

    (1)

    理想氣體守恒變量為

    Q=[ρρuρvρwρEt]T

    非平衡氣體守恒變量為

    Q=[ρjρρuρvρwρEt]T

    式中:ρ為氣體密度;ρj為氣體組分j的密度;u、v、w為直角坐標(biāo)系3個方向速度;Et為氣體總能;Re為雷諾數(shù);F、G、H和Fv、Gv、Hv分別為3個方向的無黏通量和黏性通量;W為非平衡源項(xiàng);理想氣體電磁源項(xiàng)WMHD表達(dá)式為

    WMHD=Qm·[0 (J×Β)x(J×B)y(J×B)zJ·E]T

    (2)

    其中:J為電流密度;B為磁感應(yīng)強(qiáng)度;E為電場強(qiáng)度;Qm為磁相互作用數(shù)。

    采用有限差分方法對式(1)進(jìn)行數(shù)值離散,無黏項(xiàng)采用AUSMPW+/TVD格式離散,黏性項(xiàng)采用中心差分格式離散,時間格式為LU-SGS隱式格式?;旌蠚怏w黏性系數(shù)和熱傳導(dǎo)系數(shù)由Wilke半經(jīng)驗(yàn)公式計(jì)算,輸運(yùn)系數(shù)采用Blotter曲線擬合公式和Eucken關(guān)系公式計(jì)算,擴(kuò)散系數(shù)由等效二元擴(kuò)散模型計(jì)算,具體參數(shù)及方法見文獻(xiàn)[13-17]。

    2 物理化學(xué)模型

    2.1 化學(xué)反應(yīng)模型

    數(shù)值模擬多組分電離空氣,采用Dunn-Kang模型數(shù)值模擬高溫空氣電離流動,具體化學(xué)反應(yīng)見表1。

    表1 Dunn-Kang模型Table 1 Dunn-Kang model

    注:H1=N,NO;H2=O,NO,O2;H3=O2,N2;H4=O,N,NO。

    化學(xué)反應(yīng)式的通式可寫為

    (3)

    (4)

    其中:kfi與kbi分別為第i反應(yīng)正向與逆向反應(yīng)速率常數(shù);Mj為第j組分分子量。化學(xué)模型中所有反應(yīng)產(chǎn)生的第j組分生成源項(xiàng)為

    (5)

    2.2 電動力學(xué)模型

    根據(jù)等離子體電中性假設(shè),電流密度滿足如下電流連續(xù)性方程:

    (6)

    J由廣義歐姆定律給出:

    (7)

    (8)

    耦合式(1),數(shù)值求解式(8)得到電勢函數(shù)φ和電場E,再由式(7)得到電流密度J。

    3 計(jì)算方法驗(yàn)證

    3.1 三維等離子體流場數(shù)值模擬

    等離子體的電導(dǎo)率是磁流體數(shù)值模擬最重要的參數(shù)之一,它有賴于高溫流場參數(shù)(尤其是氣體電離程度)的精確獲得。因此,在高超聲速M(fèi)HD數(shù)值模擬過程中,高溫非平衡流場等離子體分布的準(zhǔn)確模擬是關(guān)鍵。為考核非平衡流場等離子體分布數(shù)值模擬的準(zhǔn)確性,本文以RAM-C鈍錐外形作為研究對象,該外形有電子數(shù)密度飛行試驗(yàn)結(jié)果[20-21]。

    RAM-C鈍錐外形的球頭半徑Rn=0.152 4 m,全長1.295 m,半錐角為9°。計(jì)算條件為:速度U∞=7 650.0 m/s,高度H=61,71,81 km。壁面溫度設(shè)為1 500 K,飛行迎角為0°??紤]高溫氣體非平衡效應(yīng),高溫空氣化學(xué)反應(yīng)模型采用Dunn-Kang模型,壁面采用完全非催化(NCW)和完全催化(FCW)壁面條件。

    圖1為飛行高度H=61 km時RAM-C鈍錐外形流場電子數(shù)密度。可以看出,頭部駐點(diǎn)區(qū)域電子數(shù)密度較高(電離程度較高),等離子體沿流動方向向下游擴(kuò)展,形成等離子鞘套。圖2進(jìn)一步給出了RAM-C鈍錐外形的X軸向電子數(shù)密度峰值??梢钥闯觯河?jì)算得到的電子數(shù)密度分布與飛行試驗(yàn)數(shù)據(jù)符合較好,這說明本文的計(jì)算方法可以準(zhǔn)確模擬RAM-C鈍錐外形飛行器的化學(xué)非平衡流動,能為分析等離子體流動與電磁場耦合效應(yīng)提供可信度高的流場參數(shù)。此外,完全非催化條件下計(jì)算得到的電子數(shù)密度峰值高于完全催化壁面下的計(jì)算結(jié)果。

    圖1 流場電子數(shù)密度云圖(H=61 km) Fig.1 Density contours of flow field electron number at H=61 km

    圖2 RAM-C鈍錐外形電子數(shù)密度峰值Fig.2 Peak electron number density over RAM-C blunt cone shape

    由于在同等條件下,氣體電離程度越高,等離子體電導(dǎo)率越大,電磁作用越強(qiáng)。因此,高超聲速飛行器磁流體控制設(shè)計(jì)時,可以利用飛行器頭部繞流氣體電離度高的特點(diǎn),將機(jī)載磁場設(shè)置在頭部,并選取催化速率較低的表面材料,以獲得較好的磁場作用效果。

    3.2 MHD數(shù)值模擬

    采用二維鈍錐,頭部半徑為Rn=0.09 m,半錐角為15°,長度為L=0.122 4 m;計(jì)算狀態(tài)為H=40 km,來流溫度T∞=250.35 K,來流馬赫數(shù)Ma∞=15,壁面溫度Tw=800 K。外加磁偶極子磁場,模擬實(shí)際應(yīng)用的螺線管磁場:

    (9)

    式中:(r,θ)為極坐標(biāo)單位矢量,r和θ分別為對應(yīng)的標(biāo)量,偶極子中心位于坐標(biāo)原點(diǎn),即二維鈍錐頭部圓心;B0為極軸上距離偶極子中心r0處的磁感應(yīng)強(qiáng)度,B0=2.0 T,r0=0.09 m??紤]兩種不同方向的磁偶極子配置:位于坐標(biāo)原點(diǎn)指向X軸負(fù)向的偶極子,記為M1;位于坐標(biāo)原點(diǎn)指向Y軸正向的偶極子,記為M2。為了對比分析,采用與文獻(xiàn)[19]相同的處理方法:近似認(rèn)為該條件下的流動處于熱化學(xué)平衡狀態(tài);考慮理想磁流體,由近似電離模型給出所需的流場電子摩爾分?jǐn)?shù):

    (10)

    式中:ε=10-9;T0=3 000 K;D=3 000 K;X0=0.002;tanh為雙曲正切函數(shù)。

    圖3(a)和圖3(b)為磁場矢量和磁場強(qiáng)度圖,上部為M1配置,下部為M2配置。圖3(c)和圖3(d)給出了有無磁場配置的流場電子摩爾分?jǐn)?shù),上部為無磁場配置,下部為有磁場配置??梢钥闯?,附加磁場后,激波的脫體距離加寬,電子摩爾分?jǐn)?shù)較高區(qū)域顯著增大。

    圖4(a)為流向表面熱流Q分布,S為表面弧長??梢钥闯觯?種磁場配置條件,本文均與文獻(xiàn)[19]結(jié)果符合較好。雖然兩種外磁場的大小相同,即能量和重量負(fù)載相當(dāng),但配置M2的熱流比無磁場下降了約40%,而配置M1只下降了約25%。說明在熱流控制方面,配置M2優(yōu)于M1。圖4(b)進(jìn)一步給出兩種磁場配置(上部為M1,下部為M2)洛倫茲力矢量圖(背景為壓力云圖,黑色線條為洛倫茲力矢量)。由圖可以看出:洛倫茲力作用范圍集中在頭部激波以內(nèi),絕大部分區(qū)域洛倫茲力作用方向與來流方向相反,這種反向作用力使激波脫體增加;兩種磁場配置,洛倫茲力差別很大。駐點(diǎn)附近,M1配置條件下,洛倫茲力方向幾乎與物面平行,M2配置條件下,洛倫茲力為遠(yuǎn)離壁面方向,因此,M2配置條件下駐點(diǎn)壓力較低,熱流下降幅度大。在頭部和身部結(jié)合處靠近頭部區(qū)域,M1配置條件下,洛倫茲力為遠(yuǎn)離壁面方向,M2配置條件下,洛倫茲力幾乎與物面平行,因此,這一區(qū)域,M1配置條件下,壓力較低,熱流下降幅度大。

    圖3 磁場分布和流場電子摩爾分?jǐn)?shù)云圖Fig.3 Magnetic field distribution and flow field electron molar fraction contours

    圖4 表面熱流分布和流場洛倫茲力矢量Fig.4 Surface heat flux distribution and Lorentz vector of flow field

    4 RAM-C外形算例分析

    采用3.1節(jié)給出的RAM-C鈍錐外形,高度H=81, 71, 61, 51 km,Ma∞=12~25,壁面溫度Tw=1 500 K。

    外磁場由螺線管生成,其分布方式與3.2節(jié)類似,此時r0=0.152 4 m。磁偶極子方向指向直角坐標(biāo)X軸負(fù)方向,偶極子中心點(diǎn)為坐標(biāo)系原點(diǎn),即RAM-C頭部球心,B0=0~2.0 T。

    圖5 電導(dǎo)率、感應(yīng)電流強(qiáng)度和電子數(shù)密度云圖(H=61 km, NCW) Fig.5 Density contour of conductivity, current and electron number (H=61 km, NCW)

    圖6進(jìn)一步給出了H=61,71 km下不同磁場強(qiáng)度配置頭部駐點(diǎn)線溫度分布。由圖可以看出,磁場變化,溫度峰值略有變化,但變化幅值不大;磁感應(yīng)強(qiáng)度越大,激波面外推距離越大。這與文獻(xiàn)[22]結(jié)果(圖6(c)為文獻(xiàn)[22]給出的溫度分布圖,其計(jì)算條件與本文不同,實(shí)線和虛線分別對應(yīng)不同模型尺寸,X為離開駐點(diǎn)的軸向距離,Rs為頭部半徑,U∞為來流速度)有較大相似之處。值得關(guān)注的是,B0=1.0 T時,H=61 km激波脫體距離增加到無磁場時的2倍左右,而H=71 km 激波脫體距離增加到無磁場時的5倍以上。這說明,同等磁場條件下,飛行高度越高,磁場對流場影響越明顯。

    圖6 不同磁場強(qiáng)度駐點(diǎn)線溫度分布(U∞=7 650.0 m/s)Fig.6 Temperature distribution along stagnation line in different extra magnetic fields (U∞=7 650.0 m/s)

    圖7給出不同飛行高度、不同磁場強(qiáng)度條件下表面熱流沿軸向分布。可見,在壁面大部分區(qū)域,熱流隨著外加磁場強(qiáng)度增大而減??;同一飛行高度,磁場越大,對熱流的影響越大;不同飛行高度,同等磁場條件下,飛行高度越高,磁場作用越明顯,這個結(jié)論是否具有普遍性,還需要在以后的研究中進(jìn)一步分析給出。無磁場時,熱流沿軸線平滑下降;附加磁場后,熱流曲線沿軸向出現(xiàn)起伏,磁場較大時,在X=0 m附近,熱流出現(xiàn)低谷,隨后出現(xiàn)第2個熱流峰值;磁場越大,起伏幅度越大,在H=51和61 km外加磁場較大時,表面部分區(qū)域熱流甚至比無磁場還要高。究其原因,可能是這一區(qū)域在較強(qiáng)磁場作用下,洛侖茲力使高溫氣流附著造成的。

    圖7 不同條件表面熱流軸向分布(U∞=7 650.0 m/s) Fig.7 Surface heat flux axial distribution under different conditions (U∞=7 650.0 m/s)

    圖8給出了H=61 km時壓力云圖,上部為有磁場,下部為無磁場,圖中黑弧線為磁矢量??梢钥闯?,在較強(qiáng)磁場作用下,在頭錐結(jié)合處下游之后,出現(xiàn)一個壓力升高區(qū)域,該區(qū)域洛侖茲力如圖中F所示。

    圖9(a)給出了同一飛行速度、不同飛行高度有/無磁場作用飛行器駐點(diǎn)熱流及其差量。可以看出:磁場作用使駐點(diǎn)熱流下降,飛行高度越高,下降幅度越大;飛行高度H=51~81 km,駐點(diǎn)熱流下降幅度由3.4%升至63%。圖9(b)給出了同一飛行速度、不同磁場強(qiáng)度下的飛行器駐點(diǎn)熱流、頭身結(jié)合處(X=0 m)處熱流以及錐身(X≥0 m)熱流峰值。可以看出:磁場增強(qiáng),熱流呈下降趨勢;磁場增大到一定程度后,駐點(diǎn)熱流下降幅度減小(B0=0.8,0.9,1.0 T時駐點(diǎn)熱流變化幅度很小);磁場較小時,錐身處熱流起伏不大,峰值熱流在X=0 m處取得;當(dāng)磁場B0≥0.5 T時,X=0 m 處出現(xiàn)熱流低谷,第2個峰值熱流在錐身中部取得(X=0.1 m附近)。圖9(c)給出了同一飛行高度、不同飛行速度有/無磁場作用飛行器駐點(diǎn)熱流及其差量??梢钥闯觯簾o磁場時,駐點(diǎn)熱流隨馬赫數(shù)Ma∞上升明顯升高;附加1.0 T磁場,駐點(diǎn)熱流呈下降趨勢,隨馬赫數(shù)升高,下降幅度先增大后減小(Ma∞≤15時,駐點(diǎn)熱流基本與無磁場時重合;Ma∞=17~22時,駐點(diǎn)熱流基本保持在370 kW/m2左右,不隨馬赫數(shù)變化;Ma∞>22時,駐點(diǎn)熱流隨馬赫數(shù)上升而上升),最大下降幅度接近60%。

    圖8 壓力云圖(H=61 km)Fig.8 Pressure contours (H=61 km)

    圖9 不同條件下表面熱流分布 Fig.9 Surface heat flux distribution under different conditions

    5 結(jié) 論

    1) 本文考慮高溫空氣化學(xué)反應(yīng)、分子能量模態(tài)激活以及流動中的非平衡效應(yīng),建立了低磁雷諾數(shù)三維非平衡MHD數(shù)值模擬方法及程序,不僅可以模擬高超聲速飛行器理想氣體磁流體流動,而且可以模擬高溫非平衡氣體磁流體流動。

    2) 飛行器高溫氣體等離子體流場模擬結(jié)果與飛行試驗(yàn)符合良好,可為MHD數(shù)值模擬提供較為準(zhǔn)確的流場參數(shù);高超聲速飛行器的等離子體流場與磁場耦合數(shù)值模擬結(jié)果表明,氣動熱環(huán)境結(jié)果與文獻(xiàn)符合,可信度高。

    3) 對不同飛行高度、不同馬赫數(shù)、不同磁場感應(yīng)強(qiáng)度配置的高超聲速磁流體氣動熱環(huán)境,進(jìn)行了數(shù)值模擬??梢钥闯觯捍艌鍪癸w行器頭部弓形激波外推,磁場強(qiáng)度越強(qiáng),激波面外推距離越大;不同磁場強(qiáng)度,溫度峰值大小略有變化,變化幅度較?。淮艌鰧^大部分區(qū)域的熱流有減緩作用(部分條件下,局部熱流反而增大),作用的大小與飛行高度、馬赫數(shù)以及磁場的配置緊密相關(guān);當(dāng)前的計(jì)算條件下,飛行的高度越高,磁場的作用越明顯。

    本文建立的低磁雷諾數(shù)三維非平衡MHD數(shù)值模擬方法及程序,運(yùn)用了消息通訊MPI并行技術(shù),不僅可以模擬三維鈍錐外形飛行器高超聲速磁流體流動,而且可以用于各種復(fù)雜外形飛行器磁流體并行數(shù)值計(jì)算。由于高溫等離子體流場和電磁場的耦合效應(yīng),不僅受磁場配置(磁場類型、磁感強(qiáng)度、配置方向等)影響,而且與影響飛行器高溫非平衡效應(yīng)的各種因素(飛行速度、高度、邊界、飛行器幾何類型、尺寸等)直接相關(guān)。因此,在下一步的研究工作中,可以針對這些方面,展開更加細(xì)致、深入的研究。

    [1] 樂嘉陵. 再入物理[M]. 北京: 國防工業(yè)出版社, 2005: 9-21.

    LE J L. Reentry physics[M]. Beijing: National Defence Industry Press, 2005: 9-21 (in Chinese).

    [2] MACCORMACK R W. Evaluation of the low magnetic Reynolds approximation for aerodynamic flow calculations: AIAA-2005-4780[R]. Reston, VA: AIAA, 2005.

    [3] LEE J, HUERTA M A, ZHA G C. LowRem3D MHD hypersonic equilibrium flow using high order WENO schemes: AIAA-2010-229[R]. Reston, VA: AIAA, 2010.

    [4] BISEK N J, POGGIE J. Exploration of MHD flow control for a hypersonic blunt elliptic cone with an impregnated ablator: AIAA-2011-0897[R]. Reston, VA: AIAA, 2011.

    [5] CHERNYSHEV A, KURAKIN Y, SCHMIDT A. Effect of MHD system arrangement on high speed plasma flow structure around conical body: AIAA-2012-5974[R]. Reston, VA: AIAA, 2012.

    [6] FUJINO T, ISHIKAWA M. Numerical simulation of MHD flow control along super orbital reentry trajectory: AIAA-2013-3000[R]. Reston, VA: AIAA, 2013.

    [7] TAKAHASHI T, SHIMOSAWA Y, MASUDA K, et al. Numerical study of thermal protection using magneto-hydrodynamic flow control in mars entry flight: AIAA- 2015-3365[R]. Reston, VA: AIAA, 2015.

    [8] MASUDA K, SHIMOSAWA Y, FUJINO T. Three-dimensional numerical simulation of magnetohydrodynamic flow control in reentry flight: AIAA-2015-3366[R]. Reston, VA: AIAA, 2015.

    [9] 赫新, 陳堅(jiān)強(qiáng), 鄧小剛. NND格式在多維理想磁流體方程組中的應(yīng)用[J]. 空氣動力學(xué)學(xué)報(bào), 2005, 23(3): 267-273.

    HE X, CHEN J Q, DENG X G. NND scheme’s application in multi-dimensional ideal magnetohydrodynamic equations[J]. Acta Aerodynamica Sinica, 2005, 23(3): 267-273 (in Chinese).

    [10] 潘勇. 高超聲速流場磁場干擾效應(yīng)數(shù)值模擬方法研究[D]. 南京: 南京航空航天大學(xué), 2007.

    PAN Y. Numerical methods for hypersonic flowfield with magnetic interference[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2007 (in Chinese).

    [11] 孫曉輝. 新型沖壓推進(jìn)系統(tǒng)波系結(jié)構(gòu)及其MHD控制[D]. 南京: 南京理工大學(xué), 2013.

    SUN X H. Investigations of the wave structures of new concept RAM propulsion systems and their MHD control[D]. Nanjing: Nanjing University of Science and Technology, 2013 (in Chinese).

    [12] 張洪量. 高超聲速流場電磁干擾數(shù)值模擬研究[D]. 南京: 南京航空航天大學(xué), 2014.

    ZHANG H L. Numerical simulation for hypersonic flowfield with electromagnetic interference[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2014 (in Chinese).

    [13] 董維中. 熱化學(xué)非平衡效應(yīng)對高超聲速流動影響的數(shù)值計(jì)算與分析[D]. 北京: 北京航空航天大學(xué), 1996.

    DONG W Z. Numerical simulation and analysis of thermo-chemical non-equilibrium effects at hypersonic flows[D]. Beijing: Beihang University, 1996 (in Chinese).

    [14] 高鐵鎖, 李椿萱, 董維中, 等. 高超聲速電離繞流的數(shù)值模擬[J]. 空氣動力學(xué)學(xué)報(bào), 2002, 20(2): 184-191.

    GAO T S, LI C X, DONG W Z, et al. Numerical simulation of hypersonic ionized flows[J]. Acta Aerodynamica Sinica, 2002, 20(2): 184-191 (in Chinese).

    [15] 董維中, 高鐵鎖, 丁明松. 高超聲速非平衡流場多個振動溫度模型的數(shù)值研究[J]. 空氣動力學(xué)學(xué)報(bào), 2007, 25(1): 1-6.

    DONG W Z, GAO T S, DING M S. Numerical studies of the multiple vibrational temperature model in hypersonic non-equilibrium flows[J]. Acta Aerodynamica Sinica, 2007, 25(1): 1-6 (in Chinese).

    [16] 高鐵鎖,董維中,丁明松,等. 物理化學(xué)模型對高溫流場等離子體分布的影響[J]. 空氣動力學(xué)學(xué)報(bào), 2013, 31(5): 541-545.

    GAO T S, DONG W Z, DING M S, et al. The effects of physicochemical models on distribution of plasma in high-temperature flowfield[J]. Acta Aerodynamica Sinica, 2013, 31(5): 541-545 (in Chinese).

    [17] 董維中. 氣體模型對高超聲速再入鈍體氣動參數(shù)計(jì)算影響的研究[J]. 空氣動力學(xué)學(xué)報(bào), 2001, 19(2): 197-202.

    DONG W Z. Thermal and chemical model effect on the calculation of aerodynamic parameter for hypersonic reentry blunt body[J]. Acta Aerodynamica Sinica, 2001, 19(2): 197-202 (in Chinese).

    [18] SHERCLIFF J A. A text book of magnetohydrodynamics[M]. Oxford: Pergamon Press, 1995.

    [19] 陳剛, 張勁柏, 李椿萱. 磁流體流動控制中的磁場配置效率研究[J]. 力學(xué)學(xué)報(bào), 2008, 40(6): 752-759.

    CHEN G, ZHANG J B, LI C X. Efficiency analysis of magnetic field configuration in MHD flow control[J]. Chinese Journal of Theoretical and Applied Mechanics, 2008, 40(6): 752-759 (in Chinese).

    [20] DUNN M G, KANG S W. Theoretical and experimental studies of reentry plasmas: NASA-CR-2232[R]. Washington, D.C.: NASA, 1973.

    [21] CANDLER G V, MACCORMACK R W. The computation of hypersonic ionized flows in chemical and thermal nonequilibium: AIAA-1988-0511[R]. Reston, VA: AIAA, 1988.

    [22] BITYURIN V A, BOCHAROV A N. On efficiency of heat flux mitigation by the magnetic field in MHD re-entry flow: AIAA-2011-3463[R]. Reston, VA: AIAA, 2011.

    (責(zé)任編輯: 李明敏)

    URL:www.cnki.net/kcms/detail/11.1929.V.20170217.1121.002.html

    *Correspondingauthor.E-mail:dingms2008@qq.com

    Numericalsimulationof3DplasmaMHDaero-thermalenvironment

    DINGMingsong*,JIANGTao,DONGWeizhong,GAOTiesuo,LIUQingzong

    ComputationalAerodynamicsInstitute,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China

    Electromagneticflowcontroltechnique,asignificantmultidisciplinaryintersectingdirection,showswideapplicationprospectsinaerodynamiccharacteristicsoptimization,aerodynamicthermalenvironmentmitigation,boundarylayertransitionandplasmadistributionforflowcontroloverhypersonicvehicle.Inthispaper,chemicalreactions,molecularvibrationexcitationandchemicalnon-equilibriumeffectsareconsideredintheflowfieldofhypersonicvehicle,coupledwithelectromagneticfieldeffectandwiththeassumptionoflowmagneticReynoldsnumber.Bysolving3Dchemicalnon-equilibriumNavier-StokesequationsandMaxwellequations,numericalsimulationmethodandthecorrespondingcomputationalcodesaredevelopedforextramagneticfieldcoupledwithreentryplasmaflow,andarevalidatedbynumericallycalculatingtwotypicalexamples.Thesesimulationresultsareinagreementwiththoseinliteratures.Onthisbasis,theinfluenceofextramagneticfieldon3Dplasmaflowsandaero-thermalenvironmentunderdifferentflightconditionsisstudied.Theresultsshowthatextramagneticfieldcanobviouslychangethestandoffdistanceofshockwaveandreducethesurfaceheatfluxinmostsurfaceregions.Itisfoundthatthegreaterthemagneticfieldstrengthis,themoreobviousthemodificationeffectis.Theinfluencedegreeisrelevanttothefactorsofflightaltitude,velocityandextramagneticfieldconfiguration.Undercurrentcalculationconditions,theinfluencedegreeismoreobviouswhentheflightaltitudeishigher.

    MHD;plasma;chemicalnon-equilibrium;numericalsimulation;aero-thermalenvironment

    2016-12-06;Revised2016-12-30;Accepted2017-01-29;Publishedonline2017-02-171121

    2016-12-06;退修日期2016-12-30;錄用日期2017-01-29; < class="emphasis_bold">網(wǎng)絡(luò)出版時間

    時間:2017-02-171121

    www.cnki.net/kcms/detail/11.1929.V.20170217.1121.002.html

    .E-maildingms2008@qq.com

    丁明松, 江濤, 董維中, 等. 三維等離子體MHD氣動熱環(huán)境數(shù)值模擬J. 航空學(xué)報(bào),2017,38(8):121030.DINGMS,JIANGT,DONGWZ,etal.Numericalsimulationof3DplasmaMHDaero-thermalenvironmentJ.ActaAeronauticaetAstronauticaSinica,2017,38(8):121030.

    http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

    10.7527/S1000-6893.2017.121030

    V411.3; O354.7

    A

    1000-6893(2017)08-121030-10

    猜你喜歡
    磁流體駐點(diǎn)超聲速
    磁流體·吸引力
    中國寶玉石(2024年1期)2024-03-11 04:06:18
    高超聲速出版工程
    磁流體音箱
    高超聲速飛行器
    非均勻磁場下磁流體形態(tài)的研究
    電子制作(2019年9期)2019-05-30 09:42:16
    不可壓縮磁流體方程組在Besov空間中的爆破準(zhǔn)則
    基于游人游賞行為的留園駐點(diǎn)分布規(guī)律研究
    中國園林(2018年7期)2018-08-07 07:07:48
    超聲速旅行
    利用遠(yuǎn)教站點(diǎn),落實(shí)駐點(diǎn)干部帶學(xué)
    利用遠(yuǎn)教站點(diǎn),落實(shí)駐點(diǎn)干部帶學(xué)
    中国国产av一级| 亚洲欧洲国产日韩| 日韩国内少妇激情av| 秋霞伦理黄片| 特大巨黑吊av在线直播| 久久久久国产精品人妻一区二区| 成年人午夜在线观看视频| av卡一久久| av一本久久久久| 日韩三级伦理在线观看| 26uuu在线亚洲综合色| 午夜福利高清视频| 国产熟女欧美一区二区| av天堂中文字幕网| 伦理电影大哥的女人| 国产探花极品一区二区| 日韩,欧美,国产一区二区三区| 亚洲,欧美,日韩| 熟女av电影| 久久久精品免费免费高清| 国产成人精品一,二区| 欧美成人午夜免费资源| 亚洲精品国产av蜜桃| 精品99又大又爽又粗少妇毛片| 最近的中文字幕免费完整| 在线免费十八禁| 久久久久网色| 日韩免费高清中文字幕av| 18禁在线无遮挡免费观看视频| av卡一久久| 欧美xxxx黑人xx丫x性爽| 亚洲aⅴ乱码一区二区在线播放| 国产亚洲最大av| 日韩强制内射视频| 男女国产视频网站| 少妇的逼好多水| 2018国产大陆天天弄谢| 国产精品一及| 亚洲av.av天堂| 黄片wwwwww| 国产一区二区在线观看日韩| 我的老师免费观看完整版| 国产成人精品久久久久久| 在线看a的网站| 欧美xxⅹ黑人| 伦精品一区二区三区| 日日摸夜夜添夜夜添av毛片| a级一级毛片免费在线观看| 国产在线一区二区三区精| 韩国高清视频一区二区三区| 精品一区二区三区视频在线| 国产在视频线精品| 2021少妇久久久久久久久久久| 亚洲欧美清纯卡通| 亚洲av免费高清在线观看| 国产免费一区二区三区四区乱码| 一级毛片我不卡| 六月丁香七月| 好男人视频免费观看在线| 午夜福利视频1000在线观看| 最近手机中文字幕大全| 好男人在线观看高清免费视频| 国产精品久久久久久精品电影小说 | 91精品伊人久久大香线蕉| 久久99热这里只频精品6学生| 青春草亚洲视频在线观看| 嫩草影院精品99| 国产精品三级大全| 久久久久久国产a免费观看| 欧美少妇被猛烈插入视频| 激情 狠狠 欧美| 91久久精品国产一区二区三区| 国产综合精华液| 校园人妻丝袜中文字幕| av免费观看日本| 亚洲精品自拍成人| 只有这里有精品99| 午夜日本视频在线| av国产免费在线观看| 五月伊人婷婷丁香| 国产精品av视频在线免费观看| 国产一区二区在线观看日韩| 日本-黄色视频高清免费观看| 久久久成人免费电影| 一区二区三区精品91| 男人和女人高潮做爰伦理| 国内精品宾馆在线| 久久精品久久久久久久性| 欧美高清性xxxxhd video| 久久综合国产亚洲精品| 国产精品成人在线| 久久亚洲国产成人精品v| 亚洲精品一二三| 内射极品少妇av片p| 国产精品人妻久久久久久| 内射极品少妇av片p| 波野结衣二区三区在线| 国产精品一区二区在线观看99| 精品国产一区二区三区久久久樱花 | 国产黄频视频在线观看| 人妻夜夜爽99麻豆av| 日日啪夜夜撸| 国产一区二区在线观看日韩| 亚洲美女视频黄频| 免费不卡的大黄色大毛片视频在线观看| 国产男人的电影天堂91| 一级a做视频免费观看| 亚洲精品久久午夜乱码| 日本一二三区视频观看| 亚洲四区av| 51国产日韩欧美| 亚洲精品456在线播放app| 晚上一个人看的免费电影| 国产 一区 欧美 日韩| 中文资源天堂在线| 女人久久www免费人成看片| 亚洲国产高清在线一区二区三| 中文字幕av成人在线电影| 精品人妻视频免费看| 亚洲精品乱久久久久久| 亚洲国产欧美人成| 国产高清国产精品国产三级 | 日本猛色少妇xxxxx猛交久久| 国产亚洲91精品色在线| 亚洲精品乱久久久久久| 99九九线精品视频在线观看视频| 国产免费视频播放在线视频| 亚洲欧美中文字幕日韩二区| 国产亚洲精品久久久com| 国产美女午夜福利| 热re99久久精品国产66热6| 久久久久久久大尺度免费视频| 欧美成人精品欧美一级黄| 国产欧美日韩一区二区三区在线 | 日韩视频在线欧美| 80岁老熟妇乱子伦牲交| 伦理电影大哥的女人| 亚洲电影在线观看av| 国产永久视频网站| av在线观看视频网站免费| av国产免费在线观看| 五月玫瑰六月丁香| 国产精品精品国产色婷婷| 国模一区二区三区四区视频| 插阴视频在线观看视频| 黄色视频在线播放观看不卡| 国产视频内射| 亚洲美女搞黄在线观看| 中文在线观看免费www的网站| 日韩一区二区视频免费看| 日韩在线高清观看一区二区三区| 亚洲av中文字字幕乱码综合| 免费看光身美女| 欧美最新免费一区二区三区| 天堂中文最新版在线下载 | 国产精品av视频在线免费观看| 亚洲内射少妇av| 噜噜噜噜噜久久久久久91| av在线天堂中文字幕| 国产av国产精品国产| 高清欧美精品videossex| 国产 精品1| 欧美极品一区二区三区四区| 亚洲美女搞黄在线观看| 大又大粗又爽又黄少妇毛片口| 日韩成人伦理影院| 在线天堂最新版资源| 国产精品久久久久久久电影| 肉色欧美久久久久久久蜜桃 | 欧美日韩在线观看h| 国产一区二区三区av在线| 99热国产这里只有精品6| 777米奇影视久久| 久久99热这里只频精品6学生| 一级毛片我不卡| 国产精品99久久99久久久不卡 | 婷婷色麻豆天堂久久| 欧美日韩在线观看h| 亚洲av男天堂| 我要看日韩黄色一级片| 少妇猛男粗大的猛烈进出视频 | 白带黄色成豆腐渣| 纵有疾风起免费观看全集完整版| 成人黄色视频免费在线看| 少妇的逼好多水| 卡戴珊不雅视频在线播放| av女优亚洲男人天堂| 成人二区视频| 日韩电影二区| 久久久久国产网址| 久久久亚洲精品成人影院| 久久久久精品久久久久真实原创| 日本猛色少妇xxxxx猛交久久| 久久97久久精品| 亚洲精品国产av成人精品| 少妇人妻精品综合一区二区| 中文字幕亚洲精品专区| 成人午夜精彩视频在线观看| 亚洲国产精品999| 亚洲电影在线观看av| 神马国产精品三级电影在线观看| 国产伦精品一区二区三区视频9| 成人特级av手机在线观看| 午夜亚洲福利在线播放| 性色avwww在线观看| 国产亚洲一区二区精品| 欧美一级a爱片免费观看看| 九九在线视频观看精品| 亚洲真实伦在线观看| 日本猛色少妇xxxxx猛交久久| 69av精品久久久久久| 午夜爱爱视频在线播放| 伊人久久国产一区二区| 在线精品无人区一区二区三 | 在线观看av片永久免费下载| 日韩不卡一区二区三区视频在线| 国产黄a三级三级三级人| av在线播放精品| 乱系列少妇在线播放| 乱码一卡2卡4卡精品| 三级男女做爰猛烈吃奶摸视频| 亚洲精品日韩av片在线观看| 国产黄a三级三级三级人| 成人无遮挡网站| 麻豆乱淫一区二区| 国产成人精品一,二区| av在线亚洲专区| 亚洲精品视频女| 大码成人一级视频| 69av精品久久久久久| 久久精品久久久久久久性| 秋霞伦理黄片| 我的老师免费观看完整版| 久久久久久伊人网av| 91久久精品国产一区二区成人| 97精品久久久久久久久久精品| tube8黄色片| 亚洲av二区三区四区| av免费观看日本| 久久综合国产亚洲精品| 少妇 在线观看| 国产色婷婷99| 国产熟女欧美一区二区| 亚洲美女视频黄频| 18禁在线播放成人免费| 亚洲欧美一区二区三区国产| 精品99又大又爽又粗少妇毛片| 国产精品国产av在线观看| 欧美日韩综合久久久久久| 国产精品秋霞免费鲁丝片| 久久综合国产亚洲精品| 免费黄色在线免费观看| 久久久久久久亚洲中文字幕| 大陆偷拍与自拍| 日韩中字成人| 欧美日韩视频高清一区二区三区二| 丝袜脚勾引网站| 在线观看一区二区三区激情| 久久鲁丝午夜福利片| 亚洲精品乱码久久久久久按摩| 噜噜噜噜噜久久久久久91| 国产精品熟女久久久久浪| 久久久久久久久久久丰满| 夫妻性生交免费视频一级片| 哪个播放器可以免费观看大片| 久久影院123| 亚洲精品第二区| 久久精品熟女亚洲av麻豆精品| 99精国产麻豆久久婷婷| 看非洲黑人一级黄片| 欧美日韩视频高清一区二区三区二| 三级国产精品片| 人人妻人人爽人人添夜夜欢视频 | av免费在线看不卡| 久久99热这里只频精品6学生| 99热6这里只有精品| 精品久久久久久久人妻蜜臀av| 欧美精品一区二区大全| 亚洲一级一片aⅴ在线观看| 免费观看无遮挡的男女| 边亲边吃奶的免费视频| 少妇人妻 视频| 国精品久久久久久国模美| 99热6这里只有精品| 80岁老熟妇乱子伦牲交| 永久网站在线| 中国三级夫妇交换| 日韩欧美精品免费久久| 乱码一卡2卡4卡精品| 欧美日韩综合久久久久久| 日韩不卡一区二区三区视频在线| 欧美激情在线99| 超碰97精品在线观看| 国产伦理片在线播放av一区| 神马国产精品三级电影在线观看| 一级a做视频免费观看| 日本欧美国产在线视频| 老女人水多毛片| 亚洲色图av天堂| 免费观看a级毛片全部| 亚州av有码| 国语对白做爰xxxⅹ性视频网站| 久久精品国产a三级三级三级| www.av在线官网国产| 亚洲精品乱码久久久久久按摩| 一边亲一边摸免费视频| 久久久久网色| 国产日韩欧美在线精品| 久久久久久久久久成人| 国产午夜精品久久久久久一区二区三区| 国国产精品蜜臀av免费| 简卡轻食公司| 久久久精品94久久精品| 久久久色成人| 男女无遮挡免费网站观看| 国产精品一区二区三区四区免费观看| 欧美日韩在线观看h| 黄色怎么调成土黄色| 熟妇人妻不卡中文字幕| 精品国产一区二区三区久久久樱花 | 啦啦啦中文免费视频观看日本| 啦啦啦啦在线视频资源| 日韩成人伦理影院| av在线观看视频网站免费| 国产精品爽爽va在线观看网站| 国产乱来视频区| 在线精品无人区一区二区三 | 亚洲丝袜综合中文字幕| 99久久九九国产精品国产免费| 你懂的网址亚洲精品在线观看| 色吧在线观看| 中文乱码字字幕精品一区二区三区| 久久精品熟女亚洲av麻豆精品| 欧美精品人与动牲交sv欧美| 夫妻午夜视频| 国产乱来视频区| 伦精品一区二区三区| 亚洲天堂av无毛| 18禁在线无遮挡免费观看视频| 久久精品国产a三级三级三级| 三级国产精品欧美在线观看| 欧美国产精品一级二级三级 | 蜜臀久久99精品久久宅男| av在线天堂中文字幕| 美女高潮的动态| 黄色配什么色好看| 日本wwww免费看| 六月丁香七月| 成人美女网站在线观看视频| 久久热精品热| 最近中文字幕高清免费大全6| 亚洲国产精品国产精品| 免费av不卡在线播放| 亚洲av男天堂| 最近的中文字幕免费完整| 麻豆成人av视频| 中文欧美无线码| 色播亚洲综合网| 亚洲美女搞黄在线观看| 欧美三级亚洲精品| 丰满乱子伦码专区| 建设人人有责人人尽责人人享有的 | 国产亚洲91精品色在线| 亚洲欧美一区二区三区国产| 日韩在线高清观看一区二区三区| 国产精品无大码| 国产av国产精品国产| 精品人妻一区二区三区麻豆| 永久网站在线| 深爱激情五月婷婷| 久久精品国产亚洲av天美| 大话2 男鬼变身卡| 夫妻午夜视频| 建设人人有责人人尽责人人享有的 | 男女边吃奶边做爰视频| 少妇人妻精品综合一区二区| 亚洲最大成人中文| 中国三级夫妇交换| 狂野欧美白嫩少妇大欣赏| 国产一区二区三区综合在线观看 | 精品视频人人做人人爽| av播播在线观看一区| 国模一区二区三区四区视频| 国产色婷婷99| 午夜福利在线在线| 男男h啪啪无遮挡| 日韩欧美一区视频在线观看 | 久久久久久久大尺度免费视频| 欧美一级a爱片免费观看看| 欧美日韩综合久久久久久| 国产一区二区在线观看日韩| av在线天堂中文字幕| 国产91av在线免费观看| 在线天堂最新版资源| 色播亚洲综合网| 久久精品国产鲁丝片午夜精品| 男女边摸边吃奶| 国产av国产精品国产| 国产高清三级在线| 国内少妇人妻偷人精品xxx网站| 麻豆国产97在线/欧美| 少妇熟女欧美另类| 久久精品国产亚洲av涩爱| 国内揄拍国产精品人妻在线| 啦啦啦中文免费视频观看日本| 尤物成人国产欧美一区二区三区| 寂寞人妻少妇视频99o| 欧美日韩视频高清一区二区三区二| 久久久久性生活片| 人人妻人人爽人人添夜夜欢视频 | 亚洲欧美日韩卡通动漫| 亚洲久久久久久中文字幕| 九九在线视频观看精品| 欧美性猛交╳xxx乱大交人| 别揉我奶头 嗯啊视频| 天天躁夜夜躁狠狠久久av| 久久综合国产亚洲精品| 真实男女啪啪啪动态图| 男女边摸边吃奶| av卡一久久| 亚洲丝袜综合中文字幕| 亚洲成人中文字幕在线播放| 免费观看的影片在线观看| 日韩成人伦理影院| 97超视频在线观看视频| 久久久久精品性色| 亚洲欧美成人精品一区二区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 成年免费大片在线观看| 亚洲国产日韩一区二区| 777米奇影视久久| 久热久热在线精品观看| 全区人妻精品视频| 国产成人福利小说| 美女高潮的动态| 边亲边吃奶的免费视频| 亚洲av电影在线观看一区二区三区 | 精品99又大又爽又粗少妇毛片| 免费观看av网站的网址| 舔av片在线| 男人狂女人下面高潮的视频| 亚洲欧美精品自产自拍| 午夜免费鲁丝| 欧美一级a爱片免费观看看| 亚洲精品,欧美精品| videossex国产| 18+在线观看网站| 九九久久精品国产亚洲av麻豆| av女优亚洲男人天堂| 亚洲在久久综合| 黑人高潮一二区| 欧美成人午夜免费资源| 内射极品少妇av片p| 少妇的逼好多水| 日日啪夜夜爽| 日韩欧美精品v在线| 在线 av 中文字幕| 亚洲av二区三区四区| 偷拍熟女少妇极品色| 午夜免费观看性视频| 美女国产视频在线观看| 超碰97精品在线观看| 亚洲精品aⅴ在线观看| 亚洲国产精品成人综合色| 亚洲人成网站在线观看播放| 99久久中文字幕三级久久日本| 91精品国产九色| 久久久亚洲精品成人影院| 插阴视频在线观看视频| 人妻制服诱惑在线中文字幕| 麻豆国产97在线/欧美| 精品亚洲乱码少妇综合久久| 日韩欧美精品v在线| 亚洲人与动物交配视频| 亚洲国产高清在线一区二区三| 国产免费视频播放在线视频| 另类亚洲欧美激情| 亚洲国产精品专区欧美| 婷婷色av中文字幕| 日韩欧美 国产精品| 精品久久久精品久久久| 人妻一区二区av| 色哟哟·www| 真实男女啪啪啪动态图| 日韩精品有码人妻一区| 99热这里只有是精品在线观看| 纵有疾风起免费观看全集完整版| 亚洲高清免费不卡视频| 麻豆成人av视频| 欧美另类一区| 欧美一级a爱片免费观看看| 国产精品久久久久久精品电影| 特大巨黑吊av在线直播| av在线观看视频网站免费| 成年人午夜在线观看视频| av天堂中文字幕网| 国产免费视频播放在线视频| 国产成人午夜福利电影在线观看| 777米奇影视久久| videos熟女内射| 亚洲伊人久久精品综合| 国产精品伦人一区二区| 国产黄色免费在线视频| 在线免费十八禁| 亚洲第一区二区三区不卡| 大片免费播放器 马上看| 成人高潮视频无遮挡免费网站| 国产欧美另类精品又又久久亚洲欧美| 九九爱精品视频在线观看| av免费观看日本| 又爽又黄无遮挡网站| 人妻夜夜爽99麻豆av| 日本一二三区视频观看| 大陆偷拍与自拍| 最近中文字幕高清免费大全6| 国产午夜精品一二区理论片| 大片电影免费在线观看免费| 又爽又黄a免费视频| 精品久久久久久电影网| 中文字幕亚洲精品专区| 亚洲欧美精品专区久久| 国产乱来视频区| 22中文网久久字幕| 午夜爱爱视频在线播放| 麻豆久久精品国产亚洲av| 18禁裸乳无遮挡免费网站照片| 看十八女毛片水多多多| 欧美zozozo另类| 免费大片黄手机在线观看| 校园人妻丝袜中文字幕| 在线精品无人区一区二区三 | 搡老乐熟女国产| 赤兔流量卡办理| 欧美另类一区| 欧美变态另类bdsm刘玥| 精品久久久久久久久av| 夜夜爽夜夜爽视频| av又黄又爽大尺度在线免费看| 别揉我奶头 嗯啊视频| 国产精品秋霞免费鲁丝片| 在线亚洲精品国产二区图片欧美 | 久久鲁丝午夜福利片| 日本色播在线视频| 日本猛色少妇xxxxx猛交久久| 少妇人妻精品综合一区二区| 十八禁网站网址无遮挡 | 在线观看美女被高潮喷水网站| 成人美女网站在线观看视频| 免费观看a级毛片全部| 免费看av在线观看网站| 亚洲自拍偷在线| 亚洲精品中文字幕在线视频 | 国产免费视频播放在线视频| 国产精品一区二区三区四区免费观看| 一二三四中文在线观看免费高清| 日韩,欧美,国产一区二区三区| 久久97久久精品| 亚洲色图av天堂| 欧美成人一区二区免费高清观看| 欧美区成人在线视频| 男插女下体视频免费在线播放| 免费观看av网站的网址| 在现免费观看毛片| 国产免费福利视频在线观看| 国产老妇伦熟女老妇高清| 日韩欧美 国产精品| 久久久久国产精品人妻一区二区| 神马国产精品三级电影在线观看| 建设人人有责人人尽责人人享有的 | 国产视频内射| 汤姆久久久久久久影院中文字幕| 啦啦啦啦在线视频资源| 美女主播在线视频| 别揉我奶头 嗯啊视频| 成人国产av品久久久| 中文字幕av成人在线电影| 久久久久性生活片| 免费观看a级毛片全部| 国产久久久一区二区三区| 婷婷色麻豆天堂久久| 久久久久性生活片| 街头女战士在线观看网站| 国产69精品久久久久777片| 国产黄频视频在线观看| 国产欧美日韩一区二区三区在线 | 亚洲精品色激情综合| 国产欧美日韩精品一区二区| 2018国产大陆天天弄谢| 午夜视频国产福利| 在线观看国产h片| 国产探花极品一区二区| 亚洲高清免费不卡视频| 夜夜看夜夜爽夜夜摸| 狂野欧美激情性bbbbbb| 久久99精品国语久久久| 一级av片app| 亚洲美女视频黄频| 观看免费一级毛片| 日韩av免费高清视频| 国产午夜精品一二区理论片| 青青草视频在线视频观看| 成人免费观看视频高清| 男插女下体视频免费在线播放| 日本黄大片高清| 久久久久国产精品人妻一区二区| 亚洲欧美精品专区久久| 亚洲美女搞黄在线观看| 国产欧美另类精品又又久久亚洲欧美| 成人欧美大片| 超碰av人人做人人爽久久| 亚洲成色77777| 热99国产精品久久久久久7| 亚洲激情五月婷婷啪啪| 国产精品伦人一区二区| 久久久成人免费电影|