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

穩定風速和陣風交錯作用下的林木風倒分析

2014-04-27 03:50:39賈杰
浙江農林大學學報 2014年4期
關鍵詞:風速模型研究

賈杰

(東北林業大學土木工程學院,黑龍江哈爾濱 150040)

穩定風速和陣風交錯作用下的林木風倒分析

賈杰

(東北林業大學土木工程學院,黑龍江哈爾濱 150040)

風倒是形成林隙的主要原因,是發生頻率較高、波及范圍較廣、造成損失較嚴重的一種自然災害。因此,研究風倒具有重要的理論和實踐意義。針對相關理論模型的不足,利用動力學理論建立了林木風倒的非線性動態模型,分別從穩定風速作用和陣風作用2個方面對林木的風倒機理進行了研究,通過嚴格的數學推導獲得了林木發生風倒的判斷條件。結果表明:穩定風速和陣風作用下林木發生風倒的條件不同,且互為補充。基于理論分析,以西加云杉Picea sitchensis為例進行數值計算,得到了風倒的標準風速值,研究結果可以應用到其他樹種的風倒問題中。同時,穩定風速和陣風交錯作用下發生風倒的標準風速值與還與風荷載的作用時間有關。研究結果為后續的理論研究和工程應用提供了指導。圖7參18

森林生態學;動力學模型;風倒;西加云杉;風速標準值;紊流強度

在眾多自然災害中,風荷載作用所引起的林木風倒是最常見一種,同時風倒也是形成林隙的主要原因,逐漸成為森林生態的一個擾動因子。風倒現象主要表現為干折、根折、連根拔起等3種方式。全世界每年非災害性風倒造成的木材減產在15%以上[1]。近年來,人類對自然環境的破壞日益嚴重,導致全球氣候環境的急劇變化,風倒引發的森林破壞頻繁發生。因此,展開樹木風倒機制的相關研究日益受到科學研究人員的重視,研究結果對于森林生態系統和木材生產的理論研究和生產實踐都有極其重要的意義。迄今為止,國內外很多專家學者開展了林木風倒機制方面的研究,并取得了豐富的研究成果。其中,李國旗等[2]研究了位于不同風壓和高度處樹木的應力分布規律。陳少雄等[3-4]、陳士銀[5]基于統計和野外調查的方法分別研究了泡桐Paulownia tomentosa,桉樹Eucalyptus robusta,木麻黃Casuarina equisetifolia等樹種的抗風性與個體特征之間的關系。代力民等[6]從生態學角度出發研究林隙的形成原因。宋曉鶴[7],王琳[8],賴秋明[9]分別從靜力學和動力學角度建立了風倒的模型,并展開了相應的分析。吳志華等[10]進行了抗風性狀對抗風效果的影響規律的研究。段溪[11]運用層次分析法研究了杭州灣地區亞熱帶泥質海岸的68種防護樹種的抗風性能。英國、加拿大、美國等一些學者對樟子松Pinus sylvestris var. mongolica等樹種的風倒機制進行了深入廣泛的研究,以動力學理論,結合野外試驗中林木和土壤力學參數的測定,地形數字高程模型和實物模型風洞試驗,獲得了復雜地形下和不同林分密度時的風場分布,建立了一系列力學方程來描述林木的抗風性[12],對風倒的力學機制也進行了相應的闡述,形成若干關于林木風倒機制的模型,如MC2[13]和HWIND[14]等。盡管這些研究代表了風倒機制研究的國際先進水平,但主要是基于材料力學方法的線性力學模型[15]。雖然描述了風倒機制的許多方面,但受研究條件等的限制在系統設計時對關鍵環節進行了簡化,造成模型與實際情況之間的誤差。作者在已有研究成果的基礎上,建立了林木發生風倒的非線性動力學模型[16],并且利用Von-Karman風譜模型從頻率域內分析了林木風倒機制[17],驗證了模型的有效性。相比其他模型,這個模型更接近實際物理系統且考慮了實際林木發生風倒的一些不確定性影響因素。但模型也存在一些缺陷,首先,并未考慮樹冠的形狀、樹木個體差異、枝條和樹葉的作用等對風倒的影響,其次,也未考慮風荷載作用下的“流固耦合”現象,這也是后續研究的方向和內容。本研究在前期的研究基礎上,著重考慮孤立木在穩定風速和陣風作用下所引起的傾覆,通過嚴格的數學推導給出風倒發生的條件,并以西加云杉Picea sitchenrsis為例,進行了數值計算,給出樹木發生風倒的標準風速值。研究結果對森林抗風防護建設和經營提供現實的指導意義。

