999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

下擊暴流作用下超大型冷卻塔風場驅動機理與風荷載極值模型

2022-12-01 09:28:44韓光全柯世堂李文杰任賀賀
振動與沖擊 2022年22期
關鍵詞:風速

韓光全,柯世堂,楊 杰,李文杰,任賀賀

(南京航空航天大學 民航學院 土木與機場工程系,南京 211106)

下擊暴流[1]又稱雷暴沖擊風,其主要特征是強下沉氣流和近地面爆發狀強風。雷暴天氣中,下落雨滴伴隨著濕冷空氣產生局部強下行氣流,撞擊地面后沿四周迅速擴散,在近地表面形成放射狀的強風,該流動模式通常被認為反向的龍卷風。不同于大氣邊界層良態風場,下擊暴流最大風速產生在距離地面約50~100 m 高度處[2],且氣象記錄表明,其最大風速可達75 m/s[3]。同時,下擊暴流在核心射流區域和突出前緣均存在明顯的垂直分量,在與強切變的水平風速共同作用下,造成了建構筑物的嚴重損傷甚至破壞[4-5]。超大型冷卻塔作為火/核電廠重要的基礎設施,是典型的高聳旋轉薄殼建筑,風荷載是其結構設計的控制荷載之一。然而,GB/T 50102—2014《工業循環水冷卻設計規范》[6]僅給出良態風作用下的荷載,對下擊暴流等特異風作用下冷卻塔風荷載分布模式未給予考慮。因此,對下擊暴流作用下冷卻塔風場作用機理和風荷載特性進行研究具有現實意義,可為超大型冷卻塔設計和建設提供重要的參考依據。

對于下擊暴流三維風場,限于其短時突發性[7]與移動性[8],尚且難以全面即時地進行現場實測,目前國內外學者主要采用物理試驗和數值模擬等手段對其進行研究。物理試驗一般采用主動控制風洞[9-10]或沖擊射流裝置[11-13],對結構表面風壓和空氣動力學特性進行研究,但是由于難以滿足縮尺比等試驗條件,物理試驗難以對復雜結構形式的建構筑物進行下擊暴流作用下的準確模擬。隨著計算機技術的發展,數值模擬逐漸被應用到下擊暴流三維風場研究中,利用數值模擬技術,可對風場及其對結構的作用機理[14-16]進行較為全面的探究。文獻[17-18]采用沖擊射流模擬裝置,對下擊暴流作用下高層和雙坡屋面建筑風場流動特性進行了分析,結果發現在風場一定位置處建筑表面風壓顯著高于ASCE 7-05標準[19]。文獻[20-21]分別借助可移動的沖擊射流裝置和滑移網格技術對移動型下擊暴流進行試驗和數值模擬,研究表明,移動的下擊暴流作用于建筑物上會使建筑表面出現整體變化,甚至超出規范良態風風壓限值。文獻[22-23]研究發現下擊暴流會引起輸電塔間縱向拉力急劇增大或者關鍵結構構件的局部失效,進而引起失穩倒塌。對下擊暴流作用下大跨度橋梁[24-25]、風力機[26-27]、屋蓋廠房[28-30]等結構的研究也顯示相比于良態風,下擊暴流會不同程度增加結構物的潛在破壞性。而目前鮮有涉及下擊暴流三維脈動風場對超大型冷卻塔的風場作用機理與風荷載特性,亟需開展特異風作用下超大型冷卻塔風荷載取值研究。

鑒于此,本文基于沖擊射流模型與大渦模擬(large eddy simulation,LES)技術,首先對下擊暴流三維風場進行了全生命周期非定常數值模擬,將模擬結果與經驗模型和實測數據進行對比驗證,并闡釋其風場特性;在此基礎上,詳細分析了超大型冷卻塔在下擊暴流不同徑向位置處的繞流特性,升/阻力系數分布特征,以及瞬態風壓和極值風荷載分布規律,結論可為此類超大型冷卻塔抗雷暴設計提供一定的參考。

1 工程概況

該在建世界最高超大型冷卻塔位于我國內蒙地區,設計塔高突破規范190 m的限定,冷卻塔主體結構由塔頂剛性環、塔筒、支柱、環基及加勁肋構成。其中,塔筒采用雙曲線型鋼筋混凝土殼體結構,其厚度呈指數變化,最小與最大壁厚分別為0.42 m和2.45 m,塔筒外表面沿環向均勻布置120條梯形加勁肋,60榀X型支柱通過60個支墩與底部鋼筋混凝土環板型基礎承臺連接,單個支墩下布24根長30 m直徑0.8 m的鋼筋混凝土灌注樁,冷卻塔主要結構尺寸參數如表1所示。

