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

風浪聯合作用下的海上單樁基礎風力發電機組動力響應分析

2012-03-04 12:08:58賀廣零仲政
電力建設 2012年5期
關鍵詞:風速結構模型

賀廣零,仲政

(1.同濟大學力學博士后流動站,上海市,200092;

2.中國電力工程顧問集團華北電力設計院工程有限公司,北京市,100120; 3.同濟大學航空航天與力學學院,上海市,200092)

0 引言

近年來,伴隨著全球海上風電場大規模建設的趨勢,我國海上風能技術亦獲得了高度重視與迅猛發展。然而,由于根基薄弱與底蘊不足,迄今我國海上風能技術還不夠成熟。就海上風機基礎而言,現有的研究成果以擬靜力分析為主[1-2],不考慮動力效應,存在低估結構響應、荷載組合困難等難題。已有部分學者[3-4]嘗試進行海上風力發電機組結構動力響應分析,但并未給出合理的風浪耦合機制,所采用的結構模型亦顯粗糙,與工程應用相去較遠,參考價值不大。有鑒于此,本文提出一種相對合理的風浪耦合機制,結合實際工程中的海上風力發電機組結構模型,進行結構在風浪聯合作用下的動力響應分析,以正確評估結構的動力效應。

1 環境荷載模型

海上風力發電機組外部環境十分復雜,不僅承受風、海浪和地震作用,還與浮冰、海流、潮汐和沖刷等因素相關。而風力發電機組也有自己的運行狀態及多種失效模式。對于不同的運行狀態,受到不同外部激勵時所產生反應各不相同(強度破壞或疲勞失效等)。由于風力發電機組環境荷載和運行狀態自身的復雜性,加上兩者的隨機耦合,就導致應用典型規范需要考慮諸多工況。在現有的理論模型中,還很難把所有環境荷載的相關性考慮完整,只能考慮主要荷載,如風和海浪的共同作用等。

1.1 隨機Fourier風譜模型

基于隨機過程的隨機函數描述,引入截尾函數,則隨機Fourier譜[5]可定義為

式中:T為樣本持時;隨機過程X(η,t)是樣本x(t)的集合,η為影響隨機激勵發展過程且具有物理意義的隨機變量或隨機向量。根據各向同性紊流理論,可確定隨機Fourier譜的基本表達式。依據隨機建模準則,由310組實測風速數據記錄,可以給出基本隨機變量的概率分布和待定擬合參數的具體值,并最終確定隨機Fourier譜的表達式[6]為

式中:u為10 m高平均風速,服從極值I型分布;z0為地面粗糙度,服從對數正態分布。顯然,塔體上任意一點處脈動風速的隨機特性可以通過隨機Fourier譜來體現。與塔體不同的是,槳葉在運行過程中存在顯著的旋轉效應,將會對風速場產生很大的干擾。本文采用基于物理機制的旋轉Fourier譜來考慮脈動風速的旋轉效應[7-8]

式中:Fii(n)為隨機Fourier譜;~Fii(n)為旋轉Fourier譜;km(n)為銜接二者的旋轉模態。事實上,旋轉Fourier譜為作用在旋轉槳葉上的風速的相關函數經過Fourier變換所得,是一種自身蘊含了槳葉旋轉效應的紊流風速譜。

對隨機 Fourier譜與旋轉 Fourier譜進行逆Fourier變換,即可分別獲得塔筒與葉片上的風速時程,具體步驟詳見文獻[8]。海上風力發電機組葉片結構頗為復雜,且存在顯著的旋轉效應、風剪效應、尾流效應等,本文引入廣泛接受的葉素動量(blade elementmomentum)理論[9-10]以準確而簡單地確定作用在葉片上的風荷載,并考慮旋轉槳葉對風速的減緩作用。對于作用于塔筒上的風荷載,可借鑒建筑結構領域的相關研究成果確定[10]。

1.2 隨機Fourier海浪譜模型

