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

焊趾處橢圓表面裂紋的權函數與殘余應力強度因子的權函數法

2017-05-02 03:33:31徐磊黃小平
船舶力學 2017年4期
關鍵詞:裂紋有限元

徐磊,黃小平

(上海交通大學海洋工程國家重點實驗室,上海200240)

焊趾處橢圓表面裂紋的權函數與殘余應力強度因子的權函數法

徐磊,黃小平

(上海交通大學海洋工程國家重點實驗室,上海200240)

利用三維有限元計算了焊趾處半橢圓表面裂紋的應力強度因子。利用統一的權函數形式,結合得到裂紋半長比a/c=0.2;0.4;0.6;0.8,a/t=0.1~0.8的有限元數據,得到了適用于T型接頭焊趾處半橢圓表面裂紋最深點和表面點的權函數。權函數的準確性,用有限元在裂紋面施加高階載荷進行了驗證,對于表面點和最深點,半長比a/c= 0.2~0.8,a/t=0.1~0.8,權函數與有限元結果誤差在8%以下?;诘玫降臋嗪瘮?,計算了T型接頭焊趾處半橢圓表面裂紋的殘余應力強度因子Kres,并與有限元計算結果進行對比,對比誤差在10%以下,表明新的權函數能很好地預測T型接頭焊趾處的殘余應力強度因子。

T型接頭;焊趾;表面裂紋;應力強度因子;權函數;殘余應力

0 引言

工程結構中存在著大量的焊接接頭,焊接接頭焊趾處的半橢圓裂紋通常是該結構存在的最普遍的缺陷形式。計算這些結構物中出現裂紋以后,裂尖應力強度因子的大小,是損傷容限設計極其重要的一部分[1]。

對于平板含有橢圓表面裂紋應力強度因子的計算,已經有大量的學者對此進行了研究。其中認可度最廣的是Raju與Newman[2]基于大量有限元數據得到的經驗公式,英國最新的金屬結構裂紋驗收評定方法指南(BS7910)[3]便是以該公式為基礎,提出了平板半橢圓裂紋應力強度因子的計算表達式。當裂紋存在于平板對接頭,T型焊接頭的焊趾處時,焊趾處的應力集中和焊趾復雜的幾何模型就會使得計算應力強度因子變得相對困難。對于焊趾的處理,Bowness[4]在總結前人工作的基礎上,結合三維有限元分析,給出了焊趾處應力強度因子的放大系數Mk的計算公式。

對于以上的經驗公式,其使用范圍均是在已知遠場應力載荷形式的情況下得出的。實際結構中的載荷應力形式往往比已知的遠場外載荷復雜得多,尤其是在焊趾處,通常的制造過程將在焊趾處引入很大的焊接殘余應力。因此,怎么計算由殘余應力引起的殘余應力強度因子Kres,由此來評估殘余應力對疲勞壽命的影響,變成了另外一個問題。因為已有的經驗公式并不能適用于裂紋近場殘余應力導致的應力強度因子的計算,學者Bao Rui[5]運用權函數計算了幾種簡單試件的殘余應力強度因子,Labeas[6]給出了數值計算穿透裂紋的殘余應力強度因子的方法。

利用權函數來計算應力強度因子,很好地解決了殘余應力強度因子的計算問題。權函數法是求解在任意載荷分布下,尤其是裂紋處的應力場分布很復雜的情況下,計算裂紋應力強度因子具有十分高的效率。對于平板表面半橢圓裂紋的應力強度因子,Wang[7-8]通過對含有三維表面橢圓裂紋有限平板在線性載荷與拋物線載荷下進行了裂紋應力強度因子的計算,并結合Shiratori等人[9]的數據,在已有的統一權函數形式上,通過有限元計算,分別給出了平板裂紋小半長比(a/c<1)和大半長比(a/c>1)裂紋的權函數。本文將以求解平板權函數的思路,結合三維有限元裂紋分析的數據,提出一個適用于T型接頭焊趾處表面裂紋的權函數,并用二次和三次應力分布驗證了其有效性,基于得到的新權函數,計算T型接頭焊趾處半橢圓表面裂紋的殘余應力強度因子Kres,并與有限元計算結果進行對比。探討權函數在殘余應力強度因子計算中的適用性。

1 權函數的統一形式

Bueckner[10]和Rice[11]提出了任意載荷下,計算裂紋應力強度因子K的權函數法。其具體表達式為

