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

兩類沖擊桿端部的二維數值分析及實驗研究

2016-06-27 09:59:50陳東華王立權
哈爾濱工程大學學報 2016年5期

陳東華, 王立權

(哈爾濱工程大學 機電工程學院,黑龍江 哈爾濱 150001)

兩類沖擊桿端部的二維數值分析及實驗研究

陳東華, 王立權

(哈爾濱工程大學 機電工程學院,黑龍江 哈爾濱 150001)

摘要:沖擊類機械的工作載荷通常由兩實心圓柱形桿件撞擊產生,而對于沖擊桿件的設計目前主要采用忽略徑向效應的一維彈性桿波動力學方法。為了研究徑向效應的影響,采用二維Lax-Wendroff有限差分算法分別對實心圓桿和空心圓桿撞擊端附近區域的軸向應力進行了數值計算。結果表明,撞擊過程中距離實心圓桿撞擊端約1倍直徑的區域內存在應力集中現象,而對于空心圓桿撞擊端鄰近區域應力分布則比較均勻。最后通過實驗測量了撞擊過程中沖擊桿內部關鍵點的應力值,并與數值計算結果進行了對比,二者具有較好的一致性。

關鍵詞:沖擊桿;徑向慣性效應;實心桿;空心桿;有限差分法;應變測量;數值分析

網絡出版地址:http://www.cnki.net/kcms/detail/23.1390.U.20160509.1036.002.html

沖擊機械中經常采用兩圓柱形實心桿件的撞擊來產生沖擊力。目前沖擊桿件主要采用一維彈性桿的波動力學理論進行設計。在一維彈性桿的波動力學理論中,第一個基本假設就是忽略桿中質點徑向慣性作用,沿截面只有均勻分布的軸向應力[1]。文獻[2-6]的實驗研究表明,距離撞擊端2倍桿直徑區域內的軸向應力分布是不均勻的,故一維彈性桿的波動力學理論不再適用。文獻[7]通過采用考慮徑向慣性效應的數值仿真發現在實心桿的撞擊端鄰近區域的軸線上將發生應力集中,提出軸線上的應力集中點是彈性解析解的奇異點,其理論應力為無限大。文獻[8]在對釬桿斷裂的研究中發現,斷裂的位置都偏向于撞擊端,其中有些斷裂就是由內部的裂紋源引起的。由文獻[9-10]可知,沖擊機械在工作過程中,沖擊桿軸線上的應力集中將大大降低沖擊桿的強度,促進疲勞裂紋的產生。所以在進行沖擊桿的設計時有必要考慮徑向慣性作用在撞擊端附近引起的應力集中。

很多文獻對撞擊桿端部附近區域中質點的徑向慣性作用進行了研究。文獻[3]對實心桿表面和桿軸心上的軸向應變進行了實驗研究。除了在表面粘貼應變片,還在距離撞擊端不同距離的位置上嵌入應變片,實驗發現桿軸心上的應變始終比桿表面的應變大。但在試驗中采用的是圓頭沖擊桿,故不能充分說明應變的不同是由徑向慣性效應引起的。文獻[5]采用文獻[3]的方法,在實心鋁圓桿的端部嵌入應變片來測量圓桿內部的彈性應變。實驗測量了兩個3/4英寸直徑的6061-T6鋁桿撞擊的彈性應變波。實驗中是采用平頭沖擊桿,測量結果與數值解的一致性非常好。文獻[7]發展了一種有限差分算法來對應力波在線彈性固體中傳播的偏微分控制方程組進行數值積分,并用該數值算法研究了端部承受拉伸沖擊的半無限圓桿,發現了von Schmidt波在被沖擊圓桿端部軸線上的聚焦現象。文獻[9-11]通過優化沖擊桿的形狀來調整分離式霍普金森壓桿入射脈沖的形狀。文獻[12]對關于兩桿撞擊過程中,在線彈性響應范圍內,鄰近桿撞擊端部區域的應力分布情況的理論和實驗研究進行了綜述。文獻[13-14]對沖擊端部具有不同形狀及不同撞擊面積的沖擊桿與同一實心圓桿的撞擊進行了數值研究和實驗研究[13-14],其數值結果和實驗結果證明沖擊桿的端部形狀將決定被撞圓桿在沖擊端的敏感區域的長度[13]。根據理論研究的解析解敏感區域的長度大概不會超過2.2倍的直徑[14]。

