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

冰載荷仿真及其在海上風機一體化分析中的應用

2021-12-26 13:00:00王凱劉天輝萬遠琛施偉
哈爾濱工程大學學報 2021年11期
關鍵詞:海冰有限元振動

王凱, 劉天輝, 萬遠琛, 施偉

(1.中山大學 海洋工程與技術學院,廣東 珠海 519082; 2.南方海洋科學與工程廣東省實驗室(珠海),廣東 珠海 519082; 3.大連理工大學 海岸和近海工程國家重點實驗室, 遼寧 大連 116024)

面對日益嚴峻的環境污染、溫室氣體排放和氣候變化等問題,許多國家都在調整能源發展戰略,煤電投資熱度持續降溫,加快可再生能源的開發利用已成為目前世界各國的普遍共識。新世紀以來,風能發電在世界范圍內迅速發展,并在可再生能源發電領域中占據重要地位。風能發電可以按照其安裝的位置分為陸上風電和海上風電。相對陸上風電, 海上風電在風資源裕度、風電場位置、風力機裝機容量等方面具有優勢, 在風力機制造安裝復雜度、全壽命周期成本、運行維護難度、對電網調峰裕度的要求等方面存在劣勢。綜合考慮后海上風電的利大于弊,具有很大開發價值,海上風電是世界眾多國家推動能源低碳轉型的核心路徑。

海上風力發電在高緯度寒區海域,例如歐洲波羅的海、亞洲渤海灣、美洲五大湖等區域有很大的發展潛力,對于寒區海上風機而言海冰載荷是重要載荷之一。風能探測的報告顯示我國渤海及黃海北部海域風能密度大,然而冬季會有海冰影響風電設備。大量工程項目表明, 海冰會對海洋結構產生巨大的危害。1969年渤海發生特大冰封, 巨大的海冰直接將海上的“海二井”石油平臺推倒, 造成了建國以來最大的由于海冰導致的石油平臺事故。放眼全球, 阿拉斯加庫克灣的采油平臺[1]、北歐的Bothnia灣燈塔[2]等冰區海洋結構都不同程度地遭到過海冰的破壞。因此可以說, 海冰是關乎海洋結構安全的全球性問題。

1 我國冰載荷的研究現狀

目前我國關于冰載荷的研究發展尚不成熟,而且大多集中在海洋平臺領域,海冰是高緯度海域海洋結構物設計和安全運營的重要影響因素,研究表明海冰載荷往往是寒區海域影響海洋工程結構物穩定的控制載荷[3]。早在1966年,我國就開始在渤海海域建設抗冰的海上平臺,至今已成功建成了數十座具有抗冰能力的海上平臺。

冰載荷的研究只在少數的幾個重視寒區海洋能源開發的國家開展,包括美國、北歐、加拿大、俄羅斯、日本、中國等,相對于波浪載荷、風載荷、地震載荷等問題的研究,冰載荷問題的工程需要相比之下沒那么廣泛,而且各個國家海域的冰情、采用的結構形式等均有著較大的差異,研究重點也根據實際海域有所不同,造成彼此研究的相對獨立。我國的冰載荷研究主要圍繞渤海展開,目前的研究主要分布在以下幾個方面:

1)渤海的冰情調查,包括歷史冰情記錄[4-5],冰情預報方法[6]統計分析劃分海冰區域以及給定工程設計參數[7];

2)海冰的物理力學性能[8](包括溫度、鹽度、密度等物理參數和壓縮強度、彎曲強度等力學參數);

3)冰力及其檢測方法[9];

4)冰與結構的動力相互作用和結構的冰振響應[10];

5)結構的冰激疲勞分析以及抗冰設計[11]。

2 冰載荷數值仿真方法

2.1 數值仿真概況

針對冰載荷研究主要方式有現場實測、原型試驗、理論分析和數值模擬。對冰載荷的研究主要針對靜冰載荷和動冰載荷。由于海洋環境復雜,影響海冰運動的因素較多,因此如何建立合理有效的數值方法來模擬海冰動力學過程一直都是研究的重點,雖然海冰與結構物相互作用的過程至今沒有公認的數學物理模型,不過研究者們還是基于模型實驗和現場觀測為基礎,用數值模擬方法展現冰和結構相互作用的過程,特別是對海冰動態破壞的過程。

在分析冰與一般的結構物作用時,通常所研究的冰尺度較小。對于較小的尺度,若假設其為均勻介質就可以很簡單地運用目前最為常用的數值分析方法有限元法[17-21]。除此以外還有離散元法[22-23]、非連續變形分析方法[15-16]等其他數值分析方法在處理不同問題時有著各自優勢。

