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

風沙流對戈壁地區擋風墻響應規律的數值模擬分析

2016-05-13 04:00:25李曉軍蔣富強
鐵道標準設計 2016年3期

李曉軍,蔣富強

(1.蘭州交通大學土木工程學院,蘭州 730070; 2.中鐵西北科學研究院有限公司,蘭州 730000)

?

風沙流對戈壁地區擋風墻響應規律的數值模擬分析

李曉軍1,2,蔣富強2

(1.蘭州交通大學土木工程學院,蘭州730070; 2.中鐵西北科學研究院有限公司,蘭州730000)

摘要:基于FLUENT歐拉雙流體模型,對蘭新鐵路沿線既有擋風墻周圍風沙兩相流運動特性進行數值模擬,得到擋風墻背風側的流場分布特點以及積沙情況。結果表明:擋風墻背風側風速廓線變化規律呈指數增長趨勢,在0.5 m至擋風墻自身高度區間內變化較為復雜,呈先減小后增加的趨勢;擋風墻背風側近地表氣流速度反向增大后沿著初始速度的方向減小為0且繼續增大至初始速度大小,風速最大值增加的幅度保持在50%左右,風速越大,氣流的削弱作用越明顯;當初始氣流速度為較小時,線路上積沙較少,沙粒多數堆積在擋風墻背風側墻角處;隨著風速的增加,單位時間內通過擋風墻的沙粒增多,由于過流斷面減小,氣流擴散,更多沙粒沉積在線路上;在強風地區,布設擋風墻時應考察線路上風向的地表情況,沙源比較豐富時應采用工程治沙措施來減小風沙流密度,達到防沙的目的。

關鍵詞:蘭新鐵路;擋風墻;風沙兩相流;數值模擬

蘭新鐵路沿天山山前戈壁地帶穿過,絕大部分路段處在干旱荒漠風沙地帶,在寒潮天氣的影響下,形成了特殊的大風天氣,風力強勁,風向多變,大風頻繁,途經安西風區、煙墩風區、百里風區、三十里風區和達坂城風區等五大風區,尤其以百里風區風力最為強勁,風速最高可達60 m/s,局部地段每年有200 d風力在8級以上[1-2]。在強氣流運動作用下,戈壁地表嚴重失穩,風沙災害頻發,造成路基風蝕和線路積沙等危害,運輸時常迫停,每年因大風停運造成的直接經濟損失達數億多元。通過調查研究在鐵路沿線迎風側修建擋風墻是鐵路安全運營的一項重要措施[3-5]。

對風沙流的研究已經有了很長的一段歷史,近幾年來,國內外眾多學者采用風洞實驗及數值模擬等手段對路基以及擋風墻的功效進行了研究,并得到了不少成果[6-9]。文獻[10]通過對路堤周圍風沙兩相流運動特性進行數值模擬研究,發現沙粒速度與風速相互影響形成一種反饋機制并揭示了路基迎風坡積沙量大于背風坡的形成機理。文獻[11]通過對戈壁區路基周圍風沙流進行數值模擬研究,發現在路基周圍有氣流的相對高速區和相對低速區,在相對高速區產生風蝕沙害,在相對低速區產生積沙。文獻[12]通過實驗測定擋風墻的表面壓力分布、最佳疏透度、路基高度和不同軌心距的擋風墻合理高度等指標,發現擋風墻工程以高3 m緊密結構矩形體為宜。文獻[13]通過實地調查分析與現場測試,發現3 m高擋風墻在墻后距離5 m以內存在低速渦旋積沙區,當平均風速不大于41 m/s且上風區沒有新的沙源補充時,通常設置的3 m高擋風墻基本可以阻擋風沙流的運動。

