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

XFEM及其在正交異性鋼橋面板疲勞裂紋擴展中的應用

2018-04-13 08:12:09顧安邦杜柏松
關鍵詞:裂紋有限元

渠 昱,曾 勇,顧安邦,杜柏松

(1.重慶交通大學 山區橋梁與隧道工程國家重點實驗室培育基地,重慶 400074;2.重慶交通大學 山區橋梁結構與材料教育部工程研究中心,重慶 400074)

線彈性斷裂力學的主要任務是計算裂紋擴展長度、擴展率和能量釋放率。裂紋擴展計算比較困難,對于復雜結構更困難,它涉及斷裂力學、材料科學和計算方法等學科。在有限元(FE)中,裂紋計算要求有限單元邊界與裂紋對齊,裂紋擴展需要不斷地更新網格,諸如可視化后處理技術等計算方法需要對網格演化進行大量計算,因此采用普通有限元方法(FEM)計算裂紋擴展遇到極大挑戰,限制了有限元的應用。

1 擴展有限元(XFEM)

XFEM最有吸引力的是將不連續問題加入到有限元而又不要求網格與不連續函數及其導數一致,從而避免了重新劃分網格。在許多情況下,避免網格重劃分方法是可取的,可以減少大量的網格演化等計算[1]。XFEM是一個非常靈活的方法,其基本想法是將不連續位移場矩陣分解為2部分:uh(X)=ucont(X) +udisc(X),連續部分ucont(X)是標準有限元插值矩陣,不連續部分udisc(X)是根據局部單位分解定理[2,3]增加的部分(即富集),并將這些富集信息嵌入有限元插值。所謂富集是將額外的未知量濃縮到單元水平,這樣的矩陣仍然是稀疏的,因此可以用現有的有限元分析框架處理不連續問題。

當裂尖位移場已知,基于漸進場的富集函數可以被添加到近似解空間中。對于裂紋和位錯奇異的漸進場,可以作為裂尖區和位錯核的富集函數,從而大大減少了在子域細化網格的工作量,XFEM已成為斷裂分析最流行的方法。

這些方法在解決材料科學問題上的主要優點是對不連續現象建模簡化。在傳統的有限元方法中,劃分網格時必須使單元的邊緣/面與裂紋重合、在裂紋每一側布置節點,才能考慮材料沿裂紋面的分離。這樣的網格創建相當困難,特別是三維問題,因為網格還必須考慮到模型的其它特性,如晶界或夾雜物。在XFEM中,沿裂紋面不連續位移場的引入是通過添加基函數來創建近似解,擺脫了網格的約束。當XFEM與水平集方法LSM(level-set method)相結合,就可以在原來有限元節點上用節點值完整表達諸如裂紋的幾何形狀和位移場等特性。同樣,相邊界和其他不連續性的模擬,也可以用這種方法。特別是當幾何形狀演化、裂紋擴展或移動相邊界,這些優點是非常重要的。

1.1 XFEM中的裂紋的描述

裂紋采用2個水平集函數描述[3-5]:

1) 描述裂紋路徑Γα的水平集函數f(x),即有向距離函數:

(1)

2) 定義裂紋位置的水平集函數g(x),原點在裂尖,方向與裂紋垂直。該函數定義的坐標系見圖1(b),坐標系中的r和θ定義為:

圖1 富集單元、過渡單元和水平集函數Fig. 1 Enrichment unit,transition unit and horizontal set function

與一般閉合面相比,裂紋屬于開放界面,2個水平集函數f(x,t)和g(x,t)分別描述裂紋面和裂尖。點x位于裂紋面上的條件是f(x,t)=0,g(x,t) > 0;位于裂尖的條件是f(x,t)=0,g(x,t)=0。裂紋問題的另一個特點是,某位置一旦成為裂紋的一部分,則一直保持為裂紋的一部分。因此裂紋擴展時,只需要更新水平集函數f(x,t)以及g(x,t) 在裂尖的位置,而f(x,t) 在裂尖后面的位置保持不變。等于0的水平集稱為零水平集,用于描述裂紋面。裂尖后面的零水平集不必更新,從而避免了差分方程的麻煩。在界面的法向方向移動,水平集函數不需要隨界面運動速度更新,因此,不使用Hamilton-Jacobi方程更新;反之,水平集函數隨裂尖的速度運動需要更新。