2.2 有限元法

有限元分析是解決工程和數學模型時最為廣泛使用的方法。將大型的系統細分為小而簡單的元,在滿足工程要求的前提下設置元素,將求解域劃分為離散域,再根據離散域各自的方程進行求解,就可以用有限數量的未知量代數方程組去求解偏微分方程組。對于不同數學模型的問題有限元的求解方法是大致相同的,其核心思想均為“數值近似”和“離散化”,求解問題的基本步驟[12]通常是根據問題確定好計算域,然后將計算域離散化,再確定狀態變量的方程組和設定邊界條件開始進行計算,計算得出結果后再進行分析。處理步驟大概可分成3個階段,前處理、求解和后處理流程如圖1所示。前處理就是根據研究的問題建立有限元模型并且根據計算需求完成單元網格劃分;后處理則是在得到計算結果以后根據計算的設定條件分析結果,以方便其他研究者了解。有限元方法是基于連續介質理論用來描述連續結構體運動應用最廣泛的數值計算方法,適合用來分析海洋結構在冰力作用下動力響應。

2.3 離散元法

離散元法是專門用來解決不連續介質問題的數值模擬方法最初往往應用于巖石力學領域,把介質離散為元,相鄰的元之間存在某種或幾種作用力,元的運動滿足牛頓第二定律。近年來,基于離散元方法模擬海冰與船體和海洋平臺結構間的相互作用研究取得了較大進展。冰離散元方法的計算單元可采用球體、圓盤、多面體和擴展多面體等不同形態。采用離散單元模型計算更加適合模擬離散顆粒組合體在準靜態或是動態條件下的變形和破碎過程,而海冰是屬于非連續介質可以離散為具有一定質量的顆粒單元,顆粒之間允許位移重疊或分離,通過這些單元可以靈活構建不同類型的海冰,包括平整冰、浮冰、冰脊和堆積冰等。

圖1 有限元分析示意Fig.1 Finite element analysis diagram

2.4 其他數值方法

除了前述2種主流數值分析方法之外,還有許多數值方法用于冰載荷的研究,如差異元素法[13]、三維顆粒流法[13]、幾何網格法[14]等。石根華[15]于 1988 年提出了一種新型數值分析方法將數學拓撲學和工程實踐相結合的非連續變形分析(discontinuous deformation analysis,DDA)。該方法最初是為了解決巖體大變形和大位移問題時介質不連續性導致計算誤差大的問題而提出的,平行于有限元法,充分考慮結構介質的不連續性,將結構面切割而成的塊體作為分析單元,將動力學與靜力學統一起來,用最小勢能原理把塊體之間的接觸問題和塊體本身的變形問題統一到矩陣中求解,后來被學者運用到冰與結構的相互作用模擬中。張運良等[16]就基于DDA算法模擬冰與結構相互作用的過程,將冰破碎的過程可視化,計算了冰撞擊結構后每個時間步的結構動態響應。模擬出的冰破碎過程與原型觀測到的過程較為吻合,與有限元和邊界元方法相比,DDA方法對冰的破碎模擬更為接近實際。不過目前整體而言該方法在應用于冰載荷研究還相對較少。

3 數值仿真方法在抗冰研究中的應用

在研究冰與海洋相互作用時所關注的問題往往是作用在結構上的最大靜冰力和動冰力形式、海冰引起結構振動的機理以及結構在冰載荷作用下的失效形式,按照結構與冰作用部分的形狀可以分為直立結構基礎及非直立結構基礎。現在數值仿真在冰載荷方面的應用往往是通過與模型試驗數據和現場測試數據對比以驗證數值仿真方法的適用性,積累冰載荷領域數值仿真的經驗,并在此基礎上進一步模擬仿真不同壞境、不同結構時冰與結構相互作用過程,為未來的抗冰結構設計奠定基礎。

3.1 在直立結構中的應用

冰與直立結構作用時,可能出現劈裂、彎曲、屈曲和擠壓等多種破壞形式,其中以擠壓破壞為主。岳前進等[17-18]、K?rn?等[19]、李洪升等[20]對冰與直立結構的作用過程及結構的振動進行了分析,發現冰與直立結構作用主要發生屈曲破壞和擠壓破壞,并且不同速度的海冰與結構作用擠壓破壞時會形成如圖2示的3種不同形式的冰激振動:準靜態振動、穩態振動、隨機振動,海冰的擠壓破碎也因為冰速的不同而表現為3種不同的破碎模式。當冰速很慢時,出現準靜態冰力,對應冰的破碎模式為準靜態韌性破壞;隨著冰速的增加會出現導致結構穩態振動的冰力,此時的破壞模式為韌脆轉換破壞;隨著冰速的繼續增加會產生導致結構隨機振動的冰力,對應破碎模式為脆性破壞。