為了獲得考慮徑向慣性效應時,兩實心圓桿和兩空心圓桿在撞擊過程中撞擊端部附近的軸向應力分布情況。本文在文獻[9]中的數值算法基礎上,通過修改初始條件和邊界條件對兩實心圓桿及兩空心圓桿的撞擊進行了數值分析,并根據文獻[3,5]的試驗原理,測量了實心桿軸心和空心桿內表面的軸向應力。

1沖擊動力學方程和數值方法

1.1沖擊動力學方程

對于一個等截面、質量分布均勻、各向同性材料的半無限圓桿,考慮其在不發生扭轉變形,只發生彈性軸對稱變形的情況下,以圓桿的軸作為z軸,將r軸置于圓桿的端面上建立一個拉格朗日柱坐標系。此柱坐標系下的半無限圓桿彈性動力學方程為

(1)

(2)

式中:A1、A2、A3為系數矩陣。

1.2數值方法

在進行苗木育苗的過程當中,相關人員要做好前期的準備工作,技術人員在實際育苗當中,要對育苗施工現場進行實地勘察。針對施工現場的地形、土壤理化性質、肥力情況、土質等作出全面了解,充分保證該施工地區的土壤條件適合育苗。在實際的調查研究過程,技術人員要對現場雜草進行處理,并對周邊區域的場地進行適當的平整。而在實際的栽植當中,技術人員要加強對土壤的檢測,進一步分析土壤的養分、水分等,并結合實際的分析結果,制定合理的苗木種植規劃,為苗木的生長提供良好的環境。

圖1 拉格朗日計算網格Fig.1 Lagrangian computational grid

采用Lax-Wendroff有限差分格式對控制方程組(1)中的空間導數進行離散,可以得到一個兩步差分格式:

(3)

(4)

式中:μ=△t/△z=△t/△r。首先在單元上施加初始條件,然后通過式(3)及施加在邊界節點上的邊界條件計算t+1/2時刻節點的狀態,然后通過式(4)更新t+1時刻的單元狀態。

兩桿撞擊的初始條件為

2數值分析

采用6061-T6鋁作為兩實心桿和兩空心桿彈性撞擊模型的材料,材料密度為2 780 kg·m-3、縱波波速5 089 m·s-1、橫波波速3 151 m·s-1。撞擊桿的撞擊速度均為1 m·s-1。實心桿和空心桿的初始條件均為

1≤m≤M,1≤n≤N1

1≤m≤M,N1

兩實心桿的直徑均為2 cm,為了與實心桿的橫截面積相同,兩空心桿的尺寸都是外徑2.5 cm內徑1.5 cm的圓管。不考慮兩桿撞擊界面上的摩擦,設置摩擦系數k=0;設桿外壁上受到的壓應力Po=0,桿外壁上受到的切應力To=0;設空心桿內壁上受到的壓應力Pi=0,空心桿內壁上受到的切應力Ti=0。

從仿真結果中提取兩實心桿撞擊過程中,桿中出現最大應力時桿中應力的分布情況,計算結果如圖2所示。

圖2中虛線表示兩桿的碰撞端面。可以看出圖中撞擊桿和被撞擊桿的軸向應力分布基本相同。撞擊桿和被撞桿的軸向應力值都在距離碰撞端面2cm左右的軸心區域內急劇升高。這是由于在碰撞的瞬間,在整個碰撞端面上形成平面壓縮波的同時產生了一個平面膨脹波,而隨著膨脹波的傳播以及膨脹波在桿表面的卸載,將在桿的表面形成vonSchmidt波,并向桿內傳播。當vonSchmitdt波傳播到桿軸時,由于軸線上的質點不能發生徑向位移,所以vonSchmidt波將在軸線上聚焦產生應力集中。

