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

高頻諧振載荷作用下疲勞裂紋尖端變形場分析

2015-02-19 00:24:22高紅俐鄭歡斌
浙江工業大學學報 2015年2期

高紅俐,鄭歡斌,劉 歡,劉 輝

(浙江工業大學 特種裝備制造與先進加工技術教育部/浙江省重點實驗室,浙江 杭州 310014)

高頻諧振載荷作用下疲勞裂紋尖端變形場分析

高紅俐,鄭歡斌,劉歡,劉輝

(浙江工業大學 特種裝備制造與先進加工技術教育部/浙江省重點實驗室,浙江 杭州 310014)

摘要:為研究高頻諧振式疲勞裂紋擴展試驗中,帶有I型預制裂紋的CT緊湊拉伸試件裂紋尖端位移、應變場的變化規律,利用動態有限元方法,采用ANSYS和MATLAB軟件編寫程序,計算了CT試件在高頻恒幅正弦交變載荷作用下,在一個應力循環及裂紋擴展到不同長度時裂紋尖端區域的位移、應變場并分析了其變化規律.為驗證有限元計算結果的準確性,進行了高頻諧振式疲勞裂紋擴展試驗,采用動態高精度應變儀測量了CT試件疲勞裂紋擴展到5 mm時在一個應力循環內裂紋尖端點的應變.研究結果表明:在穩態裂紋擴展階段,高頻諧振載荷作用下I型疲勞裂紋尖端區域位移、應變均為和載荷同一形式的交變量;隨著裂紋的擴展,I型疲勞裂紋尖端的位移、應變幅不斷增大;裂紋尖端測量點應變有限元計算結果和實驗結果最大誤差為2.93%。

關鍵詞:高頻諧振載荷;疲勞裂紋尖端;動態有限元;位移場;應變場

Finite element analysis of displacement and strain fields at crack

tip under high frequency resonant loading

GAO Hongli, ZHENG Huanbin, LIU Huan, LIU Hui

(Key Laboratory of E&M, Ministry of Education&Zhejiang Province, Zhejiang University of

Technology, Hangzhou 310014, China)

Abstract:This work researched the variation law of displacement field and strain field at fatigue crack tip of compact tension specimen with type I pre-notch based on dynamic finite element method (FEM) in the high frequency resonant fatigue crack propagation test. The displacement field and the strain field at CT specimen fatigue crack tip in one stress cycle and at different crack lengths under constant amplitude high frequency sinusoidal alternating loading condition are calculated and the related variation laws of the displacement field and the strain field are analyzed. In order to validate the FEM result, the high frequency resonant fatigue crack propagation test was performed, and the dynamic strain gauge is used to measure the strain at crack tip of CT specimen with crack length of 5mm during one stress cycle. The research results show that during crack stable propagation stage, the displacement and the strain at type I fatigue crack tip are the same form with the high frequency resonant load; the displacement and strain increase with the crack growth. The maximum error of the strain at crack tip between the calculated result by FEM simulation and the experimental result is 2.93%。

Keywords:high frequency resonant load; fatigue crack tip; dynamic FEM; displacement field; strain field

疲勞破壞是機械零部件失效的主要形式,材料或構件在交變載荷的作用下,疲勞裂紋的產生和擴展是造成疲勞破壞的主要原因.由于目前尚不能完全通過有效的理論方法來研究,采用特定的材料疲勞裂紋擴展試驗,來探索零部件擴展斷裂過程的規律,對提高機械產品的可靠性和使用壽命有著十分重要的意義[1-3].實際工程構件中裂紋形式大多屬于I型張開裂紋,也是最危險的一種裂紋形式,疲勞裂紋擴展試驗采用具有I型預制裂紋的標準CT試件來進行試驗,通過測量試件在交變載荷作用下,疲勞裂紋擴展速率da/dN和裂紋尖端應力強度因子幅ΔK之間的關系來研究材料的疲勞裂紋擴展規律,另外,裂紋的萌生、擴展和斷裂不能完全由宏觀裂紋形態和尺寸a來表征,要系統全面研究疲勞裂紋擴展機理,除了研究宏觀裂紋的擴展規律da/dN和ΔK的關系,研究疲勞裂紋擴展過程中裂紋尖端位移、應變、應力場等力學特征參數的變化規律是非常必要的。