將海浪視為一個平穩的隨機過程,自然地可以引入功率譜的概念描述波能及其頻域分布特征,但更為合理的方式是依據風生波的物理機制發展Fourier海浪譜模型。基于擬層流風波生成機制,即給定頻率處海浪譜的能量等于具有該頻率的諧波自具有對數剖面的平行氣流中汲取而來,文獻[11-12]建立了隨機Fourier海浪譜模型。

式中:設10 m高、10 min平均風速u1為等效風速u,m/s;=μω/ωP,ωP為峰值頻率,μ為譜峰頻率調整系數;β()為能量傳遞系數;A()為波幅;γ為譜峰值調整系數;σ=σL(ω<ωP)和σ=σR(ω>ωP)為譜形參數;ρa為空氣密度,kg/m3;ρw為水密度,kg/m3;g為重力加速度,m/s2。

由于隨機Fourier海浪譜基于隨機過程的隨機函數表示建模,對其進行逆Fourier變換,并考慮隨機相位由隨機初相位和相位差譜確定,采用式(4)生成海浪時程。

式中:F(ωj)為隨機Fourier海浪譜模型在離散點處的值;φ0j為隨機初相位;Δφj為相位差譜;n為離散點數目。

根據文獻[13],如果構件直徑與波長相差很大(D/L<0.2),則可采用Morison公式來計算結構所受波浪荷載。本文所采用的模型直徑最大為6 m,為小尺度結構物,其對波長的影響不明顯,故采用Morison公式計算結構承受的波浪力是合理的。該模型將波浪力分為慣性力和拖拽力2個部分,前者與水質點運動的加速度有關,后者與水質點的速度有關。在Morison公式中,作用在單位長度直立圓柱上的波浪力為

式中:CD是拖拽力系數;CM是質量系數;ρ是海水的密度,kg/m3;ux為質點水平速度,m/s;˙ux為質點水平加速度,m/s2,由選用的波浪理論確定。

1.3 風浪耦合物理機制

風浪耦合一直是海洋工程研究的難點。文獻[13-16]對風浪耦合都有相應規定,其內容均具備復雜性、多樣性與保守性的特點。顯然,上述規范之所以存在這些共性,是因為風浪耦合的物理機制尚不明確。在無法把握物理機制的情況下,可以嘗試建立半經驗關系。Neumann和Pierson[17-18]給出了有效波高與平均風速之間的半經驗關系,Ochi[19]研究了颶風條件下二者之間的關系,以實現風浪耦合。然而,半經驗關系成立的前提條件比較苛刻,例如要求海洋狀態必須充分發育,且僅適用于深海區域,否則將會導致較大的誤差。另外,一些學者意圖純粹從數學角度來實現風浪耦合。其中,較為成功的是Turkstra法則[20-23]與多元極值理論[24]。Turkstra法則本質上是將控制荷載的最大值與其他荷載時程進行隨機組合。由于Turkstra法則是基于工程經驗提出的一種組合規則,盡管有效地簡化了荷載組合過程,但它卻很有可能低估荷載組合值[22,25]。多元極值理論是研究多個相關變量極值性質的理論,探討了確定多個相關變量極值聯合分布的方法,是考慮各種極端荷載共同作用的有效手段,在海洋工程中獲得了廣泛應用[26-27]。但由于多元極值分布中相關結構較為復雜,其相關函數多為隱式形式,只能通過繁雜的迭代求解,不利于工程應用。

