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

整體式止屈器數值模擬方法與優化研究

2020-12-16 07:48:52顏鎧陽余建星韓夢雪許偉澎張春迎
海洋工程 2020年6期
關鍵詞:模型

顏鎧陽,余建星,余 楊,韓夢雪,許偉澎,張春迎

(1. 天津大學 水利工程仿真與安全國家重點實驗室,天津 300072; 2. 天津大學 天津市港口與海洋工程重點實驗室,天津 300072)

海底管道是海洋油氣系統不可或缺的部分,被譽為海洋油氣系統中的“生命線”,海底管道運輸是一種安全性高、經濟性好以及對環境危害小的運輸方式。管道長期處于高壓環境下,所處外部環境十分復雜,即使在鋪設過程中沒有發生任何損壞,但在實際服役期間,管道可能會由于受到其他物體的碰撞,產生局部的凹坑缺陷,或是由于磨損、腐蝕產生局部壁厚的減小[1]。上述各種情況都會導致管道局部上抗屈曲性能的減弱,管道在這些局部缺陷處極易發生屈曲壓潰。這種局部上的壓潰變形會導致相鄰區域的管道也發生壓潰,使壓潰變形沿管長方向快速傳播,引發海底管道大范圍的破壞,這種現象稱為屈曲傳播[2]。對于深水海底管道,一旦發生屈曲傳播,傳播的速度可以達到數百米每秒,造成的損失無法估量,如海底管道壓潰后破裂會造成石油泄漏事故,經濟上和生態環境上都會是巨大的災難[3]。

海底管道的屈曲傳播特性會導致管道大面積的破壞,帶來巨大的經濟損失。增大管道的壁厚可以加強管道的抗壓能力,使管道不那么容易壓潰[4]。但是,增大管道的壁厚意味著需要增加原材料的使用量,增加了生產成本;從鋪設管道角度來看,管道懸跨段的重量會隨著管道壁厚的增大而增加,重量增加可能使原本使用的鋪管船的張緊器能力不足,從而需要對鋪管船的張緊器進行升級,進一步增加了成本[5]。因此,單憑增大管道的壁厚去解決屈曲傳播導致的管道破壞問題在經濟上不太可行。另一方面,若是等到屈曲傳播發生后再去修復或者更換破壞的管道顯然也是不可行的,因為修復或者更換海底管道的成本是巨大的[6]。

解決海底管道屈曲傳播問題的一個既實用又更具經濟性的方法是沿著管長方向每隔一段距離設置止屈器。止屈器的抗屈曲能力更強,如果管道發生局部壓潰并產生屈曲傳播,在屈曲傳播到止屈器處時,止屈器可以阻止屈曲傳播跨越止屈器,能使屈曲傳播僅發生到止屈器處,而不再繼續傳播下去。止屈器的局部加強作用能夠限制海底管道壓潰后發生屈曲傳播的范圍,使屈曲傳播產生的破壞止于止屈器處,保證了海底管道整體的安全性。當海底管道僅在局部發生破壞而其余部分仍完好時,管道整體的安全性不受影響。總的來說,雖然止屈器允許管道發生一定范圍內的屈曲,但極大地限制了海底管道破壞失效的范圍,很好地保持了管道的整體性。事實證明,止屈器的設置降低了維修管道的難度和所需時間,同時有效減少了經濟損失,提高了海底管道的服役年限;另一方面,裝置止屈器并不困難,成本不高,在安全性和經濟性間取得了一個比較良好的平衡[7]。

按照安裝形式的不同,海底管道止屈器的主要類型包括:扣入式止屈器、纏繞式止屈器、焊接式止屈器和整體式止屈器。整體式止屈器在工程中應用比較廣泛[8]。

整體式止屈器是在管道鋪設之前,按一定布置間距焊接在管道指定位置處的厚壁圓環,所以止屈器的內徑和管道的內徑相同,只是止屈器的壁厚比管道的壁厚大。止屈器和管道同軸排好之后在連接處進行焊接,如圖 1所示,止屈器和管道形成一個整體,黑色三角為焊接處。

圖1 整體式止屈器示意Fig. 1 Integral buckle arrestors

止屈器可以阻止屈曲的傳播,但是當外部水壓高于一定的數值時,屈曲傳播還是會“穿越”止屈器繼續向下游傳播[9],屈曲能夠穿越止屈器繼續向下游傳播時對應的外部水壓壓力值稱為止屈器的穿越壓力Px。