疲勞裂紋擴展試驗[4-5]包括基于電液式強迫振動系統的低頻疲勞試驗和基于電磁諧振式振動系統的高頻疲勞試驗,前者,試驗振動頻率一般為1~10 Hz,后者則是基于共振原理的用于測定金屬材料及其構件在高頻諧振載荷作用下疲勞特性的測試裝置,試驗振動頻率一般為80~300 Hz左右,由于其工作頻率高、能量消耗低、試驗時間短以及試驗波形好等優點被力學實驗室廣泛用來進行材料疲勞試驗.電液式強迫振動疲勞裂紋擴展試驗裂紋尖端位移、應變和應力場的計算可以采用靜力學的方法進行計算,但在采用標準CT試件進行的高頻諧振式疲勞裂紋擴展試驗中,試件在高頻諧振載荷作用下高速振動,由于要考慮慣性效應和應力波傳播效應[6-7],使得裂紋尖端位移、應變的計算問題變得更加復雜,其解析解至今還未得到.多年來,已有眾多學者對裂紋尖端位移、應變場進行了實驗和理論研究[8-10],但基本是在靜態和準靜態的情況下,筆者針對高頻諧振式疲勞裂紋擴展試驗中的標準試件CT試件進行了動態有限元計算,研究了在恒幅高頻正弦交變載荷的作用下疲勞裂紋擴展過程中I型裂紋尖端位移、應變場的變化規律,為進一步研究高頻諧振載荷作用下疲勞裂紋擴展機理和擴展參數的測量奠定了理論基礎。

1高頻諧振式疲勞裂紋擴展試驗

電磁諧振式高頻疲勞試驗機是基于共振原理的用于測定金屬材料及其構件在高頻諧振載荷作用下疲勞特性的測試裝置,目前國內主流機型結構PLG-100,如圖1所示,其為一雙自由度的諧振式振動系統[2],當電磁激振器的激振頻率與試驗機振動系統本身的固有頻率基本一致時,試驗機振動系統便發生共振,在共振狀態下工作臺的振動振幅將增加數倍,這樣作用在CT試件上的正弦交變載荷的幅值也隨之增加,使試驗在功率消耗很小的情況下進行.圖2為疲勞裂紋擴展試驗標準CT試件尺寸圖,圖3為疲勞裂紋擴展試驗CT試件和夾具安裝圖,圖4為諧振式疲勞裂紋擴展試驗正弦交變載荷,其中Fmax為最大載荷,Fmin為最小載荷,Fm為平均載荷,Fa為振幅,諧振式疲勞裂紋擴展試驗在一個應力循環內,CT試件均處于拉應力狀態,這樣在圖3的安裝方式下,可認為CT試件的上加載孔的上半圓柱面和試件的下加載孔的下半圓柱面與定位銷軸面接觸,傳遞載荷到試件上。

圖1 電磁諧振式疲勞試驗機主機結構圖Fig.1 The structure of resonant fatigue tester

圖2 標準CT試件尺寸圖Fig.2 The dimensioned drawing of CT specimen

采用圖1所示試驗裝置進行了多種材料的標準CT試件諧振式疲勞裂紋擴展試驗,試驗結果表明:在穩態裂紋擴展階段,系統的諧振頻率范圍在90~135 Hz,這樣為簡化計算過程,在研究高頻諧振載荷作用下疲勞裂紋尖端位移、應變的變化規律時采用125 Hz的正弦交變載荷進行計算。

圖3 疲勞裂紋擴展試驗CT試件和夾具安裝圖Fig.3 The installation drawing of CT specimen

圖4 諧振式疲勞裂紋擴展試驗交變載荷Fig.4 Sinusoidal alternating load

