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

        ?

        環(huán)形多孔節(jié)流空氣靜壓軸承氣膜流場(chǎng)分析*

        2022-01-17 08:20:58孫小亮盧志偉張君安劉錫堯
        潤(rùn)滑與密封 2021年11期
        關(guān)鍵詞:孔數(shù)氣膜節(jié)流

        孫小亮 盧志偉 張君安 劉 波 劉錫堯

        (西安工業(yè)大學(xué)機(jī)電工程學(xué)院 陜西西安 710021)

        空氣靜壓軸承具有無污染、功耗低和壽命長(zhǎng)等優(yōu)點(diǎn),是超精密加工裝備及航空航天測(cè)試裝備等領(lǐng)域的重要支撐部件[1-2]。由于氣體的非穩(wěn)態(tài)流動(dòng)會(huì)使得氣體在氣膜內(nèi)形成壓力波動(dòng),從而導(dǎo)致軸承產(chǎn)生振幅從幾納米至幾十納米及頻率由幾十赫茲到上千赫茲的寬頻微振動(dòng),其影響不容忽視。因此,提高氣浮支撐系統(tǒng)的工作穩(wěn)定性和減小空氣靜壓軸承產(chǎn)生的微振動(dòng),已成為空氣靜壓氣體軸承研究中普遍關(guān)注的一個(gè)熱點(diǎn)問題[3-4]。

        TALUKDER和STOWELL[5]通過物理實(shí)驗(yàn)研究了空氣靜壓軸承的氣錘振動(dòng)現(xiàn)象,為確保空氣靜壓軸承內(nèi)部湍流流場(chǎng)的穩(wěn)定工作提供了一定的理論指導(dǎo)??字锌频萚6]認(rèn)為空氣靜壓軸承氣膜間隙內(nèi)部存在的壓力波動(dòng)會(huì)引起系統(tǒng)共振,導(dǎo)致氣錘振動(dòng)現(xiàn)象的發(fā)生。CHEN等[7]采用大渦模擬方法對(duì)空氣靜壓軸承瞬態(tài)特征進(jìn)行了研究,認(rèn)為旋渦流的脫落與流場(chǎng)中氣膜壓力波動(dòng)之間存在著一定的相互關(guān)系。王偉等人[8]通過仿真研究發(fā)現(xiàn),當(dāng)空氣靜壓軸承內(nèi)部相對(duì)氣容量超過一定比例時(shí),會(huì)出現(xiàn)微振動(dòng)現(xiàn)象。LI等[9]采用大渦模擬方法分析了單孔節(jié)流空氣靜壓軸承內(nèi)的時(shí)變流場(chǎng),認(rèn)為承載面上的壓力波動(dòng)決定了微振動(dòng)的幅值大小和頻率高低;并通過數(shù)值分析和實(shí)驗(yàn)研究了軸承參數(shù)對(duì)微振動(dòng)大小的影響,指出可根據(jù)承載面氣壓波動(dòng)預(yù)測(cè)軸承的微振動(dòng)。趙曉龍等[10]運(yùn)用軟件通過加入湍流模型對(duì)單個(gè)小孔節(jié)流空氣靜壓軸承進(jìn)行了完整的Navier-Stokes方程求解,分析了節(jié)流孔出口處流場(chǎng)的特性。

        目前對(duì)于空氣靜壓軸承研究主要集中在承載力、剛度等方面且大部分都是針對(duì)單孔式的,而對(duì)于多孔式氣膜流場(chǎng)特性分析相對(duì)比較少[11-14]。因此,本文作者設(shè)計(jì)了一種環(huán)形多孔節(jié)流空氣靜壓軸承,并提出采用大渦模擬方法對(duì)其節(jié)流孔出口處附近計(jì)算區(qū)域的氣膜流場(chǎng)性能展開分析。

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

        1.1 物理模型的建立

        圖1所示為設(shè)計(jì)的一種環(huán)形多孔節(jié)流空氣靜壓軸承物理模型,L表示軸承邊長(zhǎng),H表示軸承高度,h表示氣膜厚度,d表示節(jié)流孔的直徑,b表示節(jié)流孔的孔距,l表示節(jié)流孔的長(zhǎng)度,d1表示多孔節(jié)流器直徑,D表示進(jìn)行大渦模擬數(shù)值計(jì)算區(qū)域。

        圖1 環(huán)形多孔節(jié)流空氣靜壓軸承結(jié)構(gòu)示意Fig 1 Structure of the annular multi-holerestrictor aerostatic bearing

        文中以幾何形狀為方形的環(huán)形多孔節(jié)流空氣靜壓軸承為研究對(duì)象,其具體的相關(guān)尺寸參數(shù)如表1所示。

        表1 環(huán)形多孔節(jié)流空氣靜壓軸承設(shè)計(jì)參數(shù)

        1.2 控制方程

        大渦模擬的基本思想是通過濾波將湍流中的瞬時(shí)脈動(dòng)分解為大尺度脈動(dòng)和小尺度脈動(dòng),大尺度脈動(dòng)通過可解尺度的Navier-Stokes方程直接求得,而將小尺度脈動(dòng)的動(dòng)量和能量輸運(yùn)對(duì)大尺度運(yùn)動(dòng)的影響用亞格子模型予以反映,小尺度脈動(dòng)用建立亞格子應(yīng)力模型進(jìn)行模擬。采用密度加權(quán)過濾器(Favre濾波)對(duì)可壓縮Navier-Stokes方程進(jìn)行濾波,式(1)—(3)分別為濾波后得到的連續(xù)方程、動(dòng)量方程和能量方程[15]。

        (1)

        (2)

        (3)

        2 網(wǎng)格劃分及計(jì)算方法

        文中為了提高計(jì)算精度需要對(duì)節(jié)流孔壁面和氣膜間隙壁面處進(jìn)行網(wǎng)格局部加密,如圖2所示。圖2中A處網(wǎng)格局部放大圖,如圖3所示。

        圖2 空氣靜壓軸承計(jì)算區(qū)域網(wǎng)格Fig 2 Grid of calculation area for aerostatic bearing

        圖3 A處網(wǎng)格局部放大圖Fig 3 A partial enlarged view of the grid at A

        對(duì)空氣靜壓軸承的氣膜流場(chǎng)求解區(qū)域邊界條件設(shè)為:進(jìn)口邊界條件等于供氣壓力ps;壁面、氣膜邊界條件均采用無滑移的壁面邊界條件;假設(shè)F計(jì)算區(qū)域邊界條件等于出口壓力pl,它是通過直接求解雷諾方程和流量平衡方程來確定。

        進(jìn)行數(shù)值計(jì)算時(shí),首先在穩(wěn)態(tài)下采用SIMPLE算法對(duì)壓力-速度耦合方案進(jìn)行計(jì)算,然后將其計(jì)算結(jié)果作為進(jìn)行大渦模擬LES瞬態(tài)計(jì)算的一個(gè)初始值。在進(jìn)行大渦模擬LES瞬態(tài)計(jì)算中,采用PISO算法對(duì)壓力-速度進(jìn)行耦合,動(dòng)量選擇高階的“有界中心差分”插值格式,其余能量等選取二階迎風(fēng)插值格式;設(shè)定材料為ideal-gas氣體,迭代殘差等于1×10-5時(shí)停止計(jì)算,設(shè)置空氣靜壓軸承外部環(huán)境變量參數(shù)如表2所示。

        表2 外部環(huán)境變量參數(shù)

        2.1 網(wǎng)格無關(guān)性驗(yàn)證

        選取環(huán)形多孔節(jié)流空氣靜壓軸承節(jié)流孔數(shù)n=9、節(jié)流孔的孔距b=1.25 mm、節(jié)流孔直徑d=0.2 mm、氣膜間隙h=10 μm、供氣壓力ps=0.5 MPa進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證。劃分三套網(wǎng)格:(1)粗糙:網(wǎng)格總數(shù)為270萬;(2)中等:網(wǎng)格總數(shù)為540萬;(3)精細(xì):網(wǎng)格總數(shù)為790萬。三套網(wǎng)格對(duì)空氣靜壓軸承節(jié)流孔出口處最大壓降幅度Δp(MPa)的計(jì)算結(jié)果如表3所示??芍?dāng)網(wǎng)格數(shù)量超過540萬時(shí),數(shù)值計(jì)算結(jié)果已經(jīng)不再隨著網(wǎng)格數(shù)量的增加發(fā)生明顯的變化。因此,綜合考慮數(shù)值計(jì)算精度和實(shí)際運(yùn)算設(shè)備以及數(shù)值計(jì)算時(shí)間周期,文中采用中等數(shù)量的網(wǎng)格劃分方法對(duì)物理模型進(jìn)行網(wǎng)格劃分。

        表3 網(wǎng)格無關(guān)性驗(yàn)證

        2.2 研究方法驗(yàn)證

        為驗(yàn)證文中使用的研究方法,擬采用空氣靜壓軸承節(jié)流孔數(shù)n=1、節(jié)流孔直徑d=0.2 mm、氣膜間隙h=10 μm、供氣壓力ps=0.5 MPa等參數(shù),對(duì)空氣靜壓軸承沿著徑向的壓力分布進(jìn)行數(shù)值計(jì)算,并將之與文獻(xiàn)[16]中的實(shí)驗(yàn)和理論計(jì)算結(jié)果進(jìn)行了對(duì)比,如圖4所示??梢钥闯?,由于文中選擇空氣靜壓軸承節(jié)流孔附近區(qū)域進(jìn)行理論計(jì)算而并非整個(gè)模型,且進(jìn)行數(shù)值計(jì)算所設(shè)定的工作參數(shù)有一定區(qū)別,因此與文獻(xiàn)中理論和實(shí)驗(yàn)結(jié)果有一定差異;但文中計(jì)算的空氣靜壓軸承沿著徑向的壓力分布與文獻(xiàn)中的理論十分接近,且與文獻(xiàn)實(shí)驗(yàn)結(jié)果的變化趨勢(shì)有很好的一致性,這在一定程度上驗(yàn)證了文中所采用方法的可行性與數(shù)值計(jì)算的正確性。

        圖4 文獻(xiàn)[16]實(shí)驗(yàn)和理論計(jì)算結(jié)果與文中數(shù)值計(jì)算結(jié)果的比較Fig 4 Comparison between experimental and theoretical results inreference[16] and numerical results in this paper (a)gas film pressure distribution for reference [16];(b)gas film pressure distribution curve of 1 hole bearing

        3 結(jié)果及分析

        文中分別取過空氣靜壓軸承中心節(jié)流孔與環(huán)形節(jié)流孔的平面為P1、P3,僅過中心節(jié)流孔取平面為P2;平面P2與平面P1、P3之間的夾角都為22.5°;過氣膜間隙上壁面附近位置取一個(gè)壓力平面P4,中間位置取一個(gè)速度平面P5;在氣膜間隙上分別取靠近上壁面附近位置G、中間位置H、下壁面附近位置K。選擇環(huán)形多孔節(jié)流空氣靜壓軸承節(jié)流孔數(shù)n=9、節(jié)流孔的孔距b=1.25 mm、節(jié)流孔直徑d=0.2 mm、氣膜間隙h=20 μm、供氣壓力ps=0.5 MPa,對(duì)軸承氣膜流場(chǎng)進(jìn)行了分析。

        3.1 不同截面空氣靜壓軸承氣膜流場(chǎng)的壓力分析

        在平面P1、P2、P3上,得到環(huán)形多孔節(jié)流空氣靜壓軸承計(jì)算區(qū)域氣膜流場(chǎng)的壓力云圖,如圖5所示。

        圖5 不同截面的壓力云圖Fig 5 Pressure clouds at different sections

        從圖5中可知,在平面P1、P2、P3上,空氣靜壓軸承中心節(jié)流孔出口處附近氣膜流場(chǎng)的壓力均表現(xiàn)出先減小后增加且都形成了壓力渦旋,位置靠近于氣膜上壁面,中心節(jié)流孔出口處附近B的局部放大圖,如圖6所示。

        圖6 B處的局部放大圖Fig 6 A partial enlarged view of B

        在氣膜間隙上壁面附近位置壓力平面P4上,得到環(huán)形多孔節(jié)流空氣靜壓軸承氣膜流場(chǎng)的壓力云圖,如圖7所示。

        圖7 上壁面附近氣膜壓力云圖Fig 7 Gas film pressure cloud diagram near the upper wall

        結(jié)合圖5與圖7可知,在環(huán)形節(jié)流孔附近時(shí),平面P2與平面P1、P3上氣膜的壓力變化是不同的。在平面P2上氣膜的壓力呈現(xiàn)出逐漸減小趨勢(shì);但在平面P1、P3上環(huán)形節(jié)流孔出口處附近靠近中心節(jié)流孔一側(cè)氣膜的壓力表現(xiàn)出先減小而后增加的趨勢(shì),在靠近空氣靜壓軸承出口處一側(cè)氣膜的壓力表現(xiàn)出先減小而后增加再減小的趨勢(shì),在節(jié)流孔兩側(cè)均出現(xiàn)壓力突降的現(xiàn)象,但靠近中心節(jié)流孔一側(cè)氣膜的壓力突降幅度小于靠近空氣靜壓軸承出口處一側(cè)氣膜的壓力突降幅度,且靠近中心孔一側(cè)形成了壓力渦旋。平面P1、P3上氣膜的壓力變化表現(xiàn)出一致性,環(huán)形節(jié)流孔出口處附近C的局部放大圖,如圖8所示。

        圖8 C的局部放大圖Fig 8 A partial enlarged view of C

        3.2 不同截面空氣靜壓軸承氣膜流場(chǎng)的速度分析

        在平面P1、P2、P3上,得到環(huán)形多孔節(jié)流空氣靜壓軸承計(jì)算區(qū)域氣膜流場(chǎng)的速度云圖,如圖9所示。

        從圖9中可知,在平面P1、P2、P3上,空氣靜壓軸承中心節(jié)流孔出口處附近氣膜流場(chǎng)的速度均表現(xiàn)出先增加后減小,且形成了氣旋, 中心節(jié)流孔出口處附近D的局部放大圖,如圖10所示。

        圖9 不同截面的速度云圖Fig 9 Velocity clouds at different sections

        圖10 D的局部放大圖Fig 10 A partial enlarged view of D

        在氣膜間隙中間位置的速度平面P5上,得到環(huán)形多孔節(jié)流空氣靜壓軸承氣膜流場(chǎng)的速度云圖,如圖11所示。

        結(jié)合圖9與圖11可知,在環(huán)形節(jié)流孔附近時(shí),平面P2與平面P1、P3上氣膜的速度變化是不同的。在平面P2上氣膜的速度呈現(xiàn)出先增加而后減小的趨勢(shì)并形成了氣旋,相比于平面P1、P3上氣膜的速度增加到最大值的位置向外發(fā)生偏移且其值小于環(huán)形節(jié)流孔出口處氣膜的速度;但在平面P1、P3上環(huán)形節(jié)流孔出口處附近靠近中心節(jié)流孔一側(cè)氣膜的速度表現(xiàn)出先增加而后減小的趨勢(shì),在靠近空氣靜壓軸承出口處一側(cè)氣膜的速度表現(xiàn)出先增加而后減小再增加的趨勢(shì);在節(jié)流孔兩側(cè)都形成了氣旋,且出現(xiàn)了速度突升的現(xiàn)象,但靠近中心節(jié)流孔一側(cè)氣膜的速度突升幅度小于靠近空氣靜壓軸承出口處一側(cè)氣膜的速度突升幅度,環(huán)形節(jié)流孔出口處附近E的局部放大圖,如圖12所示。

        圖11 中間位置氣膜的速度云圖Fig 11 Gas film velocity cloud diagram at the middle position

        圖12 E的局部放大圖Fig 12 A partial enlarged view of E

        在環(huán)形多孔節(jié)流空氣靜壓軸承中心節(jié)流孔出口處附近氣膜流場(chǎng)都形成了壓力渦旋、速度氣旋;環(huán)形節(jié)流孔出口處附近兩側(cè)壓力均呈現(xiàn)出突降現(xiàn)象且靠近中心孔一側(cè)形成了壓力渦旋,氣膜流場(chǎng)的速度呈現(xiàn)出突升現(xiàn)象且都形成了氣旋。

        3.3 不同位置處空氣靜壓軸承氣膜壓力、速度分布

        選擇氣膜間隙上壁面附近位置G、中間位置H、下壁面附近位置K,得到空氣靜壓軸承計(jì)算區(qū)域氣膜壓力和速度分布曲線,如圖13與圖14所示。

        圖14 氣膜速度分布曲線Fig 14 Gas film velocity distribution curves

        從圖13中可知,在節(jié)流孔出口處附近不同位置氣膜壓力均表現(xiàn)出突降,氣膜壓力在由最低點(diǎn)逐漸回升中呈現(xiàn)出波動(dòng)現(xiàn)象;在x=1.35 mm附近,G處氣膜壓力最小值為0.272 1 MPa,H處氣膜壓力最小值為0.286 4 MPa,K處氣膜壓力最小值為0.287 3 MPa,但在上壁面附近位置G處壓力突降最為顯著;空氣靜壓軸承的氣膜壓力在節(jié)流孔出口處呈現(xiàn)出分離現(xiàn)象,然而在遠(yuǎn)離節(jié)流孔出口處不同位置氣膜壓力曲線又表現(xiàn)出重合現(xiàn)象。

        從圖14中可知,在節(jié)流孔出口處不同位置氣膜的速度均表現(xiàn)出突升,氣膜的速度在由最高點(diǎn)逐漸回落中呈現(xiàn)出大幅度波動(dòng)現(xiàn)象;在x=1.35 mm附近,G處氣膜的速度最大值為241.8 m/s,H處氣膜的速度最大值為290.6 m/s,K處氣膜的速度最大值為145.8 m/s,但在中間位置H處速度突升最為明顯;空氣靜壓軸承氣膜的速度在節(jié)流孔出口處也呈現(xiàn)出分離現(xiàn)象,在遠(yuǎn)離節(jié)流孔的出口H處速度明顯高于G和K處速度,G、K處速度曲線表現(xiàn)出重合現(xiàn)象。

        從圖13與圖14中可知,在G處氣膜壓力變化明顯,H處氣膜的速度變化顯著,所以在分析空氣靜壓軸承的氣膜壓力變化時(shí)選擇氣膜上壁面附近位置G,氣膜的速度變化選取中間位置H。

        3.4 孔數(shù)對(duì)空氣靜壓軸承氣膜流場(chǎng)壓力、速度的影響

        選擇環(huán)形多孔節(jié)流空氣靜壓軸承節(jié)流孔數(shù)n=1、n=9、節(jié)流孔的孔距b=1.25 mm、節(jié)流孔直徑d=0.2 mm、氣膜間隙為h=10 μm、供氣壓力ps=0.5 MPa時(shí),得到其氣膜流場(chǎng)的壓力和速度計(jì)算結(jié)果,如圖15與圖16所示。

        圖15 不同孔數(shù)下的氣膜壓力分布曲線Fig 15 Gas film pressure distribution curves underdifferent number of holes

        圖16 不同孔數(shù)下的氣膜的速度分布曲線Fig 16 Gas film velocity distribution curves underdifferent number of holes

        從圖15中可知,在節(jié)流孔出口處氣膜壓力均表現(xiàn)出突降現(xiàn)象,隨著節(jié)流孔數(shù)增加氣膜壓力突降幅度會(huì)顯著減小;當(dāng)節(jié)流孔數(shù)n=1時(shí),在x=0.10 mm處氣膜壓力最小值為0.185 MPa;當(dāng)節(jié)流孔數(shù)n=9時(shí),在x=1.35 mm處氣膜壓力最小值為0.418 MPa;最大壓降幅度9孔為1孔的26%左右。

        從圖16中可知,在節(jié)流孔出口處氣膜的速度均表現(xiàn)出突升現(xiàn)象,隨著節(jié)流孔數(shù)增加其氣膜的速度突升幅度會(huì)顯著減小。當(dāng)節(jié)流孔數(shù)n=1時(shí),x=0.10 mm處氣膜的速度最大值為338 m/s;當(dāng)節(jié)流孔數(shù)n=9時(shí),x=1.35 mm處氣膜的速度最大值為148 m/s;最大速度突升幅度9孔為1孔的43%左右。

        從圖15與圖16中可知,增加節(jié)流孔的數(shù)量n可以顯著減小空氣靜壓軸承節(jié)流孔出口處氣膜壓力的突降幅度和速度的突升幅度。

        4 結(jié)論

        (1)對(duì)環(huán)形多孔節(jié)流空氣靜壓軸承的氣膜間隙上壁面附近位置G、中間位置H、下壁面附近位置K流場(chǎng)進(jìn)行了分析,發(fā)現(xiàn)氣膜壓力在節(jié)流孔出口處附近氣膜間隙上呈現(xiàn)出分離現(xiàn)象,但在遠(yuǎn)離節(jié)流孔出口處在不同位置氣膜壓力曲線又表現(xiàn)出一致;氣膜壓力在氣膜上壁面附近變化明顯,氣膜的速度在氣膜中間位置變化突出。該研究結(jié)果對(duì)分析空氣靜壓軸承的氣膜壓力、速度變化具有一定的參考價(jià)值。

        (2)研究表明,在節(jié)流孔出口處,節(jié)流孔數(shù)為9時(shí)的最大壓降幅度是節(jié)流孔數(shù)為1的26%左右,最大速度突升幅度是節(jié)流孔數(shù)為1的43%左右,這表明增加節(jié)流孔的孔數(shù)可以顯著減小空氣靜壓軸承節(jié)流孔出口處氣膜壓力的突降幅度和速度的突升幅度,提高空氣靜壓軸承工作的穩(wěn)定性。

        (3)在環(huán)形多孔節(jié)流空氣靜壓軸承中心節(jié)流孔出口處附近氣膜流場(chǎng)都形成了壓力渦旋、速度氣旋;環(huán)形節(jié)流孔出口處附近兩側(cè)壓力均表現(xiàn)出明顯的突降現(xiàn)象,且靠近中心孔一側(cè)形成了壓力渦旋,速度呈現(xiàn)出明顯的突升現(xiàn)象且都形成了氣旋;節(jié)流孔出口處附近壓力和速度的劇烈變化與其內(nèi)部產(chǎn)生的微振動(dòng)密切相關(guān)。該結(jié)果對(duì)研究空氣靜壓軸承的微振動(dòng)現(xiàn)象提供了一定的理論指導(dǎo)。

        猜你喜歡
        孔數(shù)氣膜節(jié)流
        T 型槽柱面氣膜密封穩(wěn)態(tài)性能數(shù)值計(jì)算研究
        天然氣井井下節(jié)流器研究現(xiàn)狀及應(yīng)用前景
        穴盤孔數(shù)對(duì)早熟菜薹生長(zhǎng)發(fā)育及產(chǎn)量的影響
        氣膜孔堵塞對(duì)葉片吸力面氣膜冷卻的影響
        超高壓氣井井下節(jié)流技術(shù)應(yīng)用和設(shè)計(jì)方法
        節(jié)流孔數(shù)對(duì)空間熱泵徑向氣體軸承承載特性影響分析
        載人航天(2020年5期)2020-10-31 06:47:38
        噴油嘴孔徑、孔數(shù)對(duì)4100QBZL柴油機(jī)PM排放的影響研究
        西部論叢(2020年7期)2020-10-20 11:41:23
        靜葉柵上游端壁雙射流氣膜冷卻特性實(shí)驗(yàn)
        PR方程模擬節(jié)流效應(yīng)的數(shù)值研究
        山東化工(2019年12期)2019-07-05 08:44:26
        躲避霧霾天氣的氣膜館
        国产精品情侣露脸av在线播放| 在线天堂www中文| 无码人妻丰满熟妇片毛片| 久久这里只精品国产2| 视频一区视频二区亚洲| 国产精品毛片无遮挡高清| 无码人妻精品一区二区在线视频| 亚洲国产福利精品一区二区| 有码中文字幕一区二区| 中文字幕一区二区三区久久网| 国产无套内射久久久国产| 国产欧美精品一区二区三区,| 国产美女高潮流的白浆久久| 日本高清一级二级三级| 欧美肥胖老妇做爰videos| 亚洲动漫成人一区二区| 日本av第一区第二区| 男人添女人囗交做爰视频| 韩国19禁无遮挡啪啪无码网站| 亚洲午夜无码视频在线播放 | 无码一区二区三区中文字幕| 青青青国产精品一区二区| 国产亚洲无码1024| 亚洲一区二区日韩精品在线| 国产白嫩护士被弄高潮| 北条麻妃在线视频观看| 日本不卡一区二区三区在线| 亚洲精品在线国产精品| 射死你天天日| 激情亚洲的在线观看| 久久一区二区三区少妇人妻| 又爽又黄又无遮挡的视频| 欧美日韩精品福利在线观看| 中文字幕人成乱码中文| 少妇久久久久久被弄高潮| 少妇精品久久久一区二区三区| 麻豆av一区二区天堂| 偷拍夫妻视频一区二区| 日韩制服国产精品一区| 国产美女裸身网站免费观看视频| 国产视频激情视频在线观看|