這些成果都是從擋風墻防止列車傾覆的角度去研究,并沒有考慮到擋風墻造成的沙害問題,事實證明,在擋風墻為機車的穩定性運營提供安全保障的同時,也會造成線路積沙。本文基于FLUENT軟件,模擬風沙兩相流途經擋風墻不同幾何條件時擋風墻背風側流場的變化情況,并對擋風墻背風側線路積沙形態進行研究,為鐵路防風防沙提供依據。

1 數值分析模型建立

蘭新線位于戈壁灘地段,沙源較為豐富,并且途經五大風區,風力強勁,大風經常卷起地面的沙子擊碎車窗玻璃,為此,蘭新線擋風墻的實際高度設定為3 m,正好位于客車車窗玻璃上方。同時,蘭新線上的風為季風,風向較為穩定,因此,蘭新線修建的是單側擋風墻[14]。具體參數設置過程如下。

1.1控制方程

由于鐵路沿線橫向風速一般小于50 m/s,馬赫數小于0.3,故流動按二維不可壓縮處理[15]。另外,本模型不考慮熱量的交換,是單純流場問題,不用包含能量方程。所以描述防風擋沙墻后的風場的控制方程主要包括連續方程、動量方程和湍流模型方程,本文采用的湍流方程式k-ε湍流模型。

動量方程

式中,ux、uy為速度分量;為速度矢量;ρ為密度; g為重力加速度;τ是因分子黏性作用而產生的作用在微元體表面上的黏性應力分量; P為流體微元體的壓強; t為時間。單位統一為國際單位。

連續方程

式中,k為湍動能;μt為湍動粘度;ε為湍動耗散率; u為速度分量; G由平均速度梯度引起的湍動能產生;σk、σε為相對應的普朗特系數; C為經驗常數。

1.2計算區域

為了讓氣流的繞流和流場發展充分,為避免路堤背風側渦旋流對出口邊界條件的影響,理論上計算區域應越大越好。應用Gambit軟件建立二維計算模型,通過試算,模型總體尺寸取為120 m×25 m。為節約計算資源,將路堤模型靠近風速入口40 m處,此時,擋風墻后面流場基本達到充分發展,計算區域更大時計算結果改變很小。模型中取路堤高度3 m,路堤邊坡坡率為1∶1.75,對在不同幾何條件下的擋風墻進行數值模擬。

1.3網格劃分及邊界設置

本文所建模型較為復雜,為節省工作量采用Pave法將區域劃分為非結構性網格,網格類型為Tri(三角形)網格形式。由風沙流的分布特點可知近地表的沙粒運動更加復雜,因此防風擋沙墻及其路堤周圍部分網格劃分較為密集,而在進口、出口和上側劃分較為稀疏。這樣劃分既考慮到模擬計算的收斂時間、計算機的計算能力,又充分考慮到了在重點部位劃分密集以保證模擬準確性的原則,從而嚴格控制了網格的數量,提高了計算精度和效率。網格整體數量超過了60萬,計算模型的網格劃分如圖1所示。入口采用速度入口邊界(VELOCITY_INLET),出口采用出口流動邊界(OUTFLOW),地面、擋風墻與其他邊界采用固體無滑移壁面邊界(WALL)。

圖1 計算模型的網格劃分

1.4FLUENT求解模型及參數選取

本次試驗采用歐拉兩相非定常模型、k-ε湍流模型(采用標準壁面函數)的標準格式。指定空氣為基本相,定義沙作為第二相,由于此次模擬戈壁地區風沙流,取沙相的體積分數為0.02[10]。風沙之間采用Schiller-naumann計算兩相之間的阻力,并考慮沙的重力。入口處的風沙流速度分別取10、20、30、40 m/s,風沙流中的沙粒粒徑一般在0.075~0.25 mm,本次計算取0.1 mm。由于采用歐拉雙流體計算模型,計算方法采用一階迎風格式。

2 數值模擬分析與討論

2.1擋風墻背風側流場特征分析

采用FLUEN軟件對給出的有關流場的基本方程數值求解后,可得到在不同條件下路堤周圍的流場特征。