表1 冷卻塔主要結構尺寸Tab.1 The main structural sizes of cooling tower

2 數值模擬計算概況

2.1 數值模型與工況設置

為全面監測冷卻塔各處風壓,在冷卻塔模型內、外表面沿子午向布置14層測點,每層沿環向逆時針間隔12°布置一個風壓系數測點,內、外表面測點一一對應,則環向共布置60個測點,塔筒內、外表面共布置840個測點。

CFD數值模擬采用三維足尺模型,沖擊射流速度入口直徑Djet=600 m,射流高度為2Djet,為保證流動的充分發展,計算域尺寸為10Djet×10Djet×3Djet。邊界條件設置如下:射流入口采用速度入口,射流初始風速Vjet=29 m/s,湍流強度為1%;計算域四周及頂面采用壓力出口邊界條件,回流湍流強度同為1%;地面采用無滑移壁面,射流口上部圓柱面設為滑移壁面。計算域及邊界條件如圖1(a)所示。

考慮到下擊暴流風場與冷卻塔結構特性,為探究冷卻塔在風場不同位置處的風荷載分布特性,共設置4個典型計算工況。其中:工況1為冷卻塔置于射流入口正下方;工況2~工況4分別為塔筒底部迎風面前端距射流中心徑向距離r=1.0Djet,r=1.5Djet,r=2.0Djet。工況設置如圖1(b)所示。

圖1 計算域邊界條件及工況設置示意圖Fig.1 Calculation domain boundary conditions and working condition settings

2.2 網格劃分與參數設置

假定沖擊射流沿高度方向均為尺寸統一的圓形,而O型拓撲適合此類兩端均為圓形截面的物體,故在射流中心區域采用O-Block二次拓撲生成雙層O型網格。同時為滿足冷卻塔結構周圍的模擬精度,采用混合網格離散形式,將整個計算域分為內、外兩個部分:核心區域采用四面體網格,并對冷卻塔加勁肋周圍網格進行局部加密,外圍區域采用規則拓撲的高質量六面體結構網格,細部及整體網格劃分如圖2所示。

圖2 細部及整體網格劃分示意圖Fig.2 Diagram of detailed and overall grid subdivision

不同網格數量下r=1.0Djet徑向位置處的風剖面,如圖3所示。由圖3可知,當網格數量達到1 200萬時,在上部區域風剖面隨網格數量的增加無明顯變化,即此時數值模型已基本滿足與網格數量無關的要求。考慮到計算效率和準確性,后續研究中使用的網格數為1 600萬。

圖3 不同網格數量下風剖面對比示意圖Fig.3 Comparison of wind field results with different gird numbers

此網格方案中近壁面網格高度為0.001 m,壁面法向網格的增長率為1.08,寬度為15 m,故近壁面最大網格高寬比為6.67×10-5。模擬采用的時間步長為Δt=0.03 s,冷卻塔表面最小網格尺寸為Δx=1.2 m,由

Courant=Δt·Δx/Vjet

(1)

計算可得Courant number為0.725,其值小于1,可判定時間步長與網格是匹配的。

基于Fluent程序進行計算,在LES模擬中,亞格子模型采用algebraic wall-modeled LES(WMLES)模型以彌補近壁面網格數量不足的問題,從而實現雷暴對流系統的高精度模擬需求。為生成隨時間變化的入口邊界條件,采用Vortex method,Number of Vortices設置為190;雷諾應力定義方法采用k or turbulent intensity。速度-壓力耦合方法采用壓力耦合方程組的半隱式方法(semi-implicit method for pressure linked equations,SIMPLE),梯度插值方案采用least squares cell based,壓力采用standard格式,對動量的計算采用second order upwind格式以實現高階離散化。質量殘差設置為10-6,x,y,z方向速度收斂標準設置為10-3。在計算時,采用Data Sampling for Time Statistics,sampling interval設置為1。

3 下擊暴流風場模擬

3.1 有效性驗證

首先對沒有冷卻塔結構的下擊暴流風場進行模擬,采用無量綱時間T=Δt·Vjet/Djet,為確保風場模擬的精確和完整性,無量綱時間步長為0.001 45,共模擬20 000步,總模擬時長T=29。