水平集方法簡單、有效和便于實現,XFEM與水平集相結合,功能更強大,因為這樣描述裂紋、位錯、晶粒邊界和位相界面可以嚴格用節點表達,而且能利用水平集這個強大的工具跟蹤這些表面的演化。

1.2 XFEM中的富集函數

在XFEM中,單元被裂紋切割有2種情況:① 被裂紋完全切割的單元稱為裂紋體單元;② 含有裂尖的單元稱為裂尖單元。裂紋體單元節點富集采用有向距離函數式(2),裂尖富集函數是基于裂尖場的漸進解。對于線性斷裂力學,取分支函數(整體函數)為富集函數[6]:

(2)

式中:r、θ為原點定義在裂尖的局部極坐標。

1.3 XFEM中的過渡單元問題

與總體富集方法相比,T. BELYTSCHKO等[3,4]介紹的富集函數局部單位分解能減少未知數量,改進數值計算條件。

J. CHESSA等[7]認為:在富集和非富集單元之間存在過渡單元產生多余項問題,近似解空間的多余項是引起誤差的根源。對于采用有向距離函數式(2)作為富集函數的過渡單元邊上,富集函數為0,因此不會產生多余項。但是對于裂尖富集,XFEM的過渡單元〔圖1(a)中交叉線所示〕 中只有2個節點被富集,形函數就不再滿足單位分解條件[8],過渡單元會產生多余項,從而降低數值算法的收斂性。

例如:如果富集節點集I*中只有一個過渡單元,這個單元的近似解為:

(3)

N. SUKUMAR等[8]發現,采用背脊型的富集函數,過渡單元會產生誤差,從而降低數值算法的收斂性。修正的背脊型富集函數能改進其收斂效率。J. CHESSA等[7]通過進一步的研究,提出假設應變法以消除過渡單元的多余項。這個方法的缺點是需要特殊的假定應變單元,而該單元又取決于富集函數和單元類型,需要逐點構建。富集函數構造困難,且不具有通用性。

對于裂尖單元,過渡單元的多余項限制了局部單位分解精度。T. P. FRIES[9]提出了斜坡函數法來解決過渡單元多余項的影響。其主要缺點是:① 裂尖過渡單元富集節點為16個,多于XFEM的4個;② 自由度大約是XFEM的3倍;③有缺秩問題,富集矩陣的秩是16,而線性無關的條件只有14。

2 不連續伽遼金XFEM(DG-XFEM)

另一個避免過渡單元多余項的方法是利用內罰函數方法[10],使富集和不富集區域之間強制連續,稱為DG-XFEM[11]。該方法比標準XFEM更精確,收斂速度得到了優化。

在基于分片區域的DG-XFEM方法中,域被分成幾個不互相重疊的區域,對這些區域分別進行富集并用內罰函數強制區域之間連續。該方法的最佳特性是富集是局部的,且所有的形函數滿足單位分解。為簡單計,筆者只論述平面二維的情況,該方法同樣適用于三維平面裂紋。

圖2 富集區域劃分示意Fig. 2 Enrichment region division

(4)

覆蓋重疊之間用內罰函數強制連續。通過弱不連續變分,得到過渡單元的剛度矩陣

(KXFEM+KDG)d=fext

(5)

式中:KXFEM為每一個覆蓋的標準XFEM剛度矩陣的集成;KDG為重疊覆蓋部分強制連續的不連續伽遼金(DG)項,是剛度的修正項;fext為節點力。

域Ω被分成np個互相不重疊的覆蓋(Ωb)P,b=1~np,b覆蓋中I節點向量為:

(6)

(7)

(8)

式中:a為罰函數;因子μ=1;dT為節點自由度向量。

采用適度的內點罰函數值,方程(5)受損不明顯,線性系統方程組的迭代方法解仍然有效:

(9)

(10)

(11)

(12)

從而,可由沒有富集的有限元、完全富集的XFEM和過渡單元構成求解域的離散方程(5)。

DG-XFEM驗證:選擇橋梁鋼板Q370qE、單邊拉伸的邊裂紋情況對DG-XFEM進行驗證。

鋼板長L=5 m,寬W=1.5 m,板厚d=0.01 m。利用XFEM計算得到沿裂紋面的應力sxx和syy分布情況,并與Westergaard應力函數[12]計算的理論結果進行比較(圖3)。可以看出,二者符合很好,在過渡單元沒有出現數值不收斂的情況。