圖2為擋風墻高度為3 m,風沙流初始速度為20 m/s時路堤周圍氣流速度矢量圖。

圖2 設置3 m高擋風墻時路堤周圍氣流速度矢量圖

按照能量守恒定理,一部分氣流通過擋風墻時速度增加,所攜帶的沙粒速度也增加,必然有一部分沙粒速度降低,沉積在擋風墻背風側,造成線路積沙,給行車運營帶來不便。從圖2中不難看出,受擋風墻和路堤影響,周圍氣流速度發生改變。當氣流途經路堤和擋風墻時,氣流速度形成分區:隨著氣流靠近擋風墻,過流斷面減小,氣流被壓縮發生繞流,在迎風側墻底形成渦旋流,攜沙風速度減小,一部分沙粒在此處堆積。當氣流繞過擋風墻以后,過流斷面增大,氣流發生擴散,攜沙風會帶走一部分沙粒直接躍過線路上方被帶向遠處,另一部分則會落下,在擋風墻背風側堆積;受逆壓梯度的影響在路基頂面上方形成氣流紊流區,并且形成了一個明顯的渦流,此處沙粒運動比較復雜,在渦流作用下,部分沙粒會堆積在擋風墻背風側墻角位置處,另一部分沙粒在反向風的作用下從路堤背風側坡腳朝著線路運動并在此處重新進行二次堆積。

2.2高度對擋風墻周圍風沙流運動特性的影響

擋風墻在為列車的安全運營提供保障的同時,不可避免地造成了線路積沙,為了研究擋風墻背風側風沙流結構與擋風墻高度的關系,分別取擋風墻高度為3.0、3.5、4.0 m,初始風沙流瞬時速度20 m/s為例進行數值模擬研究。

風沙流是指含有沙粒的運動氣流,其形成必須要具備地表較為豐富的沙物質和一定的風力。風是沙粒運動的動力來源,研究背風側流場情況主要是研究氣流速率的變化情況。如圖3所示,不同高度擋風墻背風側一定高度內的風速垂直輪廓線呈現出大致相同的走勢規律,可分為3個階段。在0~0.5 m風速增大較快,呈指數分布規律,并且隨著擋風墻高度的增加,風速增加的幅度較大,當擋風墻高度為3 m時,在這個階段內最大風速約為7 m/s,當擋風墻高度為4 m時,最大風速約為9 m/s。在0.5 m以上至擋風墻自身高度區間內風速變化較為復雜,呈先減小后增大的趨勢,并且風速較低,基本保持在2~9 m/s,擋風墻高度越高,距床面同一高度處風速越大。在擋風墻高度0.5 m以上部分區間為加速區,隨著擋風墻高度的增加風速呈指數分布規律急劇上升恢復至初始風速。擋風墻高度越高,恢復初始速度的位置距離路基頂面的高度越高,說明渦旋區范圍越大,基本都在擋風墻本身高度以上1 m位置處恢復至初始速度。

圖3 不同高度條件下風速垂直輪廓線圖(背風側距擋風墻2 m處)

風沙流是一種風吹沙動的現象,近地表風速對沙粒的運動有著顯著的作用,所以研究近地表風速對擋風墻背風側的積沙形態有著重要意義。圖4為3種高度擋風墻距床面0.5 m處的風速變化曲線。取近地表的x方向速度進行研究,負值代表該區域內存在回流,氣流速度方向與初始方向相反。由圖4可以看出,風速呈現出先增大后減小再增大的趨勢。距擋風墻10 m以內,風速沿著與初始氣流速度相反的方向增大,并且隨著擋風墻高度的增加風速增加幅度越大,當擋風墻高度為3 m時,在這個階段內最大風速約為7.5 m/s,當擋風墻高度為4 m時,最大風速約為12 m/s。距擋風墻背風側10~20 m,風速從反向的最大值減小至0并且急劇恢復至初始速度,說明此處產生了一個明顯的渦旋。距擋風墻背風側約20 m時,由于背風側路基邊坡的擾動作用,此處產生了一個較小的渦流,在此之后距擋風墻約55 m處風速恢復至初始速度,并且隨著擋風墻高度的增大,恢復至初始速度的位置距離擋風墻越遠,說明擋風墻高度越高,擋風墻背風側形成的渦流范圍越大,渦流區長度越長。