(a)極限應力準則 (b)極限動量準則 (c)極限力準則圖2 直立結構最大靜冰力的3種極限準則Fig.2 Three limit rules for calculating peak static ice force on vertical structures

在靜冰力領域,目前已有相應的工程規范,楊耀鵬等[21]以新建的海上平臺為例,通過有限元軟件分析其抗冰性能,發現結構設計靜力的安全儲備比較大,不過冰振對結構造成的潛在影響還有待進一步考究,可能會造成結構的疲勞失效。

在動冰力的研究領域,季順迎等[22]通過離散元法模擬計算了海冰與直立海洋平臺相互作用的過程,模擬海冰對不同樁徑的直立結構的動冰力過程,確定海冰的破碎規律以及冰載荷的特性,發現平均冰力與樁徑呈正比。賈賓等[23]在前者研究的基礎上,采用基于近場動力學理論的數值算法,數值模擬的主要參數參考了前者的模型,模擬了海冰與柱狀結構的作用過程,并考慮了相互作用時產生的結構振動,探究了樁徑、冰速等參數對極值靜冰力和結構振動位移的影響,冰力部分模擬結果與前者相近,可見這種新型的數值算法也是具有合理性的。

冰激振動是影響海洋結構安全以及導致其損傷和疲勞壽命的重要因素,也是冰載荷研究領域的熱點。張大勇等[24]在分析冰激振動下直立結構的疲勞壽命時,借助有限元數值分析方法,主要考慮了穩態振動和隨機振動造成的疲勞損傷影響,選取安全壽命設計方法[25]進行冰激疲勞估計,明確冰與直立結構相互作用過程,提出一套冰激直立抗冰平臺的疲勞壽命計算流程。

在分析冰載荷問題時往往把海洋平臺視為剛體結構,而在分析海洋結構的動力特性時則把冰載荷簡化成冰力函數處理,隨著冰激振動研究的深入,需要進一步研究海冰與海洋結構的動力耦合作用,不少研究者提出了建立離散元-有限元的耦合模型來分析冰激結構振動的內在機理。王帥霖[26]采用基于區域分解方法的離散元-有限元耦合模型,對單樁直立結構的冰激振動進行了仿真模擬,由梁-殼組合單元構成海洋結構的有限元模型,在耦合界面上發展了相應的接觸算法和耦合參數傳遞算法,考慮了冰速、厚度變化的影響,在冰況相近的前提下,仿真計算出冰載荷、冰激振動加速度和應力分布與現場實測數據趨勢相同,初步說明了該耦合模型在分析冰與結構相互作用過程的適用性。

3.2 在非直立結構中的應用

近年來,錐體結構在海洋工程領域不斷被采用,不少國家為冰區的海洋平臺加裝了圖3所示的錐體結構,國內外的研究人員針對冰排與固定傾斜結構作用時受到的冰載荷做了大量的理論、實驗和數值研究,通過海冰的物理力學特性發現,冰的彎曲強度遠小于擠壓強度。錐體取代直立結構被公認為降低冰力的解決辦法。

圖3 抗冰單樁基礎[24]Fig.3 Ice resistant single pile foundation[24]

對于錐體結構而言,由于冰排的彎曲強度遠小于擠壓強度,而在根據現在設計標準,即使是直立結構也有很大的靜力安全儲備,而錐體則不言而喻,如圖4所示,對于錐體而言,雖然冰排破壞的力大幅減小了,不過會產生冰與斜面的摩擦力所帶來的額外載荷,因此研究椎體結構靜冰力時不僅需要考慮錐體本身的物理參數還需要考慮冰爬坡力所帶來的額外載荷。

圖4 海冰與錐體碰撞示意[42]Fig.4 The sketch of collision between ice and cone[42]