應該強調,在一般的荷載(荷載效應)組合中,對不同荷載的效應最終通過若干組合系數或某種特殊組合方式[平方和開平方(square root of the sum of the squares,SRSS)]等來根據單個荷載(荷載效應)的極值獲得共同作用時的極值,這在本質上是擬靜力分析與設計的結果。正因為僅僅進行擬靜力分析,才出現各個荷載的極值通常不同時出現,因而不能簡單相加,由此導致荷載(荷載效應)組合的各種理論或方式。事實上,由于荷載與荷載效應是非線性的且不具有簡單的疊加性質,最大荷載組合與最不利效應之間不能建立一一對應的關系,因此上述組合方式在考慮非線性的情況下也是不適用的。更為普遍適用的方法是進行多種隨機作用下的結構隨機動力反應分析。在直接動力分析中,對不同隨機荷載按照上述方式進行加權組合既無必要、也不正確,一種簡單且有效的方法是基于物理機制來實現荷載組合。在本文所采用的隨機Fourier風速譜和隨機Fourier海浪譜模型中,它們的基本物理參數集分別為{u,z0}和{u,ωP,γ,μ,σL,σR}。可以看出,風和海浪的基本物理隨機變量有1個相同的參數u,而正是這個共有的參數給出了風與海浪之間的物理關系,揭示了風浪耦合的物理機制。這樣,風浪耦合可以通過考慮其隨機Fourier譜中基本參數的耦聯來生成風速和海浪時程。

2 一體化結構模型

所謂風力發電機組一體化結構模型,是指將槳葉、機艙、塔體和基礎同時建模,以模擬結構自身構件以及結構與其他介質之間的動力相互作用,包括不同構件(尤其是葉片與塔體)之間的動力相互作用、土-結構動力相互作用、氣彈效應等,并反映結構應力集中、局部屈曲等細部特征。總體上,風力發電機組一體化建模的目的在于構建精細化的結構模型,把握準確的結構動力特性,獲取精確的結構動力響應。海上風力發電機組結構如圖1所示。

圖1 海上風力發電機組結構Fig.1 Offshore wind turbine system

與陸上風力發電機組不同的是,海上風力發電機組囿于環境條件,需要采用獨特的基礎形式,以滿足承載力與正常使用要求。重力式基礎、單樁基礎、多腳架基礎、導管架基礎與高樁承臺基礎均為常見的海上風機基礎形式,其中又以單樁基礎應用最為廣泛。縱觀常見的海上風機基礎,大部分需要采用超大直徑鋼管樁。然而,由于土體動力弱化效應、土塞效應等現象的存在,迄今尚未出現適合超大直徑鋼管的樁土模型。由于在分析與設計過程中未全面考慮這些因素,已建成的海上風機基礎(尤其是單樁基礎)均出現了或多或少的傾斜現象。其中,部分風機基礎因傾斜過大,導致風力發電機無法正常運行。土體動力弱化效應是指在循環荷載下樁周土體剛度隨著循環次數的增加而降低。文獻[14]明確規定需考慮該效應,但未給出具體的分析方法。在考慮土體動力弱化效應的分析方法中,最常見的是文獻[15]建議的p-y曲線分析法。首先,該法主要基于循環次數少于200次的現場實驗,而循環次數大于200次時,土體動力弱化現象并未終止[28-29]。其次,該法中的p-y曲線是基于小直徑樁基礎(樁徑小于1.5 m)現場實驗所得[30-31],故而該法并不適用于超大直徑樁基礎(樁徑超過3 m),且極有可能低估其動力響應[32-33]。為此,Juirnarongrit和Ashford[34]研究了基樁尺寸對p-y曲線的影響,提出了考慮尺寸效應的p-y曲線法,本文即采用此模型進行分析。在有限元模型中,樁-土相互作用可采用彈簧單元來模擬:水平方向樁土相互作用采用水平彈簧,其剛度定義根據p-y曲線來求得;豎向樁土相互作用采用豎向彈簧,其剛度利用t-z曲線來確定;端部彈簧根據q-z曲線來確定其剛度。彈簧模型原理如圖2所示。

圖2 樁-土相互作用有限元模型Fig.2 Finite elementmodel of pile-soil interaction

3 工程實例

以我國東海海域某一具有可行性的場址為例,結合我國某知名風機廠家提供的3.6 MW海上風力發電機組,進行海上風電機組結構在風浪聯合作用下的動力響應分析。

3.1 工程地質與海洋水文