r=1.0Djet徑向位置處徑向風速豎向風剖面和z=0.2Djet高度處徑向風速徑向風剖面,分別如圖4和圖5所示。其中:圖4橫坐標為徑向風速u與最大徑向風速umax的比值,徑向風速u為下擊暴流發展成熟階段(T=8.7~14.5)瞬時風速的平均值,縱坐標為離地高度z與最大徑向風速對應高度zmax的比值;圖5橫坐標為徑向距離r與最大徑向風速對應徑向距離rmax的比值。從圖5中可知,數值模擬得出的歸一化風剖面與各經驗模型[31-34]及實測數據吻合度較高,說明本文采用沖擊射流模型和LES技術可以實現對下擊暴流風場的有效模擬。

圖4 徑向風速歸一化豎向風剖面示意圖Fig.4 The vertical profile of the radial velocity

圖5 徑向風速歸一化徑向風剖面示意圖Fig.5 The radial profile of the radial velocity

3.2 風場特性

下擊暴流部分時刻沿流域縱剖面的風速云圖,如圖6所示。T=5.8時刻為下擊暴流發生初期階段,此時沖擊風自射流入口向下方運動,氣流前端尚未到達地面,極值風速為33.51 m/s;隨著時間的推移,在T=8.7時刻,沖擊風前緣已經抵達地面,與地面接觸并開始徑向對流,觸發邊界層氣流的非定常分離,形成反向旋轉、徑向移動的第一個渦環,極值風速增大到41.00 m/s;在射流入口逐漸有穩定氣流進入的情況下,于第一個渦環之后產生一系列尺度較小的渦環,極值風速也增大到約47.30 m/s左右,如圖6(c)、圖6(d)所示。

圖6 下擊暴流發展階段不同時刻的縱剖面風速云圖Fig.6 Wind speed nephogram of longitudinal section at different times of downburst development stage

不同徑向位置處近地面風速時程曲線,如圖7所示。在r=0徑向位置處,豎向風速為主導分量,隨著沖擊風的持續輸入,近地面各測點風速從0逐漸增大,在約T=7.25時刻達到最大值,10~60 m高度處各測點極值風速分別約為1 m/s,2 m/s,3 m/s,4 m/s,5 m/s,6 m/s,而后略有減小并維持穩定波動狀態。

在r=1.0Djet,r=1.5Djet,r=2.0Djet徑向位置處,水平風速為主導分量,由于氣流撞擊地面并向四周逐漸擴散,各徑向位置處近地面測點風速首先達到首個峰值,且達到首個峰值的時間隨徑向距離的增加而延遲,首個峰值的最大值分別約為38 m/s,50 m/s,22 m/s。隨著徑向移動、反向旋轉渦環的不斷生成,風速相應呈現波動趨勢,形成一系列“波峰”與“波谷”。其中r=1.0Djet,r=1.5Djet兩處一系列風速“波峰”與“波谷”的數值較為穩定,穩定“波峰”的峰值均維持在40 m/s左右,r=1.5Djet處“波谷”風速小于r=1.0Djet處,可見r=1.5Djet處渦環運動相較于r=1.0Djet處更為劇烈。r=2.0Djet徑向位置由于遠離入流中心,湍流發展紊亂且能量削弱,一系列風速“波峰”與“波谷”呈現逐漸衰減的趨勢。

圖7 近地面風速時程曲線Fig.7 Time history curves of near-surface wind speed

4 結構流場與風荷載特性

4.1 繞流特性

由圖7可知,各工況在T=15.0之后風速基本維持穩定變化,即下擊暴流在此時刻后為發展成熟階段,故圖8給出了在T=17.4時刻沿y=0縱剖面冷卻塔周圍風壓云圖與空氣流線圖。工況1中,下沖氣流在塔筒頂部產生分離,進入塔筒內部的氣流從X型支柱處流出,與外部下行氣流產生對流,從而在塔筒外部形成大尺度渦旋,渦旋處風壓相對周圍區域較小,而塔筒內部風壓較高且分布均勻;工況2~工況4中,流動的復雜性相比工況1愈加明顯,近地面氣流一部分沿塔筒外表面爬升,與上部下行氣流在迎風面塔頂沖擊,形成渦旋;另一部分近地面氣流從X型支柱處進入塔筒內部,撞擊塔筒內壁后形成渦旋,此3個工況中塔筒迎風區外表面和背風區內表面為高風壓區。