式中:σ(x)為假想無裂紋體裂紋處應力分布,m(x,a)為權函數。

Glinka和Shen[12]指出,對于一維與二維的裂紋,權函數可有如下的統一形式

當運用到半橢圓的表面裂紋時,當x=a時表示為裂紋的最深點,當x=0時表示為裂紋的表面點,如圖1所示。統一的權函數運用到半橢圓表面裂紋最深點時的計算公式如下:

圖1 T型接頭焊趾裂紋幾何示意圖Fig.1 T-butt weld toe crack geometry

統一的權函數運用到半橢圓表面裂紋表面點時的計算公式如下:

式中:M1A,M2A,M3A和M1B,M2B,M3B為對應于T型接頭焊趾處的計算系數,可以通過在裂紋面施加常數和線性載荷,并結合推導裂紋權函數的自適應條件來求得。Glinka和Shen在推導權函數的統一形式時,指出x=a,即裂紋最深點權函數的自適應條件為:

x=0,即裂紋表面點權函數的自適應條件為

2 焊趾表面裂紋的權函數

2.1 有限元模型的建立

有限元計算模型為T型接頭,長L=100 mm;寬W=100 mm;板厚H=10 mm;焊角θ=45°;焊趾寬度t=5 mm;翼板高度h=25 mm;楊氏模量E=2.1e5MPa;泊松比ν=0.3;橢圓裂紋長度為c,深度為a。具體如圖1所示。

有限元在計算應力強度因子時,邊界條件為固定四個頂點;載荷施加方式為直接施加在兩個裂紋面上。有限元軟件為ANSYS14.0;采用的單元為20節點的solid186高階體單元,其中裂尖單元采用奇異單元處理。計算應力強度因子的方法為ANSYS自帶的位移插值法。求得的SIF值采用如下公式進行無量綱化:

圖2 T型接頭焊趾處有限元模型Fig.2 T-butt weld toe crack FE model

2.2 裂紋最深點處的權函數

為了求得焊趾處裂紋最深點的權函數,即要求出(3)式中對應于T型接頭的三個系數M1A,M2A, M3A。求解方法是,在裂紋面上分別直接施加常數載荷與線性載荷,并求得與之對應的最深點應力強度因子、,結合裂紋最深點權函數的自適應條件式(5),即可求得裂紋最深點的權函數。對于裂紋最深點處的應力強度因子,數值計算了半長比a/c=0.2;0.4;0.6;0.8,a/t=0.1~0.8之間的應力強度因子,給出了誤差在5%以下的擬合公式。

常載荷下:

表1 T-butt SIF有限元結果σ(x)=σ0Tab.1 T-butt SIF FEM result σ(x)=σ0

續表1

圖3 載荷σ(x)=σ0的SIF擬合曲面(左—表面點;右—最深點)Fig.3 SIF fitting curve when load σ(x)=σ0(left-surface point;right-deepest point)

2.3 裂紋表面點處的權函數

與求解最深點處的計算方法類似。在裂紋面上分別直接施加常數載荷σ(x)=σ0與線性載荷并求得與之對應表面點的應力強度因子和,結合裂紋表面點權函數的自適應條件,等式(6),即可求得裂紋表面點的權函數。數值計算了半長比a/c=0.2;0.4;0.6;0.8,a/t=0.1~0.8之間裂紋表面點的應力強度因子,并給出了誤差在5%以下的擬合公式。

常載荷下:

結果:

線性載荷下:

結果:

將(15)~(20)式代入(1)式,并結合裂紋表面點權函數的自適應條件式(6),求得對應于T型接頭表面半橢圓裂紋最深點處權函數的計算系數為:

表2 T-butt SIF有限元結果σ(x)=σ0(1-x/a)Tab.2 T-butt SIF FEM result σ(x)=σ0(1-x/a)

圖4 載荷σ(x)=σ0(1-x/a)的SIF擬合曲面(左—表面點;右—最深點)Fig.4 SIF fitting curve when load σ(x)=σ0(1-x/a),(left-surface point;right-deepest point)

3 權函數有效性的驗證

在Glinka和Shen提出已有的統一權函數基礎上,依據Wang求解平板權函數的思路,計算了裂紋半長比a/c=0.2;0.4;0.6;0.8,a/t=0.1~0.8之間的應力強度因子,得到了針對T型接頭焊趾處,計算