圖4 不同高度條件下背風側距床面0.5 m處風速變化曲線

2.3擋風墻背風側不同風速下風沙流運動特性

風速是沙粒運動的動力條件,也是影響沙粒運動形式的主要因素之一。為了驗證背風側氣流隨著風速的變化情況,取擋風墻高度為3.5 m,分別取風速為10、20、30 m/s和40 m/s進行數值模擬計算。

研究風速變化情況對擋風墻背風側的積沙形態以及形成機理有著重要作用。故描繪距擋沙墻背風側2 m處的風速廓線圖以及近地表的風速變化曲線,從圖5中可以看出,在0~0.5 m,風速變化趨勢基本一致,呈指數分布逐漸增大,當初始風速為10 m/s時,在這個階段內,最大風速為5 m/s,當初始風速為40 m/s時,最大風速為15 m/s,風速越大增加幅度越大。在0.5~4 m的渦旋區內,當風沙流初始速度為10 m/s時,氣流速度變化較為平緩,基本保持在5 m/s左右。隨著風速的增加,渦旋區的氣流速度變化逐漸趨于復雜化,在距床面約2 m高度處,氣流速度減至最小值,當風沙流初始風速為10 m/s時,氣流速度從5 m/s減至4 m/s,當初始風速為40 m/s時,氣流速度從15 m/s減至5 m/s,風沙流初始速度越大,氣流衰減的幅度越小,隨著遠離床面,風速逐漸增加,初始風速越大,增加的幅度越大。距床面4 m以上,風速呈線性急劇增加,恢復至初始速度。由此可以說明,當風沙流速度越大時,單位時間內通過擋風墻的沙粒越多,隨著過流斷面減小,沉積在線路上的沙粒越多更容易形成積沙。

圖5 不同風速條件下風速垂直輪廓線(背風側距擋風墻2 m處)

圖6 不同風速條件下背風側距床面0.5 m處風速變化曲線

取近地表的x方向速度進行研究,負值代表該區域內存在回流,氣流速度方向與初始方向相反。從圖6可以看出,在距擋風墻背風側0~5 m,風速變化趨勢基本一致,呈反向增大的趨勢。在這個階段內,風速最大值增加的幅度保持在50%左右,當初始風速為10 m/s時,最大風速為5 m/s,當初始風速為40 m/s時,最大風速為23 m/s,風速越大增加幅度越大。在距擋風墻背風側5~20 m,風速沿著初始風速的方向減小為0,同時保持風向不變繼續增大恢復至初始速度大小。在初始風速為40 m/s時,風速極值差最大,風速從-23 m/s變為29 m/s,說明風速越大時,氣流的削弱效果最明顯。在距離擋風墻背風側約20 m處,風速逐漸恢復至初始速度,并且風速越大,恢復至初始速度的位置距離擋風墻越遠,說明渦流區長度也隨著風速的增加而增大,表明擋風墻背風側積沙區逐漸遠離擋風墻,直接影響背風側的積沙形態和積沙位置。

2.4擋風墻背風側積沙形態分析

風是沙粒運動的直接動力,風速的大小直接影響擋風墻背風側的積沙分布。路堤對風沙流的運動起到一定的阻礙作用,使一部分攜沙風速度減小,沙粒掉落從而導致線路積沙,而布設擋風墻以后使得這種阻礙效果加劇。