關于止屈器的止屈效果和穿越壓力,國內外學者進行了諸多研究。Park等[10]進行了整體式止屈器的穿越試驗,分析止屈效果與長度、厚度等變量之間的關系。Mansour等[11]采用有限元方法對整體式止屈器進行了準靜態和動態試驗條件下穿越壓力的計算,并與試驗結果進行對比。Toscano等[12]通過試驗和數值模擬的方法,研究了整體式止屈器的穿越模式,取得了較為一致的結果,這表明有限元模型可以對止屈器性能進行較好的預測。Dyau等[13]運用模型試驗和數值模擬的方法,研究了扣入式止屈器和整體式止屈器的工作性能,分析了止屈器的長度、壁厚和材料性能等參數對止屈器穿越壓力的影響。Kamil等[14]使用有限元分析軟件ANSYS,對安裝整體式止屈器的海底管道進行仿真模擬,分析了管道壁厚、徑厚比、止屈器長度以及管道外徑4種幾何要素的變化對整體式止屈器止屈效果的影響。

進行模型試驗能夠比較真實地反映研究對象發生作用的全過程,得到的結果比較準確,但是進行模型試驗所需的材料成本和時間成本大,每次研究都進行模型試驗是不現實的,特別是整體式止屈器試驗過程復雜,無法進行大量模型試驗。而隨著科學技術的發展以及有限元分析方法的完善,現在可以使用各種有限元分析軟件進行仿真模擬,并且仿真模擬的結果具有較強的科學性和準確性,計算得出的數值比較接近實際情況,有限元數值模擬已經成為一種常用的研究方法。因此,利用ABAQUS軟件建立相應的有限元模型,并且改變管道和止屈器的重要設計參數,進行系列計算,以便進行整體式止屈器的優化。

1 計算方法

1.1 廣義弧長法

深水海底管道的屈曲傳播與止屈涉及較復雜的非線性問題,主要包括:①管道發生屈曲傳播時,鋼材已經進入塑性變形階段,這屬于材料的非線性問題;②在管道壓潰繼而發生屈曲傳播的過程中,管道的截面會發生大變形,這屬于幾何非線性的問題;③在壓力作用下,管道內壁最終會相互接觸,接觸分析屬于典型的邊界條件非線性問題。

上述3個非線性問題導致整個管道止屈器系統發生極值失穩后的荷載-位移路徑非常復雜,必須選擇合適的方式控制增量步長,才能追蹤荷載-位移的變化路徑。非線性分析中一種比較常用的方法是廣義弧長法,即通過在ABAQUS中step模塊下選擇Static,Riks類型的計算方法,設置相關參數進行求解。

這種計算方法是將施加的荷載也作為一個變量來求解,所以它的荷載增量步是變化的,但這相當于給原來的求解方程組增加了一個未知量,要求解方程組就需要再補充一個控制方程。這個控制方程形式上是一個圓的方程,它的圓心是上一個荷載步的收斂點,半徑則與上一個荷載步得到的荷載和位移有關。有了這個補充方程,就可以繼續往下求解。因為補充的方程是一個圓方程,說明使用這種計算方法的迭代路徑是一個圓弧,所以稱為弧長法;另一方面,該方法需要同時求解荷載和位移,屬于一種廣義的位移控制法,所以稱為廣義弧長法。

使用這種計算方法具有操作簡單易上手的優點,只需要在ABAQUS中step模塊下設置相關參數即可。下面結合ABAQUS操作界面和外部INP文件解釋相關參數。

INP文件中該部分相關的代碼如下。

**STEP:Step-1

**

*Step,name=Step-1,nlgeom=YES,inc=

*Static,Riks

Δlinit,lperiod, Δlmin, Δlmax,λend, node, dof,μmax

*END STEP

在ABAQUS中創建分析步,設置分析步名稱,即“Step, name=Step-1”;類型選擇“Static,Riks”;設置Nlgeom的參數為on,即“nlgeom=YES”;“inc=”后的數值表示進行分析的步數的最大值;“Static,Riks”后是數據行,并不要求填寫列示的所有參數,這些參數的定義分別為:

Δlinit和lperiod兩個參數定義分析步中與求解方程相關的初始弧長和總弧長;

Δlmin和Δlmax兩個參數組成弧長增量的上下界限;

λend、node、dof和μmax四個參數用于設置分析的終止條件,λend表示荷載超出某一數值則終止分析,node, dof,μmax表示設置一點的位移分量超出某一數值則終止分析。