表面半橢圓裂紋應力強度因子的權函數。為了驗證該權函數的有效性,將一階與二價載荷下得到的權函數擴展到高階載荷的情形,權函數計算得到的高階載荷下裂紋應力強度因子與FEA結果進行對比。驗證的思路為:若該權函數能夠在已知裂紋處近場應力σ()x的分布情況下,能準確地計算裂紋處的應力強度因子,那么T型接頭裂紋面在高階載荷的作用下FEA得到的SIF,應與得到的權函數計算的SIF的計算值有很好的吻合。

表3 T-butt SIF有限元結果_二次載荷σ(x)=σ0(1-x/a)2Tab.3 T-butt SIF FEM result_parabolic load σ(x)=σ0(1-x/a)2

表4 T-butt SIF有限元結果_三次載荷σ(x)=σ0(1-x/a)3Tab.4 T-butt SIF FEM result_Cubic load σ(x)=σ0(1-x/a)3

表5 T-butt SIF權函數結果_二次載荷σ(x)=σ0(1-x/a)2Tab.5 T-butt SIF WFM result_parabolic load σ(x)=σ0(1-x/a)2

表6 T-butt SIF權函數結果_三次載荷σ(x)=σ0(1-x/a)3Tab.6 T-butt SIF WFM result_Cubic load σ(x)=σ0(1-x/a)3

焊趾處表面裂紋高階載荷作用下FEA得到的SIF與權函數對比結果:

圖5所示,二次載荷作用下表面點SIF,在裂紋半長比a/c=0.2~0.8時,a/t= 0.1~0.8時,權函數法與FEA得到的結果,兩者之間差值的百分比最大值在a/c=0.6,a/t=0.6處取得為6.3%,在a/c=0.6,a/t=0.4處,二者差值也達到6.1%;其余點,二者之間差值的百分比均在5%以下。

圖6所示,二次載荷作用下最深點SIF,在裂紋半長比a/c=0.2~0.8時,a/t= 0.1~0.8時,權函數法與FEA得到的結果,兩者之間差值的百分比最大值在a/c=0.4,a/t=0.8處取得為14.3%,在a/c=0.2,a/t=0.3處,二者差值也達到5.2%;其余點,二者之間差值的百分比均在5%以下。

圖7所示,三次載荷作用下表面點SIF,在裂紋半長比a/c=0.2~0.8時,a/t= 0.1~0.8時,權函數法與FEA得到的結果,兩者之間差值的百分比最大值在a/c=0.4,a/t=0.6處取得為16.4%,在a/c=0.6,a/t=0.4處,二者差值也達到8%;其余點,二者之間差值的百分比均在5%以下。

圖8所示,三次載荷作用下最深點SIF,在裂紋半長比a/c=0.2~0.8時,a/t=0.1~0.8時,權函數法與FEA得到的結果,兩者之間差值的百分比最大值在a/c=0.4,a/t= 0.8處取得為10.6%,其中在a/c=0.8,a/t= 0.1,二者差值也達到9.3%;其余點,二者之間差值的百分比均在8%以下。

圖5 二次載荷作用下權函數SIF與有限元SIF結果對比(表面點)Fig.5 Comparison between weight function result and FEM result for parabolic stress distribution(Surface point)

圖6 二次載荷作用下權函數SIF與有限元SIF結果對比(最深點)Fig.6 Comparison between weight function result and FEM result for parabolic stress distribution(Deepest point)

圖7 三次載荷作用下權函數SIF與有限元SIF結果對比(表面點)Fig.7 Comparison between weight function result and FEM result for cubic stress distribution(Surface point)

4 焊趾處半橢圓表面殘余應力強度因子Kres的權函數法

圖8 三次載荷作用下權函數SIF與有限元SIF結果對比(最深點)Fig.8 Comparison between weight function result and FEM result for cubic stress distribution(Deepest point)

4.1 焊接殘余應力的分布

采用權函數法計算應力強度因子,必須知道焊趾處垂直于裂紋面上沿板厚度的方向的應力分布。在過去的研究中,對于殘余應力場中的應力強度因子計算,由于殘余應力分布的復雜,通常采用已有的殘余應力分布經驗公式。

文獻[13]提出在鋼材焊接后的焊趾表面處的殘余應力大小為

式中:σR表示殘余應力,σy表示材料的屈服極限。

在此基礎上,文獻[14]給出了厚板多道焊接后,殘余應力沿著板厚度方向上的分布,其簡化形式為如下:

但是實際過程中,焊接后殘余應力的分布是十分復雜的,焊趾處殘余應力沿板厚度方向上的分布與板厚,焊接時線熱量的輸入都有十分密切的關系。為了得到更加準確的結果,對T型接頭進行焊接過程的模擬。彈性模量E=2.1×105MPa,泊松比ν=0.3,σy=235 MPa。其它熱力學的參數參考文獻[10]。

圖9 T型接頭焊接后的典型殘余應力場Fig.9 T-butt residual stress distribution

4.2 焊趾處垂直裂紋面殘余應力的確定

由圖10可以看到T型接頭焊接后沿著焊縫方向的殘余應力場,縱向殘余應力和橫向殘余應力由板的兩端向中間的過程中,有一個上升和下降的趨勢,而在板的中間處,縱向殘余應力與橫向殘余應力隨著離起焊點距離S而變化,但其變化趨勢不大,故而為了得到簡化的沿板厚方向垂直于裂紋面的橫向殘余應力,我們可以假設:在焊趾處裂紋的位置,沿著焊縫方向裂紋面上的橫向殘余應力不變,只需要提取出焊趾處橢圓裂紋處沿著厚度方向上橫向殘余應力σ()x。

圖10 沿焊縫方向殘余應力分布Fig.10 Residual stress distribution along welding direction

圖11 焊趾處殘余應力沿板厚分布擬合曲線Fig.11 T-butt residual stress fit curve

圖11所示為T型頭底板厚度t=10 mm焊接后,提取的沿厚度方向上垂直于裂紋面橫向殘余應力的散點分布圖。對于提取出來的殘余應力離散點,采用多項式擬合后,得到了殘余應力沿厚度方向的分布函數,(24)式。沿厚度方向上的殘余應力分布確定后,結合T型接頭半橢圓裂紋的權函數,則可求出T型接頭裂紋最深點與表面點的應力強度因子。

5 Kres結果討論

對T型接頭焊趾處半橢圓表面裂紋在殘余應力場中的應力強度因子Kres,分別采用了權函數法和有限元法進行計算。其中,導出的T型接頭焊趾處半橢圓裂紋表面裂紋權函數的準確性用高次應力載荷進行了驗證,結果與有限元法相比,其求解誤差大體在8%以下。在此基礎上,驗證該權函數計算殘余應力強度因子Kres的準確性,并與有限元結果進行對比。T型接頭焊趾處裂紋在(24)式所示的殘余應力)作用下,分別用權函數法與有限元法計算在不同裂紋尺寸a/c=0.2;a/c=0.4;a/c=0.6;a/c=0.8,a/t=0.1~0.8下的殘余應力強度因子Kres,計算結果無量綱化后如表7和表8所示。

無量綱化后的結果對比如圖12所示。

圖12所示,在裂紋半長比a/c=0.2時,a/t=0.1~0.8時,權函數法與有限元法得到的最深點Kres最大差值的百分比在a/t=0.8時取得為6.2%,其余點差值均在5%以下。在a/c=0.2時,權函數法與有限元法得到的表面點Kres,兩者之間差值的百分比最大值在a/t=0.6處取得為3.8%,其余各點差值均在3%以下。

表7 半橢圓表面裂紋Kres_有限元法(a/c=0.2;a/c=0.4;a/c=0.6;a/c=0.8)Tab.7 Kresof semi-elliptical crack by FEA method

表8 半橢圓表面裂紋Kres_權函數法(a/c=0.2;a/c=0.4;a/c=0.6;a/c=0.8)Tab.8 Kresof semi-elliptical crack by weight function method

圖12 權函數法與有限元法Kres對比Fig.12 Comparison of Kresbetween FEA and WFM

在裂紋半長比a/c=0.4時,a/t=0.1~0.8時,權函數法與有限元法得到的最深點Kres最大差值的百分比在a/t=0.8時取得為9.2%,其余點差值均在3%以下。在a/c=0.2時,權函數法與有限元法得到的表面點Kres,兩者之間差值的百分比最大值在a/t=0.2處取得為2.5%,其余各點差值均在2%以下。

在裂紋半長比a/c=0.6時,a/t=0.1~0.8時,權函數法與有限元法得到的最深點Kres最大差值的百分比在a/t=0.5時取得為8.4%,其余點差值均在4%以下。在a/c=0.6時,權函數法與有限元法得到的表面點Kres,兩者之間差值的百分比最大值在a/t=0.7處取得為5.3%,其余各點差值均在5%以下。

