亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種乘波構(gòu)型邊緣鈍化方法的仿真與試驗研究

        2014-11-08 06:18:06劉建霞侯中喜張扣立馬曉偉
        空氣動力學(xué)學(xué)報 2014年2期
        關(guān)鍵詞:駐點氣動力超聲速

        劉建霞,塵 軍,侯中喜,張扣立,馬曉偉

        (1.中國空氣動力研究與發(fā)展中心,綿陽 621000;2.國防科技大學(xué) 航天科學(xué)與工程學(xué)院,長沙 410073;3.空間物理重點實驗室,北京 100076)

        0 引 言

        與常規(guī)再入飛行器不同,先進(jìn)氣動布局是確保近空間高超聲速飛行器實現(xiàn)遠(yuǎn)程達(dá)到、靈活機(jī)動、精確打擊能力的必要條件。對于近空間高超聲速滑翔飛行器,其升阻比一般要求較高。常規(guī)布局設(shè)計的高超聲速飛行器,激波從邊緣脫體,流體繞過邊緣在上下表面之間發(fā)生交互,飛行器升力性能下降;同時,在高超聲速飛行條件時,強(qiáng)激波和邊界層的存在,導(dǎo)致飛行器承受的阻力增加,由此出現(xiàn)了屈西曼[1]提出的高超聲速飛行器的“升阻比屏障”。1959年,Nonweiler教授[2]提出了新型乘波布局的設(shè)計理念,其原理是將激波后的高壓氣流限制在飛行器下表面,不允許其繞過邊緣泄漏到上表面,從而獲得了比常規(guī)布局更高的升阻比。此后的50多年間,人們通過大量數(shù)值仿真和風(fēng)洞試驗證實:基于粘性優(yōu)化設(shè)計的理想乘波構(gòu)型是突破高超聲速飛行器“升阻比屏障”的一種有效嘗試[3-10]。在近兩年里,美國的 X-51A、HTV-2等先進(jìn)高超聲速飛行器便采用乘波構(gòu)型作為基本布局開展了多次飛行演示試驗,引起世界各國的廣泛關(guān)注。

        理想乘波構(gòu)型的優(yōu)異氣動性能與其尖銳的邊緣設(shè)計存在緊密聯(lián)系,但在實際過程中,絕對尖銳的邊緣不僅從結(jié)構(gòu)的加工工藝和力學(xué)特性角度難以實現(xiàn)[11],其尖邊緣處的強(qiáng)烈氣動加熱效應(yīng)也被 Ohta[12]的風(fēng)洞試驗證實。因此,乘波構(gòu)型的絕對尖銳邊緣在實際飛行中難以保持。同時,由于乘波構(gòu)型的綜合性能對外形十分敏感,一般要求其外形在飛行過程中具有非燒蝕或低燒蝕的特點。結(jié)合以上兩方面的考慮,在開展高超聲速乘波飛行器總體設(shè)計過程中,鈍化方法的選擇顯得尤為重要。在近幾年里,關(guān)于鈍化尺度對理想乘波構(gòu)型氣動力性能和氣動加熱特點影響的研究已經(jīng)成為高超聲速飛行器設(shè)計領(lǐng)域的熱點[11,13-16]。

        從國內(nèi)外文獻(xiàn)看,鈍化造成理想乘波構(gòu)型升阻比下降的現(xiàn)象已被多位學(xué)者證實。文獻(xiàn)[11,13]的仿真結(jié)果表明,對于“Λ”型乘波構(gòu)型,其升力系數(shù)、阻力系數(shù)、升阻比等對鈍化尺度十分敏感;采用半徑不足構(gòu)型總長度1%的鈍邊緣,乘波構(gòu)型升阻比下降約50%。文獻(xiàn)[14]的風(fēng)洞試驗數(shù)據(jù)表明:對于錐導(dǎo)乘波構(gòu)型,引入鈍邊緣后,其升阻比下降幅度達(dá)到19.74%。文獻(xiàn)[15]針對相交楔錐方法生成的乘波構(gòu)型得到的仿真結(jié)果表明:當(dāng)曲率半徑為0.01m時,乘波構(gòu)型升阻比下降7.21%。文獻(xiàn)[16]分別采用0.01m、0.02m 和0.03m的半徑對優(yōu)化設(shè)計的錐導(dǎo)乘波構(gòu)型進(jìn)行鈍化,并采用數(shù)值模擬方法對三者的氣動力性能進(jìn)行了對比研究,結(jié)果表明:鈍化半徑每增加0.01m,乘波構(gòu)型的升阻比下降約6.5%。怎樣在氣動力性能和防熱設(shè)計之間尋找平衡,采用合適鈍化方法進(jìn)行修形,已成為“乘波構(gòu)型”設(shè)計概念保持其在高超聲速領(lǐng)域優(yōu)勢地位亟待解決的問題。

        分析發(fā)現(xiàn),已開展的研究工作都對乘波構(gòu)型邊緣采用了相同尺度進(jìn)行鈍化,這種鈍化方法被本文命名為“一致邊緣鈍化方法”。然而,根據(jù)文獻(xiàn)[16]提供的一致邊緣鈍化乘波構(gòu)型表面受熱特點看,其高熱流僅局限在頭部區(qū)域;絕大部分邊緣的受熱形勢遠(yuǎn)比該區(qū)域緩和。在設(shè)計過程中,若采用滿足駐點防熱需求的曲率半徑對乘波構(gòu)型所有邊緣進(jìn)行鈍化,從熱防護(hù)角度看是存在冗余的;另一方面,從氣動力性能看,邊緣的大尺度鈍化會造成該位置的氣體泄露及激波強(qiáng)度變化,進(jìn)而降低乘波飛行器的升阻比?;谝陨戏治?,本文提出了一種采用不同鈍化尺度對乘波構(gòu)型不同邊緣位置進(jìn)行修形的新方法——“非一致邊緣鈍化方法”,并采用數(shù)值模擬和風(fēng)洞試驗相結(jié)合的研究手段,對根據(jù)該方法設(shè)計的鈍邊緣乘波構(gòu)型的氣動力性能和氣動加熱特點進(jìn)行了分析,以探索該鈍化方法的可行性。

        1 計算模型及方法

        1.1 計算模型

        本文選擇錐導(dǎo)乘波構(gòu)型作為研究對象,其設(shè)計步驟如下:首先,在錐形激波底部圓上給出自由面的底部基線;其次,根據(jù)上表面平行于來流的特點向上游推進(jìn),得到上表面以及該平面與激波面的交線——飛行器的邊緣曲線;最后,根據(jù)下表面平行流線的特點向下游推進(jìn)生成下表面,同時得到該曲面與底部圓的交線——下表面底部基線。在設(shè)計馬赫數(shù)10,圓錐半錐角12°,設(shè)計長度3.5m,擴(kuò)張角45°的條件下,生成的尖邊緣乘波構(gòu)型如圖1所示。根據(jù)本文提出的鈍化方法,在駐點位置采用0.03m半徑進(jìn)行鈍化,對其靠后的80%邊緣采用0.01m半徑進(jìn)行鈍化,對兩者之間的區(qū)域采用由0.03m逐漸減小到0.01m的半徑進(jìn)行鈍化,由此得到的鈍邊緣乘波構(gòu)型如圖2所示,其長度約為3.2m。各位置的鈍化尺度根據(jù)文獻(xiàn)[17]中的經(jīng)驗公式計算得到。以駐點為例,乘波體駐點位置的對流熱流密度可以采用下式進(jìn)行估算:

        式中,K為三維效應(yīng)的修正因子;下標(biāo)∞表示來流參數(shù);下標(biāo)w表示飛行器物面參數(shù);R為駐點處鈍化半徑,H0為來流總焓。假定乘波體駐點區(qū)域采用碳化硅涂層材料進(jìn)行防熱,根據(jù)材料性能參數(shù)及輻射平衡條件,聯(lián)立方程即可求得駐點處鈍化尺度。

        圖1 尖邊緣乘波構(gòu)型Fig.1 Waverider with sharp leading edge

        圖2 鈍邊緣乘波構(gòu)型Fig.2 Waverider with blunted leading edge

        1.2 計算方法

        針對鈍邊緣乘波構(gòu)型流場,所采用的控制方程為三維N-S方程。在求解過程中,對流項采用AUSM格式進(jìn)行離散[18],粘性項采用中心差分格式進(jìn)行離散,時間項采用LU-SGS隱式迭代方法求解。邊界條件的設(shè)置如下:壁面滿足無滑移、法向壓力梯度為0等條件;入口參數(shù)由壓力遠(yuǎn)場條件得到,出口參數(shù)采用外推方法得到。

        1.3 有效性研究

        考慮到高超聲速流場熱流密度的高精度仿真對算法的要求嚴(yán)于氣動力性能仿真,采用兩組外形的試驗測量數(shù)據(jù),對本文數(shù)值算法在熱流密度模擬方面的仿真精度進(jìn)行研究。

        首先,基于鈍乘波構(gòu)型駐點區(qū)域流動與球頭繞流的相似性,進(jìn)行了半徑為0.015m的高超聲速球頭駐點熱流密度的數(shù)值仿真。該算例的試驗數(shù)據(jù)[19]由中國空氣動力研究與發(fā)展中心FD-14A高超聲速激波風(fēng)洞測得,表1列出了算例的試驗條件。鉑薄膜電阻溫度計測得的球頭駐點熱流密度為5.858×105W/m2;根據(jù)文獻(xiàn)[20]的經(jīng)驗公式,求得的駐點熱流密度為5.456×105W/m2;本文得到的仿真結(jié)果如圖3所示,駐點熱流密度為5.499×105W/m2。對比可知,仿真值與試驗測量數(shù)據(jù)的相對誤差約為6%,仿真值與經(jīng)驗公式估算值的相對誤差約為1%。

        基于鈍乘波構(gòu)型邊緣流動與圓柱繞流的相似性,針對半徑為0.038m的圓柱的表面熱流密度分布進(jìn)行數(shù)值仿真,并與Holden[21]的測量數(shù)據(jù)進(jìn)行了對比分析。表2列出了該算例的計算條件。由圖4可知,仿真結(jié)果與試驗數(shù)據(jù)在數(shù)值和趨勢上都非常吻合。試驗測得的駐點熱流密度約為6.0×105W/m2;根據(jù)文獻(xiàn)[20]的經(jīng)驗公式,求得的駐點熱流密度約為5.78×105W/m2;本文得到的仿真結(jié)果為 5.94×105W/m2。仿真值與測量值的相對誤差約為1%,仿真值與經(jīng)驗公式估算值的相對誤差約為3%。

        表1 球頭算例的試驗條件Table 1 Test conditions of sphere case

        圖3 球頭算例的網(wǎng)格及熱流密度分布Fig.3 Grid and heat flux distribution of sphere case

        表2 圓柱算例的計算條件Table 2 Computational conditions of cylinder case

        圖4 圓柱算例的網(wǎng)格及熱流密度分布Fig.4 Grid and heat flux distribution of cylinder case

        兩算例的研究結(jié)果表明,在針對類似于球頭和圓柱的高超聲速繞流問題進(jìn)行熱流密度仿真時,本文數(shù)值算法能夠確保結(jié)果達(dá)到一定精度。

        文獻(xiàn)[22]指出,高超聲速流場熱流密度預(yù)測的準(zhǔn)確性與其計算網(wǎng)格密切相關(guān),網(wǎng)格雷諾數(shù)被認(rèn)為是壁面熱流預(yù)測的重要因素。在開展數(shù)值模擬前,本文首先針對鈍邊緣乘波構(gòu)型進(jìn)行了網(wǎng)格收斂性分析。結(jié)果表明,當(dāng)網(wǎng)格雷諾數(shù)在1的量級時,其表面壓力分布及熱流密度分布都趨于收斂。鈍邊緣乘波構(gòu)型的計算網(wǎng)格如圖5所示。

        圖5 乘波構(gòu)型的計算網(wǎng)格Fig.5 Computational grid of waverider configuration

        2 試驗設(shè)備及模型

        2.1 試驗設(shè)備

        試驗在中國空氣動力研究與發(fā)展中心超高聲速所的FD-14A高超聲速激波風(fēng)洞中進(jìn)行。該風(fēng)洞為直徑0.6m的脈沖型高超聲速氣動力、氣動熱測試設(shè)備。風(fēng)洞由高壓段、低壓段、噴管、試驗段、真空箱、超高壓供氣系統(tǒng)等部分組成。風(fēng)洞的驅(qū)動氣體為氫氣和氮氣的混合物(體積比為0.85∶0.15),被驅(qū)動氣體為氮氣。試驗?zāi)P捅砻鏌崃髅芏炔捎弥睆綖?.002m的薄鉑膜電阻溫度計進(jìn)行測量。

        2.2 試驗?zāi)P?/h3>

        本文測熱模型是在圖2所示的鈍邊緣乘波構(gòu)型基礎(chǔ)上截取頭部區(qū)域獲得的。考慮到阻塞度影響,試驗?zāi)P偷某叽缡艿絿?yán)格控制。如圖6所示,在模型邊緣中心線上分布了許多直徑為0.0022m的測孔。圖6為安裝在風(fēng)洞中的模型照片。

        圖6 乘波構(gòu)型試驗?zāi)P虵ig.6 Experiment model of waverider configuration

        3 結(jié)果與討論

        3.1 數(shù)值模擬結(jié)果

        在高度40km、馬赫數(shù)10的條件下,對非一致邊緣鈍化乘波構(gòu)型氣動力性能和表面受熱情況進(jìn)行數(shù)值模擬,并與一致邊緣鈍化乘波構(gòu)型進(jìn)行了對比分析。表3給出了兩者氣動力數(shù)據(jù)的對比。由表可知,在相同來流下,與一致邊緣鈍化乘波構(gòu)型相比,非一致邊緣鈍化乘波構(gòu)型的升力增加282.4N,升阻比提高約3%。此外,由于在絕大部分邊緣采用更小的鈍化半徑,非一致邊緣鈍化乘波構(gòu)型表面的總面積增加約0.264m2,這一因素導(dǎo)致其所受總阻力出現(xiàn)小幅增加,大約20N。在非一致邊緣鈍化乘波構(gòu)型受到的阻力中,波阻仍是最主要的部分,其比例達(dá)到90%。

        考慮到兩類鈍化方法的差異主要體現(xiàn)在邊緣,圖7給出了相同橫截面上兩類鈍化乘波構(gòu)型邊緣處的壓力和密度分布對比。由圖可知,兩者的流場結(jié)構(gòu)具有相似性,但從絕對值看,采用更小鈍化尺度的非一致鈍化乘波構(gòu)型,通過邊緣向上表面泄漏的下表面高壓氣體要比一致鈍化乘波構(gòu)型少得多,這是乘波構(gòu)型氣動力性能得到改善的重要原因。

        表3 兩類鈍邊緣乘波構(gòu)型氣動力對比Table 3 Aerodynamic force of two blunted waveriders

        圖7 兩鈍邊緣乘波構(gòu)型邊緣處的流場對比Fig.7 Flow field around the leading edge of two blunted waveriders

        與一致邊緣鈍化方法相比,采用非一致邊緣鈍化方法進(jìn)行修形的乘波構(gòu)型,其下表面的熱流密度略有增加,上表面則基本不受影響。如圖8所示,下表面熱流密度增加約1.0×104W/m2,該值約為相同來流條件下駐點熱流密度的1%。由圖9可知,非一致邊緣鈍化乘波構(gòu)型的高熱流仍局限在駐點附近,該區(qū)域是乘波構(gòu)型熱防護(hù)系統(tǒng)設(shè)計的重點。對于其受熱相對緩和的絕大部分前緣,由于鈍化尺度減小,其熱流密度增加約1.0×105W/m2,該值約為駐點熱流密度的1/10。綜合兩種鈍邊緣乘波構(gòu)型的氣動力性能和表面受熱情況的對比結(jié)果看,非一致邊緣鈍化方法在基本不影響乘波構(gòu)型熱防護(hù)系統(tǒng)設(shè)計的前提下,可改善乘波構(gòu)型的氣動力性能。

        圖8 兩鈍邊緣乘波構(gòu)型表面受熱對比Fig.8 Heat flux on the surface of two blunted waveriders

        圖9 兩鈍邊緣乘波構(gòu)型邊緣處的受熱對比Fig.9 Heat flux along the leading edge of two blunted waveriders

        3.2 風(fēng)洞試驗結(jié)果

        為驗證非一致邊緣鈍化方法是否能夠確保熱流密度在很小范圍內(nèi)下降到很低水平,本文截取了圖2所示乘波構(gòu)型的頭部區(qū)域在FD-14A風(fēng)洞中進(jìn)行熱流密度的測量試驗。風(fēng)洞運行條件如表4所示。

        表4 風(fēng)洞的運行條件Table 4 Operation conditions of the wind tunnel

        圖10為試驗拍攝到的紋影照片??梢钥吹?,在試驗?zāi)P偷念^部區(qū)域出現(xiàn)一道弓形激波。弓形激波的出現(xiàn)是鈍化對理想乘波構(gòu)型流場結(jié)構(gòu)產(chǎn)生的最直接影響,也是鈍邊緣乘波構(gòu)型阻力增加的主導(dǎo)因素。

        圖10 試驗?zāi)P偷募y影圖Fig.10 Schlieren photo of the test model

        圖11給出了攻角-3°、0°、3°和6°條件下,試驗?zāi)P瓦吘壘€熱流密度的測量結(jié)果。圖中,熱流密度采用0°攻角時的測量值,長度采用圖2所示外形的總長度進(jìn)行了無量綱處理。試驗數(shù)據(jù)體現(xiàn)了非一致邊緣鈍化乘波構(gòu)型頭部區(qū)域受熱的兩個特征。第一,在駐點附近的很小范圍內(nèi),熱流密度即衰減到很低水平:如圖所示,在x/L=0.06位置,邊緣的熱流密度已減小到駐點值的1/4左右。第二,在小攻角飛行時,非一致邊緣鈍化乘波構(gòu)型頭部區(qū)域的受熱形勢基本不受攻角影響。(乘波構(gòu)型適宜采用小攻角飛行,因為其最大升阻比在0°附近獲得。)

        圖11 試驗測得的熱流密度沿模型邊緣分布Fig.11 Heat flux along the leading edge of test model obtained by experiment

        4 結(jié) 論

        (1)通過減少邊緣處的氣體泄露,非一致邊緣鈍化方法較有效地提高了乘波構(gòu)型的升力性能和升阻比,改善了傳統(tǒng)一致邊緣鈍化乘波構(gòu)型的氣動力性能。

        (2)非一致邊緣鈍化方法對乘波構(gòu)型上下表面受熱的影響很小;頭部區(qū)域是非一致邊緣鈍化乘波構(gòu)型熱防護(hù)系統(tǒng)設(shè)計的重點,該區(qū)域的氣動加熱情況基本不受攻角影響(小攻角飛行)。

        (3)對于理想乘波構(gòu)型,其邊緣的鈍化方法應(yīng)基于氣動力性能和表面受熱情況進(jìn)行優(yōu)化設(shè)計;選擇合適的鈍化方法,鈍邊緣的乘波構(gòu)型仍可作為近空間高超聲速飛行器的重要候選布局。

        致謝:感謝中國空氣動力研究與發(fā)展中心超高聲速所502室在實驗過程中給予的大力幫助。

        [1]KUCHEMANN D.The aerodynamic design of aircraft[M].Pergamon Press,1978.

        [2]NONWEILER T R F.Aerodynamic problems of manned space vehicles[J].JournaloftheRoyalAeronauticalSociety,1959,69(9):521-528.

        [3]FREDERICK F,TERRY C,STEPHEN A,et al.The development of waveriders from an axisymmetric flow field[C].AIAA 2007-847,2007.

        [4]CORDA S,ANDERSON J D Jr.Viscous optimized hypersonic waveriders designed from axisymmetric flow fields[C].AIAA-88-0369,1988.

        [5]LEWIS M J,HASTINGS D E.Bow shock matching with viscous effects on hypersonic fore bodies[C].AIAA Paper 89-2678,1989.

        [6]MANGIN B,BENAY R,CHANETZ B.Optimization of viscous waveriders derived from axisymmetric power-law blunt body flows[J].J.SpacecraftRockets,2006,43(5):990-998.

        [7]ANDERSON J D Jr,LEWIS M J.Hypersonic waveriders:where do we stand[R].AIAA 93-0399,1993.

        [8]BOWCUTT K G.Optimization of hypersonic waveriders derived from cone flows including viscous effects[D].Maryland:University of Maryland,1986.

        [9]CORDA S.Star-body waveriders with multiple design Mach numbers[J].J.SpacecraftRockets,2009,46(6):1178-1185.

        [10]YANKUI W,SHUI F Y,DONG J Z,et al.Design of waverider configuration with high lift-drag ratio[J].J.Aircraft,2007,44(1):144-148.

        [11]SANTOS W F N.Bluntness impact on lift-to-drag ratio of hypersonic wedge flow[J].JournalofSpacecraftandRockets,2009,46(2):329-339.

        [12]OHTA T,URANO A,MIYAGAWA T,et al.Flow visualization on lower surfaces of waverider configurations at Mach 5.5[R].AIAA 97-1884,1997.

        [13]SILVESTER T,MORGAN R.Computational hypervelocity aerodynamics of a caret waverider[R].AIAA 2004-3848,2004.

        [14]GILLUM M J,LEWIS M J.Experimental results on a Mach 14 waverider with blunt leading edges[J].JournalofAircraft,1997,34(3):296-303.

        [15]CAO D Y,LI C X.A numerical study on aerothermodynamic performances of a waverider based vehicle[J].JournalofAstronautics,2008,29(6):1782-1785.(in Chinese)曹德一,李椿萱.乘波飛行器氣動力、熱特性的數(shù)值模擬研究[J].宇航學(xué)報.2008,29(6):1782-1785.

        [16]CHEN X Q,HOU Z X,LIU J X,et al.Bluntness impact on performance of waverider[J].ComputersandFluids,2011,48(1):30-43.

        [17]VANMOL D O,ANDERSON J D Jr.Heat transfer characteristics of hypersonic waveriders with an emphasis on leading edge effects[R].AIAA 92-2920,1992.

        [18]LEEJ H,RHO O H.Numerical analysis of hypersonic viscous flow around a blunt body using Roe's FDS and AUSM +schemes[R].AIAA 97-2054,1997.

        [19]LI S X.Hypersonic flow characteristics of typical configurations[M].Publishing company of National Defense Industry,2007:64-65.(in Chinese)李素循.典型外形高超聲速流動特征[M].北京:國防工業(yè)出版社,2007:64-65.

        [20]DETRA R W,KEMP N N,RIDDELL F R.Heat transfer to satellite vehicles reentering the atmosphere[J].JetPropuls,1957,27(12):1256-1257.

        [21]HOLDEN M S,WIETING A R,MOSELLE J R,et al.Studies of aerothermal loads generated in regions of shock/shock interaction in hypersonic flow[R].AIAA Paper 88-0477.1998.

        [22]MENSHOV I S,NAKAMURA Y.Numerical simulation and experimental comparison for high-speed non-equilibrium air flows[J].FluidDynamicsResearch,2000,27(5):305-334.

        猜你喜歡
        駐點氣動力超聲速
        高超聲速出版工程
        高超聲速飛行器
        飛行載荷外部氣動力的二次規(guī)劃等效映射方法
        基于游人游賞行為的留園駐點分布規(guī)律研究
        中國園林(2018年7期)2018-08-07 07:07:48
        超聲速旅行
        側(cè)風(fēng)對拍動翅氣動力的影響
        利用遠(yuǎn)教站點,落實駐點干部帶學(xué)
        利用遠(yuǎn)教站點,落實駐點干部帶學(xué)
        2300名干部進(jìn)村“串戶”辦實事
        源流(2015年8期)2015-09-16 18:01:32
        高超聲速大博弈
        太空探索(2014年5期)2014-07-12 09:53:28
        国产精品理论片在线观看| 国产精品国产三级国产在线观 | 亚洲熟妇av乱码在线观看| 国产成人aa在线观看视频| 人妻少妇偷人精品视频| 国产精品久久精品第一页| 99热久久精里都是精品6| 在线观看一区二区女同| 亚洲粉嫩视频在线观看| 亚洲av永久无码精品一福利| 国产超碰人人模人人爽人人喊| 久久中文字幕久久久久| 人妻av在线一区二区三区| 国产免费久久精品99久久| 国产剧情麻豆女教师在线观看| 老熟女熟妇嗷嗷叫91| 日韩av一区二区蜜桃| 国模精品一区二区三区| 天天操夜夜操| 国产在线一区二区视频免费观看| 久久精品蜜桃亚洲av高清| 成 人 免费 在线电影| 国产日韩欧美在线| 亚洲人妻av在线播放| 亚洲最新无码中文字幕久久| 97久久精品午夜一区二区| 免费黄网站久久成人精品| 中文字幕色资源在线视频| 久久不见久久见免费影院国语| 国产精品美女| 中文字幕久久人妻av| 极品人妻被黑人中出种子| 精品无码一区二区三区亚洲桃色| 天堂AV无码AV毛片毛| 亚洲天堂久久午夜福利| 亚洲av无码精品蜜桃| 久久99久久久无码国产精品色戒| 日本女优免费一区二区三区| 欧洲熟妇色xxxx欧美老妇性| 香蕉人妻av久久久久天天| 亚洲毛片av一区二区三区|