1 模型建立

當風荷載作用于樹木的時候,空氣動力將會對樹桿和樹枝產生一個對于基礎部分的傾覆力矩,由根部所產生的力矩與其抵抗(抵抗力矩)。我們將傾覆力矩定義為mT(θ),抵抗力矩定義為mR(θ),并且假設樹為剛性體,我們可以得到如下的模型:

式(1)中:I0表示樹木相對根部的轉動慣量,θ表示樹桿在風力作用下相對于垂直位置的傾角。如圖1所示。以西加云杉為例,通過對其根部的抵抗力矩進行測量,發現抵抗力矩取決于根部鉸接部分的抵抗矩、土壤的壓力、迎風面根系的拉力、根系土壤的質量等等,如圖2所示。

由圖2可知:總體抵抗力矩是關于傾角θ的函數,是很多樹木拉力試驗中的典型曲線。并且當θ很小的時候,總體抵抗力矩急劇上升,在2°到5°左右,抵抗力矩達到最大值,在10°的時候,抵抗力矩降到其最大抵抗力矩的一半左右。基于本研究的目的,可以將抵抗力矩模擬成一個二次方程,在θ=θm的時候有最大抵抗力矩mR(θ)=Mm,在θ=0時,定義mR(θ)=M0。因此可以得到關于mR(θ)的方程:

圖1 樹木基本受力圖Figure 1 B asic force diagram of trees

圖2 總體抵抗力矩和傾角θ的關系Figure 2 R elationship between whole resistingmoment and angleθ

在0°到0.5°之間,由于其mR(θ)曲線陡然增加,近似用割線代替,從坐標軸上選取(0,M0)作為起點,得抵抗力矩模型曲線如圖3所示。

2 穩定風速計算分析

2.1 風倒原理分析

如圖1所示,空氣動力對于樹木一個微小的單元作用力為νn2d A,νn即對單元作用的風速,d A表示樹干上具有代表性的單元面積。假設風具有連續的水平速度u0,當樹木垂直的時候,風力對樹木根部產生的彎矩為Fm(u0),當樹木在風力作用下發生傾斜,傾角為θ的時候,其彎矩變為Fm(u0)cos2(θ),同時重心也將發生偏移,自身重力將產生一個傾覆力矩,定義為Wmsin(θ)。Wm表示θ=90°時,樹木自身重力對根部所產生的力矩。則方程變為:

圖3 抵抗力矩mR(θ)的模型圖Figure 3 M ode diagram of resistingmoment mR(θ)

當根部抵抗力矩較小時,可以近似地認為:sinθ≈θ,cos(θ)≈1-θ2/2,因此式變為:

且θ較小時,角加速度θ″(t)可用旋轉角速度ω=θ′(t)表示為:θ′(t)=ω′(t)=ωdω/dθ,則式(4)變為:

式(5)綜合了樹木旋轉時的動能,同時角速度ω也是θ的函數,則有:

由式(6)可知,當θ→0時,樹木在風載作用下開始傾斜,此時ω=0,C=0。則式(6)可表示為:

由以上分析可知,式(7)只有當θm<5°時才成立,此時θ可以用θm線性表示:

則式(7)可變為:

當φ取正值的時候,式(9)右邊如果也保持正值,那么角速度ω也將保持正值,這時隨著θ的增大,樹木將繼續傾斜,根據式(9)可以得到如圖4的函數圖像。從圖4中可知,當滿足條件(10)時式(9)沒有正根:

當u0足夠大能滿足不等式(10)時,樹木風倒便會發生。

針對圖1中的樹木模型,將樹干模擬成一個圓柱,其半徑為r0,高度為h,則作用于單元d y上的荷載為:ρacDr0u02d y。其

圖4 式(9)的函數圖像[a表示在條件(10)滿足的時候,b表示條件(10)不滿足的情況]Figure 4 Function diagram of the formula(a show the condition(10)is satisfied,b show the condition(10)is dissatisfied)

中:ρa表示空氣密度,cD表示無量綱的阻力系數(與樹干的半徑有關),u0表示作用于樹干的平均連續穩定風速。則風荷載相對于根部所產生的彎矩為:

式(12)中:ε=0.4ρgh3/(Er02),ρ表示樹木的密度,E表示彈性模量。式(12)相當于式(11)中的阻力系數cD被提高了1+ε倍。

2.2 數值計算

基于2.1的理論分析,以西加云杉為例進行數值計算,選取相應參數為:h=19 m,r0=0.1 m,ρ=900 kg·m-3,ρa=1.2 kg·m-3,cD=0.6。可計算得到抵抗力矩:θm=2.5°=0.044 rad,M0=8×103N·m,Mm=12.8× 103N·m。通過對根系強度的可變性進行分析,對上述大小的西加云杉樹木,其根倒彎矩值可能在10× 103N·m到50×103N·m,并且能給出圖1模型樹木對根部的轉動慣量:

實際上,樹木在風載作用下產生傾斜,樹木重心位置會發生偏移,重力矩Wm將會對彎矩Fm(u0)產生一個附加的作用,這個附加作用可以通過對式(11)進行修正來體現。

可得u0=29.9m·s-1,由于式(13)并未考慮根系和樹枝慣性矩的影響,因此計算得到風倒發生時的風速可能是被低估了的風倒標準風速。在臨界風速u0=29.9 m·s-1作用下,樹木至少需要2 s的時間讓傾斜達到θm,而需要4.5 s的時間讓樹木旋轉角速度ω達到最小值。總之,樹木如果受到連續穩定風速作用并產生彎矩Fm(u0),當滿足條件式(10),那么樹木將在連續的風速作用下發生傾覆。

除西加云杉以外,還對稻類植物、懸鈴木Platanus acerifolia也進行了相應的分析,基于各類參數計算得到稻類植物發生風倒的臨界風速u0=11.6m·s-1,懸鈴木發生風倒的臨界風速u0=22.8m·s-1,而根據所查資料[18],3種林木實際發生風倒的數值分別:26.3,8.8和19.1m·s-1。注意到,理論計算值和實際數據比較吻合,但是計算值看起來要比實際數據大一些。究其原因是:實際發生風倒與林木個體參數、土壤條件、根系條件以及枝條、樹葉與風荷載所形成的“流固耦合”效應有關。以典型的小麥Triticum aestivum來說,一般情況下都有2到6個分支,假設分支數目為n,如果每一個分支所承受的彎矩都有根系來支持,就降低了理論計算所得到的風倒臨界風速。再比如“流固耦合”效應導致流體邊界的變化,從而對流動產生反作用。不考慮“流固耦合”作用會使得計算結果是偏于安全的,也即“流固耦合”效應導致實際發生風倒的臨界風速降低。根系在風倒過程中有至關重要的地位,長期的連續風荷載作用,會降低根系的支持力和抵抗彎矩的能力,以至于林木可能在較小的風速下發生風倒。因此,理論計算值與實驗數據結果比較吻合,同時風的動態作用可以得到合理的解釋。

對于連續穩定風速引起的風倒,滿足條件式的最小風速為:

3 陣風(沖擊荷載)計算分析

實際上樹木更容易在強大的陣風作用下發生風倒,假設其作用時間為τ,在風速為u0連續風載作用下根部產生傾覆彎矩Fm(u0),則t=0時刻陣風所引起的沖量矩為最大值,樹木的角動量為M0)τ,其中Mgτ表示陣風所產生的沖擊荷載引起的沖量矩,M0τ表示陣風引起根部的抵抗力矩。則樹木的角速度為:

式(7)變為:

當t=0+時,積分常數C=I0ω02/2,陣風的作用,導致產生了力矩Mg,將Mg表達成連續穩定風速作用下彎矩Fm(u0)的函數:

顯然,式(18)右邊必須保持正值,樹木才會繼續發生傾斜,在其他所有參數已知的情況下,式(18)是關

以得到根的最大值。要保證式(18)右邊是正值,則必須保證:

成立。由式可以得到相應的判斷條件:

可以看出:式(20)是對式(10)的很好補充。其物理含義是:式(10)指出當風速足夠大的時候,在連續穩定風速的作用下,風倒便會發生。式(20)指出當連續風速不夠大,以至于不能夠引起風倒的情況下,連續的風載將和陣風共同作用下引起風倒[陣風系數GV足夠大,并且作用時間τ足夠長,能滿足式(19)]。式(19)取決于如下的無量綱量:則條件表示的破壞陣風系數必須滿足:

4 平均風速作用下的破壞規范

由無量綱化的風倒破壞式(21)可知,風倒條件取決于所定義的無量綱量,且都是風速u0和陣風作用時間τ的函數。以西加云杉為例進行數值計算,所選參數與2.2中的相同,則無量綱參數變為:

當陣風風速最大值?和作用時間τ滿足式(21),風倒便會發生,可得到風倒的風速臨界值:

式(23)中,σ/u表示紊流強度,對于獨立樹木,σ/u一般在0.2左右;而對于森林,σ/u會更高一些。依據式(23),給出不同的τ值,可以得到如圖6所示的一族直線,這些直線穿過圖5中的曲線,相應產生交點,將交點連在一起,如圖7所示,得到風倒發生的風速標準線。在這個直線下方,不可能出現陣風風速?。由圖7可知:風倒可能在較小風速17.7m·s-1就可能發生。而對于陣風來說,在風速為27.5m· s-1,作用時間5 s的情況下風倒就可能發生。盡管5 s時間給樹木以沖擊荷載看起來比較短,但于樹木來說,從受風載作用開始傾斜到其抵抗彎矩最大值位置,大約需要2 s。所以,可以尋找作用時間更短的陣風,對于陣風風速為35.1 m·s-1,作用時間τ=1 s疊加在風速為21.0 m·s-1的風速上就會引起風倒。

圖5表示出陣風風速?和連續風速u0(可以由陣風風速的平均值ū代替)在不同的時間段τ內的關系。從圖中可以看到,當u0接近于臨界風速u0=29.9m·s-1的時候,所有曲線都將終止。如果將1 h風速的作用時間按照τf時間段來劃分,可以得到這段時間內發生的最大風速。

圖5 引起風倒的陣風風速和連續風速的函數關系(在τ=0.2,0.5,1,2,5,10 s的情況下)Figure 5 Function relationship between gust wind and continuous wind about windthrow(in the case of τ=0.2,0.5,1,2,5,10 s)

圖6 紊流強度為0.2的時候破壞的平均風速圖Figure 6 Diagram of average wind speed when turbulence intensity is0.2

5 結論

林木的破壞模式有2種,一種是破壞值超過林木材料的允許值而發生的折斷現象,另外一種破壞值超過根系的支持力而發生的連根拔起現象。本研究綜合考慮在穩定風速和陣風作用下樹木的傾覆問題。將樹木風倒用一個簡單的模型來模擬根系,從動力學角度建立了風倒的動態方程,得到在穩定風速下樹木發生風倒的條件。同時將陣風作用模擬為一個沖量荷載作用在樹木上部,建立陣風作用下樹木風倒的條件。以西加云杉為例,計算了樹木在穩定風速和陣風作用下破壞的標準值,將陣風風速的作用時間與平均風速作用時間比較,得到樹木發生風倒的極限風速。但是,實際林木風倒臨界風速值低于文章所計算出的理論值,這是由于林木樹冠形狀、個體參數差異、枝條、樹葉與風荷載之間的“流固耦合”效應以及復雜的非線性動力學行為都導致實際發生風倒的風速值偏低。這也為今后的研究提供了方向和思路。總之,文章從理論方面的研究結果對森林的抗風防護以及木材的生產實踐有很重要的現實指導意義。

圖7 紊流強度為0.2,τ=1,2,5 s破壞的平均風速圖Figure 7 Average wind speed diagram of damage when turbulence intensity is 0.2,τ=1,2,5 s

[1]于振良,趙士洞.林隙(gap)模型研究進展[J].生態學雜志,1996,16(2):42-46.

YU Zhenliang,ZHAO Shidong.Research progress on gap model[J].Chin JEcol,1996,16(2):42-46.

[2]李國旗,安樹青,張紀林,等.海岸帶防護林4種樹木的風壓應力分析[J].南京林業大學學報,1999,23(4):76-80.

LIGuoqi,AN Shuqing,ZHANG Jilin,et al.The bending stress analysis of 4 species of woods caused by wind pressure in coastal shelter forest[J].JNanjing For Univ,1999,23(4):76-80.