2高頻諧振載荷作用下CT試件有限元建模

2.1疲勞裂紋尖端動態位移、應變場計算方法

對CT試件在高頻諧振載荷作用下裂紋尖端的位移、應變場進行計算,試件上所施加的交變載荷如圖4所示,它是靜載荷和正弦交變載荷共同作用的結果.當CT試件只受靜載荷作用時,可以采用ANSYS中的靜力學分析計算其位移及應變;當CT試件只受正弦載荷作用時,可以對其進行諧響應分析.但是當靜載荷和正弦載荷兩者共同作用在CT試件上,單獨采用靜力學分析或單獨采用諧響應分析都無法求解[11]。

根據振動理論,求在隨時間變化載荷作用下的響應,實際上就是求解系統動力學方程:

(1)

式中:M為系統的質量矩陣;δ為系統的總體位移列陣;C為系統的阻尼矩陣;K為系統的總體剛度矩陣;F(t)為系統的總體載荷列陣,F(t)=f(t)+f0,其中f(t)為正弦載荷,f0為靜載荷.即

(2)

將方程式(2)分解成兩個方程

(3)

(4)

Kδ=f0

(5)

假設δ1為方程式(3)的解,δ2為方程式(5)的解,由于質量矩陣M,阻尼矩陣C及剛度矩陣K均為常量矩陣,顯然δ=δ1+δ2為方程式(2)的解.所以可以利用有限元軟件ANSYS,對CT試件在正弦載荷作用和靜載荷作用下分別進行計算分析,可以得到相應的位移解δ1和δ2.ANSYS中直接給出的是解的實部δ1r,δ2r和虛部δ1i.δ1由實部δ1r和虛部δ1i組成,用極坐標形式表示復數δ1,可得δ1=δ1rcoswt+δ1isinwt,δ2僅由實部δ2r組成.將δ1和δ2進行迭加得到δ,并用時間函數形式表示為

δ=δ1+δ2=δ1rcoswt+δ1isinwt+δ2r

(6)

式中:w為正弦載荷f(t)的頻率.按方程式(6)就可以得出CT試件各節點任意時刻的位移,即動態位移.則該時刻相應的應變,表示為

ε=Bδe

(7)

式中:ε為單元內任意一點的應變;B為單元幾何矩陣;δe為單元節點位移列陣。

2.2裂紋尖端奇異性處理

由斷裂力學的理論知道,在裂紋尖端附近具有應力奇異性,即在裂紋尖端附近各應力分量都與r-1/2成正比關系,當r→0時,應力會急劇增加.在常規的有限元法中,一般均采用多項式來表示單元內部應力和位移,但在奇異點附近卻不能很好地反映出應力的變化.為了解決這個困難,一般采取以下兩種辦法:一種是將裂紋尖端附近的網格劃分得非常細,但采取這種辦法,節點會很多,計算量也會相應增大;另一種是在裂紋尖端處設置一種特殊的奇異單元,用以更好地反映出應力場的奇異性,在這些特殊奇異單元的外面,仍采用常規的單元.采用奇異等參數單元(1/4邊中點法),即把8結點等參數單元的邊中節點從正常位置移至1/4邊長處,在角點附近即出現r-1/2級的應力奇異性如圖5所示.在裂紋尖端處,布置如圖6所示的4個奇異等參數單元,即能較好地反映裂紋尖端附近的應力場。

圖5 應力奇異性Fig.5 The stress singularity

圖6 布置4個奇異等參數單元Fig.6 Four singularity isoparametric elements at crack tip

2.3疲勞裂紋尖端有限元網格劃分及約束