可以看出,初始的操作設置十分簡單,完成之后便交給計算機求解;但要使設置的弧長能夠有效地求解計算卻不容易,一個不成功的弧長設置會導致荷載-位移曲線原路返回,此時需要調整參數再次嘗試。

使用廣義弧長法進行了8次計算,其中比較典型的錯誤情況包括出現“回漂”現象以及計算到止屈器部分不收斂現象。

圖2展現了出現“回漂”現象時模型的荷載比例因子隨弧長變化的曲線。荷載比例因子是模型所受靜水壓力與模型中預設壓力值的比值,弧長與增量步數正相關。值得注意的是,曲線在出現第一個峰值點后開始下降,而在工程實際中,在產生壓潰后,屈曲開始傳播,此時壓力減小且維持穩定,直到屈曲傳播到止屈器部分,由于止屈器的結構強度更大,屈曲若想穿越,則需要更大的壓力,這個壓力便是止屈器的穿越壓力。圖 2曲線連續下降,荷載比例因子持續減小甚至出現負值,說明外部水壓力一直減小甚至從外壓變成了內壓,這顯然是不符合實際的。“回漂”現象出現在屈曲還未傳播到止屈器部分時,往前計算不能收斂,往回計算可以收斂,于是出現了屈曲往回傳播的現象,這種情況下便無法得到止屈器的穿越壓力。

另外一種情況是屈曲成功傳播到止屈器部分,計算過程中沒有產生“回漂”現象,此時管道受到的靜水壓力逐漸增大,但隨著壓力的增大,到達某個壓力值時,計算不再收斂,最后報錯,如圖3所示。查看變形云圖可知報錯的點在止屈器穿越壓力點附近。出現這種情況的原因可能是因為止屈器部分的結構比較復雜,在止屈器部分容易計算不收斂。

圖2 “回漂”現象Fig. 2 Failed cases-drift back

圖3 止屈器附近報錯Fig. 3 Failed case-aborted

總而言之,使用廣義弧長法計算管道止屈器模型,往往需要經過多次嘗試,得到適當的弧長搭配,最后才能得到需要的結果。

由于文中需要建立多組模型,使用廣義弧長法顯然會給研究帶來很大的不確定性和不便利。對此,使用另一種計算方法來解決這個麻煩的問題——靜水流體單元法。

1.2 靜水流體單元法

靜水流體單元(F3D4)又稱流體彈簧單元,是在空間上由若干個表面形成一個空腔,在指定的參考點上與受載結構外表面耦合,從而由結構變形情況確定空腔的體積變化,進而確定作用于結構上流體壓力的大小。靜水流體單元可以用來確定在結構物周圍定義的某個控制區域內的體積變化,于是將壓力控制加載轉變為體積控制加載,而壓力的大小可以通過體積的變化來計算。所以,靜水流體單元法采用的是一種特殊體積控制加載方式,這種加載方式可以很好地模擬管道模型在水壓作用下的屈曲傳播以及屈曲穿越止屈器,不會出現廣義弧長法中的“回漂”現象和計算不收斂的問題。

使用廣義弧長法計算管道止屈器模型會出現各種問題是因為分析海底管道此類結構的一個難點在于物體的變形與流體的壓力之間有耦合作用。這種耦合作用,使得結構的變形大小不只取決于外部施加的荷載,還取決于流體的壓力,而流體的壓力又會受到結構變形的影響。安裝止屈器的管道結構更加復雜,所以,不只是計算管道受到的壓力有困難,成功求得止屈器的穿越壓力更加困難。ABAQUS提供的靜水流體單元是一種能夠提供結構變形與內部流體壓力變化之間耦合作用的表面單元,它可以用來模擬被流體充滿的腔體結構,這種結構的響應不僅取決于外部的載荷,還與內部流體的壓力變化有關,而內部流體的壓力變化又會受到腔體變形的影響。靜水流體單元能提供這種結構變形與內部流體壓力變化之間的耦合作用。

靜水流體單元是一種覆蓋在腔體邊界上的表面單元,它與腔體內邊界共享節點,對所有的靜水流體單元必須定義一個共享的參考點,這個參考點是用來計算腔體體積的,從而來反映腔體內部流體的壓力,所以參考點的位置是個關鍵問題。如果腔體結構不是對稱的,參考點必須能夠完全被靜水流體單元包圍成一個封閉的區域,此時參考點的位置是任意的,不一定要在腔體內部;如果腔體結構是對稱的,通常選擇只建立一半模型以簡化建模過程,此時參考點必須位于對稱軸或對稱面上,若具有多個對稱面,參考點應位于對稱面的交點上。

