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

        ?

        基于體積力模型的1維渦輪特性評估方法

        2016-09-23 03:37:48鄒正平姚李超
        航空發(fā)動機(jī) 2016年2期
        關(guān)鍵詞:冷氣渦輪氣動

        何 杰,鄒正平,姚李超

        (1.北京航空航天大學(xué)航空發(fā)動機(jī)氣動熱力國防重點(diǎn)實(shí)驗(yàn)室;2.先進(jìn)航空發(fā)動機(jī)協(xié)同創(chuàng)新中心:北京100191)

        基于體積力模型的1維渦輪特性評估方法

        何杰1,2,鄒正平1,2,姚李超1,2

        (1.北京航空航天大學(xué)航空發(fā)動機(jī)氣動熱力國防重點(diǎn)實(shí)驗(yàn)室;2.先進(jìn)航空發(fā)動機(jī)協(xié)同創(chuàng)新中心:北京100191)

        渦輪低維氣動特性評估方法的預(yù)測精度對渦輪氣動設(shè)計有重要意義?;谶m用于軸流渦輪特性評估的1維體積力方法,分析和局部調(diào)整了A dam和Leonard的1維體積力模型中的葉片力源項(xiàng)和離心力源項(xiàng),并在特性計算中引入渦輪損失模型。通過1個低壓渦輪和1個含冷氣高壓渦輪對該方法進(jìn)行驗(yàn)證,計算結(jié)果表明:該方法對2個算例在渦輪設(shè)計工況點(diǎn)的流量效率預(yù)測偏差均在1%以內(nèi),且能夠反映渦輪氣動特性隨膨脹比的變化趨勢,具有良好的精度,可用于快速可靠地評估渦輪低維氣動特性。

        渦輪;低維設(shè)計;體積力模型;歐拉方程;特性評估

        0 引言

        渦輪部件氣動設(shè)計是1個從低維到高維逐步設(shè)計和優(yōu)化的過程,低維氣動設(shè)計是高維空間工作的基礎(chǔ),并在很大程度上決定了渦輪的性能水平,是設(shè)計流程中1個非常關(guān)鍵的環(huán)節(jié)。渦輪低維氣動特性評估方法可在渦輪設(shè)計初期預(yù)測工作特性,其預(yù)測精度對渦輪氣動設(shè)計具有重要意義。

        渦輪低維氣動特性評估常用方法為1維平均中徑法,通過分析渦輪內(nèi)部氣動損失機(jī)理并總結(jié)損失經(jīng)驗(yàn)關(guān)系式,建立計算軸流渦輪葉柵通道內(nèi)氣動損失的模型,以定量評估渦輪1維平均中徑上性能。在工程領(lǐng)域常應(yīng)用的AMDC損失模型及AMDCKO模型,可較好評估渦輪設(shè)計點(diǎn)和非設(shè)計點(diǎn)的性能[1-2]。其他比較著名的經(jīng)驗(yàn)?zāi)P腿鏣raupel模型[3]、Craig&Cox模型[4]等也能較好地預(yù)測渦輪氣動損失。平均中徑法在過去應(yīng)用廣泛,但隨著渦輪設(shè)計技術(shù)的發(fā)展和設(shè)計要求的變化,該方法被用來考慮渦輪中一些新參數(shù)的影響規(guī)律,如Coull和Hodson等[5]分別分析了流量系數(shù)、升力系數(shù)和雷諾數(shù)等參數(shù)對低壓渦輪性能的影響規(guī)律等;在中國,如姚李超等[6]發(fā)展了以載荷系數(shù)、反力度、軸向速比等多種設(shè)計變量為約束條件求取渦輪效率的多級低壓渦輪1維氣動優(yōu)化程序。平均中徑法對渦輪性能的預(yù)測能力和精度取決于所使用的損失模型精準(zhǔn)度和適用范圍,對使用者要求較高。

        除了平均中徑法,利用歐拉方程與體積力相結(jié)合對葉輪機(jī)性能評估的方法近年來也得到重視,該方法適用于1~3維和非定常等不同維數(shù)葉輪機(jī)性能的評估。高維體積力模型如Adamcyz提出的通道平均多級流動模型[7],通過體積力模型來模擬葉片對氣流的作用力以求解壓氣機(jī)中的3維流動。高維體積力模型還有其他形式,如適用于3維非定常畸變問題的Gong的周向平均體積力模型[8]和Xu的黏性體積力模型[9]等。低維體積力模型如Adam和Leonard提出的1維模型[10-12],該模型在低維歐拉方程中引入體積力源項(xiàng),采用有限體積法求解該方程組從而得到葉輪機(jī)流道平均中徑上的氣動參數(shù),用于確定葉輪機(jī)在不同工況下的氣動性能,以評估葉輪機(jī)械優(yōu)劣。在此基礎(chǔ)上,中科院工程熱物理研究所和北京航空航天大學(xué)也進(jìn)行了相應(yīng)的研究,并發(fā)展了多級軸流透平氣動特性計算程序[13-14]。

        本文在Adam和Leonard的體積力模型的基礎(chǔ)上,對1維渦輪特性評估方法進(jìn)行了研究,分析調(diào)整了模型中的葉片力源項(xiàng)和離心力源項(xiàng),在特性計算中引入了損失模型,通過雙級低壓渦輪和氣冷高壓渦輪的性能結(jié)果比較對該方法進(jìn)行驗(yàn)證。

        1 計算方法

        1.1控制方程

        渦輪流道截面如圖1所示。本文所建立體積力的源項(xiàng)在空間的分布基于均勻的軸對稱特性,在1維評估中只考慮氣動參數(shù)沿流向的變化,忽略了流道法向氣動參數(shù)變化,由此建立曲線坐標(biāo)系(m,n,θ),其中m為流向坐標(biāo),n為法向坐標(biāo),θ為周向坐標(biāo)。微元控制體如圖2所示。控制體上邊界為機(jī)匣,下邊界為輪轂,假定進(jìn)口截面為截面①,出口截面為截面②,建立守恒的歐拉微分方程組

        其中

        式中:ρ為流體的密度;S為子午流道的橫截面積;t為時間;Vm為流體的絕對速度子午流向分量;Vθ為流體的絕對速度周向流向分量;r為到渦輪中心軸線的半徑;p為流體的靜壓;e°為流體的總內(nèi)能;Qb為無黏葉片力源項(xiàng);Qf為黏性力源項(xiàng);Qg為幾何力源項(xiàng);Qc為冷氣源項(xiàng)。

        圖1 研究坐標(biāo) 

        圖2 控制體[10]

        1.2體積力模型

        基于Adam和Leonard的思路[10],本文局部修正了其體積力模型中無黏葉片力切向分量:體積力模型采取在歐拉方程上添加源項(xiàng)的方法,確定源項(xiàng)對于預(yù)測精度和能力至關(guān)重要,為了分析源項(xiàng),給定葉片平均中徑進(jìn)出口速度三角形,其參數(shù)如圖3所示。

        圖3 速度三角形

        r1、r2分別為葉片進(jìn)、出口平均中徑,對應(yīng)的周向速度分別為Vθ1、Vθ2,qm為假定通道流量,Vol為葉片排內(nèi)氣體體積。求解源項(xiàng)需要分析葉片黏性力、無黏葉片力、離心力、幾何流道形式和冷氣影響。

        黏性力源項(xiàng)Qf考慮黏性力的影響,作用于葉片區(qū)域的黏性力使得流體沿流線產(chǎn)生熵增,利用AMDCKO損失模型求解氣流通過葉柵后的熵增Δs,該模型由AMDC損失模型改進(jìn)得到,不僅適用于低亞聲渦輪葉柵,也同樣適用于高馬赫數(shù)的渦輪葉柵。由于葉柵通道幾何參數(shù)未知,故近似認(rèn)為氣流在葉片區(qū)域內(nèi)的單位流向坐標(biāo)上的熵增相等,該假設(shè)在設(shè)計點(diǎn)考慮是合理的[15-16],在變工況下的熵增分布有待進(jìn)一步優(yōu)化改進(jìn)。從該角度出發(fā),可得黏性力矢量F→f、及其流向分量和周向分量分別為

        分析動量矩方程的葉片力源項(xiàng)Qb,Adam和Leonard認(rèn)為動量矩方程中力矩起主導(dǎo)作用的為無黏葉片力力矩[10],故只考慮無黏葉片力源項(xiàng)。由于渦輪設(shè)計工況點(diǎn)流動的主控因素是位勢作用,該假設(shè)在設(shè)計點(diǎn)及設(shè)計點(diǎn)附近的狀態(tài)點(diǎn)是非常合理的,但在1維特性評估過程中,要合理評估非設(shè)計狀態(tài)的性能,在遠(yuǎn)離設(shè)計點(diǎn)的工況下,黏性的影響將越來越大,成為不應(yīng)被忽略的因素。本文在動量矩方程中將綜合考慮無黏葉片力力矩和黏性力力矩的作用,以適用于更寬廣的工作范圍。利用角動量守恒方程求解切向力總和Fθ,可得單位體積上氣體所受到的無黏葉片力切向分量Fb,θ和流向分量Fb,m分別為

        幾何力源項(xiàng)Qg包含了離心力和渦輪子午流道形式的影響。Adam和Leonard文獻(xiàn)中以軸向?yàn)榛鶞?zhǔn)坐標(biāo),將離心力分解為軸向和周向分量[10],本文以子午流線為基準(zhǔn)坐標(biāo),離心率如圖4所示。僅考慮葉片旋轉(zhuǎn)引起的離心力,沒有離心力周向分量,即離心力只作用在軸向動量方程,不影響切向動量方程。

        圖4 離心力

        渦輪子午流道形式指輪轂和機(jī)匣之間的流道面積變化,具體為垂直于平均子午流線的橫截面面積以及平均中徑的變化所帶來的影響。

        冷氣源項(xiàng)Qc主要與冷氣有關(guān),本節(jié)中所建立的1維體積力模型為簡化考慮,只考慮在質(zhì)量守恒方程和能量守恒方程中添加源項(xiàng)冷氣,即冷氣的加入會改變通道氣體流量和焓的變化,而不影響在流向的動量守恒方程和角動量守恒方程,該假設(shè)在冷氣量所占比例較小時是合理的。

        綜上所述,各體積力源項(xiàng)分別為

        1.3數(shù)值方法

        采用時間推進(jìn)的有限體積法求解歐拉控制方程,其數(shù)值格式為迎風(fēng)型矢通量分裂格式AUSM[17],無黏通量F分裂為對流項(xiàng)和壓力項(xiàng)。

        程序采用定比熱方法計算,計算時需要給出以下條件:渦輪子午面輪轂和機(jī)匣型線幾何數(shù)據(jù)、各葉片排前尾緣在子午面內(nèi)的坐標(biāo)、各葉片排個數(shù)以及進(jìn)出口幾何構(gòu)造角、冷氣參數(shù)、渦輪轉(zhuǎn)速以及邊界條件。邊界條件包括:進(jìn)口總溫和總壓、進(jìn)口切向氣流角以及出口靜壓。在計算過程中通過改變渦輪轉(zhuǎn)速和出口靜壓以得到渦輪在不同非設(shè)計工況下的工作特性。

        2 計算結(jié)果與討論

        選取2級低壓渦輪和冷氣單級高壓渦輪2個算例渦輪對程序的正確性進(jìn)行驗(yàn)證,同時與3維RANS計算結(jié)果進(jìn)行對比。

        2.12級低壓渦輪

        某型2級低壓渦輪的流道形式如圖5所示。渦輪進(jìn)口總溫為1100 K,進(jìn)口總壓為357 kPa,設(shè)計轉(zhuǎn)速為11500 r/min,軸向進(jìn)氣、總冷氣量較少,小于2%。計算得到不同轉(zhuǎn)速下渦輪流量和效率隨膨脹比變化的曲線及其與3維計算結(jié)果的對比曲線,如圖6所示。

        圖5 低壓渦輪通道

        Adam和Leonard的1維模型、修正后的1維模型和3維模型計算的110%,100%、90%、80%折合轉(zhuǎn)速下低壓渦輪出口折合流量、效率隨膨脹比的變化如圖6所示。在流量-膨脹比特性圖中,以3維RANS結(jié)果為基準(zhǔn),修正后的1維模型比原1維模型對渦輪流量的預(yù)測更為精準(zhǔn),對流量堵點(diǎn)出現(xiàn)的位置預(yù)測更加精確;在效率-膨脹比特性圖中,修正后的1維模型和原1維模型計算結(jié)果分別與3維結(jié)果對比,修正的1維模型計算結(jié)果與3維結(jié)果更加貼近。

        圖6 低壓渦輪流量、效率隨膨脹比和折合轉(zhuǎn)速的變化

        修正后的1維體積力模型計算結(jié)果與3維RANS計算結(jié)果中的渦輪流量和效率的變化趨勢基本相同,在不同的折合轉(zhuǎn)速下,利用修正后的1維體積力模型計算的渦輪的流量和效率隨著膨脹比的變化特性與3維結(jié)果基本相同。

        不同轉(zhuǎn)速和膨脹比下的流量和效率偏差如圖7所示。對比分析了各狀態(tài)點(diǎn)的預(yù)測精度。流量偏差Δm和效率偏差Δη定義如下

        圖7 低壓渦輪不同轉(zhuǎn)速和膨脹比下的流量和效率偏差

        2.2冷氣單級高壓渦輪

        高壓渦輪進(jìn)口總溫為1600 K,進(jìn)口總壓為1529 kPa,設(shè)計轉(zhuǎn)速為18500 r/min,軸向進(jìn)氣,導(dǎo)葉冷氣量為渦輪進(jìn)口質(zhì)量流量的10.4%,動葉冷氣量為3.8%,該算例中渦輪冷氣量較大,應(yīng)以考核體積力的方法對氣冷渦輪性能適應(yīng)性進(jìn)行評估。高壓渦輪的流道和冷氣摻混比例如圖8所示。

        圖8 高壓渦輪的流道和冷氣摻混比例

        Adam和Leonard的1維模型、修正后的1維模型和3維計算的100%、90%、80%折合轉(zhuǎn)速下高壓渦輪出口折合流量、效率隨膨脹比的變化如圖9所示。其中冷氣渦輪效率定義見文獻(xiàn)[18]。從圖9(a)中可見,修正后的1維模型計算結(jié)果中流量堵塞點(diǎn)出現(xiàn)的位置更符合渦輪真實(shí)工作情況,修正后的1維模型的計算結(jié)果比原模型的更貼近3維計算結(jié)果;從圖9 (b)中可見,在膨脹比為4.2的設(shè)計點(diǎn)處,修正后的1維體積力模型比原模型對效率的預(yù)測更加精確。說明了在不同的折合轉(zhuǎn)速下,采用修正后的1維模型計算的渦輪流量和效率隨著膨脹比的變化趨勢與3維結(jié)果基本相同,修正后的1維模型比原模型對流量特性預(yù)測的精度更高。

        圖9 高壓渦輪流量和效率隨膨脹比和折合轉(zhuǎn)速的變化

        高壓渦輪不同轉(zhuǎn)速和膨脹比下的流量和效率偏差如圖10所示。從圖中可見,修正后的1維結(jié)果與3維結(jié)果的數(shù)值偏差較小,設(shè)計點(diǎn)處流量偏差和效率偏差均小于1%,非設(shè)計點(diǎn)處流量偏差和效率偏差較小,修正后的1維結(jié)果可信度較高。針對冷氣量較多的高壓渦輪,修正后的1維模型仍然可以較為精準(zhǔn)地評估渦輪特性。

        圖10 高壓渦輪不同轉(zhuǎn)速和膨脹比下的流量和效率偏差

        在2個算例中,在偏離設(shè)計點(diǎn)狀態(tài)較遠(yuǎn)處都出現(xiàn)渦輪流量效率偏差較大的現(xiàn)象,說明體積力模型在非設(shè)計點(diǎn)處的預(yù)測尚需進(jìn)一步修正,在非設(shè)計點(diǎn)處氣動損失模型的精度決定著源項(xiàng)的偏差,對體積力模型的精度有一定影響。如何進(jìn)一步提高體積力模型在非設(shè)計點(diǎn)下的預(yù)測精度是未來的工作方向。

        3 結(jié)論

        本文結(jié)合修正后的1維體積力模型,改進(jìn)了計算程序,并對2個算例進(jìn)行渦輪特性評估,得出以下結(jié)論:

        (1)修正后的1維體積力模型比原模型對流量特性預(yù)測精度更高。能夠在設(shè)計點(diǎn)處比較準(zhǔn)確預(yù)測流量效率特性,在非設(shè)計點(diǎn)處能夠較為準(zhǔn)確預(yù)測流量特性,且在非設(shè)計點(diǎn)處性能偏差在工程允許的誤差范圍內(nèi),能夠預(yù)測流量效率特性隨膨脹比的變化趨勢。

        (2)低維體積力方法能夠快速準(zhǔn)確地對渦輪特性進(jìn)行評估,適用于高壓低壓、單級多級、含冷氣不含冷氣的軸流渦輪。

        (3)低維體積力方法的預(yù)測精度取決于體積力源項(xiàng)的?;?,進(jìn)一步提高體積力模型在渦輪非設(shè)計狀態(tài)下的預(yù)測精度是未來的研究方向。

        [1]Dunham J,Came P M.Improvements to the Ainley-Mathieson method of turbine performance prediction[J].Journal of Engineering for Gas Turbines and Power,1970,92(3):252-256.

        [2]Kacker S C,Okapuu U.A mean line prediction method for axial flow turbine efficiency[J].Journal of Engineering for Gas Turbines and Power,1982,104(1)1:111-119.

        [3]TraupelW.Termischeturbomaschinenzweiterbandgel?nderte getriebsbe dingungen,regelung,mechanische probleme,temperatur problem [M].New York:Springer-Verlag Berlin Heidelberg,1977:73-77.

        [4]Craig H R M,Cox H J A.Performance estimation of axial flow turbines [J].Proceedings of the Institution of Mechanical Engineers,1970,185 (1):407-424.

        [5]Coull J D,Hodson H P.Blade loading and its application in the Mean-Linedesignoflowpressureturbines[J].Journalof Turbomachinery,2013,135(2):21-32.

        [6]姚李超,鄒正平,張偉昊,等.基于粒子群算法的多級低壓渦輪一維氣動優(yōu)化設(shè)計方法[J].推進(jìn)技術(shù),2013,34(8):1042-1049.

        YAO Lichao,ZOU Zhengping,LIU Huoxing.An optimal method of one-dimensional design for multistage low pressure turbines based on particle swarm optimization[J].Journal of Propulsion Technology,2013,34(8):1042-1049.(in Chinese)

        [7]Adamczyk J J.Model equation for simulating flows in multistage turbomachinery[J].Lecture Series-van Kareman Institute for Fluid Dynamics,1996,5(1):1-28.

        [8]Gong Y,Greitzer E M,Tan C S,et al.A computational model for short-wavelength stall inception and development in multi stage compressors[J].Journal of Turbomachinery,1999,121(4):726-734.

        [9]Xu L.Assessing viscous body forces for unsteady calculations[R]. ASME 2002-GT-30359.

        [10]Adam O,Léonard O.A quasi-one dimensional model for axial compressors[R].ISABE 2005-1011.

        [11]Adam O,Léonard O.A quasi-one dimensional model for axial turbines[R].ISABE 2007-1215.

        [12]Adam O,Léonard O.A quasi-one dimensional model for multistage turbomachines[J].Journal of Thermal Science,2008,17(1):7-20.

        [13]張曉東,劉建軍,安柏濤.多級軸流透平氣動特性計算方法研究[J].工程熱物理學(xué)報,2011,32(4):569-572.

        ZHANG Xiaodong,LIU Jianjun,AN Baitao.Investigation on the computational method of axial turbine aerodynamic characteristics based on 1.5 dimensional Euler equations[J].Journal of Engineering Thermophysics,2011,32(4):569-572.(in Chinese)

        [14]姚李超.多級低壓渦輪低維氣動優(yōu)化設(shè)計技術(shù)研究[D].北京:北京航空航天大學(xué),2013.

        YAO Lichao.Studies on an optimal method of low dimensional aerodynamic design for multistage low pressure turbine[D].Beijing:Beihang University,2013.(in Chinese)

        [15]Damle S V,Dang T Q,Reddy D R.Through flow method for turbomachinesapplicableforallflow regimes[J].Journalof Turbomachinery,1997 119(2):256-262.

        [16]鄭寧,鄒正平,徐力平.風(fēng)扇進(jìn)氣畸變?nèi)S非定常數(shù)值模擬技術(shù)研究[J].航空動力學(xué)報,2007,22(1):60-65.

        ZHENG Ning,ZOU Zhengping,XU Liping.3D Unsteady numerical simulation of fan/compressor with inlet distortion[J].Journal of Aerospace Power,2007,22(1):60-65.(in Chinese)

        [17]Liou M S.A sequel to ausm:ausm+[J].Journal of Computational Physics,1996,129(2):364-382.

        [18]Hartsel J E.Prediction of effects of mass-transfer cooling on the blade-row efficiency of turbine airfoils[C]//AIAA 10th Aerospace Sciences Meeting,San Diego:AIAA Paper,1972:17-19.

        (編輯:趙明菁)

        A One-Dimensional Simulation Method to Evaluate Turbine Performance Based on Body Forced Model

        HE Jie1,2,ZOU Zheng-ping1,2,YAO Li-chao1,2
        (1.National Key Laboratory of Science and Technology on Aero-engines,Beihang University;2.Collaborative Innovation Center of Advanced Aero-engine:Beijing 100191,China)

        The precision of turbine aerodynamic performance of low dimensional evaluation method is of great significance for the turbine aerodynamic design.A one-dimensional simulation method for axial flow turbine based on body force model was studied,the blade force source term and centrifugal force source term in Adam and Leonard body force model were corrected,the turbine loss models were introduced in the computation.The method was verified by a low pressure turbine and an air cooling high pressure turbine,and the results show that efficiency deviation of two examples turbine at design point is within 1%.In addition,the results can reflect the tendency for aerodynamic performance of the turbine with the expansion ratio variation.All in all,the corrected body force model has a good precision,and can be used as a fast and reliable evaluation method of low dimensional gas turbine performance.

        turbine;low dimensional design space;body force model;Euler equation;performance evalution

        V 211.6

        A

        10.13477/j.cnki.aeroengine.2016.02.005

        2015-08-05

        何杰(1993),男,在讀碩士研究生,研究方向?yàn)槿~輪機(jī)械數(shù)值仿真;E-mail:hejie@buaa.edu.cn。

        引用格式:何杰,鄒正平,姚李超.基于體積力模型的1維渦輪特性評估方法[J].航空發(fā)動機(jī),2016,42(2):22-26.HE Jie,ZOUZhengping,YAOLichao.A one-dimensionalsimulationmethodtoevaluateturbineperformancebasedonbodyforcedmodel[J].Aeroengine,2016,42(2):22-26.

        猜你喜歡
        冷氣渦輪氣動
        中寰氣動執(zhí)行機(jī)構(gòu)
        基于NACA0030的波紋狀翼型氣動特性探索
        一種新型高分子塑料瓶成型模具
        2014款寶馬525Li渦輪增壓壓力過低
        基于反饋線性化的RLV氣動控制一體化設(shè)計
        冬天的冷氣
        渦輪增壓發(fā)動機(jī)與雙離合變速器的使用
        KJH101-127型氣動司控道岔的改造
        通用飛機(jī)冷氣加油裝置的研制
        Opel公司新型1.0L渦輪增壓直接噴射汽油機(jī)
        日韩制服国产精品一区| 久久精品国产黄片一区| 亚洲不卡高清av网站| 麻豆tv入口在线看| 四虎影视一区二区精品| 国产成人精品人人做人人爽| 国产精品日韩亚洲一区二区| 国产免费又爽又色又粗视频| 少妇被躁爽到高潮无码文| 亚洲综合伦理| 久久黄色精品内射胖女人| 在线中文字幕乱码英文字幕正常| 97精品人妻一区二区三区香蕉| 国产在线拍偷自拍偷精品| 国产精品后入内射日本在线观看| 婷婷色综合视频在线观看| 精品久久人人妻人人做精品| 国产自在自线午夜精品视频在| 国产亚洲精品一区二区在线观看 | 久久精品国产免费观看三人同眠 | 日本a级大片免费观看| 亚洲精品国产av成拍色拍| 人人妻人人做人人爽| 少妇被爽到高潮动态图| 爆乳无码AV国内| 久久日韩精品一区二区| 国产成人一区二区三区| 一级毛片不卡在线播放免费| 成人性生交大片免费看i| 无码精品一区二区三区在线| 欧美性猛交内射兽交老熟妇| 人妻无码∧V一区二区| 成人在线观看av毛片| 国产亚洲日本精品无码| 国产精品自产拍在线观看免费| 白浆高潮国产免费一区二区三区| 小说区激情另类春色| 人人做人人妻人人精| 日本一区二区三区的免费视频观看 | 国产极品喷水视频| 国产精品国三级国产a|