圖2 實心桿軸向應力分布Fig.2 Axial stress distribution in the solid bar

為進一步研究vonSchmidt波的聚焦現象,在保持桿的直徑不變的情況下,逐漸減小計算網格的邊長增大網格的密度進行數值計算。從不同網格邊長的計算結果中提取兩桿撞擊過程中,桿中出現最大軸向應力的時刻、位置及軸向應力值,如表1。

表1不同網格邊長實心桿中最大軸向應力出現的時刻、位置及軸向應力值

Table 1The axial stress, the relative positions and relative time of max axial stress with different mesh size of solid bar

網格邊長/mm時刻/μs位置/mm應力/MPa13.53(0,23)10.790.53.14(0,23.5)12.070.3332.94(0,23.7)12.920.252.89(0,24.3)13.630.22.82(0,24.8)14.250.1672.78(0,21.2)14.81

表1給出的位置是相對于桿的軸線及兩桿的撞擊端面的位置,括號中第一個值是與桿軸線的距離,第二個值是與撞擊端面的距離。表1中給出的時刻是以壓縮波的波前到達該位置的時刻為零時刻。可以看出,軸向應力最大值出現的相對位置及相對時間幾乎不隨計算網格變化而變化。根據一維初等理論兩撞擊桿的軸向應力為7.07 MPa,而從表1中可知在考慮徑向慣性效應的情況下,當網格尺寸為1 mm僅為桿半徑的十分之一時,桿中應力集中處的軸向應力最大值是一維初等理論軸向應力值的1.53倍。隨著計算網格尺寸的縮小,軸線上出現的軸向應力增大。上述現象說明實心桿軸線上的軸向應力的數值仿真結果仍未達到其極限。

從以上兩實心桿彈性撞擊的仿真結果可以看出,盡管撞擊速度僅為1 m·s-1,但在兩彈性桿軸線附近的應力集中區域也是十分明顯的。據此可知,若是兩彈塑性材料的實心桿以更大的速度撞擊時,將在距離鄰近撞擊端的軸線附近區域產生嚴重的應力集中,而該區域在撞擊過程中也將發生塑性變形,并且隨著撞擊次數的增加極易發生破壞。

2.2兩空心桿撞擊的數值分析

橫截面積相同的實心桿和空心桿的一維初等理論桿模型是完全相同的,而在二維桿模型中實心桿和空心桿在處理內邊界條件上有明顯的不同,在本文1.2節中已經給出。結合2.1節及下面兩空心桿撞擊的二維數值仿真結果可以明顯的看出由內邊界條件不同帶來的影響。同樣從仿真結果中提取出兩空心桿撞擊過程中,桿中出現最大軸向應力時,空心桿中應力的分布情況。計算結果如圖3所示。

圖3 空心桿軸向應力分布Fig.3 Axial stress distribution in the hollow bar

圖3中虛線仍表示兩空心桿的撞擊面。從圖中可以看出,在兩實心桿撞擊過程中出現的應力集中現象,在兩空心桿撞擊過程中沒再出現。在空心桿的縱截面軸向應力云圖中,壓縮波的波陣面呈明顯的帶狀分布,而且在空心桿從撞擊端到應力最大值所在的橫截面之間的區域應力分布相對于實心桿中此區域的應力分布均勻很多。

當減小計算網格的邊長時,從不同網格邊長的計算結果中提取兩桿撞擊過程中,桿中出現最大軸向應力的時刻、位置及軸向應力值,如表2所示。