在裂紋半長比a/c=0.8時,a/t=0.1~0.8時,權函數法與有限元法得到的最深點Kres最大差值的百分比在a/t=0.7時取得為3.6%,其余點差值均在3%以下。在a/c=0.6時,權函數法與有限元法得到的表面點Kres,兩者之間差值的百分比最大值在a/t=0.3處取得為8.3%,其余各點差值均在2%以下。

6 結論

本文在已有的權函數統一形式上,用三維有限元計算得到了焊趾處a/c=0.2;0.4;0.6;0.8,a/t=0.1~0.8之間裂紋最深點與表面點的應力強度因子,得到以下結果:

(1)利用已有的權函數統一形式,結合T型接頭焊趾處半橢圓裂紋最深點與表面點的三維有限元計算數據,得出了適用于焊趾表面半橢圓裂紋最深點和表面點的權函數。

(2)權函數的準確性,擴展到了二次和三次應力分布下。對于表面點,半長比a/c=0.2~0.8,a/t=0.2~0.8之間,權函數與FEA結果誤差幾乎均在8%以下;對于最深點,半長比a/c=0.2~0.8,a/t=0.1~0.8之間,權函數與FEA結果誤差均在在8%以下。權函數的準確性在這兒得到了有效驗證。

(3)基于新的權函數,對T型接頭焊趾處半橢圓表面裂紋的殘余應力強度因子Kres進行了計算,并與有限元計算結果進行對比,對比誤差在10%以下,表明新的權函數能很好地預測T型接頭焊趾處的殘余應力強度因子Kres。

[1]Servetti G,Zhang X.Predicting fatigue crack growth rate in a welded butt joint:The role of effective R ratio in accounting for residual stress effect[J].Engineering Fracture Mechanics,2009,76(11):1589-1602.

[2]Newman Jr J C,Raju I S.An empirical stress-intensity factor equation for the surface crack[J].Engineering Fracture Mechanics,1981,15(1):185-192.

[3]Standard B.7910:1999:Guide on methods for assessing the acceptability of flaws in fusion welded structures[J].British Standard Institution,2000.

[4]Bowness D,Lee M M K.Prediction of weld toe magnification factors for semi-elliptical cracks in T-butt joints[J].International Journal of Fatigue,2000,22(5):369-387.

[5]Bao R,Zhang X,Yahaya N A.Evaluating stress intensity factors due to weld residual stresses by the weight function and finite element methods[J].Engineering Fracture Mechanics,2010,77(13):2550-2566.

[6]Labeas G,Diamantakos I.Numerical investigation of through crack behavior under welding residual stresses[J].Engineering Fracture Mechanics,2009,76(11):1691-1702.

[7]Wang X,Lambert S B.Stress intensity factors for low aspect ratio semi-elliptical surface cracks in finite-thickness plates subjected to nonuniform stresses[J].Engineering Fracture Mechanics,1995,51(4):517-532.

[8]Wang X,Lambert S B.Stress intensity factors and weight functions for high aspect ratio semi-elliptical surface cracks in finite-thickness plates[J].Engineering Fracture Mechanics,1997,57(1):13-24.

[9]Shiratori M.Analysis of stress intensity factors for surface cracks subjected to arbitrarily distributed surface stresses[J].Transaction of JSME(Series A),1987,54(467):1828-1835.

[10]Bueckner H F.Novel principle for the computation of stress intensity factors[J].Zeitschrift fuer Angewandte Mathematik& Mechanik,1970,50(9).

[11]Rice J R.Some remarks on elastic crack-tip stress fields[J].International Journal of Solids and Structures,1972,8(6):751-758.

[12]Glinka G,Shen G.Universal features of weight functions for cracks in mode I[J].Engineering Fracture Mechanics,1991, 40(6):1135-1146.

[13]黃小平,賈貴磊,崔維成,祈恩榮.海洋鋼結構疲勞裂紋擴展預報單一擴展率曲線模型[J].船舶力學,2011,15(1-2):118-125. Huang Xiaoping,Jia Guilei,Cui Weicheng.Unique crack growth rate curve model for fatigue life prediction of marine steel structures[J].Journal of Ship Mechanics,2011,15(1-2):118-125.

[14]Miki C,Mori T,Tajima J.Effect of stress ratio and tensile residual stress on near threshold fatigue crack growth[C]//Proceedings of JSCE,1986(368):187-193.

[15]顏鳴皋,等.航空材料手冊[M].北京:中國標準出版社,2002:104-131.