圖3 鋼板邊裂紋計算的沿裂紋面應力分布及理論值比較Fig. 3 Comparison between the calculated crack surfaces stress distribution and theoretical values of the steel plate with side crack

3 DG-XFEM在正交異性鋼橋面板疲勞裂紋擴展中的應用

3.1 模型試驗

以某長江大橋為試驗背景。該橋為主跨長788 m的鋼箱梁懸索橋,吊索間距16 m,加勁梁為流線型扁平鋼箱梁,梁高3.5 m,主橋寬27.5 m;每一節段順橋向5個橫隔板,間距2.7 m,橫向共40個U肋,間距0.6 m。

在鋼橋面板模型設計過程中,按彈性支撐連續梁、利用有限元對實橋整體進行了仿真模擬。根據整體分析,選擇長16 m標準節段的鋼箱梁,按彈性簡支進行有限元分析。根據節段鋼箱梁的分析結果和正交異性板各疲勞細節僅對輪載敏感的局部特性[13],按照1∶2的比例尺建立簡化面板節段模型,橫向包括7個閉口U肋,縱向3個橫隔板,按縱向簡支、橫向僅限制向下的位移條件進行試驗[14]。

試驗的疲勞荷載是根據三峽庫區的實際交通狀況,按照JTG D64—2015《公路鋼結構橋梁設計規范》[13],參考AASHTO規范[15]的單車模型,考慮1.15的沖擊系數,按1∶2比例模型和等效應力脈換算得到軸載47.7 kN,在1 000 kN的MTS機上進行。二期恒載15 kN,汽車活載47.7 kN,均通過MTS機加載。考慮鋪裝效應,荷載通過分力橋施加。

在試驗中,當加載2.0×106次循環時,6 # 弧形開孔出現7.5 mm長的裂紋,2.2×106、2.4×106次循環加載時裂紋分別擴展到12、21 mm,在2.6×106次終止試驗時裂紋擴展到31 mm(圖4)。

圖4 6 # 開孔處疲勞裂紋試驗照片Fig. 4 Photo of fatigue crack in right side of 6# hole

3.2 數值模型

根據試驗模型的尺寸和約束條件,利用Abaqus平臺創建3-D不連續有限元模型,面板節段模型見圖5,劃分節點383 344個,C3D8R單元190 791個,裂紋體單元采用有向距離函數富集,裂尖單元采用分支函數富集,富集單元和過渡單元的邊界采用內點罰函數方法強制連續,采用DG-XFEM進行計算。

圖5 鋼箱梁節段面板模型Fig. 5 Segment deck model of steel box girder

在有限元分析中假定材料為理想線彈性的,14MnNbq橋梁鋼板材,彈性模量E=2.1×105MPa,屈服強度σy≥370 MPa,斷裂強度σb≥530 MPa,泊松比為ν=0.3,V型沖擊功JV≥41 J,粘性裂紋分片連續,基于損傷演化方法進行數值分析。

損傷判據為最大主應力失效準則。正交異性鋼橋面板各個疲勞細節的常幅疲勞極限不盡相同,裂紋出現在橫隔板挖孔處,疲勞細節的常幅疲勞極限為70 MPa[13],超過該值即產生裂紋。

根據損傷力學原理,裂紋在開裂和擴展過程中,在裂尖及其周圍的材料處于塑形狀態,再應用彈性材料計算就不合適了,因而采用塑形材料參數進行數值模型設置。在模型試驗中裂紋萌生和擴展采用的是應力和加載次數來體現,而在數值模擬中,加載次數不容易實現,所以采用能量原理,用法向斷裂能42.2 kN/m判斷裂紋是否擴展,擴展方向為切應力最大的方向。損傷穩定系數DS=5.0×10-5,用于判斷裂紋是否止裂。

目前裂紋擴展研究還必須將裂紋分成微觀階段和宏觀階段2個階段:

1) 疲勞是一種對微觀結構非常敏感的現象。由于微結構的易變性,裂紋萌生階段的長度與晶粒主尺寸相當,以張開(Ⅰ型裂紋) 和滑移(Ⅱ、Ⅲ型裂紋)相結合的形式生長,沿最大分切應力滑移系統擴展。最大分切應力逐點不同,擴展路徑彎曲,因此不穩定、不規則,不能用斷裂力學研究。材料的生產、制造和熱處理過程的結果,使測量的損傷趨勢具有相當大的不確定性[16],即使對于相同的材料也是如此。