為了研究風速對擋風墻背風側線路積沙的影響,分別取風沙流初始速度為10、20、30、40 m/s(風向從左向右,以下沙粒體積分數云圖相同),擋風墻高度為3 m分別進行數值模擬。圖7為t= 5 s時不同風速條件下擋風墻背風側風沙流場沙粒體積分數云圖(不同的顏色代表沙粒各種不同的體積分數,紅色最大,藍色最小),紅色表示沙粒已經沉積在此處,其他顏色代表著沙粒還在以不同的形式繼續運動。從圖7可以看出,當風速較小時,線路上基本無積沙,沙粒多數堆積在擋風墻墻角的位置以及背風側路肩的位置,隨著風速的增加,單位時間內躍過擋風墻的沙粒越多,在擋風墻背風側渦旋作用下,堆積在擋風墻背風側墻角的積沙逐漸向線路上移動,如圖7(c)、(d)所示,沙粒基本上堆積在靠近擋風墻側的線路上,在軌枕以及軌道等構筑物的影響下逐漸在線路上形成積沙,可見風速越大,沙粒越容易被風吹起,隨著過流斷面減小,沙粒越容易沉積。并且隨著風速的增大,反向風的影響也越大,沉積在背風側路基邊坡處的積沙在反向風的作用下被帶到線路上造成二次堆積。圖8為蘭新第二雙線某段鐵路風沙侵襲后的線路積沙情況,由圖8可以看出,線路積沙的位置與圖7(c)、7(d)基本一致,可見,數值模擬與現場實際調查相一致。

圖7 t=5 s時路堤頂面風沙流場沙粒體積分數云圖

圖8 蘭新第二雙線某段線路積沙

圖9為風沙流初始速度10 m/s時,不同時刻的沙粒體積分數云圖。由圖9可以看出,在風速比較低時線路上的積沙比較少,且沙粒繼續運動,沙粒大多數堆積在擋風墻背風側墻角處,也有少部分堆積在背風側路肩處。圖10為風沙流初始速度30 m/s時,不同時刻的沙粒體積分數云圖。由圖10可以看出,風速比較大時,線路上會有積沙,并且隨著時間的推移,擋風墻背風側的沙粒朝著線路上運動,一部分被線路上的軌枕軌道等構筑物攔住造成線路積沙。對比圖9和圖10可知,隨著風沙流初始速度的增大,風攜帶沙的能力增強,更容易吹起沙粒,當攜沙風通過擋風墻時,由于過流斷面減小,氣流擴散,大部分沙粒沉積下來,所以更容易引起線路積沙。另一部分沙粒在攜沙風經過擋風墻時直接被風帶走,落在擋風墻背風側路堤坡腳處,在反向風的作用下又重新被吹起來落到線路上,造成二次堆積,風速越大反向風的影響越大。故在風速較大的地區,布設擋風墻時應注意道路積沙的防治。

圖9 風沙流速度v=10 m/s時不同時刻路堤頂面沙粒體積分數云圖

圖10 風沙流速度v=30 m/s時不同時刻路堤頂面沙粒體積分數云圖

3 結論

(1)擋風墻背風側風速廓線變化呈三個階段:在0~0.5 m之間呈指數增長趨勢,并隨著擋風墻高度的增加,風速增加幅度越大;在0.5 m至擋風墻自身高度區間內變化較為復雜,呈先減小后增加的趨勢;在擋風墻0.5 m以上部分為加速區,呈指數增長趨勢急劇增加。

(2)擋風墻背風側近地表氣流變化呈3個階段:距擋風墻10 m以內,風速沿著與初始氣流速度相反的方向增大,并且隨著擋風墻高度增加風速增加幅度越大;距擋風墻背風側10~20 m,風速從反向的最大值減小至0并且急劇恢復至初始速度; 20 m之后逐漸恢復至初始速度。

(3)擋風墻高度一定時,隨著風速的增加,擋風墻背風側渦旋區的氣流豎向速度變化逐漸趨于復雜化,距床面4 m以上,風速呈線性急劇增加;擋風墻背風側近地表風速最大值增加的幅度保持在50%左右,風速越大,氣流的削弱作用越明顯。