圖2為含預制裂紋的緊湊拉伸試件,厚度為12.5mm,裂紋長度為5mm,該試件材料為45#鋼,彈性模量E為206GPa,泊松比為0.27,密度為7 800kg/m3.采用逐節點直接建模方法:先生成節點,再由節點生成單元和結構.使用ANSYS建立有限元模型時,首先選用平面8結點四邊形等參單元PLANE82建立二維模型,劃分網格時,在裂紋尖端通過ANSYS命令KSCON將其在裂紋尖端附近退化成6結點三角形等參單元,并將有關邊的邊中結點向裂紋尖端靠攏,都移置1/4邊長處,以圓周的方式布置在裂紋尖端周圍.然后再使用SOLID95單元拉伸成體單元,網格劃分如圖7所示.將CT試件上孔的上半圓柱面設置為固定狀態,在試件的下孔的下半圓柱面設置為受力面。

圖7 CT試件三維模型網格劃分Fig.7 The three-dimensional model mesh

3高頻諧振載荷作用下CT試件裂紋尖端變形場

3.1疲勞裂紋無擴展時CT試件變形場的變化規律

3.1.1CT試件裂紋尖端位移場的變化規律

研究具有一定長度疲勞裂紋的CT試件在高頻諧振載荷一個應力周期下,在裂紋沒有擴展時裂紋尖端區域位移場的變化規律.首先以疲勞裂紋長度為5mm的CT試件為例進行計算,得到裂紋尖端區域位移場的變化規律,然后對不同長度疲勞裂紋的情況也進行了相應計算,對計算結果進行了驗證.試件尺寸及裂紋尖端計算區域如圖2所示.對其分別進行靜力學分析和諧響應分析,所施加的載荷為正弦載荷f(t)=(1.35sin785t)kN,靜載荷f0=7kN.位移場的計算方法與有限元建模方法如2.1和2.3所述.為了更為直觀地觀察裂紋尖端位移場的變化,采用自編的Matlab程序進行有限元計算結果的后處理,得到裂紋尖端矩形區域位移場的三維圖及裂紋尖端處位移的二維圖。

圖8 疲勞裂紋無擴展時裂紋尖端位移場Fig.8 The crack tip region displacement fields

圖9 最大處、最小處及裂紋尖端處位移Fig.9 Maximum, minimum and crack tip displacements

圖8中九組圖直觀地表現了在一個應力循環周期8ms內,裂紋尖端區域Y向位移場的變化.位移場中最大位移值都位于如圖2所示矩形區域D點附近,最小位移值都位于如圖2所示矩形區域F點附近.造成這種現象的原因是因為矩形區域中D點離載荷施加處最近,F點離載荷施加處最遠.圖9體現的是在一個應力循環周期內,裂紋尖端區域內D處最大位移值,F處最小位移值以及裂紋尖端處位移值隨時間的變化.可以看出:裂紋尖端區域內最大位移,最小位移以及裂紋尖端處位移是正弦變化的,與施加的載荷有著相同的變化規律。

3.1.2CT試件裂紋尖端應變場的變化規律

與求解裂紋尖端區域位移場相同的方法,得到在CT試件裂紋尖端區域一個應力周期內,不同時刻的應變場以及裂紋尖端處應變的變化.圖10中九組圖直觀地表現了在一個應力循環周期內,裂紋尖端區域應變場的變化.應變場中的最大應變值都在裂紋尖端處,這是因為裂紋尖端處是應力集中位置.應變場的形狀大抵相似,在裂紋尖端矩形區域中沿裂紋擴展方向基本呈對稱分布,裂紋尖端處應變值出現明顯的突變,應變值達到最大值,沿裂紋擴展方向應變值逐漸減小.圖11體現了在一個應力循環內,在裂紋無擴展時,裂紋尖端處應變隨時間的變化.可以看出:裂紋尖端處的應變也是正弦變化的,與載荷有著相同的規律變化,計算表明:裂紋尖端區域其他點的應變值具有和裂紋尖端點同樣的變化規律。

圖10 疲勞裂紋無擴展時裂紋尖端應變場Fig.10 The crack tip region strain fields

圖11 裂紋尖端處應變Fig.11 The crack tip strain in one load cycle

3.2疲勞裂紋擴展時CT試件變形幅場的變化規律