圖8 沿y=0剖面冷卻塔周圍風壓云圖與流線圖Fig.8 Pressure cloud and streamline diagram around the cooling tower along the y=0 section

為研究在水平方向上冷卻塔結構內、外空氣流線變化,圖9給出了在T=17.4時刻沿z=40 m橫剖面冷卻塔周圍風壓云圖與空氣流線圖。由圖9可知:工況1中由于沖擊風豎直流下,故塔筒內部幾乎沒有水平方向氣流,塔筒外部氣流呈現輻射狀發展,且存在大小不一的漩渦;工況2~工況4中,高壓區同樣主要分布在塔筒迎風面外部和背風面內部,隨著冷卻塔迎風面與入流口的距離越來越遠,高壓區的范圍和大小隨之逐漸變小,并且塔筒前端氣流由集束狀態逐漸轉向水平狀態,塔筒內部和背風面尾流區域均存在漩渦。

圖9 沿z=40 m剖面冷卻塔周圍風壓云圖與流線圖Fig.9 Pressure cloud and streamline diagram around the cooling tower along the z=40 m section

4.2 風壓系數瞬態分布

定義風壓系數為

(2)

式中:P為測點的風壓;Pref為參考靜壓力;ρ為空氣密度;Vref為參考風速,在常規大氣邊界層風場中,通常取為模型頂部平均風速,而在下擊暴流風場模擬中,由于徑向風速最大值一般出現在近地面附近,且大小隨徑向位置的變化而變化,因而參考風速一般取為射流入口的初始風速Vjet。

工況1中冷卻塔位于沖擊風正下方,塔筒內、外表面風壓系數沿環向分布規律基本一致。因此,圖10給出了頂部、喉部以及底部3個典型斷面測點風壓系數時程曲線。由圖10可知:塔筒內、外表面3個不同斷面時程風壓系數變化趨勢基本相同,在初始時間范圍內,風壓系數先減小再增大,在約T=7.25時刻達到最大值,而后趨于近似正弦曲線的波動狀態,外表面底部、喉部、頂部測點風壓系數分別在0.8,0.7,0.9附近波動,波動趨勢稍顯雜亂,而內表面相應測點風壓系數分別在0.9,0.7,0.8附近波動,波動趨勢更為平穩;由于氣流到達塔筒各位置的時間不同,故不同斷面測點風壓系數在正弦曲線變化階段存在相位差,此現象在外表面更為明顯。不同于良態風作用下冷卻塔風壓分布特征,此工況下喉部斷面風壓系數整體最小。

工況2~工況4中塔筒底部、喉部、頂部典型斷面在迎風區、側風區、背風區測點的風壓系數時程曲線,分別如圖11~圖13所示。由圖11~圖13可知,各工況下相應測點時程風壓系數變化規律相似。以工況2為例分析如下:

底部外表面迎風區測點的風壓系數由0.8減小至0.3左右,隨后增加達到首個峰值,而側風區和背風區測點則是由0.8不斷減小至負壓峰值,各測點達到首個峰值的時刻與r=1.0Djet處風速達到首個峰值的時刻相吻合,此后風壓系數分別在1.0,-0.2,-0.1左右維持穩定小幅震蕩。

圖10 工況1典型斷面測點風壓系數時程曲線Fig.10 Time history curves of wind pressure coefficient at typical cross-section with working condition 1

圖11 工況2典型斷面測點風壓系數時程曲線Fig.11 Time history curves of wind pressure coefficient at typical cross-section with working condition 2

圖12 工況3典型斷面測點風壓系數時程曲線Fig.12 Time history curves of wind pressure coefficient at typical cross-section with working condition 3

圖13 工況4典型斷面測點風壓系數時程曲線Fig.13 Time history curves of wind pressure coefficient at typical cross-section with working condition 4

喉部測點風壓系數由0.7左右逐漸減小,而后在0左右維持小幅震蕩。但由于渦環主要發生在近地面,故峰值現象明顯弱于底部,同理可為頂部測點時程風壓系數變化作出解釋。由于塔內受氣流湍流擾動程度較小,故在首個峰值時刻后,各斷面內表面測點風壓系數震蕩較外表面更為平緩。

另一方面,對比不同工況測點時程風壓系數差異可知:達到首個峰值的時刻隨著徑向距離的增加而延遲,分別對應該徑向位置處風速達到首個峰值的時刻;由于r=1.5Djet徑向位置處湍流發展紊亂,故工況3中風壓系數時程曲線震蕩幅度較大,且底部外表面測點峰值絕對值較大。