[3]陳少雄,王觀明,羅建中.桉樹幼林不同株行距配置抗臺風效果的研究[J].林業科學研究,1995,8(5):582 -585.

CHEN Shaoxiong,WANG Guanming,LUO Jianzhong.Effect of different spacing on typhoon resistance of youngEucalyptusplantation[J].For Res,1995,8(5):582-585.

[4]陳少雄,楊民勝,王理平.尾葉桉造林密度與蓄積量、抗風和材性關系研究[J].林業科學研究,1998,11(4):435-438.

CHEN Shaoxiong,YANGMinsheng,WANG Liping.Effect of spacing on volume,storm-resistance and wood qualityof Eucalyptus urophylla[J].For Res,1998,11(4):435-438.

[5]陳士銀.木麻黃不同株行距配置抗臺風效果[J].廣東林業科技,1997,13(2):44-46.

CHENShiyin.Effectsofspacing on the resisitance to typhoon inCasuarina equisetifoliaforest belt[J].JGuangdong For Sci Technol,1997,13(3):44-46.

[6]代力民,徐振邦,楊麗韞.紅松闊葉林倒木儲量動態的研究[J].應用生態學報,1999,10(5):513-517.

DAILimin,XU Zhenbang,YANG Linyun,et al.Storage dynamics of fallen trees in Korean pine broad-leaved forest[J].Chin JAppl Ecol,1999,10(5):513-517.

[7]宋曉鶴.云杉風倒靜力學模型的建立及其分析[D].哈爾濱:哈爾濱工業大學,2006.

SONG Xiaohe.Builing and Analysis of Mechanics Model for W indthrow of Sp ruce[D].Harbin:Harbin Institute of Technology,2006.

[8]王琳.云杉風倒動力學模型的建立與分析[D].哈爾濱:哈爾濱工業大學,2006.

WANGLin.BuildingandAnalysisofDynamicModel forWindthrow ofSpruce[D].Harbin:Harbin Instituteof Technology,2006.

[9]賴秋明.云杉風倒動力學問題的研究[D].哈爾濱:哈爾濱工業大學,2008.

LAIQiuming.Dynamics Research forWindthrow of Spruce[D].Harbin:Harbin Institute of Technology,2008.

[10]吳志華,李天會,張華林,等.沿海防護林樹種木麻黃和相思生長和抗風性狀比較研究[J].草業學報,2010,19(4):166-175.

WU Zhihua,LI Tianhui,ZHANG Hualin,et al.Studies on growth and wind-resistance traits ofCasuarinaand Aacia stands from coastal protection forest[J].Acta Pratac Sin,2010,19(4):166-175.

[11]段溪.亞熱帶泥質海岸68個防護林樹種抗風性能比較分析[D].南京:南京林業大學,2012.

DUAN Xi.A Comparative Analysis ofWind Resistance for 68 Shelter Tree Species in Subtropical Muddy coast[D]. Nanjing:Nanjing Forestry University,2012.

[12]RUEL JC,PIN D,SPACE L,et al.The estimation of wind exposure for windthrow hazard rating:com parison between Strong blow,MC2,topex and awind tunnel study[J].Forestry,1997,70(3):254-266.

[13]ENGLAND A H,BAKER D J,SAUNDERSON S E T.A dynamic analysis of windthrow of trees[J].Forestry,2000,73(3):225-237.

[14]BENOITR,DESGAGNéM,PELLERIN P,et al.The Canadian MC2:a semi-lagrangian,semi-implicitwideband atmosphericmodel suited for fine-scale process studies and simulation[J].Mon Weather Rev,1997,125(10):2382 -2415.

[15]ACHIM A,RUEL JC,GARDINER B A,et al.Modeling the vulnerability of balsam fir forests to wind damage[J].For Ecol Manage,2005,204:35-50.

[16]賈杰,李靜輝.林木風倒動態模型的建立與分析[J].森林工程,2013,29(4):63-67.

JIA Jie,LIJinghui.Establishmentand analysisofdynamicalmodel forwindthrow of trees[J].For Eng,2013,29(4):63-67.

[17]賈杰,李洪峰,李靜輝.基于Von-Karman譜的獨立林木風倒機理的頻率域分析[J].東北林業大學學報,2013,41(8):165-169.

JIA Jie,LIHongfeng,LIJinghui.Frequency domain analysis of independent treeswindthrow based on Von-Karman spectrum[J].JNortheast For Univ,2013,41(8):165-169.

