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

軸力作用下K型焊接管節點的應力集中系數分析

2010-09-22 07:15:22邵永波
船舶力學 2010年4期
關鍵詞:焊縫有限元

邵永波

(煙臺大學土木工程學院,山東 煙臺 264005)

1 引 言

K型管節點廣泛應用于海洋平臺結構中。在支承海洋平臺的管狀結構中,節點處通常是研究者最關注的部位,因為大多數管狀結構的破壞主要是由節點處的破壞所引起的。海洋平臺中的管節點結構主要承受由海浪和風造成的循環荷載,所以節點處易發生疲勞破壞。管節點的疲勞破壞主要表現為在焊縫處產生表面裂紋,表面裂紋的擴展會產生兩方面的影響:一是造成節點結構在焊縫處的剛度降低,從而降低了節點結構的極限抗載能力,尤其是裂紋穿透管壁后,其擴展速度非常快,節點在焊接處容易產生斷裂破壞,造成很嚴重的工程事故;二是表面裂紋穿透管壁后,管道內部的石油或者天然氣會發生泄露,造成海面的污染。因此,應確保海洋平臺結構在使用年限內不發生疲勞破壞。

在研究管節點疲勞壽命的時候,國內外經常使用的方法是參照S-N曲線,即通過計算管節點在焊縫處的熱點應力大小來估計節點所能承受的疲勞載荷的循環次數。在估算焊縫周圍的熱點應力的時候,經常采用的方法是計算節點在承受不同載荷作用下的應力集中系數。國內鋼結構設計規范(GB50017)中尚無提供對管節點應力集中系數的計算方法。國際上比較通用的計算管節點應力集中系數的設計規范是CIDECT[1]。但CIDECT中有關管節點應力集中系數的計算方法并不完善,需要根據節點的幾何參數參考一定的曲線來確定節點的應力集中系數。這在工程設計中是不方便的,而且從曲線上讀取數據也容易產生比較大的誤差。國內對管節點應力集中系數的研究在80年代末就開始了,并已經取得了一些研究成果,如文獻[2-6]。但是這些結果基本上是用試驗或者有限元對某類型的節點進行測試或者計算。并沒有提出一種簡單適用的計算管節點應力集中系數的方法,所以很難應用到工程設計中去。如何提出另外一種更方便地應用于工程設計中的節點應力集中系數的計算方法是國內外很多研究者所感興趣的。本文就是通過有限元分析和試驗測試相結合,研究軸力作用下K型管節點焊接處應力集中系數的計算方法。

通常K節點的應力集中系數是由其幾何參數、載荷形式、焊縫尺寸和邊界條件等多種因素決定。工程中的K節點承受的最普通的載荷類型是軸向力(通常是一個支管受拉,另外一個支管受壓)。K節點主管兩個端部可以看成是固支的。因此,在研究K節點應力集中系數的時候,基本上集中在其幾何參數和焊縫尺寸的影響上。研究如何模擬焊縫尺寸以及調查K節點幾何參數對其應力集中系數的影響是精確估算K節點應力集中系數的主要因素。

2 KK型管節點焊縫的模擬方法

KK管節點是通過把四個支管焊接到主管外表面而形成的。支管和主管之間的交線在空間上是一條形狀非常復雜的曲線。Cao等(1997)[7]提出了這條相交曲線的方程,為模擬焊縫提供了理論上的可行性。以前對管節點的研究經常采用兩種假設:一是假定管壁的厚度與管的直徑相比是很小的,因此忽略管壁厚度方向上應力的變化而把管壁看成一個曲殼;二是忽略焊縫對管節點應力集中系數的影響。這兩個假設都過低地估計了管節點的幾何剛度,從而會過高地估計焊縫處應力集中的程度。所以以前用于計算管節點應力集中系數的方法(如CIDECT)在工程設計中都是偏于保守的。為了更精確地計算管節點的應力集中系數,就必須考慮焊縫的影響。

在不考慮焊縫的情況下,支管和主管的內部和外部交線如圖1a和1b中的實線所示。考慮焊縫的時候,支管和主管實際上的內外部交線如虛線所示。外部交線將向外移動一段距離T2,而內部交線則移動T3的距離。不考慮焊縫的時候,支管和主管初始的接觸厚度為T1,如圖1b所示。

通常T1沿著交線的大小是發生變化的,其大小取決于圖2所示的兩面角的大小和主管的曲率大小。一般情況下,支管厚度tb是要遠小于主管半徑R1的,這種情況下,T1的大小接近于tb/sin()γ 。