3.2.1CT試件裂紋尖端位移幅場的變化規律

圖12 不同裂紋長度時的位移幅場Fig.12 The displacement amplitude field with different crack length

在高頻疲勞裂紋擴展試驗中,試件施加的載荷是正弦變化的,為了方便研究CT試件在高頻諧振載荷作用下疲勞裂紋擴展過程中位移場的變化規律,引入“位移幅場”的概念,將同一應力循環下載荷最大處的位移場與載荷最小處的位移場進行相減得到裂紋尖端位移場幅的分布規律.圖12為疲勞裂紋擴展到不同長度時裂紋尖端區域位移幅場,從圖12(a)中可以看出:當裂紋長度為5mm時,裂紋尖端的位移幅場比較平緩.而從其余3組圖片中可以看出:隨著裂紋長度的增加,裂紋尖端的位移幅場越來越陡峭.這說明在疲勞裂紋擴展過程中,隨著裂紋長度的增加,裂紋尖端位移場的變化也會越來越大。

3.2.2CT試件裂紋尖端應變場幅的變化規律

同樣,為了方便研究CT試件在高頻諧振載荷作用下疲勞裂紋擴展過程中應變場的變化規律,引入“應變幅場”的概念,將同一應力循環下載荷最大處的應變場與載荷最小處的應變場進行相減得到裂紋尖端應變幅場的分布規律.圖13為疲勞裂紋擴展到不同長度時裂紋尖端區域應變幅場,從圖13(a)中可以看出:當裂紋長度為5mm時,裂紋尖端的應變幅場比較平緩.而從其余3組圖片中可以看出:隨著裂紋長度的增加,裂紋尖端的應變幅場越來越尖銳,尤其是裂紋尖端處的應變.這說明在疲勞裂紋擴展過程中,隨著裂紋長度的增加,裂紋尖端應變場的變化也會越來越大。

圖13 疲勞裂紋擴展到不同階段裂紋尖端應變幅場Fig.13 The strain amplitude field with different crack length

4疲勞裂紋尖端處應變值的實驗驗證

為驗證上述疲勞裂紋在無擴展時裂紋尖端處應變的動態有限元計算結果,在帶有5mm疲勞裂紋的CT試件表面貼上電阻應變片并進行疲勞裂紋擴展實驗,實驗是在自行研制的電磁諧振式高頻疲勞試驗平臺上進行的,實驗裝置實物如圖14所示,包括試驗載荷加載系統、裂紋尺寸在線測量系統和固有頻率跟蹤系統和載荷控制系統,試驗載荷加載系統由高頻疲勞試驗機、CT試件、電磁激振器及其放大電路組成,主要完成將電磁激振器所產生的正弦激振力作用在試驗臺上,使試驗臺產生相同頻率的正弦振動,從而使正弦試驗載荷作用于CT試件上的功能.裂紋尺寸在線測量系統:包括光源、CCD圖像傳感器、光學鏡頭、圖像采集卡及計算機,此系統主要用于在線測量疲勞裂紋擴展過程中的裂紋尺寸.載荷控制系統和固有頻率跟蹤系統,用于控制試驗載荷并跟蹤裂紋擴展過程中的固有頻率.其中疲勞試驗機采用如圖1所示的紅山PLG-100諧振式高頻疲勞試驗機,鏡頭為SONY35mm定焦鏡頭,圖像采集卡為美國NI公司所生產的PCI-1014圖像采集卡,CCD為XC-XT50CE黑白CCD攝像頭,分辨率為724X568,試件材料為退火處理后的45#鋼。

圖14 實驗平臺實物圖Fig.14 The image of the experimental facility

