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

        ?

        錯(cuò)列雙圓柱下游圓柱的升力機(jī)理

        2019-05-09 07:28:04杜曉慶張盛華王玉梁孫雅慧
        關(guān)鍵詞:旋渦流態(tài)升力

        杜曉慶, 張盛華, 王玉梁, 孫雅慧

        (1. 上海大學(xué) 土木工程系,上海 200072;2. 上海大學(xué) 風(fēng)工程和氣動(dòng)控制研究中心,上海 200072)

        大長細(xì)比圓柱型結(jié)構(gòu)在實(shí)際工程中有大量的應(yīng)用,如橋梁并列索、超高層建筑、多分裂導(dǎo)線、冷卻塔、煙囪等,并且多以柱群形式出現(xiàn).多個(gè)圓柱結(jié)構(gòu)間存在強(qiáng)烈的氣動(dòng)干擾,下游圓柱在橫風(fēng)向升力作用下常發(fā)生尾流激振現(xiàn)象,導(dǎo)致結(jié)構(gòu)發(fā)生損傷或破壞[1-4],澄清雙圓柱結(jié)構(gòu)下游圓柱的升力機(jī)理有助于理解尾流激振的發(fā)生機(jī)制.

        已有研究發(fā)現(xiàn)[5-7],對于小間距(P/D=1.2~3.5,其中P為圓心間距,D為圓柱直徑)雙圓柱,下游圓柱在小風(fēng)向角(β≈8°)下受到很大的平均升力,稱之為內(nèi)側(cè)升力;對于中等間距雙圓柱(P/D=3.0~5.0),下游圓柱在β=20°左右同樣受到較大的平均升力,稱之為外側(cè)升力.

        對于內(nèi)側(cè)升力,低雷諾數(shù)下的流跡顯示試驗(yàn)表明[6,8]:兩圓柱間的加速間隙流、下游圓柱的停滯點(diǎn)偏移以及上游圓柱的剪切流在下游圓柱表面發(fā)生再附是升力出現(xiàn)的三個(gè)原因.杜曉慶等[9]進(jìn)行了基于高雷諾數(shù)下的大渦模擬研究,對下游圓柱平均升力的出現(xiàn)給出了不同的解釋.

        對于外側(cè)升力,Zdravkovich[7]和Ting等[10]提出,兩圓柱間隙氣流速度的增大使下游圓柱間隙側(cè)表面壓力減小,導(dǎo)致平均升力的出現(xiàn),文獻(xiàn)[10-12]中也認(rèn)為下游圓柱風(fēng)壓停滯點(diǎn)的偏移對平均升力有貢獻(xiàn).Alam等[13-14]提出,受上游圓柱尾流旋渦的作用,下游圓柱的旋渦加速向下游流動(dòng),這也會(huì)引起下游圓柱間隙側(cè)表面壓力的減小,進(jìn)而導(dǎo)致平均升力的出現(xiàn).

        對以往文獻(xiàn)的分析可知,由于影響雙圓柱的因素眾多,干擾機(jī)理復(fù)雜,研究者對下游圓柱的升力機(jī)理尚未有統(tǒng)一的認(rèn)識(shí).值得注意的是,以往的機(jī)理解釋都是基于較低雷諾數(shù)下的流跡顯示試驗(yàn)得出的,但雙圓柱繞流有強(qiáng)烈的雷諾數(shù)效應(yīng)[3,15],因而在高雷諾數(shù)下很可能存在不同的干擾機(jī)理.此外,以往文獻(xiàn)主要集中于對雙圓柱平均升力的研究,受到測試條件的限制,對脈動(dòng)升力及其發(fā)生機(jī)理的研究很少.

        本文中采用大渦模擬(LES)方法,在高雷諾數(shù)(Re=1.4×105)下,研究了不同風(fēng)向角時(shí)中等間距(P/D=4)雙圓柱的氣動(dòng)力特性和流場特性,重點(diǎn)討論了下游圓柱的氣動(dòng)升力與流場結(jié)構(gòu)之間的內(nèi)在關(guān)系,提出了雙圓柱干擾流態(tài)隨風(fēng)向角的變化規(guī)律,進(jìn)一步澄清下游圓柱的平均升力和脈動(dòng)升力的流場機(jī)理.

        1 數(shù)值方法和計(jì)算模型

        1.1 控制方程和亞格子尺度模型

        在大渦模擬方法中,大尺度的旋渦由濾波后的Navier-Stokes方程直接求解,而小尺度的旋渦則采用亞格子尺度(SGS)模型模擬.與雷諾平均法(RANS)相比,大渦模擬方法可更好地模擬流場中的湍流旋渦,能捕捉到更豐富的流場脈動(dòng)信息.因此,為了準(zhǔn)確地模擬錯(cuò)列雙圓柱周圍的流場特性,采用大渦模擬方法進(jìn)行研究.經(jīng)過濾波函數(shù)的濾波,可得到大尺度旋渦的不可壓縮Navier-Stokes方程,如下所示:

        (1)

        (2)

        (3)

        (4)

        式中:μt為亞格子尺度的湍動(dòng)黏度;δij為單位張量.μt的計(jì)算式如下所示:

        (5)

        (6)

        (7)

        式中:Δi為沿坐標(biāo)軸i方向的網(wǎng)格尺寸;CS為Smagorinsky常數(shù),取0.1.

        1.2 計(jì)算模型和計(jì)算參數(shù)

        錯(cuò)列雙圓柱的布置形式如圖1所示.圖1中:CD1、CL1分別為上游圓柱的阻力和升力;CD2、CL2分別為下游圓柱的阻力和升力;θ為圓柱表面風(fēng)壓測點(diǎn)位置的角度,以圓柱迎風(fēng)點(diǎn)為零點(diǎn),順時(shí)針旋轉(zhuǎn)為正方向.兩個(gè)圓柱的直徑均為D,圓柱間距P為4D,風(fēng)向角β為自由來流與圓柱中心連線的夾角.為了掌握不同風(fēng)向角對圓柱周圍流場特性和圓柱氣動(dòng)力特性的影響規(guī)律,計(jì)算的風(fēng)向角β分別為0°、5°、10°、15°、20°、25°、30°、45°、60°、75°、90°.

        圖1 雙圓柱計(jì)算模型

        計(jì)算域網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格、O型域,徑向長度為46D,展向長度為2D,均分為20層網(wǎng)格,計(jì)算模型總網(wǎng)格數(shù)約為266萬,近壁面最小網(wǎng)格厚度為1×10-4D,最大阻塞率為4.3%(β=90°).量綱一時(shí)間步Δt*為0.005(Δt*=ΔtU0/D,其中Δt為實(shí)際計(jì)算時(shí)間步,U0為來流風(fēng)速).計(jì)算域采用均勻速度入口邊界條件,自由出口邊界條件,展向兩端采用周期性邊界條件,圓柱表面采用無滑移壁面邊界條件,如圖2所示.計(jì)算流體動(dòng)力學(xué)(CFD)數(shù)值模擬基于Fluent14.0程序?qū)崿F(xiàn),采用壓力耦合方程的半隱式(SIMPLEC)格式求解壓力速度耦合方程組,空間離散采用中心差分格式,時(shí)間離散采用二階全隱格式.

        a 計(jì)算域和邊界條件

        b 局部的計(jì)算網(wǎng)格

        2 計(jì)算結(jié)果及分析

        2.1 計(jì)算參數(shù)驗(yàn)證

        在確定雙圓柱計(jì)算模型和計(jì)算參數(shù)之前,先以單圓柱為研究對象,對六種計(jì)算模型(分別命名為Case 1~6)進(jìn)行數(shù)值模擬,研究周向網(wǎng)格數(shù)量、時(shí)間步和展向長度等參數(shù)對計(jì)算結(jié)果的影響[16].表1給出了單圓柱各計(jì)算模型的參數(shù)、平均阻力系數(shù)、脈動(dòng)升力系數(shù)和St數(shù),也列出了文獻(xiàn) [17]和[18]中的風(fēng)洞試驗(yàn)結(jié)果.通過比較可以看出,模型Case 6計(jì)算得到的氣動(dòng)力系數(shù)與文獻(xiàn)的吻合情況最好,因此采用該模型參數(shù)進(jìn)行錯(cuò)列雙圓柱的大渦模擬研究.

        2.2 氣動(dòng)力系數(shù)

        圖3是下游圓柱平均氣動(dòng)力系數(shù)曲線.由圖3可見,下游圓柱的平均升力系數(shù)和平均阻力系數(shù)隨風(fēng)向角的變化趨勢與文獻(xiàn)較為吻合.與單圓柱相比,在風(fēng)向角0°~30°內(nèi),下游圓柱的平均阻力系數(shù)均小于單圓柱,并且隨風(fēng)向角的增大而增大.從β=30°開始,下游圓柱的平均阻力系數(shù)趨于平緩,其值比單圓柱稍大.

        表1 單圓柱計(jì)算結(jié)果對比[16]

        a 阻力系數(shù)

        b 升力系數(shù)

        對于平均升力系數(shù)來說,隨著風(fēng)向角的增大,下游圓柱的平均升力系數(shù)先減小后增大,最后趨于穩(wěn)定.在β=0°~45°時(shí),下游圓柱的平均升力系數(shù)均為負(fù)值,并且在β=20°左右出現(xiàn)升力的極小值,即出現(xiàn)圓柱上側(cè)升力峰值.當(dāng)β≥45°時(shí),平均升力系數(shù)便與單圓柱相似并接近于零.圖4是下游圓柱脈動(dòng)氣動(dòng)力系數(shù)隨風(fēng)向角的變化曲線.由圖4可以看出,在β=0°~5°時(shí)下游圓柱的脈動(dòng)升力最大,其值達(dá)到單圓柱的2倍,并且在β=20°~25°時(shí)達(dá)到極小值.當(dāng)β>30°后,下游圓柱的脈動(dòng)升力系數(shù)開始緩慢減小,但依然明顯高于單圓柱;下游圓柱的脈動(dòng)阻力則在β=10°左右達(dá)到最大值,并且在β=0°~30°范圍內(nèi)都明顯大于單圓柱.

        圖4 下游圓柱脈動(dòng)氣動(dòng)力系數(shù)隨風(fēng)向角的變化

        Fig.4 Root mean square aerodynamic coefficients as a function of incidence angles for downstream cylinder

        2.3 下游圓柱表面風(fēng)壓

        圖5為下游圓柱在不同風(fēng)向角下的平均風(fēng)壓系數(shù)曲線.從圖5可以看出:當(dāng)β=0°時(shí),整個(gè)下游圓柱表面都為負(fù)壓,并具有良好的對稱性;當(dāng)β=5°時(shí),下游圓柱表面負(fù)壓遠(yuǎn)小于其他工況,并且停滯點(diǎn)偏移了20°左右;在β=15°~30°時(shí),下游圓柱表面平均風(fēng)壓呈現(xiàn)出明顯的不對稱性,圓柱下側(cè)(間隙側(cè))表面受到的負(fù)壓絕對值明顯大于上側(cè),而平均風(fēng)壓分布的不均勻則導(dǎo)致圓柱表面受到平均升力;當(dāng)風(fēng)向角大于45°時(shí),下游圓柱表面平均風(fēng)壓逐漸接近單圓柱,又重新表現(xiàn)出良好的對稱性.

        圖6為下游圓柱在不同風(fēng)向角下的脈動(dòng)風(fēng)壓系數(shù)曲線.由圖6可知:幾乎所有風(fēng)向角下下游圓柱的脈動(dòng)風(fēng)壓系數(shù)值都不同程度地大于單圓柱,這說明受上游圓柱的影響,下游圓柱表面的脈動(dòng)風(fēng)壓發(fā)生劇烈變化;在β=5°時(shí),下游圓柱表面的脈動(dòng)風(fēng)壓變化最為劇烈,也具有一定的對稱性;在β=15°~30°時(shí),下游圓柱表面的脈動(dòng)風(fēng)壓呈現(xiàn)出明顯的不對稱性,下游圓柱上側(cè)表面的脈動(dòng)風(fēng)壓變化較為平緩,而下側(cè)(間隙側(cè))表面的脈動(dòng)風(fēng)壓變化很大,脈動(dòng)風(fēng)壓的不對稱性也會(huì)促使下游圓柱受到平均升力的作用;在β=45°~90°時(shí),隨著風(fēng)向角的增大,下游圓柱表面的脈動(dòng)風(fēng)壓逐漸接近于單圓柱,并且重新具有良好的對稱性.

        圖5 下游圓柱表面平均風(fēng)壓系數(shù)分布

        圖6 下游圓柱表面脈動(dòng)風(fēng)壓系數(shù)分布

        2.4 氣動(dòng)力功率譜

        對于錯(cuò)列雙圓柱來說,不同風(fēng)向角下的St差別不大,都在0.17左右,但功率譜(PSD)峰值卻差別很大.圖7是不同風(fēng)向角下下游圓柱的升、阻力功率譜曲線.從圖7可以看出,下游圓柱升、阻力功率譜隨風(fēng)向角的變化與如圖4所示下游圓柱脈動(dòng)氣動(dòng)力系數(shù)的變化趨勢一致.

        圖7 氣動(dòng)力系數(shù)的功率譜峰值隨風(fēng)向角的變化

        Fig.7 PSD peak of aerodynamic coefficient as a function of incidence angles

        當(dāng)β=0°和5°時(shí),下游圓柱的升力功率譜峰值最大,說明此時(shí)渦脫強(qiáng)度很大,圓柱表面出現(xiàn)很強(qiáng)的脈動(dòng)風(fēng)壓,進(jìn)而產(chǎn)生很大的脈動(dòng)升力;當(dāng)β=20°時(shí),下游圓柱的升力功率譜峰值最小,表明此時(shí)的渦脫強(qiáng)度明顯減弱,相應(yīng)的下游圓柱表面受到的脈動(dòng)升力也隨之降低.

        2.5 干擾流態(tài)分類

        通過上述分析可以發(fā)現(xiàn),不同的風(fēng)向角范圍內(nèi),下游圓柱的氣動(dòng)力有著很大差異.通過觀察不同風(fēng)向角時(shí)下游圓柱周圍的流場變化,可以總結(jié)出以下四種不同的流態(tài)(見圖8):

        (1) 流態(tài)1,旋渦撞擊流態(tài).當(dāng)β=0°時(shí),上游圓柱上側(cè)脫落的旋渦撞擊到下游圓柱上側(cè),接著上游圓柱下側(cè)脫落的旋渦撞擊到下游圓柱的下側(cè)(間隙側(cè)).上游圓柱脫落的旋渦對下游圓柱的不斷撞擊是下游圓柱表面出現(xiàn)較大脈動(dòng)風(fēng)壓和脈動(dòng)升力的主要原因.

        (2) 流態(tài)2,旋渦撞擊和剪切層干擾流態(tài).當(dāng)β=5°~10°時(shí),上游圓柱上側(cè)脫落的旋渦撞擊到下游圓柱的迎風(fēng)面,上游圓柱下側(cè)脫落的旋渦并沒有撞擊到下游圓柱,而是與下游圓柱下側(cè)的剪切層相互作用.同時(shí)存在的旋渦撞擊和剪切層干擾,導(dǎo)致下游圓柱表面出現(xiàn)很大的脈動(dòng)升力.

        (3) 流態(tài)3,剪切層干擾流態(tài).當(dāng)β=15°~30°時(shí),上游圓柱上側(cè)脫落的旋渦與下游圓柱下側(cè)(間隙側(cè))的剪切層和旋渦耦合產(chǎn)生相互作用,而上游圓柱下側(cè)脫落的旋渦在向下游運(yùn)動(dòng)的過程中由于距離下游圓柱較遠(yuǎn),故沒有發(fā)生明顯的相互作用.

        a 流態(tài)1(旋渦撞擊流態(tài),β=0°)

        b 流態(tài)2(旋渦撞擊和剪切層干擾流態(tài),β=5°~10°)

        c 流態(tài)3(剪切層干擾流態(tài),β=15°~30°)

        d 流態(tài)4(尾流干擾流態(tài),β=45°~90°)

        圖8 錯(cuò)列雙圓柱在不同風(fēng)向角下的干擾流態(tài)

        Fig.8 Interaction flow patterns of two staggered circular cylinders at different incidence angels

        (4) 流態(tài)4,尾流干擾流態(tài).當(dāng)β=45°~90°時(shí),上游圓柱與下游圓柱間距較大,上游圓柱上側(cè)脫落的旋渦在向下游運(yùn)動(dòng)的過程中只有一部分尾流會(huì)與下游圓柱發(fā)生相互作用,而下側(cè)脫落的旋渦則基本不會(huì)對下游圓柱有影響.

        2.6 平均升力的機(jī)理分析

        為了進(jìn)一步澄清風(fēng)向角20°附近下游圓柱出現(xiàn)平均升力峰值(即文獻(xiàn)[15]中的外側(cè)升力)的流場機(jī)理,本節(jié)從平均流場和瞬態(tài)流場兩方面進(jìn)行分析.

        2.6.1平均流場分析

        圖9是β=20°時(shí)下游圓柱的平均風(fēng)壓場和平均流速.從圖9可以看出,下游圓柱下側(cè)(間隙側(cè))的負(fù)壓強(qiáng)度明顯大于上側(cè),并且從下游圓柱的表面平均風(fēng)壓分布可知,下游圓柱的風(fēng)壓停滯點(diǎn)相對于來流向上側(cè)偏移了6°左右,這兩個(gè)因素對下游圓柱的平均升力均有一定貢獻(xiàn).此外,以往的學(xué)者還認(rèn)為[8],加速間隙流也是導(dǎo)致下游圓柱出現(xiàn)平均升力的原因之一,但從平均風(fēng)速圖(見圖9b)和下游圓柱附近的局部風(fēng)速圖(見圖9c)可看出,兩個(gè)圓柱間隙中的風(fēng)速并未有明顯的加速現(xiàn)象,并且下游圓柱上、下側(cè)表面的最大平均風(fēng)速也未見明顯的差別.上述結(jié)果說明,在該工況下間隙流并未有明顯的加速現(xiàn)象出現(xiàn),因而對下游圓柱平均升力的貢獻(xiàn)不大.

        a 平均風(fēng)壓場

        b 平均流速

        c 下游圓柱附近的平均流速

        圖9 平均風(fēng)壓場、平均流速和下游圓柱的局部流速

        Fig.9 Mean pressure coefficient field, mean velocity ratio and local velocity ratio of downstream cylinder

        2.6.2瞬態(tài)流場分析

        從第2.6.1節(jié)分析可知,下游圓柱上、下側(cè)表面的平均風(fēng)壓存在較大差異,這是下游圓柱受到較大平均升力的主要原因之一.為了進(jìn)一步澄清該現(xiàn)象的流場機(jī)理,圖11給出了風(fēng)向角為20°時(shí)九個(gè)瞬態(tài)時(shí)刻(T1~T9)的渦量圖及下游圓柱表面的風(fēng)壓系數(shù)分布.這九個(gè)瞬態(tài)時(shí)刻在升力系數(shù)時(shí)程曲線上的具體位置如圖10所示,下游圓柱的升力系數(shù)在T1和T9時(shí)刻接近極大值,在T6時(shí)刻則接近極小值.

        由圖10還可知,上、下游圓柱的升力時(shí)程曲線均有明顯的周期性,但兩者有一定的相位差;對于下游圓柱而言,升力系數(shù)從波峰下降至波谷的時(shí)間較長,而從波谷上升至波峰的時(shí)間較短.

        圖10 β=20°時(shí)升力系數(shù)的時(shí)程曲線及瞬態(tài)時(shí)刻

        由圖11可見,上游圓柱的上、下表面交替形成分離的剪切層,剪切層受另一側(cè)旋渦的影響而發(fā)生卷起,卷起的剪切層最終脫離圓柱并在尾流中形成方向交替變化的旋渦.隨著剪切層的變化,上游圓柱表面的風(fēng)壓也發(fā)生相應(yīng)的交替變化.上游圓柱的繞流場和表面風(fēng)壓的變化與單圓柱相似.

        下游圓柱受到上游圓柱尾流旋渦的影響,導(dǎo)致其氣動(dòng)力特性與單圓柱有很大差別.由圖11可見,從上游圓柱脫落的旋渦與下游圓柱的下側(cè)剪切層(或旋渦)發(fā)生相互作用.相對來說,上游圓柱上側(cè)脫落的旋渦與下游圓柱的相互作用更為劇烈(T1、T2、T8和T9時(shí)刻),而上游圓柱下側(cè)脫落的旋渦對下游圓柱的影響則相對較小(T4、T5和T6時(shí)刻).

        特別在T8和T9時(shí)刻,當(dāng)上游圓柱上側(cè)脫落的旋渦移動(dòng)至下游圓柱的下側(cè)剪切層附近時(shí),由于旋渦與剪切層的旋轉(zhuǎn)方向不同,上游圓柱的上側(cè)旋渦吸引下游圓柱的下側(cè)剪切層從圓柱表面脫落而形成旋渦;隨后,這一對方向相反的旋渦加速向下游移動(dòng)(T9時(shí)刻).由于下游圓柱的下側(cè)旋渦迅速離開圓柱,使得下側(cè)旋渦對下游圓柱上側(cè)剪切層的影響減弱,從而導(dǎo)致下游圓柱上側(cè)剪切層的卷起現(xiàn)象不明顯.在T1和T2時(shí)刻能觀察到類似現(xiàn)象.因而,下游圓柱的上側(cè)剪切層在所有時(shí)刻中均保持較為平順的形態(tài),上側(cè)剪切層對圓柱的影響相對較小,所以下游圓柱上側(cè)表面的負(fù)壓在所有時(shí)刻中均保持在相對較低的數(shù)值,最終導(dǎo)致下游圓柱上、下側(cè)表面的平均風(fēng)壓出現(xiàn)較大差異,這是下游圓柱受到較大平均升力作用的原因之一.

        圖11 β=20°時(shí)一個(gè)周期內(nèi)的瞬時(shí)渦量及表面風(fēng)壓分布

        此外,由圖11還可見,T4~T7時(shí)刻的下游圓柱風(fēng)壓停滯點(diǎn)明顯向圓柱上側(cè)偏移,而T1和T9時(shí)刻則向下側(cè)偏移,T2、T3和T8時(shí)刻則偏移不明顯,因而總體上下游圓柱的停滯點(diǎn)是向上偏移,這是導(dǎo)致下游圓柱受到較大平均升力的另一個(gè)原因.

        3 結(jié)論

        (1)隨著風(fēng)向角的增大,兩個(gè)圓柱的干擾流態(tài)依次為旋渦撞擊流態(tài)、旋渦撞擊和剪切層干擾流態(tài)、剪切層干擾流態(tài)、尾流干擾流態(tài).在不同流態(tài)下,下游圓柱的氣動(dòng)力特性有很大差異.

        (2)在風(fēng)向角為0°~10°的旋渦撞擊流態(tài)以及旋渦撞擊和剪切層干擾流態(tài)時(shí),從上游圓柱脫落的旋渦與下游圓柱發(fā)生強(qiáng)烈的撞擊,導(dǎo)致下游圓柱的脈動(dòng)升力和迎風(fēng)面的脈動(dòng)風(fēng)壓遠(yuǎn)遠(yuǎn)大于其他流態(tài),也遠(yuǎn)大于單圓柱.

        (3)在風(fēng)向角為15°~30°的剪切層干擾流態(tài)下,下游圓柱受到顯著的平均升力作用.下游圓柱風(fēng)壓停滯點(diǎn)的偏移、上游圓柱的尾流旋渦與下游圓柱間隙側(cè)剪切層和旋渦的相互作用,是下游圓柱出現(xiàn)平均升力的兩個(gè)主要原因.

        猜你喜歡
        旋渦流態(tài)升力
        高速列車車頂–升力翼組合體氣動(dòng)特性
        側(cè)邊機(jī)組故障對泵站前池流態(tài)的影響
        無人機(jī)升力測試裝置設(shè)計(jì)及誤差因素分析
        小心,旋渦來啦
        基于自適應(yīng)偽譜法的升力式飛行器火星進(jìn)入段快速軌跡優(yōu)化
        大班科學(xué)活動(dòng):神秘的旋渦
        旋渦笑臉
        山間湖
        改進(jìn)邊界條件的非恒定流模型在城市河流橡膠壩流態(tài)模擬中的應(yīng)用
        升力式再入飛行器體襟翼姿態(tài)控制方法
        日韩精品人妻中文字幕有码在线| 国产福利酱国产一区二区| 久久国产综合精品欧美| 亚洲av天堂久久精品| 人妻少妇偷人精品一区二区三区| 日韩精品成人区中文字幕| 国产国产裸模裸模私拍视频| 亚洲日韩中文字幕一区| 精品少妇爆乳无码aⅴ区| 亚洲视频一区二区三区免费 | 亚洲av永久一区二区三区| 亚洲国产精品不卡av在线| 大胸少妇午夜三级| 人妻无码一区二区三区四区 | 精品日产卡一卡二卡国色天香 | 无码人妻精品一区二区三区免费| 九九在线精品视频xxx| 国产av在线观看91| 媚药丝袜美女高清一二区| 精品偷拍被偷拍在线观看| 亚洲精品国偷自产在线99正片| 人妻精品一区二区三区视频| 国内偷拍第一视频第一视频区| 中文字幕精品一区二区三区| 人妻夜夜爽天天爽三区麻豆av网站| 亚洲精品久久久久中文字幕二区| 无码高潮少妇毛多水多水免费| av一区二区在线免费观看| 性欧美丰满熟妇xxxx性久久久 | 99久久超碰中文字幕伊人| 狠狠色噜噜狠狠狠狠888奇禾| 91极品尤物在线观看播放| 国产亚洲专区一区二区| 国产高清在线精品一区app| 97在线观看| 天啦噜国产精品亚洲精品| 亚洲一区二区三区麻豆| 国产一区二区三区中文在线| 国产农村乱辈无码| 福利视频一二区| 亚洲中文字幕第一第二页|