表2中給出的位置是相對于桿的軸線及兩桿的撞擊端面的位置,括號中第一個值是與桿軸線的距離,第二個值是與撞擊端面的距離。表2中給出的時刻是以壓縮波的波前到達該位置的時刻為零時刻。可以看出同兩實心桿撞擊一樣,兩空心桿撞擊過程中軸向應力最大值出現的相對位置及相對時刻幾乎不隨計算網格邊長的變化而變化,相對位置穩定在距離桿撞擊端12倍空心桿壁厚左右的外壁上。橫截面積相同的實心桿和空心桿的一維初等理論模型是相同的,所以這兩個空心桿撞擊的一維軸向應力解仍為7.07 MPa。從表2可知考慮徑向效應的兩空心桿在撞擊過程中的最大軸向應力大概為一維軸向應力解的1.5倍。

根據以上兩空心桿撞擊的數值計算結果可以發現,當兩空心桿撞擊時,桿中不會出現應力集中區域。撞擊端至壓縮波波前的區域中的應力分布比較均勻。撞擊過程中的軸向應力最大值是有界的。

表2不同網格邊長空心桿中最大軸向應力出現的時刻、位置及軸向應力值

Table 2The relative positions and relative time of max stresswith different mesh size of hollow bar

網格邊長/mm時刻/μs位置/mm應力/MPa0.51.76(12.5,57)8.910.251.62(12.5,59.8)9.800.1671.54(12.5,60.5)10.230.1251.49(12.5,61.5)10.500.11.47(12.5,62.1)10.670.0831.46(12.5,62.7)10.78

3實驗驗證

為了驗證數值分析的結果,根據文獻[5,7]的實驗原理。考慮到實心桿和空心桿都是軸對稱形狀,通過在半圓桿剖截面軸線上布置應變片來近似測量撞擊過程中實心桿軸線處的軸向應力。通過在半圓管的內壁中線上布置應變片來近似測量空心桿內壁上的軸向應力。并在與之相對的表面位置上粘貼應變片,用來測得半圓桿及半圓管外表面上的軸向應變。用相同尺寸的圓管與圓桿,在距離碰撞端相同距離的外表面上粘貼應變片測量軸向應變作為對照試驗。采用同樣尺寸的圓桿與圓管作為沖擊桿。為了降低撞擊端面上的摩擦,在撞擊桿與被撞擊桿的端面上都涂有潤滑油。

為了方便實驗,僅在被撞擊桿上粘貼應變片。實驗采用的被撞擊桿如圖4所示。撞擊桿與被撞擊桿都是采用6061-T6鋁的鋁桿和鋁管,幾何尺寸都與數值仿真中采用的幾何尺寸一樣。根據數值仿真的結果,在直徑為20 mm的實心半圓桿的軸線上距離撞擊端20 mm處及對應的半圓桿外表面上都貼有應變片。作為對比的同尺寸的圓桿也在距離碰撞端20 mm處的外表面上貼有應變片。外徑為25 mm、內徑為15 mm的半圓管在距離撞擊端60 mm處的內表面中線上及對應的外表面上都貼有應變片。作為對比的同尺寸的圓管也在距離碰撞端60 mm處的外表面上貼有應變片。實驗中使用的應變片型號為BX120-1AA,敏感柵長1 mm。

圖4 被撞擊桿Fig.4 Impacted rod

圖5 距離撞擊端20 mm處實心桿軸線上的軸向應力時間歷程的實驗和仿真結果Fig.5 Experimental and numerical result at distance of  20 mm from impact on the solid rod center line

圖6 距離撞擊端20 mm處實心桿外表面上的軸向應力時間歷程的實驗和仿真結果Fig.6 Experimental and numerical result at distance of  20 mm from impact on the solid rod surface

下面各應力歷程圖中的時間都是相對時間,零時刻是指壓縮波的波前到達測量位置時的時刻。

