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

        ?

        垂直軸風(fēng)力機(jī)三維尾流模型的研究

        2023-06-13 00:00:00凌子焱趙振宙劉一格劉惠馬遠(yuǎn)卓王丁丁
        太陽能學(xué)報 2023年11期
        關(guān)鍵詞:高斯分布

        收稿日期:2022-06-12

        基金項目:國家自然科學(xué)基金(51876054;52106239;11502070);江蘇風(fēng)力發(fā)電工程技術(shù)中心開放基金(ZK22-03-01)

        通信作者:趙振宙(1982—),男,博士、教授,主要從事風(fēng)力機(jī)空氣動力學(xué)、結(jié)構(gòu)力學(xué)和尾流控制方面的研究。joephy@163.com

        DOI:10.19912/j.0254-0096.tynxb.2022-0855 文章編號:0254-0096(2023)11-0196-07

        摘 要:基于流量守恒,提出一種H型VAWT的三維解析尾流模型,該模型考慮了風(fēng)切變效應(yīng),并采用多元高斯分布。首先,通過風(fēng)洞試驗及大渦模擬數(shù)據(jù)進(jìn)行模型驗證,結(jié)果表明在[x/Dgt;3]的遠(yuǎn)尾流區(qū),橫向剖面的相對誤差小于2%,垂向剖面的相對誤差小于3%。然后,從下游4個位置([x/D=3]、6、9、12)、3種推力系數(shù)([CT=0.6]、0.7、0.8)、4種風(fēng)切變指數(shù)([α=0]、0.1、0.15、0.2)、兩種湍流強(qiáng)度([I0=5%]、8.3%)下演示了一系列預(yù)測結(jié)果,結(jié)果顯示該文模型能有效地描述尾流風(fēng)速的空間分布。由于考慮了高度影響,該模型可用于風(fēng)力機(jī)輪轂高度的優(yōu)化及風(fēng)電場的布局優(yōu)化,有利于提高風(fēng)電場的功率輸出。

        關(guān)鍵詞:垂直軸風(fēng)力機(jī);解析模型;高斯分布;三維尾流模型

        中圖分類號:TK81""""""""""" """""""""""" """""文獻(xiàn)標(biāo)志碼:A

        0 引 言

        根據(jù)風(fēng)輪旋轉(zhuǎn)軸的方向,風(fēng)力機(jī)分為水平軸風(fēng)力機(jī)(horizontal axis wind turbine,HAWT)和垂直軸風(fēng)力機(jī)(vertical axis wind turbine,VAWT)[1-2]。為了研究VAWT的尾流空氣動力學(xué)特性,研究人員開展了大量工作,包括風(fēng)洞試驗、風(fēng)場實測和數(shù)值模擬。Battisti等[3]通過風(fēng)洞試驗測量VAWT的尾流速度分布,結(jié)果表明在不同的葉尖速比和推力系數(shù)下,尾流速度剖面形狀相似。文獻(xiàn)[4]實測了全尺寸直葉片VAWT的功率性能和尾跡特性,研究結(jié)果表明,隨著風(fēng)速增加,尾流速度恢復(fù)得更慢。Shamsoddin等[5]采用大渦模擬(large eddy simulation,LES)和致動線模型(actuator line model,ALM)相結(jié)合的方法研究VAWT的尾流流場。結(jié)果表明,LES-ALM方法能精確預(yù)測尾流中的平均值和高階統(tǒng)計量。Abkar等[6]使用LES研究H型VAWT的尾流結(jié)構(gòu),發(fā)現(xiàn)尾流速度虧損具有自相似的特性。

        研究VAWT的尾流效應(yīng)最重要的目的之一是預(yù)測尾流損失。開發(fā)高效率的工程預(yù)測解析尾流模型非常有必要。文獻(xiàn)[7]基于風(fēng)洞試驗數(shù)據(jù)提出一種一維尾流模型,用于預(yù)測下游任意截面的尾流邊界和平均流向速度。實測[4]和LES數(shù)據(jù)[6]表明,尾流速度在徑向呈高斯分布,上述模型假設(shè)尾流區(qū)內(nèi)速度恒定,與實際尾流分布存在較大差異。Abkar[8]基于動量守恒和質(zhì)量守恒,建立H型VAWT的二維尾流高斯模型,通過LES和實測數(shù)據(jù)驗證了該模型的精度。二維模型在實際應(yīng)用中具有局限性。文獻(xiàn)[9]研究表明,在風(fēng)電場中使用不同輪轂高度的風(fēng)力機(jī)能提高功率輸出和降低成本,二維模型忽略了高度因素,難以應(yīng)用到此類問題。此外,隨著大容量VAWT的開發(fā),其直徑和輪轂高度可達(dá)上百米,風(fēng)輪旋轉(zhuǎn)域內(nèi)高度上的風(fēng)速差較大。二維尾流模型假設(shè)風(fēng)速在高度上均勻分布,預(yù)測輸出功率時會產(chǎn)生較大誤差,需進(jìn)一步建立三維尾流模型。

        本文考慮風(fēng)切變的影響,建立一種H型VAWT的三維尾流模型。該模型基于流量守恒和多元高斯分布提出,以精準(zhǔn)預(yù)測三維尾流速度分布。采用2個風(fēng)洞試驗和LES數(shù)據(jù)來驗證該模型的精度,并驗證該模型在不同工況、多個視圖和不同位置處對尾流風(fēng)速的預(yù)測精度。

        1 模型建立

        圖1為二維尾流模型的示意圖。二維尾流模型忽略了風(fēng)切變的影響,假設(shè)入流風(fēng)速均勻分布,尾流風(fēng)速[U(x,z)]在垂向呈標(biāo)準(zhǔn)高斯分布。

        圖2為所建立的三維尾流模型的示意圖。在該模型中,入流風(fēng)速[U0(z)]隨高度變化,尾流區(qū)風(fēng)速分布更為復(fù)雜。尾流風(fēng)速用[U(x,y,z)]表示,[x]為來流方向,[y]表示垂直來流的水平方向,[z]表示垂直方向,如圖2所示。

        切變?nèi)肓黠L(fēng)速[U0(z)]用指數(shù)率表示:

        [U0(z)=uhzzhα]""" (1)

        式中:[uh]——VAWT輪轂高度處的入流風(fēng)速;[zh]——VAWT輪轂高度,即葉片展向中間截面離地面的高度;[α]——風(fēng)切變指數(shù)。

        前人研究[4,6]發(fā)現(xiàn),H型VAWT遠(yuǎn)尾流區(qū)中的速度虧損服從正態(tài)分布的自相似特性。為此,假設(shè)下游[x]處風(fēng)速虧損滿足各項異性高斯分布,尾流速度[U(x,y,z)]表示為:

        [U(x,y,z)=A(x)e-y22σ2y+(z-zh)22σ2z+B(x)+U0(z)]"""""" (2)

        式中:[A(x)]、[B(x)]——決定高斯形狀的關(guān)鍵參數(shù);[σy]、[σz]——各維高斯標(biāo)準(zhǔn)偏差,定義為:

        [σy=k*yx+εy,σz=k*zx+εz]""" (3)

        式中:[k*y]、[k*z]——各維高斯標(biāo)準(zhǔn)偏差增長率;[εy]、[εz]——[x]=0時[σy]和[σz]的值。

        然而,這種表示在大規(guī)模工程應(yīng)用中會帶來巨大的挑戰(zhàn),原因是[σy]和[σz]的概念對工程師來說過于抽象[10]。為了解決這個問題,將[σy]和[σz]替換為更易獲得且物理上直觀的尾流邊界,其表示為:

        [σy=ryC,σz=rzCry=12D+kwx,rz=12H+kwx]"" (4)

        式中:[ry]、[rz]——[y]和[z]方向上的尾流邊界半徑;[C]——常數(shù),根據(jù)實際情況確定;[D]、[H]——VAWT的風(fēng)輪直徑和葉片長度;[kw]——尾流邊界膨脹率,采用式(5)計算:

        [kw=k0IwakeI0]""" (5)

        式中:[k0]——陸上風(fēng)力機(jī)[k0]建議值為0.075,海上風(fēng)力機(jī)[k0]建議值為0.05[11];[I0]——入流湍流強(qiáng)度;[Iwake]——尾流區(qū)的湍流強(qiáng)度,由式(6)計算[12]:

        [Iwake=0.4CT(x/D)0.5+I0.502]""" (6)

        式中:[CT]——VAWT推力系數(shù)。

        為了求解式(2)中A(x)和[B(x),]根據(jù)尾流風(fēng)速的連續(xù)性,風(fēng)速在尾流邊界處恢復(fù)到來流風(fēng)速:

        [A(x)e-y22σ2y+(z-zh)22σ2z+B(x)+U0(z)=U0(z)," y2r2y+(z-zh)2r2z=1]"""" (7)

        聯(lián)立式(4)、式(7),可得到:

        [B(x)=-A(x)e-C22] (8)

        根據(jù)動量守恒定律,任意下游界面的尾流區(qū)內(nèi)流量守恒,總流量[Q(x)]可表示為初始矩形區(qū)域[SA0]內(nèi)的流量和超出矩形區(qū)域在尾流區(qū)域[SAw]內(nèi)的流量之和:

        [Q(x)=v1A0+SAw-SA0U0(z)ds]"" (9)

        式中:[v1]——緊靠風(fēng)輪處的下游風(fēng)速,根據(jù)致動盤理論,[v1=uh1-CT];[A0]——VAWT在順風(fēng)方向的投影面積,[A0=DH];[Aw]——橢圓尾流區(qū)域的面積,[Aw=πryrz]。

        [Q(x)]也可表示為實際風(fēng)速[U(x,y,z)]在尾流區(qū)域[SAw]上的積分:

        [Q(x)=SAwU(x,y,z)ds]" (10)

        聯(lián)立式(2)、式(3)、式(9)、式(10)求解得到[A(x)]:

        [A(x)=Q(x)-SAwU0(z)dsπryrz-2C2e-C22+2C2-e-C22]"""" (11)

        式中:[Q(x)]可聯(lián)立式(1)、式(4)~式(6)、式(9)求得;[C]——一個經(jīng)驗參數(shù),受葉片數(shù)目、葉片形狀和入流風(fēng)況等多種因素影響。HAWT尾流模型研究表明[13]:[C]需根據(jù)不同的條件設(shè)為不同的值。因此,在本研究中,根據(jù)實際情況將[C]設(shè)置為不同的值。

        將式(1)、式(8)、式(11)代入式(2)即可得到VAWT三維尾流模型表達(dá)式。

        2 模型驗證

        2.1 風(fēng)洞試驗驗證

        Battisti等[3]在意大利米蘭理工大學(xué)的風(fēng)洞實驗室中進(jìn)行了一系列風(fēng)洞試驗。試驗風(fēng)力機(jī)由3個直葉片組成,翼型為NACA0021。風(fēng)輪直徑為1.030 m,葉片長度為1.457 m。入流風(fēng)速為13.14 m/s,對應(yīng)的推力系數(shù)和湍流強(qiáng)度分別為0.46和1%,與此種工況相對應(yīng)的[C]設(shè)為4.15。

        圖3給出了三維尾流模型及風(fēng)洞數(shù)據(jù)在下游x/D=1.5處的橫向速度分布。從圖3可發(fā)現(xiàn),由于馬格努斯效應(yīng)和動態(tài)失速共同作用,風(fēng)洞數(shù)據(jù)的尾流中心線向左發(fā)生了偏移。三維尾流模型假設(shè)尾流速度剖面軸對稱,因未考慮上述因素,高估了風(fēng)輪中心線左側(cè)的速度。但三維尾流模型整體上與試驗數(shù)據(jù)吻合較好,能準(zhǔn)確地預(yù)測尾流速度的最大虧損。

        2.2 LES驗證

        Abkar等[6]采用LES和ALM相結(jié)合的方法,研究大氣邊界層內(nèi)兩種葉尖速比下([λ=2.5]、3.8)單臺三葉片H型VAWT的尾流特性,這兩種葉尖速比對應(yīng)的推力系數(shù)分別為0.34和0.64。該VAWT單機(jī)容量為200 kW,風(fēng)輪直徑為26 m,葉片長度為24 m,翼型為NACA0018,輪轂高度為40 m。入流風(fēng)速為7.5 m/s,對應(yīng)的湍流強(qiáng)度和風(fēng)切變指數(shù)分別為9.1%和0.154。本文選取[λ=2.5]的數(shù)據(jù)驗證三維模型的精度,[C]設(shè)為2.3。

        圖4為下游四位置處([x/D=3]、6、9、12)三維尾流模型和LES結(jié)果的橫向速度剖面。從圖4可發(fā)現(xiàn),下游[x/D=3]處,LES結(jié)果的尾流中心線發(fā)生偏移,三維尾流模型低估了風(fēng)輪

        中心左側(cè)的速度,高估了風(fēng)輪中心右側(cè)的速度。在下游更遠(yuǎn)處,這種偏移逐漸消失,所提出的三維模型與LES結(jié)果具有較高的吻合度。

        圖5為下游四位置處三維模型與LES數(shù)據(jù)的相對誤差。從圖5可發(fā)現(xiàn),在[x/D=3]處,相對誤差波動較大,這是由于尾流中心線的偏移造成的。值得注意的是,此種情況下最大相對誤差仍不超過6%。在更遠(yuǎn)距離處,三維模型展現(xiàn)出較高的預(yù)測精度,所有數(shù)據(jù)點的相對誤差均小于2%。

        transverse profile

        圖6比較了三維模型與LES結(jié)果在下游[x/D=3]、6、9、12這4個位置的垂向速度分布??砂l(fā)現(xiàn),本文三維尾流模型能較好地預(yù)測各下游位置處垂向速度分布。僅稍稍低估了[x/D=3]處輪轂附近、高估了[x/D=3]處[z/Dgt;2.4]處的速度。這是由于[x/D=3]處位于近尾流區(qū),湍流機(jī)制復(fù)雜,本文采用的經(jīng)驗湍流強(qiáng)度公式難以精確描述。在[x/Dgt;3]處,三維模型的預(yù)測結(jié)果與LES結(jié)果保持高度一致。

        圖7為三維模型和LES結(jié)果在垂向的相對誤差。從圖7可發(fā)現(xiàn),所有位置處的相對誤差均位于(-4%,4%)的區(qū)間內(nèi),最大誤差僅為3.4%,出現(xiàn)在[x/D=3]處。這表明本文所提出的三維模型在垂向上也具有較高的預(yù)測精度。

        3 模型預(yù)測

        基于上述提出的三維尾流解析模型,從空間角度研究VAWT的尾流特性。所預(yù)測的VAWT風(fēng)輪直徑為50 m,輪轂高度和葉片長度為100 m,入流風(fēng)速為9.6 m/s,C設(shè)為2。

        圖8為[CT]=0.8、[α]=0.15、[I0]=8.3%時下游四位置處([x/D=3]、6、9、12)[yoz]平面的速度云圖。結(jié)果表明,由于本文模型采用了各項異性的高斯函數(shù),尾流速度呈現(xiàn)為橢圓分布,這與實際更為相符。從圖8還可發(fā)現(xiàn),隨著下游距離的增加,尾流影響范圍增大,風(fēng)速逐漸恢復(fù),這可由流量守恒定律解釋。同時也可發(fā)現(xiàn),[x/D=3]處的風(fēng)速仍可通過二維高斯尾流模型來估計。然而,更遠(yuǎn)下游距離處的風(fēng)速分布變得復(fù)雜,二維模型無法預(yù)測,建立三維尾流模型是必要的。

        圖9為[α]=0.15、[I0]=8.3%時不同推力系數(shù)下([CT]=0.6、0.7、0.8)風(fēng)輪中心處[xoz]平面的速度云圖??砂l(fā)現(xiàn),由于考慮了風(fēng)切變的影響,尾流速度隨高度變化。而大多數(shù)風(fēng)電場布局研究采用一維或二維尾流模型,未考慮到高度因素,將會產(chǎn)生較大誤差。從圖9還可發(fā)現(xiàn),隨著推力系數(shù)的增加,速度損失越大,尾流恢復(fù)得更慢。

        提取[CT]=0.8、[I0]=8.3%時不同風(fēng)切變指數(shù)下([α]=0、0.1、0.15、0.2)風(fēng)輪中心處[xoz]平面的速度云圖,如圖10所示。當(dāng)未考慮風(fēng)切變時,[xoz]平面的尾流速度分布沿風(fēng)輪平面軸對稱;考慮風(fēng)切變后,尾流速度分布變得不對稱:隨著風(fēng)切變指數(shù)增加,輪轂高度以上的風(fēng)速增加,輪轂高度以下的風(fēng)速降低。

        圖11為[CT]=0.8、[α]=0.15時不同湍流強(qiáng)度下([I0]=5%、8.3%)輪轂高度處[xoy]平面的速度云圖。從圖11可發(fā)現(xiàn),隨

        著湍流強(qiáng)度的增加,風(fēng)輪后低速區(qū)域的長度明顯減少。這是由于更強(qiáng)的湍流導(dǎo)致尾流區(qū)內(nèi)外流體的滲混增加,促進(jìn)了尾流快速恢復(fù)。

        4 結(jié) 論

        基于流量守恒理論,本文提出一種H型VAWT的三維尾流模型,主要結(jié)論如下:

        1)該模型采用了多元高斯分布,在每個維度使用不同的高斯標(biāo)準(zhǔn)偏差,以精確反映VAWT橫向和垂向的尾流邊界變化差異。

        2)考慮了入流的風(fēng)切變效應(yīng)。由于考慮了風(fēng)切變效應(yīng),垂向尾流速度剖面呈現(xiàn)出不對稱分布。

        3)本文模型能精準(zhǔn)預(yù)測[x/Dgt;3]的遠(yuǎn)尾流區(qū)的速度分布,垂向和橫向的最大誤差分別為2.3%和1.8%。同時,能精確預(yù)測[x/D≤3]近尾流區(qū)的最大速度虧損。

        4)演示了由本文模型預(yù)測的一系列[yoz、xoz、xoy]平面的速度云圖。從[yoz]方向看,隨著下游距離的增加,尾流影響區(qū)范圍擴(kuò)大,風(fēng)速有所恢復(fù),風(fēng)速分布更加復(fù)雜;從[xoz]方向看,風(fēng)速隨高度的變化不可忽略;隨著推力系數(shù)的增加,尾流恢復(fù)得更慢;風(fēng)切變會降低輪轂高度以下的風(fēng)速,增加輪轂高度以上的風(fēng)速;從[xoy]方向看,更高的湍流強(qiáng)度導(dǎo)致尾流恢復(fù)得更快。

        [參考文獻(xiàn)]

        [1]"""" ZHAO Z Z, WANG D D, WANG T G, et al. A review: approaches for aerodynamic performance improvement of lift-type vertical axis wind turbine[J]. Sustainable energy technologies and assessments, 2022, 49: 101789.

        [2]"""" 羅帥, 繆維跑, 李春, 等. 基于定常吸氣的垂直軸風(fēng)力機(jī)流動控制研究[J]. 中國電機(jī)工程學(xué)報, 2020, 40(21): 7078-7087.

        LUO S, MIAO W P, LI C, et al. Research on flow control of vertical axis wind turbine based on steady suction[J]. Proceedings of the CSEE, 2020, 40(21): 7078-7087.

        [3]"""" BATTISTI L, ZANNE L, DELL’ANNA S, et al. Aerodynamic measurements on a vertical axis wind turbine in"" a"" large"" scale"" wind"" tunnel[J]. Journal of energy resources technology, 2011, 133(3): 031201.

        [4]"""" LI Q A, MAEDA T, KAMADA Y, et al. Investigation of power performance and wake on a straight-bladed vertical axis" wind"" turbine"" with"" field"" experiments[J]. Energy, 2017, 141: 1113-1123.

        [5]"""" SHAMSODDIN""" S,""" PORTé-AGEL""" F.""" Large""" eddy simulation"" of"" vertical"" axis"" wind"" turbine"" wakes[J]. Energies, 2014, 7(2): 890-912.

        [6]"""" ABKAR"" M,"" DABIRI"" J"" O."" Self-similarity"" and"" flow characteristics of vertical-axis wind turbine wakes: an LES study[J]. Journal of turbulence, 2017, 18(4): 373-389.

        [7]"""" PENG H Y, LAM H F, LEE C F. Investigation into the wake aerodynamics of a five-straight-bladed vertical axis wind turbine by wind tunnel tests[J]. Journal of wind engineering and industrial aerodynamics, 2016, 155: 23-35.

        [8]"""" ABKAR M. Theoretical modeling of vertical-axis wind turbine wakes[J]. Energies, 2018, 12(1): 10.

        [9]"""" CHEN Y, LI H, JIN K, et al. Wind farm layout optimization using genetic algorithm with different hub height wind turbines[J]. Energy conversion and management, 2013, 70: 56-65.

        [10]""" GE M W, WU Y, LIU Y Q, et al. A two-dimensional model based on the expansion of physical wake boundary for wind-turbine wakes[J]. Applied energy, 2019, 233/234: 975-984.

        [11]""" BARTHELMIE R J, LARSEN G C, FRANDSEN S T, et al. Comparison of wake model simulations with offshore wind turbine wake profiles measured by sodar[J]. Journal of atmospheric and oceanic technology, 2006, 23(7): 888-901.

        [12]""" GAO X X, LI B B, WANG T Y, et al. Investigation and validation" of" 3D" wake" model" for" horizontal-axis" wind turbines based on filed measurements[J]. Applied energy, 2020, 260: 114272.

        [13]""" SUN H Y, YANG H X. Study on an innovative three-dimensional wind turbine wake model[J]. Applied energy, 2018, 226: 483-493.

        STUDY OF THREE-DIMENSIONAL WAKE MODEL FOR

        VERTICAL-AXIS WIND TURBINE

        Ling Ziyan,Zhao Zhenzhou,Liu Yige,Liu Huiwen,Ma Yuanzhuo,Wang Dingding

        (College of Energy and Electrical Engineering, Hohai University, Nanjing 211100, China)

        Abstract:Considering flow conservation, a three-dimensional analytical wake model of H-type VAWT that incorporates wind shear effects and uses a multivariate Gaussian distribution is proposed. Wind tunnel tests and large eddy simulation data are firstly used to validate the model. The results show that the relative errors are less than 2% for the transverse profile and less than 3% for the perpendicular profile in the far wake region with [x/Dgt;3]. Next, a series of results are demonstrated from four downstream locations ([x/D=3],6,9,12), three thrust coefficients ([CT=0.6],0.7,0.8), four wind shear indices ([α=0],0.1,0.15,0.2), and two turbulence intensities ([I0=5%],8.3%), and the results show that the model in this paper can effectively describe the spatial distribution of the wake wind speed. By considering the height effect, the model can be used for the optimization of wind turbine hub heights and the layout optimization of wind farms, which is beneficial to increase wind farms’power output.

        Keywords:vertical axis wind turbines; analytical models; Gaussian distribution; three-dimensional wake model

        猜你喜歡
        高斯分布
        利用Box-Cox變換對移動通信中小區(qū)級業(yè)務(wù)流量分布的研究
        2種非對稱廣義高斯分布模型的構(gòu)造
        靜態(tài)存儲器用非晶硅薄膜晶體管輸出特性研究
        科技視界(2018年20期)2018-11-21 11:47:12
        航空發(fā)動機(jī)附件機(jī)匣齒輪磨損振動統(tǒng)計分析
        今日自動化(2018年1期)2018-10-21 03:13:51
        改進(jìn)的快速遺傳算法在函數(shù)優(yōu)化中的應(yīng)用
        在航集裝箱船舶搖擺姿態(tài)的概率模型
        不同分布特性隨機(jī)噪聲的FPGA實現(xiàn)
        科技視界(2017年7期)2017-07-26 07:45:41
        改進(jìn)的自適應(yīng)高斯混合模型運(yùn)動目標(biāo)檢測算法
        一種基于改進(jìn)混合高斯模型的前景檢測
        基于音節(jié)時間長度高斯擬合的漢語音節(jié)切分方法
        人妻少妇中文字幕久久hd高清| 久久久久成人亚洲综合精品 | 亚洲精品456在线播放狼人| 亚洲日韩成人无码| 曰批免费视频播放免费直播| 亚洲AV无码一区二区二三区我| 国产麻豆成人精品av| 女优一区二区三区在线观看| 欧美a级情欲片在线观看免费| 亚洲AV无码精品呻吟| 中文字幕一区二区区免| 午夜视频在线瓜伦| 人人妻人人澡人人爽久久av| 国产精品久久中文字幕第一页| 国产黄色一区二区三区,| 免费又黄又爽又色的视频| 亚洲欧洲中文日韩久久av乱码| AV中文码一区二区三区| 一本久久精品久久综合| 国产白袜脚足j棉袜在线观看| 91av国产视频| 日本激情一区二区三区| 综合亚洲伊人午夜网| 好吊色欧美一区二区三区四区| 男人阁久久| 精品亚洲国产日韩av一二三四区| 鲁丝片一区二区三区免费| 精品欧美乱子伦一区二区三区 | 丰满多毛的大隂户毛茸茸| 亚洲丁香婷婷综合久久小说| 中文字幕色婷婷在线视频| 欧美精品欧美人与动人物牲交| 精品国产18久久久久久| 国产一区二区三区杨幂| 国产av综合网站不卡| 国产精品熟女一区二区| 国产九九在线观看播放| 日本一区二区在线免费看| 国产精品无码一本二本三本色| 国产精品女同一区二区久久| 免费看片的网站国产亚洲|