2) 宏觀裂紋擴展不受材料微結構的影響,可以用Paris定律表達。要計算裂紋擴展,首先必須確定初始裂紋長度。M.H.EL-HADDAD等[17]建議使用長裂紋應力強度因子門檻值ΔKth和虛擬裂紋長度a0表達疲勞極限σf:

(13)

式中:Y為幾何形狀系數,取決于裂紋的幾何形態,對于穿透裂紋(長度2a)無限平板Y=1.0,其它裂紋類型,Y值可以從應力強度因子求解。

式(12)中的虛擬裂紋長度a0即為相當初始裂紋長度。相當初始裂紋是指滿足Paris定律的最小宏觀裂紋,是為了避免微裂紋擴展不穩定和不規則問題而提出的。工程中常根據經驗假定初始裂紋長度,如金屬假定1~3 mm量級,或者利用無損探傷的結果來確定。對于橋梁鋼Q370qE,根據文獻[18]:

ΔKth=5.556(1-0.825R)1.147

(14)

式中:R為加載應力比。

由文中試驗應力比R=0.221 5,得到ΔKth=4.408 MPa·m0.5。根據式(11)得到初始裂紋長度為1.26 mm,以此作為文中裂紋擴展的相當初始裂紋長度。

根據模型試驗結果,利用所建立的有限元模型進行整體分析,得到的6#孔左側(PEL) 和右側(PER) 周邊單元主應力分布,如圖6。

圖6 橫隔板6 # 開孔周邊單元主應力分布Fig. 6 Principal stress distribution of perimeter elements of 6# hole at diaphragm

從圖6可以看出,橫隔板6 # 開孔最大主應力超過60 MPa,熱點應力系數取1.3進行計算[19],主應力為78 MPa,超過該疲勞細節的常幅疲勞極限70 MPa。考慮到應力集中、焊接缺陷和焊接殘余應力的影響,正交異性鋼橋面板橫隔板開孔區是疲勞裂紋多發區域,為對疲勞最敏感的區域。根據李傳習等[19]對某扁平鋼箱梁大橋裂紋統計研究,出現在橫隔板開孔區域和U肋與橫隔板焊接端頭的裂紋占全橋裂紋總數的89.4%。

根據數值分析結果,繪制5 #、6 # 開孔之間的橫隔板在沒有裂紋情況下的主應力等值線,如圖7。

圖7 橫隔板5#、6#孔之間的無裂紋主應力等值線Fig. 7 Principal stress isocline between 5# and 6# holes calculated by XFEM without crack

從圖7可以看出:

1) 開孔周邊的應力高、裂紋產生處的應力超過該疲勞細節的常幅疲勞極限,應力梯度陡,越往兩孔中間應力越低。

2) 兩開孔之間的水平區域均為高應力區,裂紋在此區域產生并沿水平方向擴展,以保證裂紋尖端獲取足夠的能量維持擴展,隨著應力強度減弱,裂紋擴展逐漸停止。

3) 裂紋在初始階段垂直于應力等值線、沿最大主應力方向向板內擴展,這與裂紋擴展的基本理論一致。

從圖4、圖8可以看出,2.6×106次循環試驗結束時裂紋擴展長度為31 mm(圖4),計算值為32.3 mm〔圖8(a)〕,二者相差很小。數值模擬中裂紋出現的位置、擴展路徑和長度與試驗結果符合很好。也與實橋橫隔板中出現的裂紋基本一致〔圖8(b)〕[20]。

圖8 橫隔板開孔處計算的裂紋擴展與實橋裂紋擴展比較Fig. 8 Comparison of crack growth calculated at diaphragm and that of real bridge

4 結 論

1) 利用平板單邊裂紋的應力函數理論解對DG-XFEM進行驗證,結果表明DG-XFEM能很好地解決過渡單元引起的數值收斂效率問題。

2) 利用XFEM得到的應力場,分析表明:弧形開孔處是高應力區,與幾何不連續造成的應力集中、焊接缺陷和殘余應力的作用疊加,形成疲勞裂紋的高發區域,是正交異性鋼橋面板對疲勞最敏感的部位。