(4)當初始氣流速度為10、20 m/s時,線路上積沙較少,沙粒多數堆積在擋風墻背風側墻角處;隨著風速的增加,初始氣流速度為30、40 m/s時,單位時間內通過擋風墻的沙粒更多,由于過流斷面減小,氣流擴散,更多沙粒沉積在線路上。

(5)在強風地區,布設擋風墻時在滿足抗傾覆力矩的條件下,不建議加高擋風墻,同時應考察線路上風向的地表情況,沙源比較豐富時應采用擋沙堤、高立式、中立式擋沙墻相結合的工程治沙措施來減小風沙流密度,達到防沙的目的。

參考文獻:

[1]高永平,錢偉平.淺析蘭新鐵路防風工程[J].資源環境與工程,2009(9):48-51.

[2]葛盛昌,蔣富強.蘭新鐵路強風地區風沙成因及擋風墻防風效果分析[J].鐵道工程學報,2009(5):1-4.

[3]董漢雄.蘭新鐵路百里風區擋風墻設計[J].路基工程,2009 (2):95-96.

[4]龐巧東,程建軍,蔣富強,等.戈壁鐵路擋風墻背風側流場特征與擋風功效研究[J].鐵道標準設計,2011(2):1-5.

[5]葛春庚,石龍,李凱崇.蘭新二線強風地區防沙措施效益評價[J].鐵道標準設計,2015,59 (9):37-40.

[6]鄭曉靜,王萍.風沙流中沙粒隨機運動的數值模擬研究[J].中國沙漠,2006,26(2):184-188.

[7]武生智,任春勇.基于歐拉雙流體模型的風沙運動模擬[J].蘭州大學學報:自然科學版,2012,48(1):104-112.

[8]王康龍,武建軍,羅生虎.風沙運動的歐拉雙流體模型參數研究[J].中國沙漠,2014,34(6):1461-1468.

[9]董治寶,慕青松,王洪濤.風沙流中風速廓線的數值模擬與實驗驗證[J].氣象學報,2008,66(2):158-166.

[10]石龍,蔣富強,韓峰.風沙兩相流對鐵路路堤相應規律的數值模擬研究[J].鐵道學報,2014,36(5):82-87.

[11]張軍平,王引生,蔣富強.蘭新鐵路戈壁地區路基周圍風沙流運動特征數值分析[J].中國鐵道科學,2011,32(4):14-18.

[12]劉賢萬,崔志剛.特大風區防翻車擋風墻工程設計的風洞實驗研究[J].中國沙漠,1994,14(3):38-46.

[13]程建軍,蔣富強,楊印海,等.戈壁鐵路沿線風沙災害特征與擋風沙措施及功效研究[J].中國鐵道科學,2010,31(5):15-20.

[14]高廣軍,段麗麗.單線路堤上擋風墻高度研究[J].中南大學學報:自然科學版,2011,42(1):254-259.

[15]朱紅均,林元華,謝龍漢.流體分析及仿真實用教程[M].北京:人民郵電出版社,2010.

Numerical Simulation Analysis of Response Law of Wind-blown Sand Flow around Wind-break Wall in Gobi Area

LI Xiao-jun1,2,JIANG Fu-qiang2
(1.School of Civil Engineering,Lanzhou Jiaotong University,Lanzhou 730070,China; 2.Northwest Research Institute Co.,C.R.E.C.Lanzhou 730000,China)