靜水流體單元的參考點可以設定一個代表壓力的自由度,參考點的壓力自由度是主要的變量。因此,可以通過定義模型的邊界條件來設定參考點的自由度,這和在結構上設定節點位移的方法是類似的。為參考點設定壓力等效于在腔體邊界上作用統一的分布載荷。如果對參考點設定了壓力的自由度,還需要定義一定的流量q來表示流體注入與排出腔體,q值為正代表向腔體內流入,q值為負代表水從腔體排出。正如流體可以根據壓力的變化自動注入或排出腔體,流體的體積也會相應的自動調整以保證充滿腔體。最后通過設置兩個變量PCAV和CVOL得到流體的壓力與流體體積的歷史輸出。

具體的操作就是,管道以及止屈器模型依然使用標準的有限元單元模擬,然后在整個模型外面建立一個流體艙,接著在流體艙與管道止屈器模型之間設置靜水流體單元,這就相當于往流體艙內“注水”,直到達到預先設置的流量值。也就是說,可以控制注入流體的量,只要注入的“水”足夠多就能夠將管道完全壓壞。所以,不會出現廣義弧長法中由于弧長設置不當導致的無法計算出穿越壓力的情況,因為可以通過注入足夠多的“水”將管道壓壞,也就可以得到需要的穿越壓力值。

2 數值模擬

海底管道的屈曲傳播與止屈涉及較復雜的非線性問題,現在還難以建立有效的力學模型求得精確的解析解,因而采用ABAQUS軟件建立三維有限元模型來模擬海底管道屈曲壓潰、屈曲傳播和屈曲穿越止屈器的全過程,并計算結構臨界失效壓力。由于焊縫的影響不是文中研究的重點,并且焊接工藝的成熟使得焊接處比較穩定,因此,模擬中將管道和止屈器設置為一個整體,不考慮它們之間的邊界約束。

在均勻外部壓力作用下管道的受力變形具有對稱性,同時模型具有幾何對稱性,因此,將x=0、y=0和z=0三個平面設置為對稱面,建立1/4的管道模型,以簡化建模過程,加快建模速度,提高計算效率。

模型材料的本構關系,可通過進行材料的力學特性試驗得到。具體的方法是在管件上裁剪試驗片來做單向拉伸試驗,并利用拉伸試驗測量得到的數據,使用Ramberg-Osgood方程來對這些數據進行擬合,最終擬合得到材料屈服應力-塑性應變曲線,求得材料硬化參數n與屈服應力σy。材料屬性如表 1所示。

表1 材料屬性Tab. 1 Material properties

采用流體加載的方式模擬海底管道的屈曲發生以及傳播的整個過程,關鍵在于引入靜水流體單元。在管道外部建立一個軸向長度相同、半徑為D的流體艙,同時需要在管道止屈器模型的外壁建立一層Skin單元,這樣才可以在流體艙與管道之間連續加入靜水流體單元,模擬注入“水”以增大管道外部壓力,如圖 4所示。注意選擇Skin單元的方向,需要使管道止屈器的外壁以及流體艙的內壁處于激活(positive)狀態,流體艙的外壁處于未激活(negative)狀態,才能在其間注入靜水流體單元,以提供管道所受外部壓力,而管道內部不存在壓力。

圖4 設置Skin單元Fig. 4 Skin manager

關于網格劃分,一方面需要盡可能減少網格的數量以提高計算速度,另一方面要保證計算的精確度,因此,管道、止屈器以及流體艙的周向上劃分為15~20份,管道以及止屈器的徑向上劃分為1~2份,管道的軸向上按每D個單位長度劃分為10份,但在管道與止屈器的接觸部分進行網格的加密處理。網格單元類型方面,管道和止屈器部分選擇八節點六面體線性非協調模式單元(C3D8I),流體艙選擇有限薄膜應變線性減縮積分殼單元(S4R),后續在INP文件中將其修改為靜水流體單元(F3D4),如圖 5所示。

圖5 模型示意和使用單元Fig. 5 Finite element model and element type

管道止屈器模型建立完畢,然后通過修改INP文件以設置靜水流體單元。同時設置兩個歷史輸出變量CVOL和PCAV,它們分別表示注入流體的體積以及模型表面受到的壓力。