本場地最大勘探揭露深度約為80 m,揭露的地基土層按地質時代、成因類型、土性的不同和物理力學性質的差異可分為7個大層,其中⑤層、⑦層各分為2個亞層,⑦1層又分為2個次亞層。各土層的埋藏條件及工程地質特性詳見表1。海底較平緩,在潮流作用下以淤積為主,灘面表層主要為淤泥,局部夾薄層粉土。本建筑場地屬于Ⅳ類,場地沿線屬于抗震不利地段。

海底高程為-10.00~-10.67 m(國家85 m高程),平均海平面高程0.23 m。平均高潮位為1.86 m,平均低潮位為-1.34 m,設計高潮位為2.55 m,設計低潮位為 -2.09 m。場址區相應歷史最高潮位為4.02 m,歷史最低潮位為-2.87 m。本海區設計潮流流速中,表層流速為315 cm/s,中層流速為257 cm/s,底層流速為148 cm/s,流速平均值為239 cm/s。平均波高為2.83 m,波周期為7.76 s,波長為74.1 m,波速為9.55 m/s,累積率1%的波高H1%為5.81 m,累積率4%的波高H4%為5.06 m,累積率5%的波高H5%為4.82 m,累積率13%的波高H13%為4.24 m。

表1 地層特性表Tab.1 Formation characteristics

3.2 結構模型

以我國某知名風機廠家提供的W 3600M-116型海上風力發電機組(裝機容量為3.6 MW)為例進行分析。風輪直徑為116 m,輪轂中心到塔低的高度為80 m,輪轂中心到平均海平面高度為90 m。組裝后風輪的質量為95 t,機艙質量為140 t,塔筒共4節,總質量為284 t,塔筒頂部直徑為3.815 m,塔頂壁厚為66 mm,底部直徑為5 m,塔底壁厚為88 mm。基礎采用單樁基礎,鋼管樁直徑6 m,樁長65 m,壁厚75 mm,入土深度約50 m,樁尖高程為-62 m,樁尖進入第⑦2層粉細砂層,樁頂高程為3 m。操作平臺高程為11 m。單樁基礎與塔筒之間設連接段鋼管過渡,連接段鋼管四周設置靠船設施、鋼爬梯及平臺等,連接段鋼管頂面通過法蘭與風機塔筒連接,連接段鋼管下端通過高強灌漿連接方式與鋼管樁連接,并設置剪力栓。整體結構尺寸如圖3(a)所示。

以大型通用有限元軟件ANSYS為建模平臺,建立海上風力發電機組一體化有限元模型,如圖3(b)所示。因葉片和塔筒1個方向的尺寸與另外2個方向的尺寸相差較大,同時葉片在工作狀態下具有顯著的應力剛化現象,故葉片和塔筒都采用了能較好體現這些特征的八節點殼體單元(SHELL 91)。在整體分析過程中,機艙及其內部構件可視為一個整體,可借助梁單元(BEAM 189)來模擬。單樁基礎依據PIPE 59單元和PIPE 16單元共同來體現。PIPE 59單元專門用于模擬浸沒在水中的桿件,可以很好地模擬海浪、海流對單樁基礎的作用。因此,采用PIPE 59單元模擬單樁基礎在水中的部分,對于泥面以下樁柱則采用PIPE 16單元來模擬。如前文所述,樁土相互作用通過非線性彈簧單元COMBIN 39來刻畫。在有限元模型中,劃分的PIPE 59單元和PIPE 16單元尺寸為0.1 m,COMBINE 39彈簧單元每隔0.5 m設置1組,遇土層分界點加設1組彈簧單元。

圖3 海上風力發電機組整體結構Fig.3 Integrated model of offshore w ind turbine system

3.3 結構動力響應分析