4.3 升/阻力系數

冷卻塔升/阻力系數由風壓系數積分計算得到[35]

(3)

(4)

式中:Ai為第i測點壓力等效覆蓋面積;θi為第i測點壓力與來流方向夾角;A為整體結構沿順風方向投影面積。由于工況1的特殊性,在此不對其進行升/阻力系數計算。

工況2塔筒底部、喉部、頂部3個典型測點層的升/阻力系數時程曲線,如圖14所示。由圖14可知,各測點層的升/阻力系數在下擊暴流初始時間范圍內均維持在0左右,隨后產生波動,波動的幅度有所差別。整體上,外表面各層的升/阻力系數波動幅度大于內表面,底部測點層的阻力系數波動幅度最大。其中:外表面均為正值,最大可達到1.5;內表面均為負值,最大可達到-0.5。工況3、工況4時程升/阻力系數變化規律與工況2相似,在此不再附圖。

對內、外表面各測點層的時程升/阻力系數進行平均,如圖15所示。由圖15可知,除工況2中冷卻塔喉部以下外表面測點層外,內外表面各層的平均升力系數基本為0,外表面阻力系數自塔頂沿塔高方向逐漸增大,在塔底達到正值最大值,內表面阻力系數變化趨勢相似,在塔底達到負值最大值。隨著冷卻塔距離射流中心越來越遠,塔底阻力系數數值由大變小,即工況2中阻力系數數值最大,外表面阻力系數為0.5,內表面阻力系數為-0.2。

圖14 典型測點層升/阻力系數時程曲線Fig.14 Time history curves of lift/drag coefficient of typical measuring point layer

4.4 極值風荷載

下擊暴流對超大型冷卻塔具有強烈的沖擊性,將冷卻塔各典型斷面測點內、外表面時程風壓系數進行疊加,并取其峰值,如圖16所示。同時給出了考慮內、外壓的規范良態風作用下風壓分布曲線,并以此計算得到考慮極值風效應的風壓系數包絡值:-0.92~2.89。由圖16分析可知:

(1) 工況1中各斷面測點總風壓系數峰值基本在-0.5上下小幅波動,處在包絡范圍內,究其原因為冷卻塔壁厚較薄,此工況中氣流沿筒壁近似豎直流下,故內、外表面瞬時風壓系數相差較小。此工況下按規范進行結構設計偏于安全。

(2) 工況2~工況4中各斷面測點總風壓系數峰值呈現明顯的對稱分布特性,沿環向分布規律與規范良態風曲線相似,均未超出包絡正值2.89,但均存在超出包絡負值-0.92的區域。其中:工況2、工況3底部背風區總風壓系數峰值分別為-1.46,-1.52,超出包絡范圍分別為58.70%,65.22%;另外工況3、工況4中塔筒中部和喉部斷面在側風區同樣存在小幅度超出包絡范圍的區域。

下擊暴流短時突發形成的渦環沖擊作用在冷卻塔內、外表面,極有可能引起瞬時極值風荷載超出規范限值,進而易引起結構的破壞,需引起設計人員的注意。

圖15 各測點層平均升/阻力系數分布曲線Fig.15 Distribution curve of average lift/drag coefficient of measuring point layers

圖16 塔筒典型斷面測點總風壓系數峰值曲線圖Fig.16 Curves of peak values of total wind pressure coefficient of typical tower sections

5 結 論

(1) 采用沖擊射流模型和大渦模擬技術模擬得出的徑向風速歸一化風剖面與經驗模型、實測數據吻合度較高;下擊暴流沖擊氣流與地面接觸碰撞后,產生一系列徑向移動,反向旋轉的渦環,在r=1.0Djet,r=1.5Djet,r=2.0Djet徑向位置處風速隨之呈現波動變化趨勢。在射流初始風速為29 m/s的情況下,r=1.5Djet處極值風速最高可達50 m/s。

(2) 在下擊暴流發展成熟階段,工況1中冷卻塔塔筒外表面附近區域存在大尺度渦旋,生成渦旋部位風壓小于周圍區域;而在工況2~工況4中渦旋主要存在于塔筒內部、迎風面塔頂區域和背風面近地尾流區,迎風區外表面和背風區內表面風壓高于其他區域。