[18]MERRENSE J,PEART D R.Effects of hurricane damage on individual growth and stand structure in a hardwood structure in a hardwood forest in New Hampshire,USA[J].JEcol,1992,80(4):787-795.

Tree windthrow with stable winds and gusts

JIA Jie
(College of Civil Engineering,Northeast Forestry University,Harbin 150040,Heilongjiang,China)

A major reason for gaps in forests is windthrow,which is a kind of natural disaster with high frequency,broad cope and much wastage making it very important to understand the theoretical and practical aspects ofwindthrow.In this study,a nonlinear dynamicmodel ofwindthrow was established using dynamics theory to improve the shortcomings of the related theoreticalmodel.The windthrow mechanism forPicea sitchensiswas studied from two aspects:stable winds and gusts.The judgment conditions of tree windthrow were obtained through strictmathematical deductions.Results showed that windthrow condition with stable winds and gusts was different,which are complementary each other.Also,the standard wind speed value was related to action wind load time for the interaction between stable winds and gusts.These research results could be applied to windthrow problems of other species and could provide guidance for subsequent theoretical research and engineering applications.[Ch,7 fig.18 ref.]

forest ecology;dynamicsmodel;windthrow;Picea sitchensis;standard wind value;turbulence intensity

S718.5;S727.2

A

2095-0756(2014)04-0604-07

2013-11-29;

2014-02-27

黑龍江省教育廳科學技術研究項目(12513033)

賈杰,講師,從事非線性動力學的分析與控制研究。E-mail:jiajiecontrol@126.com

猜你喜歡
風速模型研究
一半模型
FMS與YBT相關性的實證研究
遼代千人邑研究述論
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 中文字幕免费视频| 曰韩免费无码AV一区二区| 日本手机在线视频| 国产三区二区| 欧美精品在线看| 青青草一区| 亚洲人成网站观看在线观看| 亚洲二区视频| 精品一区二区无码av| 亚洲天堂色色人体| 激情爆乳一区二区| 看看一级毛片| 欧美第一页在线| 国产激情无码一区二区APP| 国产免费人成视频网| 色有码无码视频| 日本91在线| 亚洲一区二区精品无码久久久| 91国语视频| 亚洲国产高清精品线久久| 亚洲成人www| 日韩无码真实干出血视频| 中文字幕不卡免费高清视频| 青青热久免费精品视频6| 2021国产精品自拍| 亚洲无码高清一区| 99热这里只有精品在线观看| 久久99热这里只有精品免费看| 四虎成人精品| 日本高清在线看免费观看| 欧美综合一区二区三区| 大陆精大陆国产国语精品1024| 动漫精品啪啪一区二区三区| 午夜无码一区二区三区在线app| 国产在线拍偷自揄观看视频网站| 97国产在线视频| 三区在线视频| 久久超级碰| 日本一区二区三区精品国产| 欧美午夜久久| 一级一毛片a级毛片| 亚洲精品视频在线观看视频| 农村乱人伦一区二区| 亚洲国产欧美自拍| 天天做天天爱夜夜爽毛片毛片| 新SSS无码手机在线观看| 成人午夜网址| 久久精品一品道久久精品| 欧美精品综合视频一区二区| 日韩色图在线观看| 高潮爽到爆的喷水女主播视频 | 亚洲色图欧美在线| 就去吻亚洲精品国产欧美| 国产女人18毛片水真多1| 国产在线视频导航| 亚洲高清日韩heyzo| 亚洲区欧美区| 国产精品太粉嫩高中在线观看| 久久免费观看视频| 91九色最新地址| 麻豆精品视频在线原创| 无码人妻热线精品视频| 99久久精品久久久久久婷婷| 九九这里只有精品视频| 亚洲精品第五页| 黄片一区二区三区| 少妇精品在线| 亚洲成人网在线播放| 99久久精品免费看国产电影| 国产精品永久久久久| 精品一區二區久久久久久久網站 | 欧美成一级| 在线免费观看AV| 国产波多野结衣中文在线播放 | 亚洲成网777777国产精品| 五月综合色婷婷| 亚洲 日韩 激情 无码 中出| 亚洲中文久久精品无玛| 国产18在线播放| 欧美日韩精品一区二区视频| 国产高清国内精品福利| 一本大道无码日韩精品影视|