提交運算,結果表明該模型能夠成功地模擬管道發生局部壓潰、屈曲傳播和屈曲穿越止屈器的全過程。如圖6所示為屈曲穿越止屈器,即外界水壓等于穿越壓力Px時的Mises應力云圖。

圖6 屈曲穿越止屈器時Mises應力圖Fig. 6 Mises diagram when buckle passes through arrestors

提取歷史輸出變量CVOL和PCAV的數據,繪制得到管道表面所受壓力隨著流體艙內流體體積的變化曲線,如圖7所示。

圖7 壓力隨體積變化曲線Fig. 7 Pressure-volume curve

變化曲線可以表征整個模擬過程中管道止屈器發生的變形過程。第一個峰值表示管道在局部缺陷處,由于外界壓力作用,發生初始壓潰;隨著管道局部缺陷處發生壓潰,壓力雖然迅速回落,但壓潰導致的變形還是會不斷發展,直到管道內壁開始互相接觸;管道內壁互相接觸后,變形便停止并開始穩定向下游傳播,傳播過程中壓力穩定在屈曲傳播壓力附近,對應著曲線中間水平部分,該部分壓力的平均值可以作為管道的屈曲傳播壓力;屈曲不斷往下游傳播,直到屈曲傳播到止屈器附近,由于止屈器對屈曲傳播的阻礙作用,屈曲暫時不能通過止屈器,但“水”不斷注入,所以管道表面壓力開始逐漸增大,最后屈曲增大到能夠“穿越”止屈器并繼續向下游管道傳播,壓力又由一個峰值陡然下降,并最終穩定在其傳播壓力附近,曲線中第二個峰值的數值就是止屈器的穿越壓力。

可以看出,有限元模型能夠準確模擬屈曲穿越止屈器的全過程,通過有限元模型計算得到的曲線符合實際。另一方面,根據文獻[2]公布的數據進行了數值模擬,對比文獻中公布的試驗和數值模擬結果,對比結果如表 2所示。表中,D、t、La、h是試驗數據,分別表示管道外徑、管道壁厚、止屈器有效長度和止屈器厚度,ExperimentPx是試驗得出的穿越壓力,FEM1Px是文獻[2]中使用有限元方法得出的穿越壓力,FEM2Px是文中自行求解得出的穿越壓力,結果比較接近,可以驗證模型的準確性。

表2 對比結果Tab. 2 Comparative results

在止屈器試驗過程復雜,無法進行大量試驗的現實情況下,可以使用有限元模型對止屈器結構優化進行分析。

3 結構優化

整體式止屈器在限制深水海底管道屈曲傳播方面作用顯著,但傳統的止屈器構型上存在材料利用率低的問題,存在進行結構優化的空間。

整體式止屈器結構橫截面如圖 8所示,其中,D,t,h,Ls,La分別代表管道外徑、管道壁厚、止屈器厚度、止屈器總長度和止屈器有效長度。優化方案是在保證止屈效果的基礎上,盡可能減少止屈器材料使用量以降低生產成本。在管道外徑、管道厚度、止屈器總長、止屈器厚度不變的條件下,減小止屈器的有效長度,可以減少止屈器體積,也就減少了止屈器材料使用量。

圖8 止屈器橫截面Fig. 8 Cross section of buckle arrestors

圖9直觀展現了La=49.6 mm與La=1.0 mm兩種止屈器體積的對比,可以明顯看出,后一種止屈器可以比前者減少許多材料。

圖9 兩種止屈器體積對比Fig. 9 Volume difference between two buckle arrestors

圖10則展現了屈曲即將穿越過止屈器時,兩種止屈器的Mises應力圖。由Mises應力的分布看,有效長度為1.0 mm的止屈器,止屈器部分所受的應力基本是一致的,而有效長度為49.6 mm的止屈器,止屈器部分所受的應力并不均勻,止屈器靠近管道上游部分所受應力較大,說明這部分止屈器主要承受壓力,而它的邊緣部分所受應力較小,說明這部分止屈器并沒有充分發揮效用,實際上是一種材料的浪費。這也說明了減小止屈器的體積進行優化是有意義的。

圖10 兩種止屈器Mises應力對比Fig. 10 Mises difference between two buckle arrestors

為了更科學地分析止屈器的體積變化與穿越壓力變化之間的關系,研究該優化方法的可行性,改變止屈器的有效長度,建立了8組有限元模型,整理得到止屈器穿越壓力的變化及與初始模型的變化比、止屈器體積的變化以及與初始模型的變化比。8組止屈器的穿越壓力與體積的數值及變化值如表 3所示。