(3) 塔筒內、外表面風壓系數時程分布三維效應明顯,工況1中典型區域測點風壓系數隨時間先減小后增大,隨后在下擊暴流發展成熟階段維持穩定的正弦分布狀態;而工況2~工況4風壓系數脈動趨勢更為顯著,由大幅震蕩發展成為維持在固定數值附近的小幅震蕩,其中底部區域測點受渦環影響震蕩尤為明顯。

(4) 層平均升力系數較小,工況2中塔底層平均阻力系數最大,在內、外表面分別可達到-0.2,0.5;工況1處在規范極值風壓包絡范圍內,工況2~工況4均有區域超出包絡范圍,其中工況2、工況3的底部區域超出幅度分別達到58.70%,65.22%。

綜上所述,下擊暴流對超大型冷卻塔的風場驅動主要表現為近地面渦環的沖擊作用,其瞬態極值風壓發生位置從規范良態風作用下的喉部側風區轉移到底部背風區,同時最大數值較規范良態風作用增大65.22%。主要結論可為此類特異風作用下超大型冷卻塔結構抗風設計提供科學依據。

猜你喜歡
風速
邯鄲市近46年風向風速特征分析
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
基于時間相關性的風速威布爾分布優化方法
陜西黃土高原地區日極大風速的統計推算方法
陜西氣象(2020年2期)2020-06-08 00:54:38
基于GARCH的短時風速預測方法
快速評估風電場50年一遇最大風速的算法
風能(2016年11期)2016-03-04 05:24:00
考慮風切和塔影效應的風力機風速模型
電測與儀表(2015年8期)2015-04-09 11:50:06
GE在中國發布2.3-116低風速智能風機
考慮風速分布與日非平穩性的風速數據預處理方法研究
主站蜘蛛池模板: 999在线免费视频| 国产日韩AV高潮在线| 伊人福利视频| 久久精品丝袜| 国产亚洲一区二区三区在线| 热99精品视频| 国产高潮流白浆视频| 亚洲国产成人久久精品软件| 99视频只有精品| 免费A级毛片无码免费视频| 在线观看国产精美视频| 777国产精品永久免费观看| 男女男免费视频网站国产| 特级毛片免费视频| 日韩中文无码av超清| 尤物国产在线| 亚洲熟女中文字幕男人总站| 亚洲欧美日韩高清综合678| 日本三级精品| 亚洲欧美不卡视频| 国产精品视频a| 国产在线观看一区精品| 四虎国产精品永久在线网址| 任我操在线视频| 国产永久无码观看在线| 午夜a级毛片| 1769国产精品免费视频| 国产精品自在线天天看片| 国产特一级毛片| 国产一区二区丝袜高跟鞋| 啊嗯不日本网站| 五月天久久综合| 亚洲三级a| h视频在线播放| 久久精品波多野结衣| 全部免费特黄特色大片视频| 中文字幕在线日韩91| 亚洲欧美日韩色图| 久久精品视频一| 国产精品视频系列专区| 久久这里只精品国产99热8| 免费人成又黄又爽的视频网站| 国产成人精品优优av| 国产麻豆福利av在线播放| 欧美国产综合视频| 日韩在线欧美在线| 日韩亚洲高清一区二区| 影音先锋丝袜制服| 99久久精品美女高潮喷水| 91区国产福利在线观看午夜 | 亚洲成人黄色网址| 国产拍揄自揄精品视频网站| 宅男噜噜噜66国产在线观看| 蜜芽国产尤物av尤物在线看| 无码精品一区二区久久久| 久久semm亚洲国产| 亚洲成人一区二区| 国产精品va免费视频| 九色91在线视频| 在线视频亚洲欧美| 中国一级特黄视频| 美美女高清毛片视频免费观看| 激情无码视频在线看| 久久精品女人天堂aaa| 国产精品毛片一区视频播 | 国产精品丝袜在线| 国产在线一区视频| 熟女成人国产精品视频| 久久香蕉国产线| 伊人激情久久综合中文字幕| lhav亚洲精品| 思思99思思久久最新精品| 麻豆精品在线| 亚洲AV无码乱码在线观看代蜜桃| 国产亚洲精久久久久久无码AV| 激情六月丁香婷婷| 尤物在线观看乱码| 免费 国产 无码久久久| 久久狠狠色噜噜狠狠狠狠97视色 | 国产玖玖玖精品视频| 视频一本大道香蕉久在线播放 | 日韩无码真实干出血视频|