為模擬焊趾W0(如圖1所示),可以從兩管交線上一點A0向外延伸一段距離T2。沿著主管和支管的交線,T2隨二面角γ0的不同而變化。角度γ0的變化范圍是從假設主管和支管之間夾角θ的最小值(工程中的K節點中,此夾角的最小值大約為30°)到180°。T2的大小可以通過下面的公式計算:

其中,T2是外部焊縫厚度,k2是外部交線修正因子,FOSouter是比例因子,m是一個常數,θs是主管和支管之間的最小夾角。

經過修正以后,焊趾部分的外部交線方程可以寫成如下形式:

其中,φ0為焊縫法線方向與主管軸向之間的夾角。

在公式(3)中,點 A0(XA0,YA0,ZA0)是外部交線上一點,點B0是由點A0沿支管方向移動(T1+T4)而得到的,如圖3a和3b所示。

按照相同的方法,內部相交線上的點Ai到焊跟Wi的距離為T3,如圖3所示。二面角γi的范圍是從30°到 90°。 當 γi=30°時 T3取極大值,而當 γi=90°時 T3=0。 T3可以表示為如下形式:

其中,T3是外部焊縫厚度,k3是外部交線修正因子,FOSinner是比例因子,n是一個常數,θs是主管和支管之間的最小夾角。

經過修正以后,焊跟部分的內部交線方程可以寫成如下形式:

管節點的焊接厚度Tw是由T1,T2和T3共同決定的,且必須滿足有關規范(如AWS,2000[8])規定的焊縫厚度的最小要求。Tw可表示為:

其中kaws是由AWS(2000)[3]規范規定的焊縫厚度參數。

3 KK管節點有限元網格的產生

3.1 單元類型的選擇

以往對管節點進行數值模擬時,基本上采用的是殼單元。這是由于大部分管節點的壁厚大小遠小于其直徑的原因。而且先前的研究工作中往往忽略了焊縫尺寸對管節點應力集中系數的影響。但是,殼單元在模擬管節點時有兩個很難克服的缺點:一是在管壁厚度與直徑相比不是非常大的情況下,應力沿著壁厚方向的變化難以反應出來;二是難以比較精確地模擬焊縫形狀。Lee和Bowness[9]通過對管節點的研究發現管壁內外表面的應力大小并不相同,這種差異是由于管壁內的彎曲應力造成的。為了克服殼單元在模擬管節點上的缺陷,Herion等[10]通過研究提出了三維二次實體單元是最適合模擬帶焊縫節點的單元。在本研究中,這種20結點的二次單元也被采用來模擬帶焊縫的K節點。

3.2 有限元網格的產生原理和方法

采用有限單元法對結構進行分析時,所得到的數值結果的精度在很大程度上取決于所采用有限元網格中單元的質量和數量。所謂單元的質量是指在有限元網格中劃分的單元不能過度扭曲,尤其是在應力梯度比較大的區域單元的邊長比例不能太大。單元的數量是用來保證有限元結果的收斂性。在對一個結構進行有限元分析時,往往關心的是應力梯度高的區域,因為這個區域一般是結構發生破壞的地方。因此,在高應力梯度的區域,有限元網格可以劃分得精密一些,而且應該控制這個區域的單元質量。反之,低應力區一般是結構分析中不關心的區域,為了提高計算效率,這個區域的網格劃分可以相對粗糙一些,而且對單元質量的要求也不必太高。

根據以上原則,可以根據計算精度的要求,將整個K節點結構分成不同的區域,在每個區域單獨進行網格劃分。由于熱點應力位于焊縫位置處,這個位置應力集中很嚴重,所以焊縫位置的網格劃分需要精密,而且單元質量要求高。遠離焊縫的區域對K節點應力集中的影響比較小,在這些區域可以采用比較粗糙的網格以減少計算時間。按照這個標準,整個K節點可以劃分為幾個不同的區域,如圖4a所示,這幾個區域的網格均單獨生成。當每個區域網格生成以后,K節點的整體有限元網格可以通過合并各個區域的網格而得到,如圖4b所示。

采用分區法產生K節點有限元網格可以保證焊縫處的網格的高質量,同時也保證這個區域網格密度大。在焊縫區域,沿著管壁厚度方向上共劃分了三層單元。為了驗證有限元結果的收斂性,本研究中提出了一種K節點網格自動加密方案。在這個方案中,每個單元的每條邊都被平均分成兩份,這樣一個母單元被平均分成8個子單元。經過這樣加密后的網格中單元數量是初始網格中單元數量的8倍。圖5a中顯示了焊縫區域的網格的加密過程,圖5b則顯示了對圖4b中的K節點網格的加密結果。

