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

        ?

        帶灘槽地形的連續(xù)彎道中縱向流速橫向分布解析1)

        2021-03-10 09:46:18劉玉嬌余明輝田浩永
        力學(xué)學(xué)報(bào) 2021年2期
        關(guān)鍵詞:床面垂線水深

        劉玉嬌 余明輝 田浩永

        (武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,武漢 430072)

        引言

        天然河道的演變是水流攜帶泥沙并與河道邊界相互交換相互影響的辯證統(tǒng)一體.其中,縱向流速分布規(guī)律與河道過流能力[1]和泥沙輸移[2]的密切相關(guān),也與崩岸、切灘等突發(fā)現(xiàn)象相關(guān).對(duì)于河流流速的分布規(guī)律,研究主要從模型試驗(yàn)[3-5]、數(shù)值模擬[6-9]和解析解計(jì)算[10-13]進(jìn)行.解析解由于效率高,方法簡單,在實(shí)際應(yīng)用方面被廣泛采用[14].在順直復(fù)式河道中,Shiono 和Knight[15]提出了SKM 方法,認(rèn)為二次流導(dǎo)致的切應(yīng)力沿橫斷面為線性分布,對(duì)動(dòng)量方程沿水深積分,推導(dǎo)出垂線平均縱向流速(以下簡稱為Ud,下標(biāo)d 表示為垂線平均量)的解析解;Devi 等[16]進(jìn)一步根據(jù)實(shí)測數(shù)據(jù)量化二次流中的參數(shù),發(fā)現(xiàn)參數(shù)與定水深和變化的水深之間有不同的關(guān)系;復(fù)式河槽漫灘區(qū)域常為順直的,但主河槽則傾向于彎曲的,存在彎道環(huán)流的影響,Ervine 等[17]考慮了復(fù)式河道彎曲主槽中水流的強(qiáng)三維性,認(rèn)為二次流項(xiàng)與Ud的二次方成正比關(guān)系,對(duì)SKM 方法進(jìn)行了改進(jìn);Knight等[18]分析了環(huán)流結(jié)構(gòu)與二次流項(xiàng)取值正負(fù)之間的關(guān)系,認(rèn)為環(huán)流中心左側(cè)和右側(cè)二次流項(xiàng)符號(hào)相反;Liu 等[19]和Huai 等[20]則進(jìn)一步在動(dòng)量方程中增加了灘地植被影響,在解析解公式中考慮了植被對(duì)流速的影響;Zhong 等[21]考慮了在兩岸及水面結(jié)冰條件下,解析解的形式,并分析了床面阻力、冰蓋阻力和二次流等因素的影響;楊中華和高偉[22]比較了三種不同的積分假定,將計(jì)算值與實(shí)測值比較,得出最為吻合的計(jì)算公式;許唯臨[11]將二次流項(xiàng)歸并入雷諾切應(yīng)力項(xiàng),對(duì)復(fù)式順直河道和復(fù)式游蕩型河道中漫灘流速進(jìn)行了研究,復(fù)式游蕩型河道流速計(jì)算其采用了分層計(jì)算的方法.哈岸英等[23]利用SKM 方法研究了沖積河流漫灘水流水力因子隨水深的變化特性,提出了定量劃分復(fù)式斷面灘與槽的判別指標(biāo).Rezaei 和Knight[24]探究了SKM 方法在非對(duì)稱復(fù)合水槽中的應(yīng)用.Sharifi 和Sterling[2]利用SKM 方法計(jì)算床面切應(yīng)力的橫向分布,并將其與泥沙輸移公式進(jìn)行的對(duì)比,以獲得對(duì)泥沙輸移的預(yù)測.

        彎曲型河道是常見的天然河道形態(tài),水流經(jīng)過彎道時(shí)受到離心力的作用,產(chǎn)生橫向壓力梯度,形成彎道橫向環(huán)流,橫向環(huán)流與縱向流速結(jié)合形成彎道螺旋流[25].縱向水流在彎道進(jìn)口處,主流靠近凸岸側(cè),沿著彎道,主流逐漸向彎道凹岸側(cè)移動(dòng),在出口處,主流貼近凹岸.Tang 和Knight[26]將SKM 方法應(yīng)用于彎道中,基于彎道水流充分發(fā)展的條件下,認(rèn)為二次流橫向變化率與橫向位置為線性關(guān)系,推導(dǎo)出彎道Ud的解析式.天然彎曲河流常為多個(gè)彎道相連,床面形態(tài)存在著深槽和邊灘分布[27],本文嘗試以SKM 方法為基礎(chǔ),推導(dǎo)出彎道中Ud橫向分布計(jì)算公式,運(yùn)用實(shí)測資料率定計(jì)算參數(shù),將其應(yīng)用于帶灘槽地形的反向連續(xù)彎道中Ud橫向分布解析計(jì)算研究,為下一步的彎道縱向流速快速預(yù)測奠定基礎(chǔ).

        1 帶灘槽地形彎道垂線平均縱向流速解析解推導(dǎo)

        1.1 理論基礎(chǔ)

        結(jié)合水流連續(xù)方程,曲線坐標(biāo)系下,單位體積上的水流縱向動(dòng)量方程可簡化為

        式中,s,r,z分別表示水流的切向、徑向和垂向三個(gè)方向,U,V,W和u′,v′,w′分別為s,r,z方向的時(shí)間平均流速分量和脈動(dòng)流速分量,上劃線表示時(shí)間平均值,p為壓強(qiáng),fs為s方向的單位體積力,ρ 為水流密度,g為重力加速度

        式中,S0為床面縱比降,S0=??zb/?s,zb為床面高程,H為水深,τrs和τzs分別表示垂直于r和z平面上的切向雷諾應(yīng)力.

        將式(2)~式(4)代入式(1)中,動(dòng)量方程可表達(dá)為

        對(duì)式(5)沿垂線進(jìn)行積分可得

        式中,τb為床面切應(yīng)力;1/Sr為床面橫比降;腳標(biāo)d表示垂線平均值;Se為總水頭能量的比降,用于表示??(zb+H+U2/g)/?s[24];Ua為斷面平均流速;n為糙率,若實(shí)驗(yàn)水槽為水泥面,查表得糙率值為0.014[28];R為水力半徑.

        1.2 解析解推導(dǎo)

        Tang 和Knight[26]根據(jù)van Balen 等[29]180?彎道試驗(yàn)數(shù)據(jù),分析了二次流項(xiàng)和二次流附加項(xiàng)的橫斷面分布規(guī)律,發(fā)現(xiàn)其沿橫斷面為線性分布,本文對(duì)式(6)右邊二次流項(xiàng)和彎道附加應(yīng)力項(xiàng)也采用線性分布假設(shè),即

        式中,m和k分別為線性分布規(guī)律常數(shù)項(xiàng)和一次項(xiàng)系數(shù),根據(jù)等式左邊可知,m的單位為kg/(m·s2),k的單位為kg/(m2·s2),隨彎道平面形態(tài)及研究斷面所處位置的不同而變化,可采用實(shí)測資料率定獲得.

        式(6)可寫為

        式中,λ 為無量綱渦黏性系數(shù),其取值在河道內(nèi)變化不大,取值為0.07[10,19];f為阻力系數(shù).

        將式(13)和式(14)代入式(12)中,可得

        本文將帶灘槽地形的橫斷面形態(tài)概化成如圖1所示(圖中灰色區(qū)域?yàn)闉┑?,河槽區(qū)床面高程不變,河灘區(qū)床面線性抬高.對(duì)不同的床面形態(tài),由式(16)可求得Ud橫向分布的解析解,解析解的表達(dá)式如下.

        (a)對(duì)河槽區(qū),即水深恒定的區(qū)域,Ud的解析解為

        圖1 帶灘槽地形斷面形態(tài)概化圖Fig.1 A sketch of a cross-section with a pool-point bar

        (b)對(duì)河灘區(qū),即床面線性變化的區(qū)域(床面橫比降為1/Sr),Ud的解析解為

        式中,A1~A4為未知的系數(shù),ξ 表示床面上的水深,ξ=H?r/Sr.

        1.3 邊界條件

        對(duì)于Ud的解析式,存在兩種邊界條件,一是無滑移邊界,二是交界面處要滿足連續(xù)性,即:

        (1)兩岸(左岸和右岸)邊界處流速為0;

        (2)灘槽交界處,主槽和灘地流速相等Ud,i=Ud,i+1(i和i+1 表示參數(shù)取值不同的兩個(gè)相鄰區(qū)域);

        (3)灘槽交界處,主槽和灘地流速梯度相等,即?Ud,i/?r=?Ud,i+1/?r.

        根據(jù)上述邊界條件,解析式中系數(shù)A1,A2,A3,A4均可求解.

        2 解析解模型的應(yīng)用

        2.1 帶灘槽地形的連續(xù)彎道縱向流速斷面橫向分布試驗(yàn)介紹

        試驗(yàn)在武漢大學(xué)泥沙重點(diǎn)實(shí)驗(yàn)室的反向連續(xù)彎道水槽中進(jìn)行,水槽進(jìn)出口均有一段長7 m 的順直段,上下游兩彎道曲率保持恒定,上游彎道的中心角為90?,下游彎道的中心角為135?,水槽寬度為1.2 m,彎道中心線曲率半徑為2.4 m,兩彎道間的順直連接段長2 m,水槽橫斷面形狀如圖2(b)所示,過渡段中間斷面和進(jìn)出口5 m 內(nèi)順直段均為矩形斷面,底部高程為0.1 m,彎道橫斷面形狀為凸岸高凹岸低,由凸岸到凹岸,床面高程以1:8 的橫比降由0.2 m 減小為0.1 m(水深減小0.1 m),接著以1:2 的橫比降床面高程減小到0 m (水深減小0.1 m),之后床面高程和水深保持不變,順直段內(nèi)矩形斷面與彎道內(nèi)灘槽地形斷面通過漸變的方式進(jìn)行連接.

        試驗(yàn)采用三維超聲多普勒流速儀(ADV)進(jìn)行三維流速測量,儀器采樣頻率為200 Hz,斷面的測點(diǎn)布置橫向間距為5 cm,垂向間距為2 cm,靠近床面區(qū)域加密測量,測點(diǎn)垂向間距減小為1 cm,由于測量探頭的大小及儀器采樣體積為探頭下方5 cm,因此在邊壁和水面5 cm 以內(nèi)無法進(jìn)行測量,縱向流速的垂線平均值根據(jù)各測線的流速數(shù)據(jù)取垂線平均獲得.

        圖2 連續(xù)彎道試驗(yàn)布置Fig.2 Layout of the consecutive bends experiment

        三組工況的水力學(xué)條件如表1 所示,典型天然彎曲河流的Fr范圍為0.07~0.26[30],三種試驗(yàn)工況的Fr取值在此范圍內(nèi),并保持流量不變,研究水深變化的影響,采用試驗(yàn)測量的五個(gè)斷面(斷面A~E)進(jìn)行驗(yàn)證.表中,Q為進(jìn)口流量,H為出口斷面水深,U為進(jìn)口斷面平均流速,Fr為佛汝德數(shù),Re為雷諾數(shù).

        表1 試驗(yàn)條件Table 1 Experimental conditions

        2.2 結(jié)果與分析

        圖3 為解析解的擬合曲線與試驗(yàn)結(jié)果的比較(圖中灰色區(qū)域?yàn)闉┑?,圖中圓點(diǎn)為根據(jù)實(shí)測的三維流速數(shù)據(jù)計(jì)算出的Ud,由圖可知實(shí)測的Ud沿橫斷面的分布規(guī)律為深槽流速較小,流速沿橫向變化率較小,灘地上流速較大,流速沿橫向變化率增大,Ud的最大值在靠近凸岸側(cè),下游彎道中Ud橫向分布較上游彎道更不均勻,由工況一到工況三,出口水深逐漸減小即水流流速增大,在下游彎道中Ud的最大值逐漸向凹岸偏移.

        圖3 流速計(jì)算值與實(shí)測值對(duì)比Fig.3 Comparison of modelled depth-averaged velocity distribution with experimental data

        圖3 中曲線為計(jì)算出的擬合曲線,水深恒定不變的河槽區(qū)(距凹岸0~0.2 m 范圍內(nèi))采用式(17)進(jìn)行擬合,床面線性變化的邊坡區(qū)(距凹岸0.2~1.2 m 范圍內(nèi))由式(20)擬合,參數(shù)m的取值根據(jù)地形分為三個(gè)區(qū)域,分別為水平段m1、橫比降為1/2 的斜坡段m2和橫比降為1/8 的斜坡段m3,參數(shù)k的取值分為兩個(gè)區(qū)域(靠近凹岸側(cè)為k1,靠近凸岸側(cè)為k2),分區(qū)位置在主流附近,解析解表達(dá)式中的系數(shù)A1~A4可根據(jù)邊界條件求解.由圖可知擬合曲線與實(shí)測數(shù)據(jù)吻合良好,在床面地形變化處,Ud沿程變化率較大時(shí)以及Ud峰值的大小與位置,計(jì)算值對(duì)實(shí)測值均能較好擬合,隨著出口水深減小,Ud最大值的位置變化,計(jì)算結(jié)果也能較好地反映這一規(guī)律,僅在靠近邊壁處擬合曲線與實(shí)測數(shù)據(jù)稍有偏差,凹岸側(cè)偏差可能是凹岸次生反向小環(huán)流影響導(dǎo)致,凸岸邊壁附近偏差則可能是由于水深較淺,水流結(jié)構(gòu)較為復(fù)雜導(dǎo)致,各斷面的均方根誤差均小于0.015 m/s,平均相對(duì)誤差均在5%以內(nèi),實(shí)測值與計(jì)算值之間誤差較小.綜上所述,在帶有灘槽的反向連續(xù)彎道中,此解析式能較好擬合不同水流條件下沿程各斷面垂線平均縱向流速的橫向分布規(guī)律.

        3 計(jì)算參數(shù)敏感性分析及取值

        對(duì)參數(shù)m和k的取值進(jìn)行敏感性分析,此處以工況二的斷面A擬合結(jié)果為基礎(chǔ)進(jìn)行分析.圖4對(duì)k的數(shù)值大小進(jìn)行了敏感性分析,在0~1.1 m 和1.1~1.2 m 取值不同(分別表示為k1和k2),將k1和k2的取值分別在原基礎(chǔ)上增大和減小0.05 (每次僅變化一個(gè)參數(shù)的取值),由圖可知k1數(shù)值變化對(duì)凹岸深槽區(qū)域流速計(jì)算值基本無影響,主要變化發(fā)生在靠近凸岸的邊灘附近,且計(jì)算值變化較小,k2的變化僅影響1~1.2 m 內(nèi)的局部區(qū)域,流速計(jì)算值變化較小,流速最大值的大小和位置未發(fā)生明顯變化;圖5 對(duì)k1和k2分區(qū)位置的敏感性進(jìn)行了分析,虛線Lk表示分界位置,將初始分界位置分別向左和向右移動(dòng)0.1 m,由圖可以發(fā)現(xiàn),分界位置的變化對(duì)流速的峰值影響較大,流速最大值變化明顯,且位置也發(fā)生變化,由此可將主流線的位置作為k分區(qū)的初步依據(jù).

        圖4 k 變化時(shí)垂線平均流速計(jì)算值Fig.4 Modelled depth-averaged velocity for various k values

        圖5 k 分區(qū)調(diào)整時(shí)垂線平均流速計(jì)算值Fig.5 Modelled depth-averaged velocity for various k values in different panels

        參數(shù)m分為三個(gè)區(qū)域進(jìn)行取值,為m1(0~0.2 m),m2(0.2~0.4 m)和m3(0.4~1.2 m),圖6 為對(duì)參數(shù)m取值的敏感性分析結(jié)果,分別對(duì)m1,m2和m3的初始值增大和減小0.05 (每次僅變化一個(gè)參數(shù)的取值),由圖可知,m1的取值增大或減小對(duì)流速計(jì)算值基本無影響,m2的取值變化僅在0.2~0.4 m 附近的流速計(jì)算值產(chǎn)生影響,流速計(jì)算值變化較小,m3取值變化對(duì)整個(gè)斷面的流速計(jì)算值基本均有影響,在0.4~1.2 m 區(qū)域內(nèi),流速計(jì)算值發(fā)生明顯變化;圖7為m分區(qū)位置變化的敏感性分析結(jié)果,圖中虛線Lma和Lmb分別為m1與m2和m2與m3分區(qū)的位置,將初始分界位置分別向左和向右移動(dòng)0.1 m(每次僅變化一個(gè)分界位置),由圖可知m1和m2分界位置向斜坡段移動(dòng)對(duì)流速的計(jì)算值影響較為明顯,m2和m3分界位置的變化對(duì)流速計(jì)算值影響較小.

        圖6 m 變化時(shí)垂線平均流速計(jì)算值Fig.6 Modelled depth-averaged velocity for various m values in different panels

        圖7 m 分區(qū)位置變化時(shí)垂線平均流速計(jì)算值Fig.7 Modelled depth-averaged velocity for various panels of m

        參數(shù)m和k的取值如表2 所示,Lk1和Lk2分別為k1和k2取值區(qū)域的長度,由表中數(shù)據(jù)可知,m和k的值與Tang 和Knight[26]文中的取值是相近的.沿?cái)嗝嫠蕉蝝1的取值小于斜坡段m2和m3的取值,可能與凹岸主槽側(cè)流速較小有關(guān),參數(shù)m的取值整體上隨水深的減小逐漸減小,水深減小到一定程度,河灘區(qū)m2和m3的取值由正變?yōu)樨?fù);k1和k2取值分區(qū)的界限與主流位置較為接近,分區(qū)的位置變化較小,可能與主流位置變化較小相關(guān),k1和k2取值的相對(duì)大小整體上由工況一的k1大于k2變?yōu)楣r三的k1小于k2.根據(jù)參數(shù)的取值敏感性分析可知m1,m2和k1,k2的取值變化對(duì)計(jì)算結(jié)果的影響較小,計(jì)算結(jié)果對(duì)m3的值較為敏感,但m3的取值沿程變化較小,因此給出各參數(shù)沿水槽的均值作為初步取值參考.

        4 討論

        本文對(duì)二次流項(xiàng)和二次流附加項(xiàng)采用與Tang 和Knight[26]相同的線性假設(shè),其僅分析彎道出口附近135?斷面的數(shù)據(jù),本文將解析解用于彎道中的多個(gè)斷面,此處增加繪制van Balen 等[29]的180?彎道進(jìn)口附近29?斷面的二次流項(xiàng)ρH(UV)d和二次流附加項(xiàng)的斷面分布(如圖8 和圖9所示),由圖可知其在彎道進(jìn)口斷面也符合線性假設(shè).進(jìn)一步分析二次流項(xiàng)和二次流附加項(xiàng)在本文試驗(yàn)水槽中的斷面分布規(guī)律,選取工況二csC斷面實(shí)測數(shù)據(jù)進(jìn)行討論,由圖10 和圖11 可知二次流項(xiàng)和二次流附加項(xiàng)的斷面分布符合線性分布規(guī)律,由此可知線性假設(shè)在本試驗(yàn)中也是符合的.由此可知,對(duì)二次流項(xiàng)和二次流附加項(xiàng)沿橫斷面的線性分布假設(shè)在彎道沿程多個(gè)斷面均適用,并在本文連續(xù)彎道中也是符合的,表明本文采取的線性假設(shè)在本文的研究問題中是合理的.

        表2 參數(shù)取值Table 2 The values of parameters

        圖8 二次流項(xiàng)ρH(UV)d29?橫斷面分布[29]Fig.8 Lateral distribution of secondary flow term,ρH(UV)dat cross-section of 29?[29]

        圖9 二次流附加項(xiàng)dz 29?橫斷分布[29]Fig.9 Lateral distribution of additional secondary flow term,dz at cross-section of 29?[29]

        圖10 二次流項(xiàng)ρH(UV)d橫斷面分布(csC)Fig.10 Lateral distribution of secondary flow term,ρH(UV)d(csC)

        本文推導(dǎo)出彎道中的Ud沿橫斷面解析解表達(dá)式,與帶灘槽地形的反向連續(xù)彎道實(shí)測數(shù)據(jù)符合良好,根據(jù)實(shí)測數(shù)據(jù)率定的參數(shù)m和k的取值與Tang和Knight[26]的取值是相近的,表明參數(shù)的取值對(duì)與不同的彎道存在相近的范圍,并分析了參數(shù)取值隨水位變化的規(guī)律,但具體的參數(shù)取值與水流條件和彎道形態(tài)之間的關(guān)系仍需更加深入的研究,以更好地推廣應(yīng)用解析解公式.

        圖11 二次流附加項(xiàng)橫斷分布(csC)Fig.11 Lateral distribution of additional secondary flow term,

        5 結(jié)論

        本文基于SKM 方法,對(duì)動(dòng)量方程進(jìn)行垂線積分,推導(dǎo)出了彎道中Ud的橫向分布解析解公式,將其應(yīng)用于帶灘槽地形的反向連續(xù)彎道中.

        (1)將解析解公式應(yīng)用于有灘槽的反向連續(xù)彎道,調(diào)整線性假設(shè)中參數(shù)的值,可以使計(jì)算值與實(shí)測值在不同的出口水流條件下,對(duì)沿程各斷面流速值橫向分布規(guī)律都能吻合良好,表明本文的線性假設(shè)具有一定的合理性.

        (2)參數(shù)取值,一次項(xiàng)系數(shù)k的大小對(duì)流速橫向分布趨勢有較大影響,k的分區(qū)對(duì)主流有較大的影響,調(diào)整分區(qū)可使計(jì)算值更好擬合灘地區(qū)流速的峰值大小和位置.常數(shù)項(xiàng)m整體上隨水深的減小而減小,可對(duì)m按照斜率變化分區(qū)對(duì)局部流速計(jì)算值進(jìn)行調(diào)整,河槽區(qū)m的取值小于河灘區(qū),流速計(jì)算值對(duì)m在河槽區(qū)與河灘區(qū)間分界線的變化較為敏感,河灘區(qū)內(nèi)m的分區(qū)則對(duì)流速計(jì)算值影響較小.

        (3)本文推導(dǎo)的解析解公式應(yīng)用范圍較廣,可用于地形變化的反向連續(xù)彎道中,但對(duì)于該模式的應(yīng)用需要進(jìn)行更加詳細(xì)的研究:①增加試驗(yàn)工況,進(jìn)一步分析參數(shù)m和k取值與進(jìn)出口水流條件(流量、水位和斷面平均流速等)和彎道形態(tài)(曲率、比降和過渡段長度等)之間的規(guī)律,總結(jié)參數(shù)分區(qū)依據(jù)和合理的取值范圍;②分析m和k的物理意義,利用更加詳細(xì)精確的數(shù)據(jù)計(jì)算其表示的各項(xiàng)具體分布規(guī)律,以更深入理解m和k取值變化規(guī)律的具體原因.

        猜你喜歡
        床面垂線水深
        魚鱗狀床面粗糙特性
        書法靜水深流
        河北水利(2022年10期)2022-12-29 11:48:12
        基于水深分段選擇因子的多光譜影像反演水深
        多角度思維實(shí)現(xiàn)平面與立體的轉(zhuǎn)化——學(xué)習(xí)微專題《明修棧道(作垂線)、暗度陳倉(找垂足)》有感
        畫垂線的方法
        對(duì)瓦里安碳纖維治療床面模型的評(píng)估
        近岸懸沙垂線分布多元線性回歸分析
        淹沒植物明渠床面沖淤及其對(duì)水流運(yùn)動(dòng)的影響
        改進(jìn)的投影覆蓋方法對(duì)遼河河道粗糙床面分維量化研究
        GPS RTK技術(shù)在水深測量中的應(yīng)用
        东京道一本热码加勒比小泽| 亚洲旡码a∨一区二区三区| 99视频一区| 亚洲国产一区二区三区视频在线| 亚洲女av中文字幕一区二区| 国产国产人免费人成免费视频| 久久综合精品国产丝袜长腿| 亚洲VR永久无码一区| 美女被黑人巨大入侵的的视频| 免费无码又爽又高潮视频| 吸咬奶头狂揉60分钟视频| 综合久久久久6亚洲综合| 日本一区二区三区精品免费| 国产亚洲美女精品久久久2020| 拍摄av现场失控高潮数次| 国产AV高清精品久久| 亚洲一区二区三区重口另类| 亚洲人成绝费网站色www| 国产精品久久婷婷免费观看| 在线观看的网站| 婷婷亚洲综合五月天小说| 少妇被粗大猛进进出出| 免费av日韩一区二区| 天天躁日日躁狠狠久久| 成人午夜免费无码视频在线观看 | 无码丰满少妇2在线观看| 资源在线观看视频一区二区| 久久久中文字幕日韩精品| 四虎影视免费永久在线观看| 亚洲熟妇乱子伦在线| 国产理论亚洲天堂av| 亚洲一区二区三区99| 最新精品国偷自产在线| 精品无码av不卡一区二区三区| 国产熟妇一区二区三区网站| 少妇被又大又粗又爽毛片久久黑人| 欧美老熟妇欲乱高清视频| 国产激情久久久久久熟女老人| 久久日韩精品一区二区| 亚洲精品无码国产| 337p日本欧洲亚洲大胆色噜噜|