圖5和圖6中的數值仿真結果都是網格邊長為0.167 mm的仿真結果。從圖5和圖6中可以看出,數值計算得到的應力歷程與實驗測得的應力歷程是相似的。由圖5可知,與仿真結果一樣,實驗過程中在桿的軸線上確實存在應力升高的現象。由圖6可以看出,在距離碰撞端相同位置處,半圓桿外表面與圓桿外表面上的應力歷程基本相符,說明半圓桿實驗并不會對桿表面上的應力產生顯著影響,因此可知在半圓桿軸線上測得的應力歷程與圓桿軸線上的應力歷程近似相符。

對比圖5和圖6可以看出,當壓縮波的波前到達測量位置后的2~6 μs這段時間內,各時刻桿軸線與桿表面上的軸向應力相差較大。可以說明這段時間內桿中的軸向應力分布極不均勻,結合圖5中測量到的應力峰值,據此可知這段時間內圓桿此處的軸線附近區域將發生應力集中。由數值結果及試驗結果可知,一維初等理論不適用于兩實心沖擊桿的撞擊端附近區域。

圖7 距離撞擊端60 mm處管內壁中線上的軸向應力 時間歷程的實驗和仿真結果Fig.7 Experimental and numerical result at distance of  60 mm from impact on the tube’s internal surface  center line

圖8 距離撞擊端60 mm處管外壁中線上的軸向應力 時間歷程的實驗和仿真結果Fig.8 Experimental and numerical result at distance of  60 mm from impact on the tube’s external surface center line

圖7和圖8中的數值仿真結果都是網格邊長為0.083 mm的仿真結果。從圖7和圖8中可以看出,數值計算得到的應力歷程與實驗測得的應力歷程基本相符。由圖8可以看出,在距離碰撞端相同位置處,半圓管外表面與圓管外表面上的軸向應力歷程非常相似,說明采用半圓管進行實驗并不會明顯影響管外表面上的軸向應力,因此可知在半圓管內表面中線上測得的應力歷程與圓管內表面上的應力歷程近似相符。對比圖7和圖8可以看出,半圓管內表面和半圓管外表面上應力變化基本一致,幾乎同時達到應力最大值,據此可知碰撞過程中圓管壁中的應力分布是比較均勻的,不會發生應力集中現象。

4結論

利用二維Lax-Wendroff有限差分格式對兩半無限圓桿和兩半無限圓管的碰撞進行了數值模擬,并采用半圓桿和半圓管進行了碰撞實驗并近似測量了內部的軸向應變。經對比仿真結果與實驗數據有較好的一致性。

1)數值仿真表明,對于兩實心桿撞擊,將在靠近撞擊端1倍直徑左右的區域發生明顯的應力集中現象,當計算網格為0.167 mm時,應力集中區域的最大軸向應力為14.81 MPa,是一維初等理論軸向應力的2.1倍,并且數值結果還有繼續升高的趨勢。實驗測得的半圓桿軸線上的應力歷程及半圓桿外表面上的應力歷程同樣表明此區域確實發生了應力集中。

2)從對兩圓管撞擊的數值仿真及半圓管的實驗可知,圓管內表面和外表面上的軸向應力歷程幾乎相同。管壁內的軸向應力分布比較均勻。

3)半圓桿軸線上測得的軸向應力歷程可以近似圓桿軸線上的軸向應力歷程。半圓管內表面中線上測得的軸向應力歷程可以近似圓管內壁上的軸向應力歷程。

綜上,等橫截面積的實心桿和空心桿在彈性撞擊過程中,空心桿的軸向應力分布均勻,受力情況要優于實心桿。故在沖擊機械中采用空心桿作為沖擊桿件將能改善沖擊桿的受力情況,提高設備的使用壽命。

參考文獻:

[1]王禮立. 應力波基礎[M]. 2版. 北京: 國防工業出版社, 2005: 7-9.

WANG Lili. Foundation of stress waves[M]. 2nd ed. Beijing: National Defence Industry Press, 2005: 7-9.

[2]DAVIES R M. Stress waves in solids[J]. British journal of applied physics, 1956, 7(6): 203-209.

[3]BAKER W E, DOVE R C. Measurement of internal strains in a bar subjected to longitudinal impact[J]. Experimental mechanics, 1962, 2(10): 307-311.