4 K節點應力集中系數的計算方法

4.1 K節點應力集中系數的定義

管節點在焊縫處的應力分布不均勻,如圖6所示。在軸向力作用下,支管在遠離加力點和焊縫處截面上的應力分布是均勻的,大小為σn,稱為名義應力。在焊縫周圍,由于存在著應力集中,應力要比名義應力大得多,這個區域內的應力稱為熱點應力σh。熱點應力和名義應力的比值稱為應力集中系數,它反應了主管和支管交線處由于曲率不連續而造成的應力“突變”程度。

按照上述定義,K型管節點的應力集中系數(SCF)可以由以下公式計算得到:

公式(8)中的名義應力可以由下面的公式計算得到:

公式(9)中,d和t分別為支管的直徑和厚度。

公式(8)定義了管節點應力集中系數的計算方法。但是,以前的研究工作中對如何計算熱點應力一直存在著不同的看法。目前比較普遍采用的對熱點應力的定義有兩種:第一種是認為熱點應力區內某點的熱點應力就是該點處的最大主應力;第二種觀點認為焊縫處某點的熱點應力是垂直于焊縫的法線方向的應力。由于熱點應力是用來估算管節點的疲勞壽命的,管節點的疲勞破壞主要是沿著焊縫處表面裂紋的萌生和擴展,因此垂直于焊縫方向的應力應該對裂紋的擴展起主要作用。所以,此處采用第二種規定作為熱點應力的定義。

4.2 K節點應力集中系數的試驗測試

為了驗證有限元方法在計算K型管節點應力集中系數的準確性,測試了兩個大型K節點在軸力作用下焊縫周圍應力集中系數的分布。兩個K節點試件的幾何構造如圖7a和7b所示。K節點試件由管狀空心鋼管制作而成,所用的鋼管符合API 5L B級規范要求。支管與主管的相交部位采用輪廓線切割技術加工得到,這樣可以保證支管和主管在相交部位吻合得很好,已便于焊接。焊接方式采用的穿透焊縫。 K型管節點的中幾何參數為:α=43.3(α=2L/D,L為主管長度,D為主管直徑),β=0.516(β=d/D,d 為支管直徑),τ=0.75(τ=t/T,t為支管厚度,T 為主管厚度)和 γ=5.376(γ=D/2T)。K 節點試件的 γ 取值小于有關規范(如CIDECT)規定的范圍,這是為了驗證焊縫尺寸對應力集中系數的影響。以往對大型K節點試件的測試也曾有過報道,如Soh和Soh[11]曾對薄壁的K節點的應力集中系數進行過測試。由于管壁很薄,所以γ值很大。在數值模擬時,焊縫尺寸就被忽略了。當γ很小時,焊縫尺寸相對比較大,因而就不能忽略它對管節點應力集中系數的影響。

在進行試驗測試的時候,K節點試件主管的兩個端部通過螺栓固定在墻體上,可以視為固支邊界條件。兩個支管中的一個在端部也用螺栓固定在墻體上,在另外一個支管的端部施加軸向力,如圖8所示。

在試驗測試中,采用應變片測量熱點應力區附近的應力。因為應變片有其自身的尺寸要求,所以它不能直接放置在焊趾部位來測量這個部位的應力,而是通常在一個規定的插值區域內布置應變片。在管節點中,用線性外推插值法來確定熱點應力。這個插值區距離焊縫有一定的距離限制,其范圍可參考有關規范 (如CIDECT)。如圖9所示,在規范規定的插值區內可以粘貼兩片(線性插值)或者三片(二次插值)應變片。焊縫處的應力可以通過將插值區內的應力外推插值得到。

在試驗中,支管和主管上在焊縫附近的插值區內每隔15°布置應變片,如圖10所示。在每個點上,用兩個應變片測量垂直于焊接路徑上的應變,兩個應變片布置在CIDECT規范所規定的外推插值范圍內。試驗中在主管上應變片離焊縫的最短距離Lr,min和最遠距離Lr,max分別為10mm和18mm,支管上這兩個距離分別為8mm和20mm。

4.3 K節點應力集中系數的試驗和數值結果