Sand等[27]和Derradji-aouat[28],利用有限元軟件ANSYS計算了冰排與斜坡結構作用時結構受到的冰力,將冰排理想化成作用在彈塑性基礎上的線彈性板,將得到的結果與前人基于實驗數據的理論公式進行比對。龔榆峰等[29]在各項同性海冰失效準則的基礎上結合Abaqus有限元軟件,開發用于模擬海冰失效的新型單元并且結合了Ralston[30]給出的塑性極限分析解進行比較,直觀地給出冰損傷過程及對應結構受力的變化過程。Bjrnar等[27]使用非線性有限元方法模擬冰對固定圓錐體的作用,并分析錐角及摩擦力對冰荷載的影響。利用一種特殊的接觸算法模擬了冰原與錐狀結構之間的相互作用,使跟蹤冰原與結構之間逐漸發展的接觸成為可能。在相互作用過程中,由于冰原形態的變化,浮力也隨之變化。通過引入一個連續的非線性地基模型來跟蹤這一過程,該模型包括了浮力和冰的比重的影響。用2種不同的本構模型來近似冰的力學行為。將材料非線性和冰與結構摩擦的影響也納入有限元計算當中,精確追蹤冰蓋與斜坡結構之間的接觸,包括摩擦,計算部分或全部被淹沒的冰蓋上的浮力。探究了冰與結構之間錐角和摩擦系數的變化對結構的影響,并將數值計算的結果與基于塑性極限分析的分析結果進行了比較,發現兩者非常接近,可見該有限元模型的可靠性。

錐體結構對減小動冰力以及減小冰激振動方面相比于直立結構表現出色,前文已經提及,冰激振動是目前寒區海洋工程抗冰結構設計的面臨的重大課題,不過目前仍沒有給出計算冰激結構振動的冰力函數,對海冰引起結構振動的力學機制的解釋還沒統一,研究者們也在不斷探究適合計算海冰動力學的數值方法。王剛等[31]利用有限元軟件,對不同接觸寬度的冰與錐體相互作用過程進行了分析,對徑向開裂型和環向開裂型2種冰排彎曲斷裂形式進行數值仿真計算,研究了接觸寬度與冰排彎曲斷裂形式的關系。李海等[32]綜合考慮了計算的效率和精確性,發展了一種將拉格朗日坐標與歐拉坐標相耦合的方法,將有限差分法和質點流體力學方法相結合,考慮了海冰的堆積以及渦旋風場的作用,對遼東灣海冰動力過程進行了72 h的數值模擬計算,并且與衛星的遙感圖像進行比較,發現模擬的海冰的厚度、密集度分布以及冰速與衛星遙感圖像有良好的一致性,可以較好地模擬出渤海海冰的動力學過程。

3.3 數值仿真結果的驗證

目前對錐體海洋結構與海冰的相互作用主要通過現場實測、室內試驗以及數值模擬展開研究,由于冰-錐相互作用過程的復雜性,很難通過現場實測和試驗分析不同參數對該過程的影響,而數值模擬方法可以準確修改參數,再現出實測情況的冰-錐相互作用過程,并在確定數值仿真方法可行性的基礎上進一步改變參數模擬各種情況,從而達到冰激振動預警分析、疲勞分析以及動力優化設計的目的。以上數值仿真工作大多基于已有的現場實測或室內試驗數據仿真研究以確定數值仿真工作的合理性。關湃等[33]基于模型試驗結果,借助有限元軟件對冰錐結構進行優化設計,模擬不同冰況下結構的響應,不斷改變結構的各種形狀參數以達到結構優化的目的。也有不少研究者專門設計試驗以驗證數值分析的結果是否合理,汪震宇[34]用有限元方法分析斜坡與冰的相互過程,并設計了相應的模型試驗進行誤差分析以及模型的校正。

研究冰載荷的各種研究方法都是相互驗證,不斷改進的過程,而隨著數值分析方法不斷發展成熟,而且能夠精確地改變各種參數以適應需求,對理解冰與結構相互作用機制具有成本低效率高的優勢。

4 在海上風機一體化數值分析的應用

4.1 海上風機一體化分析需求

海上風機也作為一種海洋結構,在冰載荷方面的研究與前述的直立結構和錐面結構有諸多相似之處,在這便不再贅述了。海上風機相比于其他海洋結構,最大的特點在于風機為典型的柔性結構,海冰會引起柔性結構明顯的動力放大,冰振響應十分顯著,對結構及上部設備造成一定影響[35]。張大勇等[36]基于有限元分析法研究動冰載荷對結構的影響,選取單立柱式三樁形式的風電基礎并與相同結構形式的導管架平臺進行動力特性對比,文章中歸納總結了冰載荷作用下風電基礎的各種失效模式,包括靜冰力作用下的結構失效和動冰力作用下的疲勞失效、振動失效和土壤失效,發現現有的風電基礎設計雖能滿足極值冰載荷的要求,但在動冰力作用下對結構造成的影響還有待考量。通過有限元法仿真模擬可以得到驗證風電基礎結構主要是塔架在載荷下的響應,但并沒有模擬出動冰力振動下對內部設備如電機、限速安全機構等的影響,但可以通過結構的冰振位移、速度等指標來間接地考慮動冰力下結構振動對內部設備的影響程度,為寒區風電基礎的抗冰設計的改良、相關的實驗設計及安全保障提供合理依據。