3) 利用DG-XFEM模擬裂紋擴展,分析表明:裂紋擴展的初始階段與主應力等值線垂直,且沿水平方向的高應力區向板內擴展,以保證裂尖獲取足夠的能量維持裂紋擴展,隨著應力強度減弱,裂紋擴展逐漸停止,裂紋擴展路徑與試驗和實橋出現的情況一致,也與裂紋擴展的基本理論一致,證明了利用DG-XFEM對正交異性鋼橋面板疲勞裂紋擴展分析的有效性和合理性。

參考文獻(References):

[1] RABCZUK T.Computational methods for fracture in brittle and quasi-brittle solids:State-of-the-art review and future perspectives[J].ISRNAppliedMathematics,2013,2013:1-38.

[2] ALVES M K,ROSSI R.An extension of the partition of unity finite element method[J].JournaloftheBrazilianSocietyofMechanicalSciencesandEngineering,2005,27(3):209-216.

[3] BELYTSCHKO T,GRACIE R,VENTURA G.A review of extended/generalized finite element methods for material modeling[J/OL].ModellingandSimulationinMaterialsScienceandEngineering,2009,17(4):043001.doi:10.1088/0965-0393/17/4/043001.[2017-05-12] http://iopscience.iop.org/article/10.1088/0965-0393/17/4/043001/meta.

[4] BELYTSCHKO T,MOЁS N,USUI S,et al.Arbitrary discontinuities in finite elements[J].InternationalJournalforNumericalMethodsinEngineering,2001,50(4):993-1013.

[5] MOЁS N,GRAVOUIL A,BELYTSCHKO T.Non-planar 3D crack growth by the extended finite element and level sets-Part I:Mechanical model[J].InternationalJournalforNumericalMethodsinEngineering,2002,53(11):2549-2568.

[6] ABDELAZIZ Y,HAMOUINE A.A survey of the extended finite element[J].Computers&Structures,2008,86(11):1141-1151.

[7] CHESSA J,WANG Hongwu,BELYTSCHKO T.On the construction of blending elements for local partition of unity enriched finite elements[J/OL].InternationalJournalforNumericalMethodsinEngineering,2003,57(7):1015-1038.http://onlinelibrary.wiley.com/doi/10.1002/nme.777/full.

[8] SUKUMAR N,CHOPP D L,MOЁS N,et al.Modeling holes and inclusions by level sets in the extended finite element method[J].ComputerMethodsinAppliedMechanicsandEngineering,2001,190(46):6183-6200.

[9] FRIES T P.A corrected XFEM approximation without problems in blending elements[J/OL].InternationalJournalforNumericalMethodsinEngineering,2008,75(5):503-532.http://onlinelibrary.wiley.com/doi/10.1002/nme.2259/full.

[10] ARNOLD D N.An interior penalty finite element method with discontinuous elements[J].SIAMJournalonNumericalAnalysis,1982,19(4):742-760.

[11] GRACIE R,WANG Hongwu,BELYTSCHKO T.Blending in the extended finite element method by discontinuous Galerkin and assumed strain methods[J/OL].InternationalJournalforNumericalMethodsinEngineering,2008,74(11):1645-1669.http://onlinelibrary.wiley.com/doi/10.1002/nme.2217/full.

[12] Sanford R J.Principles of Fracture Mechanics[J].Upper Saddle River Nj,2002,12:29-44.

[13] 中交規劃設計研究院有限公司.公路鋼結構橋梁設計規范:JTG D64—2015[S].北京:人民交通出版社股份有限公司,2015.

CCCC Highway Consultants CO.,Ltd.SpecificationforDesignofHighwaySteelBridge:JTGD64-2015[S].Beijing:China Communications Press Co.,Ltd,2015.

[14] 曾勇,渠昱,譚紅梅,等.大跨懸索橋主纜-吊桿-鋼箱梁一體模型疲勞試驗研究[J].四川大學學報(工程科學版),2016,48(4):78-84.

ZENG Yong,QU Yu,TAN Hongmei,et al.Fatigue experimental study on integration model including main cable,hangers and box main girder of long-span suspension bridge[J].JournalofSichuanUniversity:EngineeringScienceEdition,2016,48(4):78-84.

[15] The American Association of State Highway and Transportation Officials(AASHTO).LRFDBridgeDesignSpecifications[S].4thEd.Washington,D C:AASHTO,2007.