[4]BERTHOLF L D, KARNES C H. Axisymmetric elastic-plastic wave propagation in 6061-T6 Aluminum bars of finite length[J]. Journal of applied mechanics, 1969, 36(3): 533-541.

[5]HABBERSTAD J L, HOGE K G, FOSTER J E. An experimental and numerical study of elastic strain waves on the center line of a 6061-T6 Aluminum bar[J]. Journal of applied mechanics, 1972, 39(2): 367-371.

[6]HE L, MA G W, KARP B, et al. Investigation of dynamic saint-venant’s principle in a cylindrical waveguide-experimental and numerical results[J]. Experimental mechanics, 2015, 55(3): 623-634.

[7]LIN Xiao, BALLMANN J. A numerical scheme for axisymmetric elastic waves in solids[J]. Wave motion, 1995, 21(2): 115-126.

[8]黃呂權, 夏琴香, 唐耀輝. 錘桿的疲勞強度研究[J]. 機械強度, 1992, 14(1): 33-39, 63.

HUANG Lvquan, XIA Qinxiang, TANG Yaohui. Study on the fatigue strength of hammer rod[J]. Journal of mechanical strength, 1992, 14(1): 33-39, 63.

[9]李夕兵, 古德生, 賴海輝. 沖擊載荷下巖石動態應力-應變全圖測試中的合理加載波形[J]. 爆炸與沖擊, 1993, 13(2): 125-130.

LI Xibing, GU Desheng, LAI Haihui. On the reasonable loading stress waveforms determined by dynamic stress-strain curves of rocks by SHPB[J]. Explosion and shock waves, 1993, 13(2): 125-130.

[10]李夕兵, 劉德順, 古德生. 消除巖石動態實驗曲線振蕩的有效途徑[J]. 中南工業大學學報, 1995, 26(4): 457-460.

LI Xibing, LIU Deshun, GU Desheng. Effective method of elimilating the oscillation of rock dynamic stress-strain-strain rate curves[J]. Journal of central south university, 1995, 26(4): 457-460.

[11]鄒廣平, 王新征, 陳劍杰, 等. 空心錐撞擊桿對分離式霍普金森壓桿入射脈沖彌散的抑制效果分析[J]. 兵工學報, 2012, 33(11): 1346-1351.

ZOU Guangping, WANG Xinzheng, CHEN Jianjie, et al. Restrain effect of hollow cone striker bar on incident pulse dispersion of split Hopkinson pressure bar[J]. Acta armamentarii, 2012, 33(11): 1346-1351.

[12]KARP B, DURBAN D. Saint-Venant’s principle in dynamics of structures[J]. Applied mechanics reviews, 2011, 64(2): 020801.

[13]KARP B, DOROGOY A, WANG Zonggang. Non-uniform impact excitation of a cylindrical bar[J]. Journal of sound and vibration, 2009, 323(3/4/5): 757-771.

[14]HE L, MA G W, KARP B, et al. Investigation of dynamic Saint-Venant's principle in a cylindrical waveguide-Analytical results[J]. International journal of impact engineering, 2014, 73: 135-144.

本文引用格式:

陳東華, 王立權. 兩類沖擊桿端部的二維數值分析及實驗研究[J]. 哈爾濱工程大學學報, 2016, 37(5): 707-712.

CHEN Donghua,WANG Liquan. Two-dimensional numerical analysis and experimental study of two types of impact bars[J]. Journal of Harbin Engineering University, 2016, 37(5): 707-712.

Two-dimensional numerical analysis and experimental study of two types of impact bars

CHEN Donghua,WANG Liquan

(College of Mechanical and Electrical Engineering, Harbin Engineering University, Harbin 150001, China)