Abstract:Numerical Simulation of the motion characteristics of wind-blown sand flow around the existing wind-break wall along Lanzhou-Urumqi Railway based on FLUENT eulerian two-fluid model is conducted to identify the distribution features of the flow field and accumulated sand at the leeward side of wind-break.The results show that,wind profile change law shows an exponential increasing tendency at the leeward side of wind-break wall,the changes at the height of 0.5 m to wind-break wall itself are complex and tend to decrease and then increase; the airflow velocity on near-surface at the leeward side of wind-break wall is increased in the reserve direction,reduced to 0 along the direction of the initial velocity,continued to increase to the initial velocity,and the maximum wind speed increase rate remains at about 50%,the higher the wind speed,the more obvious the airflow of weakening effect; when the initial velocity is small,less sand is accumulated on the line,and sand accumulation is concentrated at the corner of the leeward side of wind-break wall; with the increase of wind speed,more sand passes through the windbreak in unit time and even more sand is deposited along the line due to the reduction of flow cross-section and diffusion of air; in strong wind areas,the surface conditions affecting wind direction should be investigated before wind-break wall is located,engineering measures should be taken where sand source is rich to reduce sand flow density and prevent sand invasion.

Key words:Lanzhou-Uremqi Railway; Windbreak wall; Wind-blown-sand of two phase flow; Numerical simulation

作者簡介:李曉軍(1989—),男,碩士研究生,E-mail:472775921@ qq.com。

收稿日期:2015-08-13;修回日期:2015-08-19

文章編號:1004-2954(2016) 03-0047-05

中圖分類號:U216.41+3

文獻標識碼:A

DOI:10.13238/j.issn.1004-2954.2016.03.011

主站蜘蛛池模板: 秋霞国产在线| 国产99热| 亚洲中文字幕23页在线| 国内毛片视频| 国产色伊人| 一级毛片高清| 国产福利拍拍拍| 亚洲精品成人福利在线电影| 日韩av电影一区二区三区四区| 国产第一页第二页| 亚洲国产成人精品无码区性色| 国产91无码福利在线| 亚洲日韩精品无码专区97| 亚洲中文字幕无码爆乳| a亚洲视频| 蜜桃视频一区| 波多野结衣一级毛片| 欧美在线综合视频| 九色视频最新网址| 喷潮白浆直流在线播放| 自拍偷拍一区| 欧洲熟妇精品视频| 男人天堂亚洲天堂| 国产国产人免费视频成18| 国产亚洲一区二区三区在线| 亚洲综合久久成人AV| 国产第一页免费浮力影院| 久草青青在线视频| 久久精品人人做人人爽电影蜜月 | 91九色视频网| 九九热精品免费视频| 伊人AV天堂| 浮力影院国产第一页| 人妻精品久久久无码区色视| 免费在线色| 国模私拍一区二区三区| 亚洲午夜18| 亚洲第一极品精品无码| 亚洲v日韩v欧美在线观看| 51国产偷自视频区视频手机观看| 国产精品yjizz视频网一二区| 草逼视频国产| 国产精品福利一区二区久久| 国产在线精品99一区不卡| 成·人免费午夜无码视频在线观看| 粗大猛烈进出高潮视频无码| 国产18在线| 久久99国产乱子伦精品免| 五月婷婷激情四射| 五月婷婷中文字幕| 日韩美毛片| 日韩精品一区二区三区中文无码| 亚洲国产成人久久77| 国产精品偷伦视频免费观看国产 | 67194在线午夜亚洲| 日本高清在线看免费观看| 亚洲成人高清在线观看| 91精品综合| 久久精品aⅴ无码中文字幕| 色欲色欲久久综合网| 一区二区三区在线不卡免费| 国产亚洲精品无码专| 99尹人香蕉国产免费天天拍| 国产福利影院在线观看| av大片在线无码免费| 九九视频免费在线观看| 国产极品美女在线观看| 欧美精品1区| 国产浮力第一页永久地址| 亚洲综合片| 国产日韩欧美视频| 亚洲水蜜桃久久综合网站| 亚洲综合狠狠| 欧美高清三区| 亚洲精品天堂在线观看| 狠狠久久综合伊人不卡| 亚洲丝袜中文字幕| 日韩天堂视频| 欧美三级自拍| 亚洲第一视频免费在线| 男人天堂亚洲天堂| 国产成人亚洲综合a∨婷婷|