4.2 常用海上風機一體化分析軟件及應用

目前用于海上風機一體化的軟件主要為FAST軟件,作為風機耦合的軟件,不僅能夠建立基礎結構的模型還能搭建旋轉葉片-機艙-塔架組件的詳細模型,更加貼近實際的模擬仿真風機在各種載荷作用下的響應。FAST軟件[37]是美國可再生能源實驗室(National Renewable Energy Laboratory,NREL)開發,用于計算氣動彈性及結構動力學仿真的開源軟件。該軟件的仿真流程圖如圖5所示,FAST作為一款開源的軟件,為拓展研究、增加功能提供了基礎,部分研究者在軟件基礎上進行二次開發以適應不同的研究需求。

現代風力機作為最大的旋轉機械,扇葉長度甚至可以超過最大的飛機翼展,風載荷對風機的影響是其他海洋結構無法比擬的,所以分析風和海冰的聯合載荷與風機作用的過程十分重要。許子非等[38]考慮了海冰與湍流風載荷的聯合作用,利用FAST軟件,以近海單樁風力機為研究對象,數值計算在不同厚度及海冰移速的海冰碰撞條件下塔架的動力學響應特性,并對比了有無抗冰錐時響應的不同,發現抗冰結構能夠明顯減小結構的疲勞載荷。

除了FAST以外,目前用上風機一體化數值分析主流的軟件還有HAWC2[39](horizontal axis wind turbine simulation code 2nd generation)和Bladed,都是被一些風力機廠商所認證的仿真軟件,作為風機仿真軟件應用于很多研究項目和工程當中,能夠模擬許多現有的風力渦輪機。

HAWC2能夠綜合考慮空氣動力學載荷、水動力載荷和土力學載荷,在時域內使用線性和非線性的方法結合以計算結構動力響應,基于多體公式的非線性氣動彈性模型,可以處理復雜的結構。Shi等[40]利用了HAWC2代碼,建立了半經驗的全耦合冰載荷模型,研究冰載荷對安裝了抗冰結構的海上風電機的影響,同時考慮了風載荷和冰載荷的基礎下進行了耦合分析并且與非耦合分析的結果進行了比較,討論了耦合分析的重要性,并且研究了風速、平均冰厚、漂移速度不同時,冰的物理性質和破冰錐幾何形狀的對海上風機響應的影響。

圖5 FAST軟件仿真流程圖[38]Fig.5 Software simulation flow chart[38]

除了冰與結構的相互作用之外,海上風力機仍需考慮結冰對風力機的影響,結冰會引起額外的載荷和振動,不同葉片上冰載荷的不平衡也會在一定程度上增加機組部件的疲勞。易賢等[41]采用數值計算的方法進行結冰計算,模擬結冰形成的動態過程,分析結冰的主要部位以及厚度,還有對風機氣動特性的影響,為設計結冰探測方法和除冰方案奠定基礎。

風力機作為高聳結構,相比于其他海洋結構柔性特性更為明顯,也意味著風力機受冰激振動的影響更為顯著,而目前對于冰激振動的研究仍處于探索階段,離徹底揭示冰激振動機理和抗冰振結構設計的完善還有很長的路要走,隨著數值仿真技術的成熟,數值計算方法相比于原型觀測和室內試驗研究的進程更短,花費小效率高,能規避很多不必要的風險。

5 結束語

目前冰載荷研究的靜冰力領域發展相對成熟,各個研究地區均有明確的規范,不過仍未形成統一計算公式或是數學模型,靜冰力分析理論還有待進一步發展。新型的海洋結構或是抗冰結構對于靜冰力均有較大的安全儲備,應當探索降低靜冰力工程設計值以提高經濟效益的可能性;在動冰力和冰激振動領域,仍需探索可以用于工程的動冰力函數,進一步優化抗冰結構的減震性能,形成相關的工程規范,并且發展數值分析方法以適應動冰力和冰激振動數值計算需求。