依據擬建風電場址的風場以及水文地質資料,可確定其基本風速 u為31.2 m/s(對應風壓為0.6 kPa),地面粗糙度z0可取為1 mm[13],波浪峰值頻率ωP為0.566,譜峰值調整系數γ為8.274,譜峰頻率調整系數μ為1.728,譜形參數σL和σR分別為0.161和0.764。依據上述物理方法,模擬特定環境條件下的風速時程(圖4)、風壓時程(圖5)、水質點速度時程(圖6)、水質點加速度時程(圖7)以及波浪力時程(圖8),將風壓時程、波浪力時程作用于海上風力發電機組一體化結構模型上,利用有限元動力分析方法,進行在風浪聯合作用下的結構動力響應分析,可獲得泥面處的水平位移動力時程響應、平臺處的水平位移時程響應、樁身彎曲應力時程響應分別如圖9~11所示,并進而獲得結構動力響應最值(表2)。在分析過程中,采用瑞利阻尼[35]。考慮到風力發電機組的底阻尼特性,結合美國學者Prowell等的實驗結果[36-37],本文取阻尼比0.6%。

根據西歐各地區海上風機基礎的設計與施工經驗,單樁基礎適用于水深約為15m、地基土為中等密實土體及土層的內摩擦角在30°以上的環境條件。事實上,西歐各地區能滿足上述條件的海域頗為寬廣,加上單樁基礎的經濟性與簡潔性,單樁基礎在諸多海上風機基礎中應用最為廣泛。與之不同的是,我國各海域的地質條件較為薄弱,例如,本海域海床上部土層為粉質粘土及粉砂,地基剛度相對較小。為了驗證對單樁基礎進行地基處理的必要性,本文就地基未處理以及地基處理2種情況進行了對比分析(表2):(1)地基未處理,按原有地質條件進行分析;(2)地基處理,對海床5 m深度范圍內的土體進行注漿碎石加固,圍繞樁身位置填筑出一個高約為5 m、頂邊半徑為15 m、坡度為1∶1的圓臺體,碎石頂表面鋪1層高強土工網裝碎石(直徑大于150 mm的碎石占80%以上)以防水流沖刷。研究表明,若未進行地基處理,泥面處水平位移以及平臺處水平位移均偏大,從而引起單樁基礎的轉角偏大,最終導致海上風力發電機組傾斜過度而無法正常運行。值得說明的是,在實際工程中,這種由單樁基礎轉角過大而導致海上風力發電機組無法正常運行的現象屢見不鮮。究其原因,在于對海上風力發電機組在風浪聯合作用下的動力效應與時間效應認識不足,尚未明確2種效應導致的地基弱化現象。為此,有必要研究海上風力發電機組的精細化荷載模型與結構模型(尤其是樁土模型),并輔以必要的地基處理方法,方可確保結構可靠穩定。

4 結論

(1)基于隨機過程的隨機函數描述,給出了風速與海浪隨機物理模型,即隨機Fourier風速譜與隨機Fourier海浪譜。依據隨機物理模型中共有隨機變量的耦聯關系,提出了基于物理機制的風浪耦合方法。

(2)鑒于海上風電機組單樁基礎采用超大直徑管樁,引入了適用于超大直徑管樁的樁土相互作用機制,構建了海上風力發電機組一體化結構模型。

(3)以我國東海海域某具有可行性的場址為例,結合我國某知名風機廠家提供的3.6 MW海上風力發電機組,進行海上風電機組結構在風浪聯合作用下的動力響應分析,并對工程設計給出了若干建議。

[1]黃維平,劉建軍,趙戰華.海上風電基礎結構研究現狀及發展趨勢[J].海洋工程,2009,27(2):130-134.

[2]李靜,陳健云.海上風力發電結構動力研究進展[J].海洋工程,2009,27(2):124-129.

[3]李德源,劉勝祥,張湘偉.海上風機塔架在風波聯合作用下的動力響應數值分析[J].機械工程學報,2009,45(12):46-52.

[4]李德源,劉勝祥,張湘偉.風波聯合作用下的風力機塔架疲勞特性分析[J].太陽能學報,2009,30(10):1250-1256.

[5]Li J,Chen JB.Stochastic dynamics of structures[M].Singapore: John Wiley&Sons,2009:1-216.

[6]李杰,張琳琳.實測風速資料的隨機Fourier譜研究[J].振動工程學報,2007,20(1):66-72.

