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

    亞大氣壓六相交流電弧放電等離子體射流特性數(shù)值模擬?

    2018-03-27 06:12:30郭恒張曉寧聶秋月李和平曾實李志輝
    物理學報 2018年5期
    關鍵詞:射流等離子體飛行器

    郭恒 張曉寧 聶秋月 李和平 曾實 李志輝

    1)(清華大學工程物理系,北京 100084)

    2)(哈爾濱工業(yè)大學空間基礎科學研究中心,哈爾濱 150001)

    3)(哈爾濱工業(yè)大學電氣工程及自動化學院,哈爾濱 150001)

    4)(中國空氣動力研究與發(fā)展中心超高速所,綿陽 621000)

    5)(國家計算流體力學實驗室,北京 100191)

    (2017年11月29日收到;2018年1月4日收到修改稿)

    1 引 言

    在距離海平面20—100 km的空間區(qū)域,由于稀薄的大氣環(huán)境有效地減小了飛行器在飛行過程中所受的空氣阻力,從而較易實現(xiàn)高超聲速飛行,成為人類航空與航天活動的重要區(qū)域之一,具有廣闊的商業(yè)和軍事應用價值,日益受到各國研究者的廣泛關注.臨近空間飛行器在高超聲速飛行以及飛行器再入大氣過程中,其表面激波熱致產(chǎn)生高密度、強碰撞、非均勻以及熱力學和化學非平衡態(tài)等離子體[1,2],從而使得飛行器面臨著復雜的介質(zhì)環(huán)境.這一方面表現(xiàn)為能量的注入和等離子體的存在所導致的飛行器周圍局部區(qū)域當?shù)芈曀僭黾?、邊界層特性的顯著變化以及由此所引起的飛行器飛行阻力的變化[3];另一方面則表現(xiàn)為包覆于飛行器表面的等離子體鞘套環(huán)境對外來電磁波的吸收、散射和折射,引發(fā)“黑障”效應[4?6],進而影響飛行器的通信和控制系統(tǒng).因此,開展高超聲速飛行器與其周圍復雜介質(zhì)環(huán)境間相互作用機制的理論和實驗研究對于指導飛行器的設計具有重要的意義.

    從該領域目前的研究方法來看,主要包括了飛行實驗、地面風洞實驗和數(shù)值模擬.其中飛行實驗[7?11]和地面風洞實驗[12?14]盡管運行成本昂貴,且所能獲得的實驗數(shù)據(jù)有限,但卻為科學研究積累了寶貴的實驗數(shù)據(jù).相比而言,采用螺旋波等離子體源、電弧等離子體源、輝光放電、電感/電容耦合放電等離子體源等[15?17]在地面開展小規(guī)模的實驗室模擬以及借助于高性能電子計算機開展數(shù)值模擬研究,則是相對經(jīng)濟高效的科學研究方法.這樣的研究方法盡管無法完全再現(xiàn)臨近空間飛行器所面臨的復雜介質(zhì)環(huán)境,但從模擬飛行器與其周圍復雜介質(zhì)環(huán)境在某些方面的相互作用過程,如“黑障”問題中等離子體-波相互作用過程的研究以及高超聲速飛行器在飛行過程中所面臨的氣動力和氣動熱環(huán)境來講依然是十分有意義的.

    以研究臨近空間高超聲速飛行器以及航天器再入大氣環(huán)境飛行過程中“黑障”問題為背景,一方面,需要在地面模擬實驗中產(chǎn)生具有足夠大體積(如有效實驗區(qū)長度L3=30—50 cm)、足夠厚(如鞘套厚度δ=5—10 cm)以及足夠高的等離子體密度(如電子數(shù)密度ne=108—1012cm?3)的包覆飛行器表面的等離子體鞘套環(huán)境(如圖1所示),進而能夠研究等離子體鞘套特性及其與入射電磁波間的相互作用規(guī)律;另一方面,飛行器飛行過程中所面臨的高超聲速流動環(huán)境將與飛行器壁面附近所產(chǎn)生的等離子體強烈耦合,從而呈現(xiàn)出氣動力、氣動熱與非平衡態(tài)等離子體間復雜的耦合效應.在文獻[15]中,作者著重討論了將雙射流直流電弧放電等離子體發(fā)生器的設計理念[16]與多相交流電弧放電等離子體發(fā)生器的設計理念[17,18]相結合,所建立的能夠在大氣壓和亞大氣壓條件下放電產(chǎn)生大體積等離子體射流的六相交流電弧放電等離子體實驗平臺(multiphase gas discharge plasma experimental platform-2015,MPX-2015),并給出了不同放電參數(shù)(包括相鄰電極之間的放電間距、弧電流、等離子體工作氣體流量、真空腔壓強以及水冷約束管直徑等參數(shù))下等離子體射流區(qū)氣體溫度、電子溫度、電子數(shù)密度以及射流長度和直徑等參數(shù)的變化規(guī)律.本文以MPX-2015為研究對象,著重采用數(shù)值模擬的手段研究非平衡態(tài)等離子體的特性.

    圖1 包覆飛行器的等離子體鞘套示意圖Fig.1.Schematic of the plasma layer surrounding an aerocraft.

    開展臨近空間高超聲速飛行器與其周圍復雜介質(zhì)環(huán)境間傳熱、流動以及電磁相互作用的研究,需要建立起一套能夠自洽描述包覆飛行器表面的、具有可壓縮、非平衡、熱-電-磁-流動-化學反應強耦合特性的等離子體環(huán)境的物理數(shù)學模型,包括基本控制方程、輸運系數(shù)和定解條件以及高效率的數(shù)值求解方法.在過去的幾十年間,一方面,以計算流體力學模型為基礎,國內(nèi)外眾多研究者針對高超聲速飛行器飛行過程中的流動特性、熱力學-化學非平衡特性開展了大量的研究工作,極大地推動了高超聲速計算流體力學的發(fā)展[19?23];另一方面,在氣動力/熱與作為物質(zhì)第四態(tài)的等離子體間復雜相互作用的研究中,盡管已經(jīng)開展了不少有關描述非平衡態(tài)等離子體輸運特性的物理數(shù)學模型(主要包括基本控制方程和輸運系數(shù))的研究(如文獻[24—29]),但目前不同研究者所采用的物理數(shù)學模型依然不盡相同,特別是在基本控制方程與等離子體輸運系數(shù)的定義和計算方法方面有時還不完全自洽,這在很大程度上影響了數(shù)值模擬所預測的臨近空間飛行器飛行特性的準確性.因此,需要針對非平衡態(tài)等離子體中復雜的物理化學過程,建立和完善自洽描述非平衡態(tài)等離子體與環(huán)境氣體間質(zhì)量、動量和能量協(xié)同輸運特性的物理數(shù)學模型,包括基本控制方程、初邊值條件和輸運性質(zhì)參數(shù)數(shù)據(jù)庫等.本研究組在前期的研究中應用漸進擾動的方法求解Boltzmann方程得到了描述粒子輸運特性的守恒方程,進一步通過粒子間的相互作用計算得到了守恒方程中的輸運系數(shù)[30,31],從而保證了等離子體中組分輸運方程與其輸運系數(shù)的自洽性,建立了自洽描述非平衡態(tài)等離子體輸運過程的物理數(shù)學模型和以氬等離子體為模式體系的非平衡態(tài)等離子體輸運系數(shù)的數(shù)據(jù)庫[32].本文在此基礎上針對MPX-2015實驗系統(tǒng)開展非平衡態(tài)氬等離子體射流特性及其與環(huán)境冷氣體和下游鈍體間相互作用特性的數(shù)值模擬研究.其中,MPX-2015實驗系統(tǒng)的6個電極均勻布置在同一圓周上,等離子體工作氣體除了沿6個電極表面注入外,還可從電極平面上方注入;放電過程中,相鄰兩個電極間所加交流電壓的相位差為60°;每相電流的波形均為標準正弦波,當其中1相的電流過零時,其他5相的電流并不會同時經(jīng)過零點,從而保證了電弧的持續(xù)和一定程度上的穩(wěn)定;每相的電流在70—110 A的范圍內(nèi)變化,兩相之間的工作電壓則維持在24 V左右;氣體放電所產(chǎn)生的等離子體經(jīng)過水冷約束管后進入真空腔.實驗過程中可以采用發(fā)射光譜法測量等離子體射流區(qū)的氣體溫度、電子溫度和電子數(shù)密度等參數(shù).有關實驗平臺的具體描述可參見文獻[15].本文第2部分討論非平衡態(tài)等離子體特性數(shù)值模擬研究的物理數(shù)學模型;第3部分給出亞音速條件下自由射流特性計算結果以及與MPX-2015平臺實驗測量數(shù)據(jù)的對比,從而在一定程度上驗證模型的正確性;第4部分給出從亞音速到超音速的非平衡態(tài)等離子體射流沖擊鈍體的數(shù)值模擬結果;第5部分給出本論文研究得到的主要結論,簡要討論后續(xù)研究需要解決的主要問題.

    2 物理數(shù)學模型

    2.1 基本控制方程

    針對亞大氣壓條件下的非平衡態(tài)氬等離子體開展數(shù)值模擬研究,所采用的基本假定包括[30,31,33]:1)作為初步的研究工作,僅就六相交流放電等離子體射流的準穩(wěn)態(tài)特性開展數(shù)值模擬研究,相應的實驗對比亦采用時間平均的測量結果;2)盡管由于在物理上電子質(zhì)量遠遠小于重粒子質(zhì)量,從而導致電子與重粒子間彈性碰撞的能量交換效率極低,但由于整個等離子體體系中化學反應過程(如激發(fā)、電離及其相應的逆過程等)的存在,電子子體系和重粒子子體系并非兩個孤立的子系統(tǒng),而是相互間存在質(zhì)量、動量和能量的交換過程;3)等離子體中粒子的自弛豫時間遠小于能量交換的弛豫時間,于是,等離子體中電子和重粒子兩個子體系的速度分布函數(shù)分別滿足Maxwell分布,可分別定義電子溫度(Te)和重粒子溫度(Th);4)粒子間動量交換的弛豫時間遠小于能量交換的弛豫時間,于是,宏觀上可以將等離子體看作一種流體,流體速度為體系中所有粒子的質(zhì)量平均速度;5)粒子的自弛豫時間遠小于體系化學反應的弛豫時間,即等離子體的輸運性質(zhì)參數(shù)將不受化學反應的影響;6)i等離子體為光學薄,即等離子體對輻射能的吸收與全波長范圍內(nèi)的輻射能相比可忽略不計;7)等離子體工作氣體及環(huán)境氣體均為氬氣,且處于層流流動狀態(tài);8)由于等離子體發(fā)生器水冷約束管的型線結構以及發(fā)生器進氣口處的壓強和真空腔內(nèi)的環(huán)境壓強均可調(diào),從而使得等離子體射流可處于亞音速或超音速流動狀態(tài);9)不考慮鈍體表面材料的燒蝕對其附近等離子體射流性質(zhì),特別是輻射特性的影響.

    基于上述基本假定,在二維圓柱坐標系(z,r)下需要求解的基本控制方程包括連續(xù)性方程、組分質(zhì)量平均的動量守恒方程以及重粒子與電子能量守恒方程,即:

    在上述方程中,p,v,vr,vz,Th,Te,θ,jr,jz,Bθ,E,φ分別為等離子體的壓力、速度矢量和其在徑向、軸向的分量、重粒子溫度、電子溫度、電子溫度與重粒子溫度比(即Te/Th)、電流密度在徑向和軸向的分量、自感應磁場的切向分量、電場強度矢量和電位;分別為氬等離子體的質(zhì)量密度、黏性系數(shù)、重粒子平動熱導率和非平衡熱導率、電子平動熱導率和非平衡熱導率以及電導率.在局域化學平衡條件下,上述等離子體輸運性質(zhì)參數(shù)均為電子溫度、重粒子溫度以及壓力的函數(shù),可根據(jù)文獻[30,31]中的定義與求解方法計算得到;na,ni,ne分別為原子、離子和電子數(shù)密度;τh:?v為黏性耗散項;Urad和分別為等離子體單位體積的輻射熱損失和電子-重粒子間彈性碰撞能量交換項;分別表示電子-原子以及電子-離子間的動量傳輸碰撞頻率;me和mh分別為電子和重粒子的質(zhì)量.由于在本文的數(shù)值模擬中僅對等離子體射流區(qū)特性進行計算,因此無需求解電磁場,即電位方程(7)無需求解;相應地,方程(2)—(5)中與電磁場有關的源項亦均忽略.

    2.2 氬等離子體輸運性質(zhì)參數(shù)

    開展非平衡態(tài)等離子體特性數(shù)值模擬的另一前提條件是與上述2.1節(jié)基本控制方程自洽的等離子體輸運性質(zhì)參數(shù)數(shù)據(jù)庫的建立.數(shù)值模擬采用氬氣為等離子體工作氣體,基于文獻[30,34]的方法開發(fā)了非平衡態(tài)氬等離子體組分、黏性系數(shù)、電子和重粒子平動熱導率、組分擴散系數(shù)等輸運性質(zhì)參數(shù)計算的程序代碼.在局域化學平衡條件下,上述等離子體輸運性質(zhì)參數(shù)均為電子溫度、重粒子溫度以及壓力的函數(shù),通過與求解2.1節(jié)偏微分方程組[方程(1)—(5)]的計算程序進行耦合迭代求解,最終獲得不同工況下非平衡態(tài)氬等離子體射流特性的計算結果.

    圖2 氬等離子體重粒子平動熱導率(λh,實線)、電子平動熱導率(λe,虛線)和黏性系數(shù)(μ,點線)隨電子溫度和環(huán)境壓強的變化規(guī)律(θ=2)Fig.2.Variations of the translational thermal conductivity of heavy particles(λh,solid line)and that of electrons(λe,dashed line),and viscosity(μ,dotted line)with the electron temperature and pressure(θ=2).

    圖2給出了環(huán)境壓強p=100,1000和10000 Pa,電子與重粒子溫度比θ=2的條件下氬等離子體黏性系數(shù)(μ)、電子平動熱導率(λe)和重粒子平動熱導率(λh)隨電子溫度(Te)的變化規(guī)律.為了便于比較,圖2中上述3個輸運系數(shù)均進行了歸一化處理,歸一化處理的參考值分別為10000 Pa下的μmax(0.0015 Pa·s@9200 K),λe(8.1600 W/(m·K)@30000 K)和λh(0.0114 W/(m·K)@9000 K). 從圖2可以看到,在環(huán)境壓強(p)和溫度比(θ)一定的條件下,隨著電子溫度的升高,1)電子數(shù)密度和熱運動速度均增大,從而導致電子平動熱導率的單調(diào)增大;2)黏性系數(shù)和重粒子平動熱導率的變化則相對比較復雜,在300—8000 K的溫度區(qū)間內(nèi),隨著電子溫度的升高,二者單調(diào)增大,之后隨著氣體電離反應的增強,帶電粒子間的長程的庫侖作用勢越來越重要,二者迅速降低,在15000 K后其值變得很小,并基本保持不變[32].在電子溫度(Te)及溫度比(θ)保持不變的情況下,環(huán)境壓強(p)亦對等離子體的輸運性質(zhì)參數(shù)有著顯著的影響,隨著壓力從100 Pa增加到10000 Pa,等離子體的黏性系數(shù)、電子和重粒子平動熱導率均增大.此處需要說明的是,限于篇幅,本小節(jié)僅給出了非常有限的參數(shù)條件下的氬等離子體黏性系數(shù)以及電子和重粒子平動熱導率隨電子溫度的變化曲線,而事實上,數(shù)值模擬需要廣闊溫度范圍和壓力范圍內(nèi)完備的等離子體輸運性質(zhì)參數(shù).有關非平衡態(tài)等離子體輸運性質(zhì)參數(shù)的理論推導和計算方法可參見文獻[34].

    3 亞音速氬等離子體自由射流非平衡特性

    3.1 計算域和邊界條件的確定

    為了驗證本文數(shù)值模擬所建立的物理數(shù)學模型的正確性,首先在亞音速條件下對氬等離子體自由射流的非平衡特性進行了二維數(shù)值模擬研究,并與文獻[15]中所得到的相同工況下的實驗測量結果進行了對比分析.本節(jié)數(shù)值模擬采用第2節(jié)所給出的基本控制方程和輸運性質(zhì)參數(shù),不同之處僅是由于較低的流速,忽略了重粒子能量方程中與可壓縮效應相關的項,即方程(4)中的τh:?v項.氬等離子體自由射流數(shù)值模擬所采用的計算域如圖3所示,相應的邊界條件列于表1.需要說明的是:本節(jié)數(shù)值模擬以文獻[15]3.1節(jié)中典型的等離子體工作參數(shù)下(氬氣總流量Qtot=16.0 slpm(standard liter per minute),其中從每個電極環(huán)縫注入的氣體流量為Qe=1.0 slpm,電極平面上游的冷氣體注入量為Qc=10.0 slpm;單相弧電流峰值I=80 A,相鄰兩電極間距l(xiāng)e=40 cm,水冷約束管出口內(nèi)徑r1=2.5 cm,真空腔壓強p=1000 Pa)測量得到的水冷約束管出口軸線上時間平均的等離子體電子溫度(8000 K)和重粒子溫度(4000 K)作為射流入口的溫度邊界條件,而速度邊界條件則在保證總工作氣體質(zhì)量流量(對應于標準狀況下的體積流量為Qtot=6×Qe+Qc=16.0 slpm)守恒的條件下折算得到入口處的速度邊界條件,即vz=ρs·Qtot/(ρin·Sin)≈180 m/s,其中ρs和ρin分別代表標準狀況下的氣體密度和入口溫度條件下的等離子體密度,Sin為水冷約束管噴嘴出口截面積.

    圖3 亞音速自由射流數(shù)值模擬計算域示意圖Fig.3.Schematic of the calculation domain for the modeling of subsonic free plasma jet.

    表1 亞音速自由射流數(shù)值模擬用邊界條件Table 1.Boundary conditions used for the modeling of subsonic free plasma jet.

    3.2 數(shù)值模擬結果與實驗測量結果的對比分析

    圖4給出了Qe=1.0 slpm,Qc=10.0 slpm,I=80 A,le=40 mm,r1=2.5 cm,p=1000 Pa的典型工況下二維數(shù)值模擬得到的等離子體射流區(qū)氣體溫度分布與等離子體射流實驗圖像的對比.可以看到,數(shù)值模擬所得到的等離子體射流區(qū)重粒子溫度的分布與實驗拍攝到的射流強度定性一致.另外,計算得到的以電子數(shù)密度ne=1010cm?3的等值線度量的等離子體射流直徑和長度分別為5.0 cm和65.0 cm;而實驗測量得到的等離子體射流直徑和長度則分別為5.5 cm和60.0 cm,二者亦符合較好.

    圖4 計算得到的等離子體射流區(qū)氣體溫度分布與實驗照片的比較Fig.4.Comparison between the calculated gas temperature distribution and the discharge image.

    圖5則比較了相同工況下計算和實驗得到的沿發(fā)生器幾何軸線等離子體射流區(qū)重粒子溫度(Th)、電子溫度(Te)、電子數(shù)密度(ne)以及電子-中性粒子碰撞頻率的變化規(guī)律.可以看到:1)上述四個參數(shù)的數(shù)值模擬結果與實驗測量結果在定性的變化趨勢上是完全一致的;2)由實驗測量得到的在z=50.0 cm處的電子溫度和重粒子溫度均較數(shù)值模擬結果分別低22%和26%,而實驗測得的電子數(shù)密度則與計算結果相差約18%,計算得到的電子-中性粒子碰撞頻率則與實驗值相差10%.但從總體上來看,本節(jié)二維非平衡數(shù)值模擬結果與實驗測量結果較好地定性符合,這在一定程度上驗證了本文所建立的物理數(shù)學模型用于亞音速流動條件下等離子體特性模擬的正確性.在此基礎上,考慮可壓縮性的影響,在第4部分開展跨音速和超音速流動條件下等離子體射流特性的數(shù)值模擬,這將為我們進行超音速等離子體射流特性的實驗研究提供理論上的指導,而有關數(shù)值模擬和實驗測量結果的對比分析,則還需要后續(xù)開展進一步的實驗測量.

    圖5 計算得到的射流區(qū)等離子體參數(shù)沿其幾何軸線的變化規(guī)律及其與實驗測量結果的比較Fig.5.Comparisons of the calculated and measured plasma parameters along the jet axis.

    4 超音速氬等離子體沖擊射流非平衡特性

    4.1 計算域和邊界條件的確定

    為了更加真實地模擬臨近空間包覆高超聲速飛行器的等離子體鞘套環(huán)境,不僅需要產(chǎn)生具有足夠大尺寸和足夠高電子數(shù)密度的等離子體,而且需要等離子體處于超音速流動狀態(tài)下.在本節(jié)的研究中,針對MPX-2015實驗系統(tǒng)首先從理論上對產(chǎn)生特定尺寸、超聲速等離子體射流環(huán)境進行了數(shù)值模擬預測.為此,將文獻[15]中約束電弧等離子體的水冷約束管內(nèi)部型線從直筒型改為縮放噴管結構,通過對高溫電弧等離子體流動的約束以及調(diào)整發(fā)生器進氣口處與水冷約束管外真空腔邊界的壓差Δp實現(xiàn)等離子體流動從亞音速向超音速的跨聲速流動.本節(jié)數(shù)值模擬所選取的計算域如圖6所示,在放電區(qū)所產(chǎn)生的電弧等離子體經(jīng)過縮放型噴管的入口(OA)后進行噴管內(nèi)部,并從其出口(BH)流出后沖擊其下游的鈍體,進而在鈍體表面(EFG)形成包覆鈍體的等離子體鞘套.在本節(jié)數(shù)值模擬中,由于希望達到的最大馬赫數(shù)為3.5左右,根據(jù)理想工作條件下噴管出口與喉部截面積之比與出口馬赫數(shù)的關系[35],可以得到所采用的縮放型噴管的出口半徑與其喉部半徑之比應為2.0.等離子體發(fā)生器噴管下游鈍體的尺寸根據(jù)RAM實際尺寸確定[36],其表面型線的表達式為

    其中x和y分別代表距離坐標原點O的軸向距離和徑向距離.圖6中紅色虛弧線代表了距鈍體頭部垂直距離為L的空間位置,而點K則為軸線上距離鈍體頭部5 cm處的位置.

    本節(jié)數(shù)值模擬所采用的邊界條件列于表2,其中邊界CD和DE給定出口壓力值,OA邊界給定入口壓力值,邊界ABC和EFG則采用固壁邊界條件.計算中參考壓力面選在計算域入口處,即OA邊界,并將壓力值固定為5000 Pa,因此計算域的入口與出口之間的壓力差為Δp,計算中通過改變Δp研究其對射流特性的影響.

    圖6 等離子體沖擊射流數(shù)值模擬計算域示意圖Fig.6.Schematic of the calculation domain for the modeling of plasma impinging jet.

    圖7給出了縮放型噴管進出口壓差Δp=900,2700,3500和4500 Pa條件下二維非平衡數(shù)值模擬計算得到的等離子體射流區(qū)馬赫數(shù)(Ma)、電子數(shù)密度(ne),電子溫度(Te)和重粒子溫度(Th)的空間分布,圖中的虛線表示距離鈍體表面5 cm的空間位置.從圖7可以看到:1)隨著Δp的增大,等離子體流經(jīng)縮放型噴管后逐步從亞音速流動發(fā)展成為超音速射流,噴管喉部下游最大馬赫數(shù)從0.56增大到3.00;2)隨著Δp的增大,包覆鈍體的等離子體參數(shù)分布的空間均勻性顯著降低,這一方面表現(xiàn)為鈍體壁面附近高于臨界電子數(shù)密度(ncrt=1010cm?3)的區(qū)域沿射流幾何軸線明顯縮小;而另一方面,Δp的增大亦使得在流動過程中射流區(qū)等離子體的內(nèi)能與動能間隨著空間位置不斷發(fā)生轉(zhuǎn)化,從而使得在鈍體前緣附近隨著動能向內(nèi)能轉(zhuǎn)化而導致電子數(shù)密度的顯著增高(如圖8所示),而在鈍體后緣處則呈現(xiàn)出高電子密度區(qū)域收縮的變化規(guī)律.也就是說,隨著噴管進出口壓差的增大,鈍體表面附近的電子數(shù)密度平均值增大,高電子數(shù)密度區(qū)域更加向鈍體頭部集中.

    圖7 不同Δp下等離子體射流區(qū)馬赫數(shù)、電子數(shù)密度、電子溫度和重粒子溫度的等值線分布 (a)900 Pa;(b)2700 Pa;(c)3500 Pa;(d)4500 PaFig.7.Distributions of Mach number,electron number density,electron temperature and heavy-particle temperature in the plasma jet region with different values of Δp:(a)900 Pa;(b)2700 Pa;(c)3500 Pa;(d)4500 Pa.

    表2 跨音速沖擊射流數(shù)值模擬用邊界條件Table 2.Boundary conditions used for the modeling of plasma impinging jet.

    圖8 不同Δp下鈍體上游等離子體射流區(qū)馬赫數(shù)、電子數(shù)密度、電子溫度和重粒子溫度沿軸線的變化規(guī)律 (a)900 Pa;(b)2700 Pa;(c)3500 Pa;(d)4500 PaFig.8.Axial profiles of Mach number,electron number density,electron temperature and heavy-particle temperature in the upstream of the bluff body with different values of Δp:(a)900 Pa;(b)2700 Pa;(c)3500 Pa;(d)4500 Pa.

    圖9給出了不同Δp下包覆鈍體頭部(FG)的等離子體區(qū)電子數(shù)密度高于臨界電子數(shù)密度(ne,crt=1010cm?3)的最小鞘套厚度(δmin)以及鈍體表面的總熱流密度(Γtot)的變化規(guī)律.可以看到:1)隨著Δp從300 Pa增大到2700 Pa,覆蓋鈍體頭部的等離子體鞘套最小厚度從7.0 cm單調(diào)增加到15.0 cm,這與實驗結果得到的趨勢一致,而由等離子體注入鈍體表面的總熱流密度則變化較小,在5.0—6.0 kW之間波動;2)當Δp超過3000 Pa后,由于等離子體內(nèi)能的損失,導致射流溫度在噴管出口與鈍體表面之間的區(qū)域相比于亞音速工況有所降低,因此δmin迅速降低到約3.0 cm,但是之后隨著Δp的持續(xù)增加而緩慢增加到6.0 cm;而與此同時,Γtot則從5.5 kW單調(diào)降低到2.0 kW.

    圖9 不同Δp下鈍體表面總熱流密度和包覆鈍體頭部的等離子體鞘套最小厚度的變化規(guī)律Fig.9.Variations of the total heat flux to the bluffbody surface from plasmas and the minimum thickness of the plasma layer covering the head of the bluffbody under different values of Δp.

    圖10 不同Δp下鈍體表面局部熱流密度沿射流幾何軸線的分布Fig.10.Profiles of the local heat flux to the bluffbody surface from plasmas along the geometrical axis of the plasma jet under different values of Δp.

    圖10則給出了不同噴嘴進出口壓差條件下鈍體表面熱流密度沿射流幾何軸線的分布規(guī)律.可以看到,在Δp保持不變的條件下,鈍體頭部(FG)區(qū)域的局部熱流密度(Γ)較其后部區(qū)域(EF)要大得多,大約是后者的2—4倍;而且隨著Δp的增加,沿鈍體表面的熱流密度分布變得愈加不均勻.由此我們可以得出這樣的結論:在維持等離子體發(fā)生器進口氣體壓強不變的條件下,減小真空腔壓強有利于提高等離子體射流的流動速度、降低等離子體射流向鈍體表面的傳熱量,但覆蓋鈍體頭部的等離子體鞘套的均勻性則顯著降低,且其厚度呈現(xiàn)出先降低再升高的非單調(diào)性變化趨勢.

    5 結論與展望

    本文以臨近空間高超聲速飛行器與其周圍復雜介質(zhì)環(huán)境間傳熱、流動以及電磁相互作用的研究為背景,進行了不同真空腔環(huán)境壓強下等離子體自由射流和沖擊射流的二維非平衡數(shù)值模擬,得到的主要結論包括:

    1)在亞音速流動條件下,數(shù)值模擬結果與MPX-2015實驗平臺上的實驗測量結果符合良好,這在一定程度上驗證了本文數(shù)值模擬所建立的物理數(shù)學模型的正確性;

    2)理論上通過調(diào)節(jié)縮放噴管進出口壓力差可以實現(xiàn)高溫等離子體流動從亞音速向超音速的發(fā)展.隨著噴管進出口壓差的增大,流經(jīng)噴管的等離子體流動狀態(tài)從亞音速發(fā)展為超音速;覆蓋鈍體表面的等離子體鞘套的厚度先減小而后又增加,鞘套的空間均勻性顯著降低;等離子體射流向鈍體表面的總傳熱量隨之降低,而鈍體頭部壁面附近的電子數(shù)密度平均值則增大,即包覆鈍體的等離子體鞘套的高電子數(shù)密度區(qū)域更加向鈍體頭部集中.

    盡管數(shù)值模擬表明我們所建立的數(shù)學物理模型能夠自洽描述包覆飛行器表面的、具有可壓縮、非平衡、熱-電-磁-流動-化學反應強耦合特性的等離子體環(huán)境,但本文僅對準穩(wěn)態(tài)氬等離子體射流區(qū)進行了初步的數(shù)值模擬,并未開展電弧區(qū)更加復雜的包括電磁效應的非平衡態(tài)等離子體特性的數(shù)值模擬研究;而放電區(qū)的計算結果無疑會對射流區(qū)特性產(chǎn)生重要的影響.我們將在今后的工作中進一步深入開展不同工作氣體,特別是空氣條件下的跨音速非平衡態(tài)等離子體電弧射流特性的數(shù)值模擬研究,并通過測量不同流動狀態(tài)下等離子體射流的速度分布[37,38]來與數(shù)值模擬結果進行更加仔細的對比分析.

    [1]Keidar M,Kim M,Boyd I D 2008J.Spacecr.Rockets45 445

    [2]Morris R A,Bench P M,Golden K E,Sutton E A 199937th Aerospace Sciences Meeting and ExhibitReno,NV,USA,January 11–14,1999 AIAA-99-0630

    [3]Evans J S,Schexnayder Jr C J,Huber P W 1973NASA TN D-7332

    [4]Gillman E D,Foster J E,Blankson I M 2010NASA/TM-2010-216220

    [5]Rybak J P,Churchill R J 1971IEEE Trans.Aerosp.Electron.Syst.7 879

    [6]Mather D E,Pasqual J M,Sillence J P,Lewis P 2005AIAA/CIRA 13th International Space Planes and Hypersonics Systems and Technologies ConferenceCapua,Italy,May 16–20,2005 AIAA-2005-3443

    [7]Akey N D 1971NASA Special Publication252 19

    [8]Watillon P,Berthe P,Chavagnac C 2003AIAA/ICAS International Air and Space Symposium and Exposition:The Next 100 yearsDayton,OH,USA,July 14–17,2003 AIAA-2003-2913

    [9]Shirouzu M,Yamamoto M,Shirouzu M,Yamamoto M 1996Space Plane and Hypersonic Systems and Technology ConferenceNorfolk,VA,USA,November 18–22,1996 AIAA-96-4524-CP

    [10]Yanagihara M,Munenaga T 200424th International CongressoftheAeronauticalSciencesYokohama,Japan,August 29–September 3,2004 p2004-7

    [11]Sakurai H,Kobayasi M,Yamazaki I,Shirouzu M,Yamamoto M 1997Acta Astronaut.40 105

    [12]Auweter Kurtz M,Kurtz H L,Laure S 1996J.Propul.Power12 1053

    [13]Zhao L,Liu X X,Su H S 2015J.Telemetry,Tracking and Command36 28(in Chinese)[趙良,劉秀祥,蘇漢生2015遙測遙控36 28]

    [14]Hermann T,L?hle S,Zander F,Fulge H,Fasoulas S 2016J.Thermophys.Heat Transf.30 673

    [15]Guo H,Su Y B,Li H P,Zeng S,Nie Q Y,Li Z X,Li Z H 2018Acta Phys.Sin.67 045201(in Chinese)[郭恒,蘇運波,李和平,曾實,聶秋月,李占賢,李志輝 2018物理學報67 045201]

    [16]Wang Z,Wu G Q,Ge N,Li H P,Bao C Y 2010IEEE Trans.Plasma Sci.38 2906

    [17]Yao Y,Hossain M M,Watanabe T,Matsuura T,Funabiki F,Yano T 2008Chem.Eng.J.139 390

    [18]Watanabe T,Liu Y,Tanaka M 2014Plasma Chem.Plasma Process.34 443

    [19]Savino R,Paterna D,Fumo M S,D’Elia M 2010Open Aero.Eng.J.3 76

    [20]Scalabrin L C,Boyd I D 20069th AIAA/ASME Joint Thermophysics and Heat Transfer ConferenceSan Francisco,CA,USA,June 5–8,2006 AIAA-2006-3773

    [21]Gnoffo P A,Gupta R N,Shinn J L 1989NASA-TP-2867

    [22]Candler G V,MacCormack R W 1991J.Thermophys.Heat Transf.5 266

    [23]Huang H,Qu Z H 1999Acta Aerodyn.Sin.17 462(in Chinese)[黃華,瞿章華1999空氣動力學學報17 462]

    [24]Sun S R,Wang H X 2015Acta Phys.Sin.64 143401(in Chinese)[孫素蓉,王海興 2015物理學報 64 143401]

    [25]Wang H X,Sun S R,Chen S Q 2012Acta Phys.Sin.61 195203(in Chinese)[王海興,孫素蓉,陳士強2012物理學報61 195203]

    [26]Wang H X,Sun W P,Sun S R,Murphy A B,Ju Y 2014Plasma Chem.Plasma Process.34 559

    [27]Niu C,Chen Z,Rong M,Wang C,Wu Y,Yang F,Wang X,Pang Q 2016J.Phys.D:Appl.Phys.49 405204

    [28]Wang C,Wu Y,Chen Z,Yang F,Feng Y,Rong M,Zhang H 2016Plasma Sci.Technol.18 732

    [29]Rat V,Murphy A B,Aubreton J,Elchinger M F,Fauchais P 2008J.Phys.D:Appl.Phys.41 183001

    [30]Zhang X N,Li H P,Murphy A B,Xia W D 2013Phys.Plasmas20 033508

    [31]Li H P,Zhang X N,Xia W D 2013Phys.Plasmas20 033509

    [32]Zhang X N,Li H P,Murphy A B,Xia W D 2015Plasma Sources Sci.Technol.24 035011

    [33]Zhang X N,Li H P,Xia W D 2013High Voltage Engineering39 1640(in Chinese)[張曉寧,李和平,夏維東2013高電壓技術39 1640]

    [34]Zhang X N 2015Ph.D.Dissertation(Hefei:University of Science and Technology of China)(in Chinese)[張曉寧2015博士學位論文(合肥:中國科學技術大學)]

    [35]Tong B G,Kong X Y,Deng G H 1990Aerodynamics(Beijing:Higher Education Press)p73(in Chinese)[童秉綱,孔祥言,鄧國華1990氣體動力學 (北京:高等教育出版社)第73頁]

    [36]Fang M,Li Z H,Li Z H,Tian Y 2017Acta Aerodyn.Sin.35 39(in Chinese)[方明,李志輝,李中華,田穎2017空氣動力學學報35 39]

    [37]Li G,Xu Y J,Mu K J,Nie C Q,Zhu J Q,Zhang Y,Li H M 2008Acta Phys.Sin.57 6444(in Chinese)[李鋼,徐燕驥,穆克進,聶超群,朱俊強,張翼,李漢明 2008物理學報57 6444]

    [38]Meng X,Pan W X,Wu C K 2004J.Eng.Thermophys.25 490(in Chinese)[孟顯,潘文霞,吳承康2004工程熱物理學報25 490]

    猜你喜歡
    射流等離子體飛行器
    深海逃逸艙射流注水均壓過程仿真分析
    低壓天然氣泄漏射流擴散特性研究
    煤氣與熱力(2022年4期)2022-05-23 12:45:00
    高超聲速飛行器
    連續(xù)磁活動對等離子體層演化的影響
    基于低溫等離子體修飾的PET/PVC浮選分離
    復雜飛行器的容錯控制
    電子制作(2018年2期)2018-04-18 07:13:25
    等離子體種子處理技術介紹
    神秘的飛行器
    射流齒形噴嘴射流流場與氣動聲學分析
    地鐵站臺活塞風附壁射流起始段的實測和實驗驗證
    制冷學報(2014年3期)2014-03-01 03:07:17
    久久久久久久久中文| 午夜福利在线在线| 一进一出好大好爽视频| 欧美日韩乱码在线| 国产视频内射| 国产亚洲精品av在线| 成人特级黄色片久久久久久久| 一a级毛片在线观看| 免费看美女性在线毛片视频| 国内精品久久久久久久电影| 日日爽夜夜爽网站| 亚洲 欧美 日韩 在线 免费| 国产精品二区激情视频| 久久久国产成人免费| 九色国产91popny在线| 亚洲专区国产一区二区| 国内揄拍国产精品人妻在线 | 91大片在线观看| 不卡av一区二区三区| 99国产极品粉嫩在线观看| 久久久水蜜桃国产精品网| 午夜免费激情av| 成人18禁在线播放| 婷婷亚洲欧美| 又紧又爽又黄一区二区| 国产欧美日韩一区二区三| 亚洲国产精品成人综合色| 性欧美人与动物交配| netflix在线观看网站| 国产精品综合久久久久久久免费| 亚洲av日韩精品久久久久久密| 欧美色视频一区免费| 亚洲在线自拍视频| 国产成人欧美| 欧美性猛交黑人性爽| av福利片在线| 精品国产国语对白av| 天天躁夜夜躁狠狠躁躁| 欧美性猛交╳xxx乱大交人| 久久这里只有精品19| 日日摸夜夜添夜夜添小说| 黄色 视频免费看| 午夜免费观看网址| 久久午夜亚洲精品久久| 久久国产亚洲av麻豆专区| 国产爱豆传媒在线观看 | 99在线视频只有这里精品首页| 久久国产精品人妻蜜桃| 日本成人三级电影网站| 国产精品久久久久久精品电影 | 麻豆久久精品国产亚洲av| 国产成年人精品一区二区| 十八禁人妻一区二区| 日本三级黄在线观看| 欧美日韩精品网址| 18禁国产床啪视频网站| 最近在线观看免费完整版| 精品久久久久久久久久免费视频| 国产精品综合久久久久久久免费| 99国产综合亚洲精品| 在线看三级毛片| 热re99久久国产66热| 亚洲全国av大片| 琪琪午夜伦伦电影理论片6080| 久99久视频精品免费| 国产1区2区3区精品| 久久久久久九九精品二区国产 | 国产国语露脸激情在线看| 国产精品野战在线观看| 波多野结衣av一区二区av| 国产黄色小视频在线观看| 女人被狂操c到高潮| 国产伦一二天堂av在线观看| av欧美777| 亚洲国产欧美一区二区综合| 黄色a级毛片大全视频| 999久久久精品免费观看国产| 88av欧美| 一级作爱视频免费观看| 美国免费a级毛片| 波多野结衣av一区二区av| 亚洲av第一区精品v没综合| 人妻久久中文字幕网| 亚洲午夜精品一区,二区,三区| 亚洲国产毛片av蜜桃av| 亚洲aⅴ乱码一区二区在线播放 | 欧美绝顶高潮抽搐喷水| 久久精品国产亚洲av香蕉五月| 老鸭窝网址在线观看| 丰满人妻熟妇乱又伦精品不卡| 国产亚洲精品久久久久久毛片| 99精品欧美一区二区三区四区| 午夜精品在线福利| 亚洲一区二区三区不卡视频| 少妇被粗大的猛进出69影院| 麻豆国产av国片精品| 国产激情偷乱视频一区二区| 99国产极品粉嫩在线观看| 免费观看人在逋| 欧美另类亚洲清纯唯美| 黄色 视频免费看| 欧美日韩中文字幕国产精品一区二区三区| 狂野欧美激情性xxxx| 亚洲国产精品合色在线| 男女床上黄色一级片免费看| 国产区一区二久久| 色婷婷久久久亚洲欧美| 美女高潮喷水抽搐中文字幕| 亚洲成a人片在线一区二区| 亚洲自偷自拍图片 自拍| 中文字幕人成人乱码亚洲影| 国产av又大| 看片在线看免费视频| 99国产极品粉嫩在线观看| 精品免费久久久久久久清纯| 黑人欧美特级aaaaaa片| 此物有八面人人有两片| 色精品久久人妻99蜜桃| 国产精品一区二区免费欧美| av中文乱码字幕在线| 黄频高清免费视频| 天堂影院成人在线观看| 好男人电影高清在线观看| 亚洲九九香蕉| 美女扒开内裤让男人捅视频| 国产成人欧美| 狠狠狠狠99中文字幕| 欧美国产日韩亚洲一区| 亚洲全国av大片| 久久人妻福利社区极品人妻图片| 天堂影院成人在线观看| 日韩欧美国产一区二区入口| 国产黄色小视频在线观看| 亚洲国产精品久久男人天堂| 一进一出抽搐gif免费好疼| 亚洲人成77777在线视频| 久久久久国产一级毛片高清牌| 白带黄色成豆腐渣| 亚洲国产欧洲综合997久久, | 中文字幕人成人乱码亚洲影| 欧美人与性动交α欧美精品济南到| 久久久久久久久中文| 丁香欧美五月| 午夜精品在线福利| 色老头精品视频在线观看| 国产蜜桃级精品一区二区三区| 成人手机av| 成人精品一区二区免费| 在线永久观看黄色视频| 中文字幕人妻熟女乱码| 激情在线观看视频在线高清| 久久久久久久久久黄片| 可以免费在线观看a视频的电影网站| 亚洲色图av天堂| 在线免费观看的www视频| 国产极品粉嫩免费观看在线| 午夜激情av网站| 美女高潮喷水抽搐中文字幕| 久久久久国产一级毛片高清牌| 亚洲人成网站在线播放欧美日韩| 变态另类丝袜制服| 欧美中文综合在线视频| 一级毛片女人18水好多| 国产又爽黄色视频| 真人一进一出gif抽搐免费| 久久精品国产综合久久久| 亚洲成人精品中文字幕电影| 欧美乱妇无乱码| 黑人操中国人逼视频| 巨乳人妻的诱惑在线观看| 亚洲黑人精品在线| 精品卡一卡二卡四卡免费| 国产真人三级小视频在线观看| 久久久久国内视频| 好看av亚洲va欧美ⅴa在| 婷婷六月久久综合丁香| 精品一区二区三区四区五区乱码| 曰老女人黄片| 久久精品影院6| 国产97色在线日韩免费| 999久久久精品免费观看国产| 国产精品自产拍在线观看55亚洲| 久久人妻av系列| 亚洲精品在线观看二区| 国产精华一区二区三区| 欧美日韩瑟瑟在线播放| 亚洲中文av在线| 欧美成人免费av一区二区三区| 国产午夜福利久久久久久| 免费高清在线观看日韩| 精品第一国产精品| 麻豆成人av在线观看| 成人一区二区视频在线观看| 成年女人毛片免费观看观看9| 亚洲欧美日韩高清在线视频| 黄网站色视频无遮挡免费观看| 90打野战视频偷拍视频| 久久热在线av| 欧美在线黄色| 无人区码免费观看不卡| 日韩中文字幕欧美一区二区| 欧美成人性av电影在线观看| 看黄色毛片网站| 中文字幕人成人乱码亚洲影| 中文字幕人妻熟女乱码| 国产亚洲av高清不卡| 国产精品爽爽va在线观看网站 | 国产亚洲精品久久久久久毛片| 国产一区二区三区在线臀色熟女| 淫妇啪啪啪对白视频| 久久人人精品亚洲av| 这个男人来自地球电影免费观看| 日韩国内少妇激情av| 美女高潮喷水抽搐中文字幕| svipshipincom国产片| 亚洲 国产 在线| 一级片免费观看大全| 人人妻人人看人人澡| 国产蜜桃级精品一区二区三区| 中文字幕精品亚洲无线码一区 | 亚洲欧美一区二区三区黑人| 亚洲中文日韩欧美视频| 1024香蕉在线观看| 久久精品国产99精品国产亚洲性色| 亚洲国产高清在线一区二区三 | 丁香六月欧美| 美女高潮喷水抽搐中文字幕| 欧美一级毛片孕妇| 伊人久久大香线蕉亚洲五| 久久亚洲真实| 热re99久久国产66热| 久久精品国产清高在天天线| 人妻久久中文字幕网| 999久久久国产精品视频| 一级片免费观看大全| 脱女人内裤的视频| 嫩草影视91久久| АⅤ资源中文在线天堂| 日韩欧美免费精品| 免费在线观看成人毛片| 黄片播放在线免费| 久久精品国产综合久久久| 国产爱豆传媒在线观看 | 欧美成人性av电影在线观看| 久久久久免费精品人妻一区二区 | 亚洲男人的天堂狠狠| 国产一级毛片七仙女欲春2 | 国产v大片淫在线免费观看| 亚洲av成人一区二区三| 人人妻,人人澡人人爽秒播| 露出奶头的视频| 久久久国产成人精品二区| 两性午夜刺激爽爽歪歪视频在线观看 | 欧美成狂野欧美在线观看| 精品久久久久久久久久久久久 | 国产精品久久久久久精品电影 | 国产视频一区二区在线看| 色综合亚洲欧美另类图片| www.www免费av| 丰满人妻熟妇乱又伦精品不卡| 色婷婷久久久亚洲欧美| 国产高清激情床上av| 欧美在线一区亚洲| 亚洲国产中文字幕在线视频| 一边摸一边做爽爽视频免费| 国产精品久久久久久精品电影 | 日韩中文字幕欧美一区二区| 深夜精品福利| 欧美日韩福利视频一区二区| 久久久久久久久中文| 色在线成人网| 黄色片一级片一级黄色片| 亚洲五月色婷婷综合| 搡老妇女老女人老熟妇| 精品国产国语对白av| 亚洲久久久国产精品| 美女扒开内裤让男人捅视频| 波多野结衣av一区二区av| 日韩有码中文字幕| 老汉色∧v一级毛片| 一级片免费观看大全| 国产成+人综合+亚洲专区| 亚洲精品美女久久av网站| 国产久久久一区二区三区| 女人高潮潮喷娇喘18禁视频| 国产国语露脸激情在线看| 日本黄色视频三级网站网址| 老鸭窝网址在线观看| 国产成+人综合+亚洲专区| 久久香蕉国产精品| 国产99久久九九免费精品| 亚洲片人在线观看| 精品午夜福利视频在线观看一区| 人人妻人人澡欧美一区二区| 精品无人区乱码1区二区| 国产午夜福利久久久久久| 最近最新免费中文字幕在线| 欧美不卡视频在线免费观看 | 亚洲欧美激情综合另类| 男人舔奶头视频| 欧美人与性动交α欧美精品济南到| 黑人巨大精品欧美一区二区mp4| 日本a在线网址| 看黄色毛片网站| 亚洲国产欧洲综合997久久, | 中文字幕精品免费在线观看视频| 老汉色av国产亚洲站长工具| 啦啦啦观看免费观看视频高清| 他把我摸到了高潮在线观看| 丰满人妻熟妇乱又伦精品不卡| 精华霜和精华液先用哪个| 非洲黑人性xxxx精品又粗又长| 亚洲精品国产区一区二| 两性夫妻黄色片| 在线观看免费视频日本深夜| av视频在线观看入口| 在线观看舔阴道视频| 哪里可以看免费的av片| 91老司机精品| 国产成年人精品一区二区| 国产高清videossex| 在线观看一区二区三区| 亚洲av电影不卡..在线观看| 国产又色又爽无遮挡免费看| 欧美成人免费av一区二区三区| 村上凉子中文字幕在线| 夜夜夜夜夜久久久久| 亚洲全国av大片| 男人操女人黄网站| 丝袜在线中文字幕| 麻豆国产av国片精品| 欧美三级亚洲精品| 大香蕉久久成人网| 亚洲国产高清在线一区二区三 | 免费在线观看视频国产中文字幕亚洲| 亚洲人成77777在线视频| 麻豆一二三区av精品| 国产av不卡久久| 国产在线精品亚洲第一网站| 在线免费观看的www视频| 欧美国产日韩亚洲一区| 一区福利在线观看| 法律面前人人平等表现在哪些方面| 搡老岳熟女国产| 成人国语在线视频| 丁香欧美五月| 人妻久久中文字幕网| 97人妻精品一区二区三区麻豆 | 18禁观看日本| 99热6这里只有精品| 国产精品综合久久久久久久免费| 一进一出抽搐动态| 自线自在国产av| av福利片在线| 动漫黄色视频在线观看| 少妇裸体淫交视频免费看高清 | 中文资源天堂在线| 久久人妻福利社区极品人妻图片| 国产亚洲精品综合一区在线观看 | 国产精品av久久久久免费| 黑人欧美特级aaaaaa片| 欧美在线黄色| 成年免费大片在线观看| 日韩成人在线观看一区二区三区| 免费看日本二区| 欧美av亚洲av综合av国产av| 老司机福利观看| 欧美在线黄色| 欧美中文综合在线视频| 一进一出好大好爽视频| 亚洲人成网站高清观看| 久久精品影院6| 精品国内亚洲2022精品成人| 国产精品日韩av在线免费观看| 成人一区二区视频在线观看| 亚洲无线在线观看| 黑人操中国人逼视频| 国产一区二区三区在线臀色熟女| 国产伦在线观看视频一区| 精品国产乱码久久久久久男人| 欧美日本视频| 日韩欧美免费精品| 成人国产综合亚洲| 久久精品成人免费网站| 中文亚洲av片在线观看爽| 国产成人影院久久av| 在线观看免费视频日本深夜| 欧美成狂野欧美在线观看| 不卡一级毛片| 日本 欧美在线| 欧美乱码精品一区二区三区| 国产精品久久久人人做人人爽| 这个男人来自地球电影免费观看| 中文资源天堂在线| 在线免费观看的www视频| www国产在线视频色| 成人av一区二区三区在线看| 性色av乱码一区二区三区2| 欧美大码av| 亚洲国产毛片av蜜桃av| 午夜精品久久久久久毛片777| 欧美日本视频| 2021天堂中文幕一二区在线观 | 91九色精品人成在线观看| 一区福利在线观看| 国产麻豆成人av免费视频| 国产成+人综合+亚洲专区| 99热只有精品国产| 国产日本99.免费观看| 亚洲国产高清在线一区二区三 | 国产精品永久免费网站| 91成人精品电影| 操出白浆在线播放| 亚洲专区国产一区二区| 成人亚洲精品一区在线观看| 亚洲男人的天堂狠狠| 妹子高潮喷水视频| 性色av乱码一区二区三区2| 中文字幕高清在线视频| 日本精品一区二区三区蜜桃| 午夜久久久在线观看| 精品电影一区二区在线| 18禁裸乳无遮挡免费网站照片 | 最好的美女福利视频网| 精品欧美国产一区二区三| 777久久人妻少妇嫩草av网站| 国产亚洲精品久久久久5区| 久久青草综合色| 他把我摸到了高潮在线观看| 国产精品亚洲av一区麻豆| 每晚都被弄得嗷嗷叫到高潮| 国产精品 国内视频| 亚洲国产欧美网| 亚洲成av片中文字幕在线观看| 搡老熟女国产l中国老女人| 午夜免费成人在线视频| 成人午夜高清在线视频 | 黄色女人牲交| 老司机午夜十八禁免费视频| 精品国产亚洲在线| 91av网站免费观看| 国产精品野战在线观看| 国产成年人精品一区二区| 精品久久久久久成人av| 波多野结衣av一区二区av| 日本一区二区免费在线视频| 亚洲专区国产一区二区| 欧美日韩亚洲综合一区二区三区_| 亚洲欧美日韩高清在线视频| 国产av一区二区精品久久| 桃色一区二区三区在线观看| 久久婷婷成人综合色麻豆| 一个人观看的视频www高清免费观看 | 精品福利观看| av在线天堂中文字幕| 色播在线永久视频| 成人免费观看视频高清| 国产人伦9x9x在线观看| 亚洲片人在线观看| 精品久久蜜臀av无| 黄色视频不卡| 亚洲第一欧美日韩一区二区三区| 老司机靠b影院| 免费看十八禁软件| 曰老女人黄片| 听说在线观看完整版免费高清| 一级毛片高清免费大全| 19禁男女啪啪无遮挡网站| 国产精品一区二区免费欧美| 亚洲激情在线av| 免费看日本二区| 精品人妻1区二区| 高清毛片免费观看视频网站| 婷婷亚洲欧美| 两个人视频免费观看高清| 91在线观看av| 国产精品九九99| 国产又黄又爽又无遮挡在线| 亚洲成人国产一区在线观看| 日韩欧美在线二视频| 亚洲国产精品sss在线观看| 夜夜夜夜夜久久久久| 日本免费一区二区三区高清不卡| 人人妻人人澡人人看| www.自偷自拍.com| 久久这里只有精品19| 日韩欧美国产一区二区入口| 很黄的视频免费| 午夜视频精品福利| 真人一进一出gif抽搐免费| 在线观看66精品国产| 日日干狠狠操夜夜爽| 久久 成人 亚洲| 亚洲第一青青草原| 亚洲第一欧美日韩一区二区三区| 久久久久亚洲av毛片大全| 女警被强在线播放| 精品不卡国产一区二区三区| 亚洲欧美一区二区三区黑人| 精品久久久久久久人妻蜜臀av| 午夜a级毛片| 最新在线观看一区二区三区| 欧美av亚洲av综合av国产av| 国产麻豆成人av免费视频| 国产精品综合久久久久久久免费| 一区二区三区精品91| 久热爱精品视频在线9| 亚洲性夜色夜夜综合| 亚洲在线自拍视频| 午夜福利一区二区在线看| 午夜福利成人在线免费观看| 国产精品久久久久久亚洲av鲁大| 久久久国产成人免费| 亚洲成人久久性| 久久天躁狠狠躁夜夜2o2o| 白带黄色成豆腐渣| 一夜夜www| 欧美 亚洲 国产 日韩一| 国产三级黄色录像| 免费观看精品视频网站| 欧美激情极品国产一区二区三区| 国产激情欧美一区二区| 国产黄a三级三级三级人| 亚洲欧美日韩高清在线视频| 午夜免费观看网址| 国产v大片淫在线免费观看| 波多野结衣高清无吗| 中文字幕精品亚洲无线码一区 | 美女免费视频网站| 中文字幕精品亚洲无线码一区 | 亚洲精品久久成人aⅴ小说| 国产成人一区二区三区免费视频网站| 久久香蕉精品热| 啪啪无遮挡十八禁网站| 丁香六月欧美| 久久久久久久久久黄片| 午夜久久久久精精品| 巨乳人妻的诱惑在线观看| 午夜老司机福利片| 国产精品综合久久久久久久免费| 久久性视频一级片| 侵犯人妻中文字幕一二三四区| 叶爱在线成人免费视频播放| 日日摸夜夜添夜夜添小说| 1024香蕉在线观看| 国产激情久久老熟女| 亚洲国产精品999在线| 日本撒尿小便嘘嘘汇集6| 老司机靠b影院| 精品欧美国产一区二区三| 中文亚洲av片在线观看爽| 日本 欧美在线| 成年版毛片免费区| 国产精品久久久久久人妻精品电影| 美女国产高潮福利片在线看| 亚洲欧美一区二区三区黑人| 亚洲午夜理论影院| 变态另类丝袜制服| 50天的宝宝边吃奶边哭怎么回事| 99热只有精品国产| 50天的宝宝边吃奶边哭怎么回事| 亚洲久久久国产精品| 啦啦啦韩国在线观看视频| 男女床上黄色一级片免费看| 精品高清国产在线一区| 精品国产乱码久久久久久男人| 激情在线观看视频在线高清| 久久久久九九精品影院| 每晚都被弄得嗷嗷叫到高潮| aaaaa片日本免费| 美女免费视频网站| 一本精品99久久精品77| 黄色视频不卡| 一本综合久久免费| 国产精品久久视频播放| 免费在线观看影片大全网站| 在线看三级毛片| 亚洲片人在线观看| 我的亚洲天堂| a级毛片a级免费在线| www日本在线高清视频| 亚洲五月婷婷丁香| 大型黄色视频在线免费观看| www国产在线视频色| 欧美午夜高清在线| 国产激情欧美一区二区| 国产亚洲精品久久久久5区| 少妇 在线观看| 久久香蕉激情| 日韩成人在线观看一区二区三区| 黄色丝袜av网址大全| 十八禁网站免费在线| 9191精品国产免费久久| 久久久久久久午夜电影| 免费在线观看日本一区| 波多野结衣高清作品| 国产在线精品亚洲第一网站| 国产午夜精品久久久久久| 国产1区2区3区精品| 制服丝袜大香蕉在线| www日本在线高清视频| 免费搜索国产男女视频| 亚洲精品在线美女| 亚洲色图 男人天堂 中文字幕| ponron亚洲| 在线av久久热| 黄色女人牲交| 脱女人内裤的视频| 成在线人永久免费视频| www日本在线高清视频| 久久久久国内视频| 亚洲成人久久性| e午夜精品久久久久久久| 午夜a级毛片| 91老司机精品|