采用提出的有限元模型,可以分析得到承受軸力的K節點試件焊縫處的應力集中系數分布。數值結果的精確程度可以通過和試驗測試得到的結果進行比較,以驗證有限元模型的適用性。圖11a和11b中對兩個K節點試件的有限元結果和試驗測試結果進行了比較。從圖中可以看出,有限元結果和試驗結果吻合得比較好,它比較準確地反映了應力集中程度沿著焊縫的分布情況。從圖11中也可以發現:對于承受軸力的K節點試件,當管壁厚度相對于直徑而言比較大的時候(即γ比較?。?,主管上沿著焊縫周圍的最大應力位于冠點處。支管上沿著焊縫應力變化不是特別明顯,它受焊縫尺寸的影響比較明顯。

通過對K節點的有限元分析結果和試驗測試結果的比較,可以得到如下結論:本文提出的對焊接K節點的幾何和數值模擬方法可以比較準確地模擬其熱點應力區域內沿著焊縫周圍的應力分布情況,因此這種提出的模型可以可靠地用來分析軸力作用下K節點應力集中系數的大小。

5 計算K節點應力集中系數的參數公式

5.1 幾何參數對K節點應力集中系數的影響

K節點焊縫周圍的應力分布情況是受許多因素影響的,如幾何參數、載荷類型、邊界條件、焊接質量和初始缺陷等。在本研究中,將不考慮焊接質量和初始缺陷的影響,而且只考慮K節點承受軸向平衡力作用力的作用(如圖12所示),所以載荷類型是確定的。在對K節點進行分析時,主管兩個端部是固支或簡支的,這里采用固支的邊界條件??紤]了以上情況后,承受軸力的K節點的應力集中系數就由節點的幾何參數決定。國內外對管節點進行分析的時候,都是用幾個無量綱的幾何參數來表示不同的管節點類型。對于K節點,采用的幾何參數有:α(主管長度L/主管半徑R)、β(支管半徑r/主管半徑R)、γ(主管半徑R/主管厚度T)、τ(支管厚度t/主管厚度T)和主管與支管之間的夾角θ。在α大于一定的數值的時候,邊界條件(固支或者簡支)對應力分布基本沒有影響,所以,在模型分析中,將α的值固定在15上,從而不需要再考慮α和邊界條件的影響?;谝陨峡紤],在模型分析中,只需要研究幾何參數β、γ、τ和θ對K節點在承受軸向力作用下沿著焊縫周圍應力集中系數的大小。在模型分析中,如果兩個K節點具有相同的幾何參數β、γ、τ和θ,則這兩個K節點具有相同大小的應力集中系數,而且沿著焊縫周圍應力集中系數的分布情況也是相同的。

在K節點的參數分析中,模型是根據四個幾何參數(β、γ、τ和θ)的不同而定義的。K節點幾何參數范圍如下:

· 0.30≤β≤0.60

· 12≤γ≤30

· 0.25≤τ≤0.95

· 30°≤θ≤60°

·兩個支管尺寸相同,且與主管有相同大小的夾角

·支管之間有間隙,不考慮交疊K節點。兩個支管和主管的軸線交于一點,即無偏心。

在模型的參數分析中,每個幾何參數的取值如表1所示。根據表1中所列幾何參數的范圍和步長,總共分析了1008個K節點模型。

表1 模型分析中幾何參數的范圍Tab.1 Geometrical range of the analyzed models

根據參數分析的結果,可以總結出K節點在軸向荷載作用下的應力集中系數隨著四個幾何參數的變化而改變的規律。在研究某個幾何參數對K節點應力集中系數的影響時,把其它三個參數的大小固定。γ和β對K節點應力集中系數(SCF)大小的影響如圖13a和13b所示。在θ,γ和τ一定的情況下,主管上的應力集中系數隨著β的增大而增大,而支管上的應力集中系數卻隨著β的增大而減小。在θ,β和τ固定的情況下,主管和支管的應力集中系數均隨γ的增大而增大。

其它兩個幾何參數(τ和θ)對K節點應力集中的影響如圖14a和14b所示。K節點主管和支管上的應力集中系數均隨著τ的增大而增大。同樣,主管和支管上的應力集中也都將隨著θ的增大而增大。

5.2 計算K節點應力集中系數的參數方程

目前國內關于K型管節點應力集中系數的計算公式尚無報道。國際上比較流行使用的CIDECT規范中對軸力作用下K節點應力集中系數的計算曾提供了一種計算方法。這種方法是把K節點的應力集中系數表示為幾何參數的函數,但是需要從有關圖中確定應力集中系數的大小。這樣在從圖中讀取數據時必然存在誤差。為了簡化對K節點應力集中系數的計算方法,這里把K節點應力集中系數表示為如下形式:

指數a和b的變化范圍為0到1.1。SCF0是一個關于參數β的方程,SCF0可以由二次多項式的曲線擬合得到。μ是一個修正系數,它是由夾角θ決定。

在擬合用來估算K節點應力集中系數參數公式時,所提出的參數方程必須滿足如下要求:在參數適用范圍內,所提方程必須能夠準確反映各幾何參數對應力集中系數的影響規律,且所提參數方程必須簡單適用,可以方便地用于工程設計中。方程(10)便是一個典型的可用來計算K型管節點應力集中系數的參數公式。它是基于曲線擬合法得到的,該方程形式滿足以上要求。在參數分析中,分別對θ=30°,45°和60°時K節點的應力集中系數進行分析,幾何參數對節點SCF的影響都可以反映在方程中。

使用軟件DataFit6.0[12]對1008個K節點模型的應力集中系數結果進行了非線性曲線回歸分析,得出了離散度在允許范圍內的K節點應力集中系數的參數方程,共有六個參數方程,如表2和表3所示。表2為K節點在軸向荷載作用下的支管的應力集中系數參數方程,表3為主管的應力集中系數參數方程。當K節點的主管和支管間夾角θ不等于30°,45°或60°時,其應力集中系數值可以由相同幾何參數(β,τ和γ)的K節點在θ為30°,45°或60°時的應力集中系數值通過內部插值法得到。

表2 軸力下K節點支管上應力集中系數的參數方程Tab.2 Parametric equations of SCFs on the brace for K-joint under axial loads

表3 軸力下K節點主管上應力集中系數的參數方程Tab.3 Parametric equations of SCFs on the chord for K-joint under axial loads

5.3 K節點應力集中系數參數方程精確性的驗證

擬合得到的用于計算K節點應力集中系數的參數方程需要驗證其精確性。表2和表3所列的參數方程是基于對1008個K節點模型結果的基礎上。這1008個節點模型的幾何范圍可以從表1中得到。如果所提出的K節點應力集中系數的參數方程的精確性得到驗證,那么它可以用來比較精確地計算表1所規定幾何范圍內的任何一個K節點的應力集中系數。

在驗證提出的參數公式的精確性時,可以通過將從參數公式計算得到的結果與有限元計算得到的結果進行比較,觀察其相對誤差的大小。相對誤差Er可以通過下式表示:

其中,SCFE是由參數方程得到的K節點應力集中系數值,SCFN是有限元分析計算得到的K節點應力集中系數值。

對參數分析中的1008個K節點模型的應力集中系數值的相對誤差按照公式(11)進行計算,所得到的結果如圖15所示。圖15中,橫坐標N表示模型個數。從圖15中很容易看出:對于1008個K節點模型,由參數公式計算得到的K節點應力集中系數的大小和有限元結果比較,其相對誤差絕大多數都在±20%以內。所以,所提出的參數方程在計算軸力作用下K節點應力集中系數大小方面是精確可靠的。

6 結 論

本文提出了對K型管節點的幾何和數值模擬,尤其是對焊縫的模擬。通過對兩個大型K節點模型的試驗測試,驗證了所提出的有限元模型在計算軸力作用下K節點沿著焊縫處的應力集中系數時是準確可靠的。在此基礎上,通過對1008個K節點模型的參數分析,研究了幾何參數對K節點應力集中系數的影響。最后通過曲線擬合提出了K節點在軸向力作用下支管和主管上應力集中系數的參數方程。通過將參數方程計算所得的應力集中系數的結果與有限元分析得到的應力集中系數的結果來進行比較,發現所提出的參數方程具有相當的精確度,從而可以安全可靠地用于工程設計中。

[1]Zhao X L,Herion S,Packer J A,Puthli R,Sedlacek G,Wardenier J,Weynand K,van Wingerde A,Yeomans N.Design guide for circular and rectangular hollow section joints under fatigue loading[K],CIDECT,TUV,2000.

[2]陳鐵云.近海鉆井平臺管狀接頭應力分析的進展[J].力學進展,1985,4(25):425-433.

[3]陳伯真.海洋平臺管節點應力分析研究工作展望[J].力學進展,1987,17(3):331-336.

[4]李華明.海洋平臺管節點的熱點應力測量[J].振動、測試與診斷,1991,11(3):23-29.

[5]陳鐵云,朱正宏.海洋平臺T、Y、K及塔接K型管接頭的彈塑性有限元分析[J].海洋工程,1992,10(2):8-15.