[7]賀廣零,李杰.風力發電高塔系統基于物理機制的旋轉樣本功率譜研究[J].中國電機工程學報,2009,29(26):85-91.

[8]賀廣零,李杰.風力發電高塔系統風場模擬[J].同濟大學學報,2010,38(7):976-981.

[9]Burton T,Sharpe D,Jenkins N,et al.Wind energy handbook[M].Chichester:John Wiley&Sons,2001:41-172.

[10]賀德馨.風工程與工業空氣動力學[M].北京:國防工業出版社,2006:80-91.

[11]徐亞洲,李杰.風浪相互作用的Stokes模型[J].水科學進展,2009,20(2):281-286.

[12]徐亞洲,李杰.近海風力發電高塔波浪隨機動力響應分析[J].振動工程學報,2011,24(3):315-322.

[13]Det Norske Veritas.DNV-OS-J101 Design of offshore wind turbine structures offshore standard[S].1st Edition.Oslo,Norway:Det Norske Veritas,2004.

[14]Germanischer Lloyd.Rules and guidelines IV—Industrial services,part2:guideline for the certification of offshore wind turbines[S].Hamburg:Germanischer Lloyd,2005.

[15]American Petroleum Institute.RP2A-WSD recommended practice for planning,designing and constructing fixed off-shore platformsworking stress design[S].API Recommended Practice 2A-WSD,21st edition,Dallas:American Petroleum Institute,2000.

[16]JTS 144-1—2010港口工程荷載規范[S].北京:人民交通出版社,2010.

[17]Neumann G,Pierson W J.Principle of physical oceanography[M].Englewood Cliffs:Prentice Hall,1963:1-24.

[18]Kuhn M J.Dynamics and design optimization of offshorewind energy conversion systems[D].Delft:Delft University of Technology,2001.

[19]Ochi M K.On hurricane-generated seas[C]//Proceedings of the Second International Symposium on Ocean Wave Measurement and Analysis.New Orleans,USA:ASCE,1993:374-87.

[20]Turkstra C J,Madsen H O.Load combinations in codified structural design[J].Journal of the Structural Division,1980,106(12): 2527-2543.

[21]金偉良.工程荷載組合理論與應用[M].北京:機械工業出版社,2006:10-30.

[22]Cheng P W.A reliability based design methodology for extreme response of offshore wind turbine[D].Delft:Delft University of Technology,2002.

[23]陳為飛.近海風機塔架風浪荷載分析[D].杭州:浙江大學,2010.

[24]Coles SG,Tawn JA.Modelling extreme multivariate events[J].J R Statist Soc B,1991,53(2):377-392.

[25]Melchers R E.Structural reliability,analysis and prediction[M].Chichester:John Wiley&Sons,1999:199-215.

[26]劉德輔,董勝.隨機工程海洋學[M].青島:中國海洋大學出版社,2004:89-118.

[27]歐進萍,段忠東,肖儀清.海洋平臺結構安全評定:理論、方法與應用[M].北京:科學出版社,2003:12-49.

[28]Hettler A.Deformation responses of rigid and elastic foundations in sand under static and cyclic loadings[D].Karlsruhe:University of Karlsruhe,1981.

[29]Achmus M,Abdel-Rahman K,Kuo Y S.Numericalmodeling of large diameter steel piles under monotonic and cyclic horizontal loading[C]// International Symposium on Numerical Models in Geomechanics(NUMOGX).Rhodes,Greece:Taylor& Francis,2007:453-459.

[30]Cox W R,Reese L C,Grubbs B R.Field testing of laterally loaded piles in sand[C]//Proceedings of the 6th Annual Offshore Technology Conference.Houston,USA:OTC Committee,1974: 459-472.

[31]Reese L C,CoxW R,Koop FD.Analysis of laterally loaded piles in sand[C]//Proceedings of the 6th Annual Offshore Technology Conference.Houston,USA:OTC Committee,1974:473-483.

