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

        ?

        基于流固耦合的雙向流道軸流泵裝置振動(dòng)特性研究

        2016-03-22 04:57:15李彥軍袁壽其江蘇大學(xué)國家水泵及系統(tǒng)工程技術(shù)研究中心江蘇鎮(zhèn)江212013
        中國農(nóng)村水利水電 2016年12期
        關(guān)鍵詞:軸流泵雙向葉輪

        孟 凡,裴 吉,李彥軍,袁壽其,陳 佳(江蘇大學(xué)國家水泵及系統(tǒng)工程技術(shù)研究中心,江蘇 鎮(zhèn)江 212013)

        0 引 言

        雙向流道軸流泵裝置由于具有占地面積小,運(yùn)行管理方便等特點(diǎn),廣泛應(yīng)用于我國排澇灌溉領(lǐng)域中。近年來,雙向軸流泵裝置的安全穩(wěn)定運(yùn)行越來越受到重視,由于雙向流道軸流泵裝置的可靠性很大程度上依賴于規(guī)范的結(jié)構(gòu)設(shè)計(jì),而結(jié)構(gòu)設(shè)計(jì)的基礎(chǔ)是正確掌握泵各部件應(yīng)力分布及變形的情況。所以在研究泵裝置時(shí)需要將流場計(jì)算與結(jié)構(gòu)計(jì)算綜合考慮。

        流固耦合作用是自然界中客觀存在的一種特殊現(xiàn)象,是指流體與固體之間的相互作用。在水力機(jī)械領(lǐng)域中,流固耦合理論首先應(yīng)用于水輪機(jī)領(lǐng)域[1-4]。隨著水泵向大型化發(fā)展,流固耦合理論也逐漸應(yīng)用于泵領(lǐng)域。袁壽其等[5]及其他學(xué)者[6-9]分析了葉輪流固耦合作用對(duì)離心泵內(nèi)部流場的影響。裴吉等[10]及其他學(xué)者[11-13]利用流固耦合理論對(duì)離心泵葉輪進(jìn)行應(yīng)力及應(yīng)變分析,施衛(wèi)東等[14]及其他學(xué)者[15,16]通過流固耦合模型對(duì)軸流泵葉輪的應(yīng)力應(yīng)變進(jìn)行了數(shù)值研究。吳廣寬等[17]基于瞬態(tài)流固耦合技術(shù)分析了混流式轉(zhuǎn)輪葉片裂紋產(chǎn)生的原因。張德勝[18]采用雙向耦合方法對(duì)蝸殼結(jié)構(gòu)產(chǎn)生的振動(dòng)位移和振動(dòng)速度進(jìn)行了數(shù)值模擬。但是,目前為止,在軸流泵領(lǐng)域中流固耦合方法主要應(yīng)用于泵段的數(shù)值模擬,并沒有考慮進(jìn)出水流道的影響。而雙向流道軸流泵裝置中,雙向進(jìn)水流道由于其形狀特征容易產(chǎn)生渦帶,進(jìn)水流態(tài)不良將直接導(dǎo)致葉輪室內(nèi)的流態(tài)不良,發(fā)生脫流、漩渦等現(xiàn)象,使得葉輪產(chǎn)生強(qiáng)烈振動(dòng)[19],嚴(yán)重時(shí)會(huì)導(dǎo)致葉輪葉片產(chǎn)生裂紋,影響泵裝置的安全運(yùn)行。預(yù)防葉輪葉片出現(xiàn)裂紋的一種措施是對(duì)葉片不同區(qū)域進(jìn)行不同程度的加厚,而葉片加厚是以葉片的應(yīng)力與變形分布為理論基礎(chǔ)。

        為了獲得葉輪應(yīng)力及應(yīng)變的分布情況,本文首次采用雙向同步求解的方法,在考慮雙向進(jìn)出水流道的情況下,對(duì)軸流泵內(nèi)流場和葉輪結(jié)構(gòu)響應(yīng)進(jìn)行聯(lián)合求解,通過數(shù)值模擬可以獲得在不同工況下葉輪葉片的應(yīng)力及變形規(guī)律,不僅可以為葉輪葉片加厚提供理論指導(dǎo),也可以為雙向流道泵裝置的優(yōu)化設(shè)計(jì)與安全運(yùn)行提供參考。

        1 數(shù)值模擬

        1.1 計(jì)算模型

        本文以國內(nèi)某大型雙向流道泵站為研究對(duì)象,該泵站在運(yùn)行過程中會(huì)產(chǎn)生一定程度的水力振動(dòng),從而產(chǎn)生一定程度的機(jī)械振動(dòng),使葉輪根部產(chǎn)生不同程度的裂紋。以下將以該泵站的模型泵裝置為研究對(duì)象進(jìn)行計(jì)算。

        雙向流道軸流泵裝置主要包括雙向進(jìn)出水流道,軸流泵葉輪,軸流泵導(dǎo)葉。其中,葉輪葉片數(shù)為3,導(dǎo)葉葉片數(shù)為5,其結(jié)構(gòu)如圖1所示。模型泵裝置設(shè)計(jì)參數(shù)為:流量Q=0.33 m3/s,轉(zhuǎn)速n=1 550 r/min,揚(yáng)程H=2.76 m。

        圖1 雙向軸流泵裝置三維圖Fig.1 3D model of pump station with two-way passages

        1.2 網(wǎng)格劃分

        利用ICEM對(duì)流體域進(jìn)行網(wǎng)格劃分,其中葉輪水體與導(dǎo)葉水體采用結(jié)構(gòu)網(wǎng)格,進(jìn)出流道采用非結(jié)構(gòu)網(wǎng)格。利用ANSYS對(duì)固體域進(jìn)行網(wǎng)格劃分,固體域只考慮葉輪結(jié)構(gòu)模型。流體域網(wǎng)格與固體域網(wǎng)格如圖2所示。水體域網(wǎng)格總數(shù)為4 030 944,固體域網(wǎng)格數(shù)為18 540。

        圖2 流體域及固體域網(wǎng)格Fig.2 Mesh of fluid and structure

        1.3 邊界設(shè)置

        流體域計(jì)算中,進(jìn)水流道進(jìn)口采用質(zhì)量流量,前池的表面設(shè)置為自由水面,自由水面對(duì)速度和湍動(dòng)能均采用對(duì)稱平面處理。出水流道出口采用固定總壓,總壓設(shè)定為101 325 Pa,采用自由出流邊界條件。所有壁面為光滑壁面,采用無滑移邊界條件。將葉輪與進(jìn)水流道,葉輪與導(dǎo)葉的交界面設(shè)置為瞬態(tài)凍結(jié)轉(zhuǎn)子。將導(dǎo)葉與出水流道的交界面設(shè)置為普通連接。將葉輪與流體的接觸面設(shè)置為動(dòng)網(wǎng)格。葉輪每旋轉(zhuǎn)3°作為一個(gè)時(shí)間步,時(shí)間步長為0.000 32 s。

        葉輪材料選用結(jié)構(gòu)鋼,葉輪材料特性參數(shù)如表1所示。定義葉輪輪轂圓柱面為固定約束,設(shè)置葉輪葉片為流固耦合作用面。固體域的時(shí)間步長與流體域的時(shí)間步長保持一致。

        表1 葉輪材料特性參數(shù)Tab.1 Parameters of impeller material

        雙向耦合計(jì)算的順序是流體域計(jì)算與固體域計(jì)算同步進(jìn)行,并通過耦合作用面進(jìn)行數(shù)據(jù)傳遞,當(dāng)固體域計(jì)算與流體域計(jì)算同時(shí)收斂時(shí),進(jìn)入下一步的耦合計(jì)算。具體計(jì)算步驟如圖3所示。

        圖3 雙向同步耦合計(jì)算過程Fig.3 Solving process of FSI two-way coupling

        為了保證雙向流固耦合的精度,雙向耦合的初始值為流體域非定常計(jì)算10個(gè)周期后的數(shù)據(jù)。

        2 結(jié)果與分析

        2.1 外特性結(jié)果分析

        由于低揚(yáng)程泵裝置實(shí)驗(yàn)測試誤差較大,本文中實(shí)驗(yàn)值為原型泵裝置運(yùn)行時(shí)的單機(jī)揚(yáng)程值,流固耦合計(jì)算值是指模型泵進(jìn)行雙向流固耦合計(jì)算后的第二個(gè)周期的時(shí)均揚(yáng)程。數(shù)值模擬中采用實(shí)際尺寸縮小10倍的模型泵,轉(zhuǎn)速利用等揚(yáng)程幾何相似換算定律所得:

        (1)

        式中:H,HM為實(shí)際泵揚(yáng)程,模型泵揚(yáng)程;D2,D2M為實(shí)際泵葉輪出口直徑,模型泵葉輪出口直徑;n,nM為實(shí)際泵轉(zhuǎn)速,模型泵轉(zhuǎn)速;ηh,ηhM為實(shí)際泵水力效率,模型泵水力效率。

        如圖4所示,流固耦合計(jì)算值在小流量工況點(diǎn)略低于實(shí)驗(yàn)值,在其他工況點(diǎn)皆高于實(shí)驗(yàn)值。兩條曲線變化趨勢基本相同,在最高效率點(diǎn)(Q=0.396 m3/s),流固耦合計(jì)算值與實(shí)驗(yàn)值相對(duì)誤差為6.2%,說明雙向同步耦合計(jì)算結(jié)果是可靠的,也驗(yàn)證了所選湍流模型與網(wǎng)格類型的適用性。

        圖4 泵裝置揚(yáng)程對(duì)比圖Fig.4 Comparison of head curves

        2.2 葉輪等效應(yīng)力及變形分布圖

        圖5(a)和圖5(b)分別為葉輪吸力面與葉輪壓力面的變形分布。如圖5所示,隨著流量增加葉輪葉片的變形量逐漸減小,且葉輪吸力面的變形量大于葉輪壓力面的變形量。葉輪壓力面與吸力面的最大變形量都分布在靠近輪緣的進(jìn)水邊處,且葉輪壓力面與吸力面的變形量由輪緣向輪轂方向逐漸減小。分析其原因主要是葉輪葉片厚度由輪緣到輪轂方向逐漸增大,葉片輪緣厚度較小,容易受應(yīng)力影響產(chǎn)生變形。

        圖6(a)和圖6(b)分別為葉輪吸力面與葉輪壓力面的等效應(yīng)力分布。如圖所示,葉片吸力面的等效應(yīng)力大于葉片壓力面,且隨著流量增大,葉片的等效應(yīng)力逐漸減小。葉片等效應(yīng)力由葉緣向輪轂方向逐漸增大,主要原因是由于葉片形狀,越接近輪轂葉輪吸力面越容易出現(xiàn)脫流現(xiàn)象導(dǎo)致越接近輪轂葉片壓力面與吸力面的壓差越大。

        圖5 葉輪變形分布圖(單位:mm)Fig.5 Total deformation of impeller

        圖6 葉輪應(yīng)力分布圖(單位:MPa)Fig.6 Equivalent stress of impeller

        2.3 等效應(yīng)力時(shí)域圖分析

        表2為最大等效應(yīng)力與最大變形量值隨流量的變化。如表2所示,在Q=0.8Qopt處,出現(xiàn)最大值分別為137.6 MPa與0.086 771 mm。隨著流量增加,最大等效應(yīng)力與最大變形量逐漸減小。

        由于在不同工況下的最大等效應(yīng)力值相差較大,很難在一張圖內(nèi)表現(xiàn)出3個(gè)工況下的最大等效應(yīng)力的變化趨勢,故采用最大等效應(yīng)力方差的形式來表現(xiàn)變化趨勢:

        D(X)=[X-E(X)]2

        (2)

        式中:X為常數(shù);E(X)為X的平均值;D(X)為X的方差。

        圖7為不同工況下,最大等效應(yīng)力的方差隨時(shí)間的變化。如圖所示,在Q=0.8Qopt時(shí),最大等效應(yīng)力在一個(gè)葉輪周期內(nèi)出現(xiàn)了2個(gè)極大值。在Q=1.0Qopt時(shí),最大等效應(yīng)力在一個(gè)葉輪旋轉(zhuǎn)周期內(nèi)只出現(xiàn)了1個(gè)極小值。在Q=1.1Qopt時(shí),最大等效應(yīng)力在一個(gè)葉輪旋轉(zhuǎn)周期內(nèi)并沒有出現(xiàn)極值??梢钥闯?隨著流量的增加,最大等效應(yīng)力的變化頻率逐漸減弱,3個(gè)工況下最大等效應(yīng)力的相位互不相同。

        表2 最大等效應(yīng)力及最大總變形量比較Tab.2 Comparision of maximum equivalent stress and maximum total deformation

        圖7 最大等效應(yīng)力方差時(shí)域圖Fig.7 Maximum Equivalent stress of impeller

        由于在3種工況下,葉輪的最大等效應(yīng)力都分布在輪轂附近,而輪轂與葉輪的交界面處容易產(chǎn)生應(yīng)力集中。所以對(duì)葉輪前緣根部特征點(diǎn)P1,P2,以及葉輪尾緣根部特征點(diǎn)P3,P4進(jìn)行分析研究,其中P1,P3位于葉片吸力面上,P2,P4位于葉片壓力面上。特征點(diǎn)分布圖如圖8所示。

        圖8 特征點(diǎn)分布圖Fig.8 Distribution of feature points

        圖9-圖12為P1,P2,P3,P4處最大等效應(yīng)力方差隨時(shí)間的變化。如圖所示,P1,P2,P3,P4處的等效應(yīng)力脈動(dòng)變化趨勢基本相同。在Q=0.33 m3/s時(shí),4個(gè)特征點(diǎn)處均有比較明顯的等效應(yīng)力脈動(dòng),但隨著流量增加,4個(gè)特征點(diǎn)處的等效應(yīng)力脈動(dòng)都逐漸減小。4個(gè)特征點(diǎn)中,P1,P2處的等效應(yīng)力脈動(dòng)幅值遠(yuǎn)遠(yuǎn)大于P3,P4處的等效應(yīng)力脈動(dòng)幅值。P2處的等效應(yīng)力脈動(dòng)范圍最大,為0~3 MPa,P4處的壓力脈動(dòng)范圍最小僅為0~0.015 MPa。

        圖9 P1點(diǎn)處的等效應(yīng)力方差時(shí)域圖Fig.9 Time domain of equivalent stress on P1

        圖10 P2點(diǎn)處的等效應(yīng)力方差時(shí)域圖Fig.10 Time domain of equivalent stress on P2

        圖11 P3點(diǎn)處的等效應(yīng)力方差時(shí)域圖Fig.11 Time domain of equivalent stress on P3

        圖12 P4點(diǎn)處的等效應(yīng)力方差時(shí)域圖Fig.12 Time domain of equivalent stress on P4

        圖13 不同工況下4個(gè)特征點(diǎn)處的最大等效應(yīng)力Fig.13 The maximum equivalent stress on 4 feature point under different conditions

        圖13為P1,P2,P3,P4處的最大等效應(yīng)力隨流量的變化。如圖所示,隨著流量增大,除了P4處的最大等效應(yīng)力出現(xiàn)了先降低后上升的趨勢,其余特征點(diǎn)處的最大等效應(yīng)力皆逐漸下降。4個(gè)特征點(diǎn)中,P1,P2處的最大等效應(yīng)力值遠(yuǎn)大于P3,P4處的最大等效應(yīng)力值。在3個(gè)工況下,P2處的最大等效應(yīng)力均大于其余特征點(diǎn)處,P4處的最大等效應(yīng)力均小于其他特征點(diǎn)處。

        通過以上比較可以看出,由于在葉輪旋轉(zhuǎn)時(shí),入流造成葉片前緣處的內(nèi)部流動(dòng)不穩(wěn)定且葉片的最大等效應(yīng)力區(qū)分布在葉輪輪轂處。所以P1,P2處的等效應(yīng)力脈動(dòng)和等效應(yīng)力值遠(yuǎn)大于P3,P4處。以上計(jì)算分析結(jié)果對(duì)泵站在長期運(yùn)行后葉片根部產(chǎn)生一定程度裂紋起到了很好的驗(yàn)證作用。

        3 結(jié) 語

        (1)雙向軸流泵裝置中,葉輪所受軸向力隨葉輪旋轉(zhuǎn)呈周期性變化,其幅值受流量變化影響明顯。

        (2)雙向軸流泵裝置中葉輪葉片厚度對(duì)葉片的變形量起主導(dǎo)作用,最大變形區(qū)集中分布在葉輪前緣與輪緣的夾角處,且葉片變形量由輪緣到輪轂方向逐漸降低。葉輪形狀特征導(dǎo)致葉片所受等效應(yīng)力由輪緣到輪轂方向逐漸增大,最大應(yīng)力區(qū)出現(xiàn)在輪轂附近。

        (3)本文中,最大等效應(yīng)力脈動(dòng)隨流量逐漸降低,在泵裝置最高效率點(diǎn)Q=1.0Qopt時(shí),并不是葉輪的最高效率點(diǎn)。進(jìn)一步證明了對(duì)泵裝置進(jìn)行模擬計(jì)算時(shí)需要考慮進(jìn)出水流道對(duì)泵裝置的影響,才能準(zhǔn)確預(yù)測泵裝置的性能。

        (4)在雙向軸流泵裝置中,葉輪進(jìn)水邊所受應(yīng)力遠(yuǎn)遠(yuǎn)大于葉輪出水邊所受應(yīng)力。原因可能是雙向進(jìn)水流道形狀導(dǎo)致葉輪進(jìn)口處流態(tài)紊亂,入流造成葉輪前緣處的內(nèi)部流道不穩(wěn)定。

        綜上所述可以看出葉輪根部和葉輪進(jìn)水邊所受應(yīng)力較大,需對(duì)該區(qū)域進(jìn)行特別加厚。葉輪所受應(yīng)力隨流量增加有明顯下降,在大流量工況下,為了提高葉輪做功效率可以適當(dāng)降低葉輪葉片厚度。本文所得結(jié)論可以為雙向軸流泵裝置的結(jié)構(gòu)優(yōu)化與性能預(yù)測提供有效參考。

        [1] 張 亮, 何環(huán)宇, 張學(xué)偉, 等. 垂直軸水輪機(jī)單向流固耦合數(shù)值研究[J]. 華中科技大學(xué)學(xué)報(bào): 自然科學(xué)版, 2014,42(5):80-84.

        [2] 方 兵, 金連根, 張仁貢, 等. 基于Ansys-CFX的混流式水輪機(jī)轉(zhuǎn)輪雙向流固耦合數(shù)值模擬方法[J]. 水電站機(jī)電技術(shù), 2015,38(6):1-3.

        [3] 劉德民, 劉小兵, 李 娟. 基于流固耦合的水輪機(jī)振動(dòng)分析[J]. 流體傳動(dòng)與控制, 2008,(1):21-25.

        [4] 肖若富, 朱文若, 楊 魏, 等. 基于雙向流固耦合水輪機(jī)轉(zhuǎn)輪應(yīng)力特性分析[J]. 排灌機(jī)械工程學(xué)報(bào), 2013,31(10):862-866.

        [5] 袁壽其, 徐宇平, 張金鳳, 等. 流固耦合作用對(duì)螺旋離心泵流場影響的數(shù)值分析[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2013,44(1):38-42.

        [6] 江 偉, 郭 濤, 李國君, 等. 離心泵流場流固耦合數(shù)值模擬 [J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2012,43(9):53-56.

        [7] 吳賢芳, 談明高, 劉厚林, 等. 流固耦合作用對(duì)離心泵關(guān)死點(diǎn)內(nèi)流的影響[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào), 2015,23(1):172-181.

        [8] 劉厚林, 徐 歡, 吳賢芳, 等. 流固耦合作用對(duì)離心泵內(nèi)外特性的影響[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2012,28(13):82-87.

        [9] 裴 吉, 袁壽其, 袁建平. 流固耦合作用對(duì)離心泵內(nèi)部流場影響的數(shù)值計(jì)算[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2009,(12):107-112.

        [10] 裴 吉, 袁壽其, 袁建平. 基于雙向耦合的污水離心泵動(dòng)應(yīng)力分析[J]. 機(jī)械工程學(xué)報(bào), 2014,(5):45.

        [11] 黃浩欽, 劉厚林, 王 勇, 等. 基于流固耦合的船用離心泵轉(zhuǎn)子應(yīng)力應(yīng)變及模態(tài)研究[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2014,30(15):98-105.

        [12] Benra F K, Dohmen H J. Comparison of pump impeller orbit curves obtained by measurement and FSI simulation[C]∥ASME 2007 Pressure Vessels and Piping Conference. American Society of Mechanical Engineers, 2007:41-48.

        [13] Schneider A, Will B C, B?hle M. Numerical Evaluation of Deformation and Stress in Impellers of Multistage Pumps by Means of Fluid Structure Interaction[C]∥ ASME 2013 Fluids Engineering Division Summer Meeting. American Society of Mechanical Engineers, 2013.

        [14] 施衛(wèi)東, 徐 燕, 張啟華. 基于流固耦合的多級(jí)潛水泵葉輪結(jié)構(gòu)強(qiáng)度分析[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2013,44(5):70-73.

        [15] 朱 利, 楊昌明, 鄭 軍, 等. 基于流固耦合的軸流泵葉輪結(jié)構(gòu)分析[J]. 流體機(jī)械, 2013,41(3):20-23.

        [16] 唐學(xué)林, 王秀葉, 賈玉霞. 基于流固耦合的燈泡貫流泵葉輪強(qiáng)度分析[J]. 排灌機(jī)械工程學(xué)報(bào), 2014,32(11):921-926.

        [17] 吳廣寬, 羅興锜, 馮建軍, 等. 基于瞬態(tài)流固耦合的混流式轉(zhuǎn)輪葉片裂紋成因分析[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2015,(8).

        [18] 張德勝, 張 磊, 施衛(wèi)東, 等. 基于流固耦合的離心泵蝸殼振動(dòng)特性優(yōu)化[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2013,44(91):40-45.

        [19] 劉為民. 泵站進(jìn)水流道對(duì)水泵性能影響的數(shù)值模擬研究[D]. 江蘇揚(yáng)州:揚(yáng)州大學(xué), 2005.

        猜你喜歡
        軸流泵雙向葉輪
        雙向度的成長與自我實(shí)現(xiàn)
        出版人(2022年11期)2022-11-15 04:30:18
        潛水軸流泵運(yùn)行故障分析與排除研究
        潛水軸流泵電機(jī)運(yùn)行工況的特點(diǎn)及可靠性探討
        1.4317 QT2鋼在高能泵葉輪上的應(yīng)用
        應(yīng)用石膏型快速精密鑄造技術(shù)制造葉輪
        離心泵葉輪切割方法
        一種軟開關(guān)的交錯(cuò)并聯(lián)Buck/Boost雙向DC/DC變換器
        濃縮軸流泵干氣密封改造
        一種工作頻率可變的雙向DC-DC變換器
        基于CFD/CSD耦合的葉輪機(jī)葉片失速顫振計(jì)算
        日韩av毛片在线观看| 亚洲另类自拍丝袜第五页| 欧美亚洲国产人妖系列视| 亚洲国产剧情一区在线观看| 中文字幕一区二区在线看| 男女啪啪视频高清视频| 中国午夜伦理片| 午夜不卡av免费| 91啦视频在线观看| 无码少妇a片一区二区三区| 日本丰满妇人成熟免费中文字幕| 亚洲国产综合精品久久av| 男人天堂亚洲一区二区| 男人国产av天堂www麻豆| 国产探花在线精品一区二区| 亚洲最大天堂无码精品区| 福利一区二区三区视频在线| 国产麻豆剧传媒精品国产av| 久久国产成人精品av| 亚洲国产精品福利片在线观看| 国产男女插插一级| 亚洲第一页在线观看视频网站| 国产精品亚洲三级一区二区三区| 久久久久免费看成人影片| 少妇人妻偷人精品一区二区| 亚洲红杏AV无码专区首页 | 人妻系列无码专区久久五月天| 日韩精品视频中文字幕播放| 日本在线精品一区二区三区| 性欧美丰满熟妇xxxx性久久久| 亚洲成人色区| 色综合久久久久综合999| 国产三级在线观看不卡| 青青草成人在线免费视频| 久久精品国产99国产精品亚洲| 成人欧美一区二区三区的电影| 2021国产最新无码视频| 亚洲一区二区三区av天堂| 无色码中文字幕一本久道久| 久久精品国产亚洲av久| 天躁夜夜躁狼狠躁|