首先將帶有預制裂紋的CT試件安裝在疲勞試驗機上,此時系統的諧振頻率為126.4Hz,施加試驗載荷:Fmax=11.65kN,Fmin=6.35kN,Fm=9kN,頻率為126.4Hz,進行疲勞試驗,由裂紋尺寸在線測量系統實時測量擴展疲勞裂紋長度,固有頻率跟蹤系統和載荷控制系統實時跟蹤固有頻率和控制實驗載荷,當裂紋擴展至5mm時停機,取下試件,在試件表面粘帖電阻應變片,由于裂紋尖端存在三維效應及應變梯度,所以電阻應變片不能過于接近裂紋尖端,但為了保證裂紋尖端應變值的測量精度,電阻應變片也不能過于遠離裂紋尖端[12].采用的貼片方式如圖15所示,根據現有經驗,θ=54.27°,r=12.5mm,φ=68.01°.按儀器要求組橋連線,為了測量CT試件裂紋尖端處隨時間變化的動態應變,實驗采用XL2102A型動態電阻應變儀,其工作頻率在DC~100kHz,電磁諧振式高頻疲勞試驗系統的工作頻率為50~300Hz,可見,此儀器完全滿足實驗條件.將帶有5mm疲勞裂紋并貼有電阻應變片的試件安裝在諧振式疲勞試驗機上,此時系統諧振頻率為125Hz,為防止裂紋快速擴展將試驗載荷降至Fmax=8.35kN,Fmin=5.65kN,Fm=7kN,系統在此參數下起振進行試驗,采用動態電阻應變儀測得一個載荷周期內,裂紋尖端處的應變值,并與有限元法得出的應變值做比較,結果如表1所示。

圖15 應變片貼片圖Fig.15 Arrangement of one strain gauge

時間/ms012345678有限元法/μm1.4401.6431.7121.6431.4401.2461.1671.2461.440應變片法/μm1.4071.5991.6741.6071.4061.2151.1351.2211.399誤差/%2.352.752.272.242.412.552.822.052.93

從表1可知:運用動態有限元法在求解CT試件裂紋尖端處應變時,具有較高的精度.與用貼應變的實驗方法做比較,最大的誤差是2.93%.有限元計算結果高于實際測量結果,這是由于裂紋尖端附近應變梯度非常大,而應變片法測量的是測試范圍內的平均應變。

5結論

采用ANSYS靜力學分析與諧響應分析結合的方法,研究了CT試件在高頻恒幅正弦交變載荷下,在一個應力循環周期內及裂紋擴展到不同長度時裂紋尖端區域的位移及應變場的變化規律;通過在CT試件上貼應變片的實驗方法驗證了動態有限元法在求解裂紋尖端處應變值的準確性.通過比較分析結果表明:對于同一疲勞裂紋,I型裂紋尖端位移及應變場均為和載荷同一形式的交變量;隨著裂紋的擴展,I型裂紋尖端的位移及應變場不斷增大;利用動態有限元法計算CT試件裂紋尖端處的應變值,有較高的精度。

參考文獻:

[1]熊纓,陳冰冰,鄭三龍,等.16MnR鋼在不同條件下的疲勞裂紋擴展規律[J].金屬學報,2009,45(7):849-855。

[2]蘇彬,邢海軍,趙穎娣,等.循環載荷應力比對微動疲勞特性的影響[J].浙江工業大學學報,2011,39(4):445-451。

[3]邢海軍,趙穎娣,俞新宇,等.45號鋼微動疲勞門坎值特性研究[J].浙江工業大學學報,2011,39(6):661-665。

[4]張妮.高頻疲勞試驗機動態特性的研究[D].杭州:浙江工業大學,2009。

[5]高紅俐,張立彬,姜偉,等.電磁諧振式疲勞裂紋擴展試驗固有頻率跟蹤系統[J].兵工學報,2013,34(7):896-903。

[6]WALCKER A, WEYGAND D, KRAFT O. Inertial effects on dislocation damping under cyclic loading with ultra-high frequencies[J].Materials Science and Engineering,2010,31(8):4023-4028。

[7]AMIRKHIZI A V, TEHRANIAN A, NASSER S N. Stress-wave energy management through material anisotropy[J].Wave Motion,2010,47(8):519-536。

[8]許楊中,金杰.數值模擬在棒材熱軋過程及缺陷預測中的應用[J].浙江工業大學學報,2011,39(3):308-311。