[32]Achmus M,Abdel-Rahmen K,Peralta P.Analysis of load-bearing behavior of offshore wind turbine foundations[C]//Symposium of Piles 2005.Hannover,German:University of Hannover,2005: 137-158.

[33]Lesny K,Richwien W,Hinz P.Assessment of the offshore wind turbine foundations[C]//Symposium of offshore wind energy,the aspects on the construction and environments.Hannover German: University of Hannover,2007:175-183.

[34]Juirnarongrit T,Ashford SA.Effect of pile diameter on themodulus of subgrade reaction[D].San Diego:University of California,2005.

[35]Clough R W,Penzien J.Dynamics of structures[M].New York: McGraw-Hill,1993:15-17.

[36]Prowell I,Veletzos M,Elgamal A,et al.Experimental and numerical seismic response of a65kW Wind turbine[J].Journal of Earthquake Engineering,2009,13(8):1172-1190.

[37]范洪軍,金全州,劉鐵英,等.風力機地震響應分析的研究現狀與展望[J].結構工程師,2010,26(6):155-163.

猜你喜歡
風速結構模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
基于Kmeans-VMD-LSTM的短期風速預測
基于最優TS評分和頻率匹配的江蘇近海風速訂正
海洋通報(2020年5期)2021-01-14 09:26:54
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
3D打印中的模型分割與打包
基于GARCH的短時風速預測方法
主站蜘蛛池模板: 一本大道AV人久久综合| 国产精品性| 色综合天天视频在线观看| 午夜日韩久久影院| 国产女人水多毛片18| 第一区免费在线观看| 亚洲首页国产精品丝袜| 狠狠操夜夜爽| 亚洲AV无码乱码在线观看代蜜桃| 午夜在线不卡| 好紧好深好大乳无码中文字幕| 免费一级大毛片a一观看不卡| 午夜视频在线观看免费网站 | 亚洲专区一区二区在线观看| 成人午夜福利视频| 香蕉国产精品视频| 亚洲人成在线精品| 免费啪啪网址| 国产高清精品在线91| 国产区在线看| 国产无吗一区二区三区在线欢| 99在线视频免费| 国产哺乳奶水91在线播放| 欧美天堂在线| 波多野结衣一区二区三区AV| 久久精品人人做人人爽97| 欧美第一页在线| 久久久噜噜噜久久中文字幕色伊伊| 青青青伊人色综合久久| 2021无码专区人妻系列日韩| 一区二区自拍| www.youjizz.com久久| 久久99精品国产麻豆宅宅| 一级福利视频| 黄色网站在线观看无码| 亚洲精品不卡午夜精品| 18禁色诱爆乳网站| 日韩精品少妇无码受不了| 欧美一道本| 亚洲国产成人超福利久久精品| 又污又黄又无遮挡网站| 久久不卡国产精品无码| 国产草草影院18成年视频| 欧美一级高清免费a| 99激情网| 成人无码一区二区三区视频在线观看 | 一级片一区| 欧美天堂在线| 亚洲V日韩V无码一区二区| 第一区免费在线观看| 亚洲欧美极品| 国产精品视屏| 久久精品中文字幕少妇| 毛片网站免费在线观看| 国产00高中生在线播放| 国模沟沟一区二区三区| 亚洲va欧美ⅴa国产va影院| 97se亚洲综合在线韩国专区福利| 欧美精品v欧洲精品| 99草精品视频| 欧美va亚洲va香蕉在线| 国产主播一区二区三区| 亚洲三级成人| 国产丝袜精品| 激情无码字幕综合| 国产成人毛片| 免费黄色国产视频| 亚洲成A人V欧美综合天堂| 亚洲精品国偷自产在线91正片| 欧美成人综合在线| 国产毛片基地| 成人福利免费在线观看| 国产美女无遮挡免费视频| 无码精油按摩潮喷在线播放| 在线观看国产精品一区| 亚洲人成网7777777国产| 在线观看视频一区二区| 天天干伊人| 人人爽人人爽人人片| 丰满人妻一区二区三区视频| AV无码一区二区三区四区| 日韩精品高清自在线|