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

        ?

        磁控?zé)岱雷o(hù)系統(tǒng)高溫流場與電磁場耦合計算方法

        2017-06-15 14:33:34劉偉強
        宇航學(xué)報 2017年5期
        關(guān)鍵詞:效率

        李 開,柳 軍,劉偉強

        (國防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院,長沙410073)

        磁控?zé)岱雷o(hù)系統(tǒng)高溫流場與電磁場耦合計算方法

        李 開,柳 軍,劉偉強

        (國防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院,長沙410073)

        為了研究真實氣體條件下霍爾效應(yīng)對磁流體(MHD)控制熱防護(hù)效果的影響,建立了熱化學(xué)非平衡流場、外加磁場、霍爾電場的耦合計算方法。分析了耦合計算效率與電場更新間隔、電勢殘差收斂極限的關(guān)系。給出了采用非平衡霍爾系數(shù)模型時電場更新間隔和電勢收斂極限的較優(yōu)取值。研究表明,當(dāng)霍爾系數(shù)較小(β=1.0)時,電場更新間隔大于100流場計算時間步時耦合計算效率較高,且導(dǎo)電壁面和絕緣壁面結(jié)論一致。當(dāng)霍爾系數(shù)較大時,耦合計算時間過長,可適當(dāng)增加電場迭代間隔和電勢收斂極限以提高耦合計算效率。

        磁流體控制;熱防護(hù);霍爾電場;耦合計算

        0 引 言

        臨近飛行器在高超聲速飛行時氣動熱問題日趨嚴(yán)重,其“熱障”問題已成為制約飛行器設(shè)計的瓶頸[1-3]。最近十年間高超聲速領(lǐng)域的新概念熱防護(hù)技術(shù)層出不窮。得益于近年來電磁流動控制技術(shù)、超導(dǎo)磁鐵技術(shù)的進(jìn)步[4-5],磁控?zé)岱雷o(hù)技術(shù)的現(xiàn)實可行性逐漸增強,應(yīng)用價值逐年提高,得到了各國研究者的普遍關(guān)注[6-7]。但值得注意的是,在典型軌道再入及超軌道再入飛行條件下,霍爾效應(yīng)明顯,感應(yīng)產(chǎn)生的霍爾電場會改變前緣流場電流分布,進(jìn)而改變洛倫茲力,影響磁控?zé)岱雷o(hù)系統(tǒng)的效率。磁控?zé)岱雷o(hù)系統(tǒng)的霍爾效應(yīng)研究涉及外加磁場、感應(yīng)電場和非平衡流場三場的耦合計算,是該領(lǐng)域研究的難點[8]。

        由于高超聲速飛行器前緣激波后電離度較低,低磁雷諾數(shù)假設(shè)往往可以得到滿足[9]。此時,外加磁場、感應(yīng)電場與非平衡流場的耦合可以通過耦合求解電勢泊松方程和非平衡磁流體方程實現(xiàn)。由于大霍爾系數(shù)條件下電場求解剛性嚴(yán)重、收斂性差[10],電磁場和熱化學(xué)非平衡流場的耦合計算效率較低。在現(xiàn)有文獻(xiàn)的研究中,耦合計算的研究可以分為兩大類:一類側(cè)重于構(gòu)建精度較高的霍爾電場求解方法[11-13];另一類側(cè)重于分析霍爾效應(yīng)對磁控系統(tǒng)的影響效果[10,14-15]。

        盡管上述研究大都實現(xiàn)了流場和霍爾電場的耦合計算,但現(xiàn)有方法對于拓?fù)浣Y(jié)構(gòu)復(fù)雜、網(wǎng)格規(guī)律大的流體域耦合計算的效率仍較低。因此,非平衡流場和電場耦合計算效率還需仔細(xì)研究。在該松耦合過程中有兩個重要的控制參數(shù),一是電場更新間隔S,即非平衡流場計算S時間步電場開始并完成一次迭代;二是電勢殘差收斂極限ε,即當(dāng)電勢殘差小于該值時判定電場收斂。不同霍爾系數(shù)下電場的收斂性差別較大,給兩參數(shù)取值帶來困難。文獻(xiàn)[16]采用的電場更新間隔在10~100之間,具體數(shù)值取決于算例的磁場強度。文獻(xiàn)[17]在進(jìn)行參數(shù)霍爾效應(yīng)研究時,對于不同的霍爾系數(shù)采用了同樣的電場更新間隔(S=200)和同樣的電勢收斂極限(ε=1×10-4)。可見,現(xiàn)有文獻(xiàn)尚未考慮不同霍爾系數(shù)下電場收斂性對耦合效率的影響,也并未確定不同霍爾系數(shù)下電場更新間隔和收斂極限的合理取值方法。本文將首先建立三維熱化學(xué)非平衡流場和霍爾電場的數(shù)值方法以及流場-電磁場之間的耦合計算方法,而后分析在不同霍爾系數(shù)下電場更新間隔和收斂極限對耦合計算效率的影響,以便為高超聲速飛行器前緣磁控?zé)岱雷o(hù)系統(tǒng)的霍爾效應(yīng)研究提供參考。

        1 數(shù)學(xué)模型

        1.1 控制方程

        針對高超聲速飛行器前緣流動特點,電磁場和流場耦合分析采用低磁雷諾數(shù)假設(shè)下的三維熱化學(xué)非平衡磁流體(Magnetohydrodynamic,MHD)控制方程,如式(1)所示。其中,熱化學(xué)非平衡模型采用Park雙溫度模型,化學(xué)反應(yīng)動力學(xué)模型采用11組元20反應(yīng)的Gupta模型。

        (1)

        式中:U為守恒變量,F(xiàn)、G、H分別為x、y、z方向的對流通量矢量,F(xiàn)v、Gv、Hv為三方向上的黏性通量項,W為化學(xué)反應(yīng)和振動能量源項矢量,具體表達(dá)式見文獻(xiàn)[18]。相對于熱化學(xué)非平衡黏流,上述方程增加了一個電磁源項通量SMFD,見下式

        SMFD=

        (2)

        式中:J=(Jx,Jy,Jz)為電流密度矢量,E為電場強度矢量,B=(Bx,By,Bz)為磁感應(yīng)強度矢量。γ表征不同非平衡模式間的電磁能量分配比,介于0和1之間,這里取γ=0.5[13]。磁感應(yīng)強度B在給定磁場形態(tài)后即可確定,這里采用磁偶極子磁場[14]??紤]霍爾效應(yīng)后,結(jié)合低磁雷諾數(shù)假設(shè)可以將對感應(yīng)場強矢量E的求解轉(zhuǎn)化為對標(biāo)量電勢φ的求解。由廣義歐姆定律(3)和電流守恒方程(4)

        (3)

        ▽·J=0

        (4)

        可以得到關(guān)于φ的電勢泊松方程

        (5)

        (6)

        式中:B為磁感應(yīng)強度矢量B的模。σ為標(biāo)量電導(dǎo)率,采用式(7)非平衡流電導(dǎo)率模型進(jìn)行計算。β為霍爾系數(shù)。為全面分析不同霍爾系數(shù)下耦合效率,采用兩種方法確定β:1)式(7)非平衡模型[14];2)均勻分布。

        (7)

        1.2 數(shù)值方法

        磁流體流動方程(1)離散時,對流項差分采用基于MUSCL插值的AUSMPW格式,隱式求解采用了LU-SGS格式,并且對化學(xué)反應(yīng)、振動以及電磁源項采用了隱式處理以削弱源項過大帶來的剛性從而提高收斂性。

        電勢泊松方程(5)的離散基于單元中心有限體積法。采用交替隱式近似因子分解法(ADI-AF),轉(zhuǎn)化為下式迭代求解

        (8)

        其中,系數(shù)a1i-1,b1i,c1i+1如下所示

        (9)

        此外,a2j-1,b2j,c2j+1,a3k-1,b3k,c3k+1形式類似,不再贅述。其中,

        (10)

        2 耦合方法與邊界條件

        流場邊界條件:流動入口為自由來流;出口采用超聲速外推條件;壁面采用全催化等溫壁,壁溫取決于試驗工況。霍爾電場邊界條件:絕緣壁面,J·n=0; 導(dǎo)電壁面,φ=0 V; 其余邊界為Neumann邊界,▽φ·n=0。

        3 仿真校驗

        3.1 霍爾電場校驗

        算例2為間隔電極霍爾效應(yīng)算例,如圖3所示[13]。±y的通道邊界內(nèi)通入+x向的導(dǎo)電流體,電導(dǎo)率σ=1Ω-1m-1。有限寬度的平行電極沿周期性重復(fù)布置在通道兩側(cè)壁面上。由于通道無限長,圖中的進(jìn)口和出口邊界設(shè)置為周期性邊界。為使通道內(nèi)流體速度場對電勢場結(jié)果無影響,則必須滿足▽×(u×B)=0,因此可假定u=f(y),B=f(z),且通道內(nèi)流動為充分發(fā)展的Poiseuille流動。其中,通道中心線上的流體速度umax=1m/s,通道半寬h=0.5m。

        當(dāng)B=0 T時,兩電極間只存在靜電場。圖4為無霍爾效應(yīng)時的本文和文獻(xiàn)[13]的電勢等值線對比圖。圖5為有霍爾效應(yīng)時(β=1.0、B=1T)時的本文與文獻(xiàn)的電流流線對比圖。從圖4~5可以看出,兩種情況下,本文計算結(jié)果與文獻(xiàn)結(jié)果吻合良好。

        3.2 磁流體氣動熱校驗

        氣動熱模擬的準(zhǔn)確性是檢驗耦合非平衡磁流體計算方法正確性的重要標(biāo)準(zhǔn)。選用日本1996年發(fā)射的軌道再入試驗飛行器(Orbital reentry experiment, OREX)返回艙前體作為計算模型,如圖6所示。選取再入飛行時間在7471.5s(Ma=17.61、H=59.6km)的飛行工況,并采用了有限催化壁面模型(γ=5.0×10-3)進(jìn)行了氣動熱數(shù)值模擬。圖7分別給出了三種外加磁場情況(B0=0.0、0.3、0.5 T)下的平動溫度和壁面熱流分布。通過與Fujino等[20]計算結(jié)果的比較可以看出,本文的激波脫體距離、壁面熱流計算結(jié)果和文獻(xiàn)[20]吻合良好。驗證了本文外加磁場作用下的非平衡流場及氣動熱數(shù)值模擬的準(zhǔn)確性。值得一提的是,圖7(a)中平動溫度峰值的差異以及圖7(b)中當(dāng)B0=0.3 T時在y/L=1.3(L為參考長度,L=1.0139m)附近的壁面熱流差異,可能與本文和文獻(xiàn)采用了不同的平動-振動松弛模型有關(guān)。

        4 計算工況

        為了提高耦合計算效率,需要確定電磁場參數(shù)相對于流動參數(shù)的最優(yōu)的更新方式,即研究電場更新間隔S、電場殘差收斂極限ε對耦合計算效率的影響。其中耦合計算效率以平均單步迭代時間ta,總收斂時長tc表示。

        計算模型網(wǎng)格和第2.2節(jié)相同。采用OREX飛行器7461.5s飛行工況(Ma=20.04,H=63.60km)。駐點磁感應(yīng)強度B0=0.2 T。流場計算采用當(dāng)?shù)貢r間步長,庫朗特數(shù)取0.2。鑒于霍爾系數(shù)不同對電場收斂性影響巨大,從而會對耦合計算效率產(chǎn)生重要影響。這里根據(jù)霍爾系數(shù)不同設(shè)計了兩組算例:1)均布霍爾系數(shù)β=1.0,步進(jìn)因子ap=0.002;2)Fujino霍爾系數(shù)模型,步進(jìn)因子ap=0.006。前者霍爾系數(shù)相對較小,電場收斂快;后者霍爾系數(shù)相對較大(β≈5.0~10.0),電場收斂較慢。耦合計算以忽略霍爾效應(yīng)的流場收斂解作為初場。

        5 結(jié)果與分析

        表1給出了均布霍爾系數(shù)β=1.0時不同電場更新間隔下的單步平均迭代時間和收斂時的總迭代時間。從表1可以看出,和S=10相比,電場更新間隔S≥100時,單步迭代時間和收斂總時間大幅減小,而后隨著電場更新間隔的增加,單步平均迭代時間和收斂總時間雖略微減少但整體變化不大。導(dǎo)電壁面和絕緣壁面規(guī)律一致。

        表1 不同電場更新間隔下的迭代時間對比Table 1 Iteration time under different updating intervals of electric field (β=1.0)

        圖8(a)和(b)分別為導(dǎo)電壁面不同電場更新間隔下流場和電場的收斂曲線。從圖8可以看出,電場更新間隔對流場收斂曲線形狀影響較小,隨電場更新間隔的增加收斂步數(shù)略微減少。電場更新間隔對電勢場收斂影響較大,并且電場更新間隔越大,相同電場迭代步數(shù)流場收斂程序越高,相應(yīng)的電勢收斂所需的電場虛擬時間步越少。

        表2給出了采用非平衡霍爾系數(shù)模型時不同電場更新間隔下的單步平均迭代時間和收斂總時間。表2還對比了兩種電勢收斂殘差極限下的結(jié)果??梢钥闯?,電勢殘差收斂極限為10-4時,電場收斂耗時較多,單步平均迭代時間較長(S=10時約為58s),耦合計算時間呈現(xiàn)先降后升的趨勢,最優(yōu)的電場更新間隔為1000。而當(dāng)電勢殘差收斂極限為10-3時,在S=10~5000范圍內(nèi),電場更新間隔越大,單步平均迭代時間和總收斂時間越小,同β=1.0時規(guī)律一致。另外,當(dāng)S=1000時,兩個收斂極限(ε=10-3、10-4)下的總收斂時間相近,但電場間隔越小二者收斂時間相差越大,S=10時總收斂時間10-4可達(dá)10-3的110倍,此時采用相對較大的收斂殘差極限(10-3)可以極大地提高耦合計算效率。

        表2 不同電場更新間隔下的迭代時間對比Table 2 Iteration time under different updating intervals of electric field (nonequilibrium Hall parameter model)

        5 結(jié) 論

        以高速飛行器磁控?zé)岱雷o(hù)系統(tǒng)為研究對象,針對其在真實氣體條件下高溫流場和電磁場的耦合求解問題,建立了非平衡流場、外加磁場和霍爾電場的松耦合計算方法,在此基礎(chǔ)上分析了耦合計算效率與電場更新間隔、電勢殘差收斂極限的關(guān)系。研究表明,霍爾系數(shù)較小(β=1.0)時,電場更新間隔S=100相比S=10單步迭代時間和收斂總時間大幅減少。S>100后隨電場更新間隔的增加略微減少,計算效率變化不大,并且導(dǎo)電壁面和絕緣壁面規(guī)律一致?;魻栂禂?shù)較大時,電場收斂緩慢,耦合計算時間過長,可以通過適當(dāng)增加電場迭代間隔以及提高電勢收斂極限的方法有效提高耦合計算效率。采用霍爾系數(shù)非平衡模型時,電場更新間隔的建議取值為1000,電勢收斂極限可取10-3。

        [1] 孟松鶴,楊強,霍施宇,等. 一體化熱防護(hù)技術(shù)現(xiàn)狀和發(fā)展趨勢[J]. 宇航學(xué)報, 2013, 34(10): 1295-1302. [Meng Song-he, Yang Qiang, Huo Shi-yu, et al. State of arts and trend of integrated thermal protection systems [J]. Journal of Astronautics, 2013, 34(10): 1295-1302.]

        [2] 李鋒,艾邦成,姜貴慶. 一種熱平衡等溫機制的新型熱防護(hù)及相關(guān)技術(shù)[J]. 宇航學(xué)報, 2013, 34(12): 1644-1650. [Li Feng, Ai Bang-cheng, Jiang Gui-qing. A new thermal protection technology based on heat balance isothermal mechanism [J]. Journal of Astronautics, 2013, 34(12): 1644-1650.]

        [3] 潘勇,王江峰,伍貽兆. 非結(jié)構(gòu)網(wǎng)格高超聲速MHD流場逆風(fēng)格式數(shù)值模擬[J]. 宇航學(xué)報, 2008, 29(1): 104-109. [Pan Yong, Wang Jiang-feng, Wu Yi-zhao. Numerical simulation of hypersonic MHD flows using upwind scheme on unstructured grids [J]. Journal of Astronautics, 2008, 29(1): 104-109.]

        [4] 張康平,丁國昊,田正雨,等.磁流體動力學(xué)控制二維擴(kuò)壓器流場數(shù)值模擬研究[J].國防科技大學(xué)學(xué)報, 2009, 31(6):39-41. [Zhang Kang-ping, Ding Guo-hao, Tian Zheng-yu, et al. Numerical simulation of magnetohydrodynamic (MHD) control on 2D diffuser′s flow field [J]. Journal of National University of Defense Technology, 2009, 31(6):39-41.]

        [5] 田正雨,張康平,潘沙,等.磁流體動力學(xué)斜激波控制數(shù)值模擬分析[J]. 力學(xué)季刊, 2008, 29(1):72-77. [Tian Zheng-yu, Zhang Kang-ping, Pan Sha, et al. Numerical investigation and analysis for MHD oblique shock control[J]. Chinese Quarterly of Mechanics, 2008, 29(1):72-77. ]

        [6] Bityurin V A, Bocharov A N. Study of catalytic effects at reentry vehicle [C].The 52nd Aerospace Sciences Meeting, National Harbor, Maryland, USA, January 13-17, 2014.

        [7] Cristofolini A, Borghi C A, Neretti G, et al. MHD interaction around a blunt body in a hypersonic unseeded air flow[C].The 18th AIAA/3AF International Space Planes and Hypersonic Systems and Technologies Conference, Tours, France, September 24-28, 2012.

        [8] 劉深深,桂業(yè)偉,唐偉,等.一種多場耦合數(shù)據(jù)傳遞新方法[J]. 宇航學(xué)報, 2016, 37(1):61-67. [Liu Shen-shen, Gui Ye-wei, Tang Wei, et al. A new data transfer method in fluid-thermal-structrure coupling problems [J]. Journal of Astronautics, 2016, 37(1): 61-67.]

        [9] 呂浩宇, 李椿萱.Hall效應(yīng)對磁流體壓縮管道電磁流動特性的影響[J].科學(xué)通報, 2010, 55(12):1182-1188. [Lv Hao-yu, Li Chun-xuan. Influence of Hall effects on characteristics of magnetohydrodynamic converging channel [J]. Chin. Sci. Bull, 2010, 55(12):1182-1188.]

        [10] 胡海洋, 楊云軍, 周偉江. 大霍爾系數(shù)下電離氣體與磁場相互作用規(guī)律數(shù)值研究[J]. 力學(xué)學(xué)報,2011, 43(3): 453-460. [Hu Hai-yang, Yang Yun-jun, Zhou Wei-jiang. Numerical research on the interaction between ionized gas and magnetic field under high Hall parameter [J]. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(3):453-460.]

        [11] Gaitonde D V, Poggie J. Elements of a numerical procedure for 3-D MGD flow control analysis [C]. 40th AIAA Aerospace Sciences Meeting and Exhibit, Reno, USA, January 14-17, 2002.

        [12] Wan T, Candler G V, Macheret S O, et al. Three-dimensional simulation of the electric field and magnetohydrodynamic power generation during reentry [J]. AIAA Journal, 2009, 47(6):1327-1336.

        [13] Bisek N J. Numerical study of plasma-assisted aerodynamic control for hypersonic vehicles [D]. Michigan: The University of Michigan, 2010.

        [14] Fujino T, Matsumoto Y, Kasahara J, et al. Numerical studies of magnetohydrodynamic flow control considering real wall electrical conductivity [J]. Journal of Spacecraft and Rockets, 2007, 44(3):625-632.

        [15] 呂浩宇,李椿萱,董海濤. 三維超聲速磁流體發(fā)生器的流動特性[J]. 中國科學(xué)G輯,2009,39(3):435-445. [Lv Hao-yu, Li Chun-xuan, Dong Hai-tao. Characterization of the three-dimensional supersonic flow for the MHD generator [J]. Science in China Series G-Physics Mechan. Astron. , 2009, 39(3):435-445.]

        [16] Bisek N J, Gosse R, Poggie J. Computational study of impregnated ablator for improved magnetohydrodynamic heat shield [J]. Journal of Spacecraft and Rockets, 2013, 50(5): 927-935.

        [17] Otsu H, Konigorski D, Abe T. Influence of Hall effect on electrodynamic heat shield system for reentry vehicles [J]. AIAA Journal, 2010, 48(10):2177-2186.

        [18] 柳軍.熱化學(xué)非平衡流及其輻射現(xiàn)象的實驗和數(shù)值計算研究[D].長沙:國防科技大學(xué),2004. [Liu Jun. Experimental and numerical research on thermo-chemical nonequilibrium flow with radiation phenomenon [D]. Changsha: National University of Defense Technology, 2004.]

        [19] Gnoffo P A, Gupta R N, Shinn J L. Conservation equations and physical models for hypersonic air flows in thermal and chemical nonequilibrium [R]. NASA, TP-2867, 1989.

        [20] Fujino T, Ishikawa M. Numerical simulation of control of plasma flow with magnetic field for thermal protection in earth reentry flight [J]. IEEE Transactions on Plasma Science, 2006,34(2):409-420.

        通信地址:湖南省長沙市國防科技大學(xué)航天科學(xué)與工程學(xué)院(410073)

        電話:15616246843

        E-mail:likai898989@126.com

        (編輯:牛苗苗)

        Numerical Methods of Coupling High Temperature Flow Field with Electro Magnetic Field for Magnetohydrodynamic Heat Shield System

        LI Kai, LIU Jun, LIU Wei-qiang

        (College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China)

        Under the real flight conditions, the Hall effect usually cannot be neglected which may affect the performance of the magnetohydrodynamics(MHD) thermal protection. Therefore, the coupled numerical simulations are conducted for the thermochemical nonequilibrium flow field, the externally applied magnetic field and the induced electric field around a typical hypersonic vehicle. After that, the relations between the coupled calculation efficiency and the updating interval of the electric field as well as the converging residual limit of the electric potential are analyzed. Results show that the coupled calculation efficiency is higher under the condition that the updating interval is greater than 100 flow steps during a relatively smaller Hall parameter. However, during a larger Hall parameter, the coupling efficiency is lower but can be improved by increasing the updating interval of the electric field and potential residual limit. The suggested updating interval and potential residual limit are also given while the nonequilibrium Hall parameter model is utilized.

        MHD flow control; Thermal protection; Hall electric field; Coupled computation

        2016-09-26;

        2017-02-27

        國家自然科學(xué)基金(90916018);湖南省自然科學(xué)基金重點資助項目(13JJ2002)

        V211.1

        A

        1000-1328(2017)05-0474-07

        10.3873/j.issn.1000-1328.2017.05.005

        李 開(1989-),男,博士生,主要從事磁流體流動控制,高溫氣體動力學(xué)以及飛行器熱防護(hù)系統(tǒng)設(shè)計方面的研究。

        猜你喜歡
        效率
        你在咖啡館學(xué)習(xí)會更有創(chuàng)意和效率嗎?
        提升朗讀教學(xué)效率的幾點思考
        甘肅教育(2020年14期)2020-09-11 07:57:42
        注意實驗拓展,提高復(fù)習(xí)效率
        效率的價值
        商周刊(2017年9期)2017-08-22 02:57:49
        引入“倒逼機制”提高治霾效率
        質(zhì)量與效率的爭論
        跟蹤導(dǎo)練(一)2
        提高食品行業(yè)清潔操作的效率
        OptiMOSTM 300V提高硬開關(guān)應(yīng)用的效率,支持新型設(shè)計
        “錢”、“事”脫節(jié)效率低
        AV无码中文字幕不卡一二三区| 一色桃子中文字幕人妻熟女作品| 精品少妇人妻av无码专区| 亚洲精品中文字幕无乱码麻豆| 亚洲va精品va国产va| 日本一区二区三区光视频| 久久久www成人免费毛片| 亚洲妓女综合网99| 日本最新一区二区三区免费看| 中文字幕女同人妖熟女| 久久97久久97精品免视看| 成全视频高清免费| 国产自产拍精品视频免费看| 国产精品老熟女乱一区二区| 中文字幕乱码熟妇五十中出| 成人无码午夜在线观看| 中文字幕一二区中文字幕| 中文字幕影片免费人妻少妇 | 国产网红一区二区三区| 在线人成视频播放午夜| 军人粗大的内捧猛烈进出视频| 亚洲 欧美 激情 小说 另类| 亚洲av高清一区二区| 西西午夜无码大胆啪啪国模| 日韩精品无码一区二区三区免费| 资源在线观看视频一区二区| 日本系列中文字幕99| 樱桃视频影视在线观看免费| 国产精品丝袜在线不卡| av国产免费在线播放| 亚洲av综合av国产av中文| 国产精品久久久久免费a∨| 99精品国产av一区二区| 国产极品少妇一区二区| 亚洲精品无码久久久久| 午夜无码熟熟妇丰满人妻| 色播视频在线观看麻豆 | 第九色区Aⅴ天堂| 中文字幕免费人成在线网站| а√资源新版在线天堂| 永久免费看免费无码视频|