理論分析和數值模型的發展必然離不開工程實例與試驗,目前我國僅天津大學有冰池實驗室以及哈爾濱工程大學建設的室外冰水池實驗室,但與德國漢堡水池相比仍有著顯著的差距,有待進一步改善實驗條件,開展更多冰力實驗。

當然,冰載荷領域的研究要想應用到海上風機還需要結合其他領域知識,諸如空氣動力學、水動力學、土力學等領域的研究,并結合現代的數值仿真技術,不斷發展一體化數值分析技術以貼近實際。

相信隨著數值模擬技術與計算機的蓬勃發展,以及目前各國對冰載荷愈加重視,冰載荷的數值仿真方法也會越來越成熟,同時也會有更多的文獻方便學習。我國渤海海域能源豐富,不論是傳統能源還是新能源的開采都避不開樁式結構,而渤海易受海冰影響,因此我國在抗冰方面的步伐從未停止,特別是近年來各地興起了海上風電的熱潮。通將過數值方法與模型試驗、現場觀測等手段結合不斷發展冰力研究,并在工程實踐中不斷積累相關經驗,相信在不久之后便會形成一套完善海上風機的冰力規范。

猜你喜歡
海冰有限元振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
末次盛冰期以來巴倫支海-喀拉海古海洋環境及海冰研究進展
海洋通報(2021年3期)2021-08-14 02:20:38
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
中立型Emden-Fowler微分方程的振動性
基于SIFT-SVM的北冰洋海冰識別研究
磨削淬硬殘余應力的有限元分析
應用MODIS數據監測河北省近海海域海冰
河北遙感(2014年4期)2014-07-10 13:54:59
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
基于SolidWorks的吸嘴支撐臂有限元分析
基于TerraSAR-X全極化數據的北極地區海冰信息提取
主站蜘蛛池模板: 日韩精品一区二区三区视频免费看| 老司国产精品视频91| 国产av一码二码三码无码| 欧美视频在线播放观看免费福利资源| 22sihu国产精品视频影视资讯| 99精品伊人久久久大香线蕉| 国产美女主播一级成人毛片| 91视频区| 二级特黄绝大片免费视频大片| 欧美一区二区福利视频| 手机精品福利在线观看| 国产系列在线| 国产对白刺激真实精品91| 国产欧美日韩专区发布| 国产无人区一区二区三区| 日韩第九页| 中文字幕无码中文字幕有码在线| 三级毛片在线播放| 99在线免费播放| 亚洲精品视频在线观看视频| 久久永久视频| 日韩成人在线视频| 亚洲欧美另类日本| 蝴蝶伊人久久中文娱乐网| 午夜高清国产拍精品| 国产va在线| 在线播放精品一区二区啪视频| 国产成在线观看免费视频| 另类欧美日韩| 中文字幕免费在线视频| 日韩专区第一页| 成人在线不卡| 真实国产乱子伦高清| 亚洲熟妇AV日韩熟妇在线| 69av在线| 亚洲人成成无码网WWW| 亚洲色无码专线精品观看| 天天综合网色中文字幕| 国产日韩欧美成人| 欧美在线视频a| 91亚洲国产视频| 日韩一区二区三免费高清| 亚洲永久免费网站| 色婷婷色丁香| 亚洲中文字幕23页在线| 欧美中文字幕无线码视频| 99视频在线免费| 99视频国产精品| 亚洲美女视频一区| 不卡无码h在线观看| 午夜a级毛片| 国产精品自拍露脸视频 | 午夜国产精品视频黄| 亚洲国产一成久久精品国产成人综合| 国产精品爆乳99久久| 久久精品国产在热久久2019 | 天堂成人在线视频| 九九久久99精品| 成人夜夜嗨| av在线无码浏览| 日韩福利视频导航| 亚洲久悠悠色悠在线播放| 欲色天天综合网| 国产午夜人做人免费视频中文| 亚洲区视频在线观看| 中文字幕1区2区| 玖玖精品在线| 亚洲国产成人精品青青草原| 2021国产乱人伦在线播放 | 天天激情综合| 国产欧美成人不卡视频| 69精品在线观看| 伊人天堂网| 国产精品爽爽va在线无码观看 | 中文字幕资源站| 亚洲a级在线观看| 亚洲欧美日韩成人高清在线一区| 欧美丝袜高跟鞋一区二区| 伊人久久影视| 色噜噜狠狠狠综合曰曰曰| 91视频日本| 丰满人妻久久中文字幕|