[16] CHOWDHURY P,SEHITOGLU H.Mechanisms of fatigue crack growth-a critical digest of theoretical developments[J].Fatigue&FractureofEngineeringMaterials&Structures,2016,39(6):652-674.

[17] EL-HADDAD M H,TOPPER T H,SMITH K N.Prediction of non-propagating cracks[J].EngineeringFractureMechanics,1979,11(3):573-584.

[18] 劉艷萍,陳傳堯,李建兵,等.14MnNbq焊接橋梁鋼的疲勞裂紋擴展行為研究[J].工程力學,2008,25(4):209-213.

LIU Yanping,CHEN Chuanyao,LI Jianbing,et al.Fatigue crack growth behavior for the weld heat-affected zone of 14MnNBQ bridge steel[J].EngineeringMechanics,2008,25(4):209-213.

[19] NIEMI E,FRICKE W,MADDOX S J,et al.FatigueAnalysisofWeldedComponents——Designer’sGuidetotheStructuralHot-SpotStressApproach:IIW-1430-00[S].Cambridge:Woodhead Publishing Limited,2006.

[20] 李傳習,李游,陳卓異,等.鋼箱梁橫隔板疲勞開裂原因及補強細節研究[J].中國公路學報,2017,30(3):121-131.

LI Chuanxi,LI You,CHEN Zhuoyi,et al.Fatigue cracking reason and detail dimension of reinforcement about transverse diaphragm of steel box bridge[J].ChinaJournalofHighwayandTransport,2017,30(3):121-131.

猜你喜歡
裂紋有限元
裂紋長度對焊接接頭裂紋擴展驅動力的影響
一種基于微帶天線的金屬表面裂紋的檢測
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
微裂紋區對主裂紋擴展的影響
磨削淬硬殘余應力的有限元分析
預裂紋混凝土拉壓疲勞荷載下裂紋擴展速率
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 亚洲国内精品自在自线官| 制服丝袜一区| 国外欧美一区另类中文字幕| 伊人网址在线| 一级香蕉视频在线观看| 在线看AV天堂| 国产亚洲欧美在线人成aaaa| 久久综合九色综合97婷婷| 中文无码毛片又爽又刺激| 日韩欧美国产区| 不卡无码h在线观看| 国产欧美在线观看精品一区污| 中文字幕免费播放| 成人一级黄色毛片| 国产欧美视频综合二区| 亚洲视频在线青青| 福利小视频在线播放| 福利在线不卡一区| 日韩午夜福利在线观看| 午夜不卡视频| 国产人成午夜免费看| 久久视精品| 91蝌蚪视频在线观看| 国产精品视频免费网站| 成人在线天堂| 国产精品露脸视频| 欧美自慰一级看片免费| 久久中文电影| 性欧美精品xxxx| 91小视频在线观看| 无码AV日韩一二三区| 国产激爽大片高清在线观看| 欧美a级完整在线观看| 国产91小视频| 成人在线综合| 日本成人福利视频| 国产精品9| 不卡午夜视频| 91成人在线观看| 99精品免费欧美成人小视频| 国产成人凹凸视频在线| 中文成人在线视频| 精品无码人妻一区二区| 亚洲天堂网视频| 精品久久国产综合精麻豆| 日韩资源站| 亚洲欧美综合在线观看| 97久久人人超碰国产精品| 九九久久99精品| 毛片免费高清免费| 婷婷六月激情综合一区| 欧美日韩导航| 亚洲精品日产AⅤ| 911亚洲精品| 色视频国产| 亚洲天堂网在线播放| 黄色一及毛片| 午夜福利无码一区二区| 亚洲大尺度在线| 国产情精品嫩草影院88av| 88av在线| 91精品国产自产在线观看| 波多野结衣在线一区二区| 露脸一二三区国语对白| 亚洲人成网站色7777| 久久国产黑丝袜视频| 免费看久久精品99| 欧美精品1区2区| 亚洲中文无码av永久伊人| 特级aaaaaaaaa毛片免费视频| 天天综合网站| 综合人妻久久一区二区精品| 不卡网亚洲无码| 99激情网| 自拍偷拍欧美日韩| 国产日韩丝袜一二三区| 好久久免费视频高清| 一级片一区| 精品人妻AV区| 国产农村1级毛片| 国产人人射| 亚洲人成影院午夜网站|