Weight function for weld toe semi-elliptical surface crack and calculating residual stress intensity factors by weight function method

Xü Lei,HUANG Xiao-ping
(State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

The stress intensity factors(SIFs)of semi-elliptical surface cracks at the weld toe were analyzed by using finite element analysis(FEA).Based on the obtained stress intensity factor data of surface cracks with aspect ratio from a/c=0.2~0.8 and a/t=0.1~0.8 by FEA,new weight functions for the calculation of T-butt weld toe surface cracks SIFs at both surface and deepest point were derived under constant and linear distributed stresses.The new weight functions were extended for higher order stress distribution situation and validated through the results obtained by FEA under parabolic and cubic stress loading.The difference between weight function and FEA is less than 10%for cracks with aspect ratio a/c=0.2~0.8,a/t= 0.1~0.8 at both surface point and deepest point.The derived weight function is used to calculate the SIFs of weld toe surface crack due to residual stress,results comparison with FEA method shows new weight function can make good predication for the weld toe surface crack SIFs.

T-butt;weld toe;surface crack;stress intensity factor;weight function;residual stress

U661.4

A

10.3969/j.issn.1007-7294.2017.04.009

1007-7294(2017)04-0443-12

2016-12-27

徐磊(1989-),男,碩士生;黃小平(1963-),男,副教授,通訊作者,E-mail:xphuang@sjtu.edu.cn。

猜你喜歡
裂紋有限元
裂紋長度對焊接接頭裂紋擴展驅動力的影響
一種基于微帶天線的金屬表面裂紋的檢測
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
微裂紋區對主裂紋擴展的影響
磨削淬硬殘余應力的有限元分析
預裂紋混凝土拉壓疲勞荷載下裂紋擴展速率
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 免费又黄又爽又猛大片午夜| 秘书高跟黑色丝袜国产91在线| 久久无码av三级| 亚洲自偷自拍另类小说| 大陆精大陆国产国语精品1024| 欧美精品H在线播放| 亚洲欧美综合精品久久成人网| 国产福利一区视频| 国产一区二区三区在线观看免费| 91精品国产自产91精品资源| 青草娱乐极品免费视频| 国产手机在线ΑⅤ片无码观看| 不卡无码h在线观看| 久久美女精品国产精品亚洲| 毛片手机在线看| 久久久久无码精品国产免费| 欧美区一区二区三| 99久久国产精品无码| 亚洲国产一区在线观看| 国产va在线观看| 亚洲欧洲免费视频| 激情六月丁香婷婷| 欧美成人免费午夜全| 国产人成午夜免费看| 国产永久免费视频m3u8| 色婷婷在线影院| 国产爽歪歪免费视频在线观看| 国产精品久久久久鬼色| 99热这里只有精品免费| 精品无码日韩国产不卡av| 国产高清在线丝袜精品一区| 国产乱视频网站| 老色鬼欧美精品| 国产精品v欧美| 国产精品深爱在线| 亚洲中文字幕av无码区| 久久久久久久蜜桃| 国产女人18水真多毛片18精品| 麻豆精品久久久久久久99蜜桃| 在线va视频| 污污网站在线观看| 国产噜噜噜| 69视频国产| 欧美性精品| 免费人成在线观看视频色| 夜精品a一区二区三区| 久久久四虎成人永久免费网站| 萌白酱国产一区二区| 在线国产综合一区二区三区| 91精品aⅴ无码中文字字幕蜜桃| 国产制服丝袜91在线| 依依成人精品无v国产| 尤物在线观看乱码| 亚洲国产精品成人久久综合影院| 成人在线不卡| 99热这里只有精品免费国产| 99热这里只有免费国产精品| 少妇被粗大的猛烈进出免费视频| 国产成人盗摄精品| av一区二区无码在线| 99视频在线免费| 呦女亚洲一区精品| yy6080理论大片一级久久| 亚洲伊人天堂| 亚洲高清在线播放| 人妻丰满熟妇AV无码区| 国产精品美女网站| 啊嗯不日本网站| 看国产毛片| 亚洲综合第一区| 制服无码网站| 午夜啪啪福利| 国产99视频精品免费视频7| 一区二区三区精品视频在线观看| 无码国产伊人| 国产XXXX做受性欧美88| 久久大香香蕉国产免费网站| 91在线免费公开视频| 久久大香香蕉国产免费网站| 日韩黄色在线| 在线色国产| 青青青视频91在线 |