[6]王春光,李少甫,石永久.焊接鋼管節點熱點應力集中系數參數方程的適用性研究[J].工程力學,1999,16(6):44-53.

[7]Cao J,Yang G,Packer J A.FE mesh generation for circular joints with or without cracks[C]//In:Proceedings of the 7th International Offshore and Polar Engineering Conference.USA,1997,IV:98-105.

[8]American Welding Society.Structural welding code-steel[S].ANSI/AWS,D1.1-2000.Miami,USA,2000.

[9]Lee M M K,Bowness D.Prediction of stress intensity factors in semi-elliptical weld toe cracks in offshore tubular joints[C]//Proceedings of the 9th International Symposium and Euroconference on Tubular Structures.Dusseldorf,Germany,2001:299-308.

[10]Herion S,Mang F,Puthli R.Parametric study on multiplanar K-joints with gap made of circular hollow sections by means of the finite element method[C]//In:Proceedings of the 6th International Offshore and Polar Engineering Conference.Los Angeles,USA,1996,IV:68-73.

[11]Soh A K,Soh C K.Stress concentrations of K-tubular joints subjected to basic and combined loadings[C].Proceedings of the Institute of Civil Engineering-Structures and Buildings,1996,116:19-28.

[12]DataFit 6.0.Oakdale Engineering[DB].23 Tomey Road,Oakdale,Pa,USA,1998.

猜你喜歡
焊縫有限元
基于焊縫余高對超聲波探傷的影響分析
焊縫符號在機械設計圖上的標注
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
TP347制氫轉油線焊縫裂紋返修
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
焊縫跟蹤遺傳算法優化PID控制仿真研究
機器人在輪輞焊縫打磨工藝中的應用
光譜分析在檢驗焊縫缺陷中的應用
磨削淬硬殘余應力的有限元分析
主站蜘蛛池模板: 久久不卡精品| 日韩高清中文字幕| 尤物午夜福利视频| 久久精品丝袜| 无码日韩人妻精品久久蜜桃| 国产v精品成人免费视频71pao| 无码丝袜人妻| 高潮毛片免费观看| 五月天婷婷网亚洲综合在线| 免费在线不卡视频| 中文字幕无码av专区久久| 久久久久88色偷偷| 久久国产精品影院| 日韩精品一区二区三区swag| 国产成人欧美| 欧美无专区| 精品国产Ⅴ无码大片在线观看81| 日本在线免费网站| 久久精品国产精品青草app| 国产微拍一区| 粗大猛烈进出高潮视频无码| 亚洲黄网视频| 毛片一级在线| 国产主播在线一区| 国产精品99一区不卡| 中文字幕乱码中文乱码51精品| 无码免费的亚洲视频| 国产精品美乳| 美女被躁出白浆视频播放| 永久免费AⅤ无码网站在线观看| 成年人视频一区二区| 国产成人高清精品免费软件| 日韩一区精品视频一区二区| 日本亚洲欧美在线| 五月天婷婷网亚洲综合在线| 亚洲精品视频免费| 精品一区二区三区中文字幕| 久久精品66| 伊人成色综合网| 激情综合图区| 五月天久久婷婷| 黄色三级毛片网站| 国内视频精品| 99热最新网址| 精品自窥自偷在线看| 欧美成人二区| 亚洲性色永久网址| 这里只有精品在线播放| 中文字幕第1页在线播| 亚洲视频一区| 国产精品久久久久久搜索| 亚洲日韩精品欧美中文字幕| 国产丝袜91| 老司国产精品视频91| 日韩毛片免费视频| 欧美精品v欧洲精品| 欧美在线伊人| 国产欧美精品一区aⅴ影院| 自偷自拍三级全三级视频| 玖玖精品在线| 日本不卡免费高清视频| 午夜国产小视频| 午夜一区二区三区| 久久国产精品麻豆系列| 日本一区二区三区精品视频| 国产91成人| 中文字幕欧美日韩| 91网红精品在线观看| 少妇极品熟妇人妻专区视频| 天堂va亚洲va欧美va国产 | 亚洲福利一区二区三区| 国产自在线播放| 中文字幕亚洲第一| 免费中文字幕一级毛片| a天堂视频| 国产精品成人久久| 九九久久精品免费观看| 亚洲AⅤ永久无码精品毛片| 亚洲欧美不卡中文字幕| 国产成年无码AⅤ片在线 | 尤物成AV人片在线观看| 国产91线观看|