[9]SAKAUE K, UCHIYAMA Y, TANAKA H, et al. Evaluation of crack tip stress and strain fields under nonproportional loading in a viscoelastic material[J].Engineering Fracture Mechanics,2008,75(14):4140-4150。

[10]楊楠,佟景偉,舒慶璉,等.裂紋紋尖彈塑性應力應變場的實驗研究和數值分析[J].實驗力學,2003,18(3):338-341。

[11]JIANG F C, ASHISH R, KENNETH S V, et al. Crack length calculation for bend specimens under static and dynamic loading[J].Engineering Fracture Mechanics,2004,71(13/14):1971-1985。

[12]CHAKRABORTY D, MURTHY K S R K, CHAKRABORTY D. A new single strain gage technique for the accurate determination of mode I stress intensity factor in orthotropic composite materials[J].Engineering Fracture Mechanics,2014,124:142-154。

(責任編輯:陳石平)

中圖分類號:TP394.1;TH691.9

文獻標志碼:A

文章編號:1006-4303(2015)02-0190-07

作者簡介:高紅俐(1968—),女,河北承德人,副教授,博士研究生,主要從事系統動態特性分析與控制方面的研究,E-mail:ghlzjut@126.com。

收稿日期:2014-12-10

主站蜘蛛池模板: 色欲国产一区二区日韩欧美| 日韩免费毛片| a级免费视频| 亚洲精品人成网线在线| 国产h视频免费观看| 国产亚洲精品自在久久不卡| 国产成人毛片| 欧美福利在线| 精品国产美女福到在线直播| 欧美亚洲网| 五月婷婷伊人网| 亚洲视频二| 亚洲国产日韩在线成人蜜芽| 久久综合干| 亚洲国产欧美目韩成人综合| 久操线在视频在线观看| 99国产在线视频| 国产理论精品| 在线观看国产黄色| 视频一区视频二区中文精品| 干中文字幕| 久久一色本道亚洲| 九九视频在线免费观看| 久久99精品久久久久久不卡| 精品天海翼一区二区| 国产成人精品免费视频大全五级| 国产成人精品一区二区不卡| 在线精品亚洲一区二区古装| 啪啪啪亚洲无码| 亚洲欧洲日韩久久狠狠爱 | 无套av在线| 日韩精品一区二区三区大桥未久| 亚洲熟妇AV日韩熟妇在线| 青青草原偷拍视频| 91精品国产91久无码网站| 九九热这里只有国产精品| a级免费视频| 久久黄色免费电影| 精品国产香蕉伊思人在线| 亚欧成人无码AV在线播放| 美女内射视频WWW网站午夜| 日韩 欧美 国产 精品 综合| 久久亚洲精少妇毛片午夜无码| 色综合激情网| 国产日韩欧美在线视频免费观看| 亚洲精品无码不卡在线播放| 刘亦菲一区二区在线观看| 久久频这里精品99香蕉久网址| 欧美成人一级| 婷婷六月天激情| 日韩av电影一区二区三区四区| 中文字幕在线不卡视频| 全部免费毛片免费播放| 久久青草热| 国产一区二区丝袜高跟鞋| 91在线日韩在线播放| 伊人蕉久影院| 欧美激情网址| 久久精品电影| 久久综合九九亚洲一区| 日韩a级片视频| 欧美一区福利| 香蕉久久国产超碰青草| www欧美在线观看| 亚洲天堂区| 99久视频| 国产激情第一页| 毛片大全免费观看| 四虎影视国产精品| 啪啪永久免费av| 91小视频在线播放| 国产丝袜91| 国产在线观看91精品亚瑟| 日韩午夜伦| 精品福利国产| 午夜福利无码一区二区| 波多野结衣一二三| 久久久久亚洲精品无码网站| 久久毛片网| 亚洲国产成人麻豆精品| 国产福利小视频在线播放观看| 美女内射视频WWW网站午夜|