Abstract:Impact machines often generate impact loads by two solid cylindrical rods that are known as impact bars. The current design of striker bars mainly adopts one-dimensional elementary theory that ignores the radial-inertia effects. To research the radial-inertia effects, the two-dimensional Lax-Wendroff finite difference method was applied to calculate axial stress adjacent to the impact ends of a solid round bar and a hollow round bar, respectively. The numerical simulation shows that the stress concentration occurs in the radial area at one diameter away from the impact end of a solid round bar, whereas stress uniformly distributes over a longitudinal section near the end of a hollow round bar. Finally, internal strain experimental data of key points was acquired through strain gages placed on the center lines of the semi-circle rod and the semi-circle tube used in experiments. Numerical results showed good agreement with the experimental data.

Keywords:impact bar; radial-inertia effects; solid round bar; hollow round bar; finite difference method; strain measurement ;numerical analysis

收稿日期:2015-04-24.

基金項目:國家自然科學基金資助項目(51279042).

作者簡介:陳東華(1986-), 男, 博士研究生; 通信作者:王立權,E-mail:wangliquan@hrbeu.edu.cn.

DOI:10.11990/jheu.201504055

中圖分類號:TB125

文獻標志碼:A

文章編號:1006-7043(2016)05-0707-06

網絡出版時間:2016-05-09.

王立權(1957-), 男, 博士生導師.

主站蜘蛛池模板: 91成人在线免费观看| 日韩av高清无码一区二区三区| 97国产精品视频人人做人人爱| 日韩黄色大片免费看| 日韩精品一区二区三区免费在线观看| 亚洲日韩精品无码专区| 55夜色66夜色国产精品视频| 久久亚洲精少妇毛片午夜无码| 午夜性刺激在线观看免费| 女人av社区男人的天堂| 真实国产乱子伦视频| 精品91在线| 欧美激情视频一区| 精品一区二区三区自慰喷水| 国产亚洲精| 久久这里只有精品免费| 亚洲一区二区黄色| 国产女人在线视频| 国产成人精品高清不卡在线| 亚洲无线观看| 欧美日韩一区二区三| 9966国产精品视频| 欧美日韩国产成人在线观看| 女同久久精品国产99国| 日韩中文欧美| 国产精品30p| 亚洲永久色| 最新日韩AV网址在线观看| 亚洲嫩模喷白浆| 国产又黄又硬又粗| 国产亚洲精久久久久久无码AV| 精品国产aⅴ一区二区三区| 成人精品视频一区二区在线| 久久a级片| 国内精品视频区在线2021| 日韩在线永久免费播放| 亚洲成人免费在线| 天堂在线亚洲| 国产成人精品一区二区三区| 伊人激情综合网| 国产在线视频欧美亚综合| 无码免费的亚洲视频| 她的性爱视频| 香蕉视频国产精品人| 精品乱码久久久久久久| 亚洲综合香蕉| 六月婷婷精品视频在线观看 | 欧美性久久久久| 国模视频一区二区| 在线观看亚洲精品福利片| 久久性视频| 5388国产亚洲欧美在线观看| 欧美特级AAAAAA视频免费观看| 日日拍夜夜操| 国产日韩精品欧美一区灰| 91娇喘视频| 国产精品亚洲专区一区| 国产成人艳妇AA视频在线| 99久久精品免费观看国产| 欧美色图久久| 精品国产亚洲人成在线| 亚洲成在人线av品善网好看| 高清无码一本到东京热| 青青操国产视频| 无码专区在线观看| 久久精品国产91久久综合麻豆自制| 国产69囗曝护士吞精在线视频| 精品国产成人a在线观看| 欧美全免费aaaaaa特黄在线| 国产麻豆精品在线观看| 成年女人a毛片免费视频| 日韩麻豆小视频| 国产精品成人久久| 亚洲中文字幕日产无码2021| 毛片网站免费在线观看| 国产极品美女在线播放| 97视频在线观看免费视频| 99精品这里只有精品高清视频| 国产精品伦视频观看免费| 免费三A级毛片视频| 亚洲国产综合第一精品小说| 欧美成人午夜视频免看|