表3 各模型穿越壓力與體積對比Tab. 3 Px and volume difference among models

由對比結果可知,相比于有效長度較長的初始止屈器,即第⑧組有效長度為49.6 mm的止屈器,將止屈器有效長度減小后的止屈器具備幾乎相同的止屈能力,即使是對于第①組有效長度最小的止屈器,與初始模型相比,其穿越壓力的損失也不超過6%,卻能大大減小止屈器的體積,從而節省材料重量,最多能夠節省超過30%的材料重量,使得成本大大降低。也就是說,減少止屈器體積既可以保證止屈效果,同時節省材料用量,降低生產成本,體積的減少還使得安裝難度降低,有助于降低整個生產安裝過程的成本,因此這種減少止屈器體積的結構優化方案具有可行性。

4 結 語

以深水海底管道系統中常用的整體式止屈器為研究對象,利用有限元軟件進行數值模擬,分析了發生壓潰、屈曲傳播以及屈曲穿越止屈器的全過程,并以減少體積降低成本為目標,對止屈器進行了結構優化,得到以下結論:

1) 對比了廣義弧長法與靜水流體單元法在計算止屈器模型時的優缺點,選擇使用靜水流體單元法,并詳細介紹靜水流體單元相關內容。

2) 使用ABAQUS軟件建立數值模型,對帶有整體式止屈器的海底管道進行了從發生壓潰、屈曲傳播、屈曲穿越止屈器、屈曲向下游傳播全過程的模擬,提取相關數據繪制變化曲線,其各階段變化趨勢符合實際,驗證了有限元模型計算穿越壓力的準確性。

3) 減少止屈器體積的結構優化方案具有可行性,能夠在保證止屈能力的前提下,減少材料用量,有利于焊接、鋪設等工作的開展,具有一定的應用價值。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 伊人福利视频| 91精品在线视频观看| 午夜电影在线观看国产1区| 九九久久99精品| 美女一级毛片无遮挡内谢| 看你懂的巨臀中文字幕一区二区 | 人妻中文久热无码丝袜| 凹凸国产熟女精品视频| 国产美女一级毛片| 国产视频资源在线观看| 久久精品欧美一区二区| 97青青青国产在线播放| 久久精品国产999大香线焦| 大陆精大陆国产国语精品1024| 国产精品主播| 亚洲欧美一区二区三区图片| 亚洲精品无码av中文字幕| 久久久久88色偷偷| 成人福利一区二区视频在线| 国产黑丝一区| 国产亚洲精品自在线| 色婷婷亚洲综合五月| 美女国内精品自产拍在线播放| 亚洲国产系列| 亚洲成a人片| 国产91九色在线播放| 91成人在线免费视频| 国产高清在线观看91精品| 成人久久精品一区二区三区| 欧美成a人片在线观看| 亚洲一区精品视频在线| 亚洲成人精品在线| 永久免费精品视频| 全部免费特黄特色大片视频| 国产成人综合久久| 国产剧情伊人| 国产理论最新国产精品视频| 一级全免费视频播放| 国产香蕉国产精品偷在线观看| a级毛片视频免费观看| 亚洲天堂高清| 婷婷亚洲综合五月天在线| 国产又粗又猛又爽视频| 黄色网址手机国内免费在线观看| 中文字幕在线欧美| 天堂成人在线| 精品在线免费播放| Jizz国产色系免费| 国产黄在线免费观看| 国产三级视频网站| 亚洲国产中文欧美在线人成大黄瓜 | 久久semm亚洲国产| 456亚洲人成高清在线| 国产成人一级| 丝袜高跟美脚国产1区| 国产欧美日韩在线一区| 99视频国产精品| 国产女人在线| 亚洲精品va| 青青青国产在线播放| 无码一区中文字幕| 一级毛片在线播放| 亚洲视频色图| 国产精品永久免费嫩草研究院| AV不卡在线永久免费观看| 日韩在线播放中文字幕| 一本久道久久综合多人| 91香蕉国产亚洲一二三区 | 久久久久无码国产精品不卡| 国产亚洲现在一区二区中文| 国产精品香蕉在线观看不卡| 国产一级裸网站| 中文字幕乱码中文乱码51精品| 国产欧美日本在线观看| 久久综合国产乱子免费| 1024你懂的国产精品| 亚洲av无码成人专区| 国产第八页| 九九这里只有精品视频| 亚洲熟女中文字幕男人总站| 无码精品一区二区久久久| 999精品免费视频|