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

一般宏觀應力狀態下凹角蜂窩結構的屈曲性能分析*

2023-03-10 08:07:16周世奇侯秀慧鄧子辰
應用數學和力學 2023年1期
關鍵詞:模態方向有限元

周世奇,侯秀慧,2,鄧子辰,2

(1.西北工業大學 力學與土木建筑學院,西安 710129;2.復雜系統動力學與控制工業和信息化部重點實驗室,西安 710072)

引言

蜂窩是典型的多孔材料,這種材料以平面內單元的二維排列和面外方向平行堆疊的周期性拓撲分布為特征.單胞相互連接形成的網絡結構使得蜂窩相較于其他材料具有更高的孔隙率和更低的相對密度,這也導致了更高的比剛度強度和比能量吸收.單胞的拓撲結構可以顯著地影響這種超輕材料的力學性能.因此,合理設計的微觀結構使蜂窩具有一些超結構屬性,比如負Poisson 比[1]、負熱膨脹[2]、負剛度[3]等.這些特殊的性質源于蜂窩不同的微觀幾何拓撲,而不是構成蜂窩的基體材料.例如將傳統正六邊形蜂窩結構的凸形胞壁替換為凹形胞壁,可以得到相應的負Poisson 比特性,即結構在張力作用下縱向和橫向雙向膨脹而在壓力作用下縱向和橫向雙向收縮.其他拉脹蜂窩結構還包括雙V 形[4]、手性[5]和星形[6-7]胞元等,其中凹角蜂窩結構作為最簡單也是最經典的拉脹蜂窩,受到了廣泛的關注.以往的研究大多是基于凹角蜂窩基本的物理力學行為,如壓縮[8]、沖擊[9-10]、熱學[11]、光學[12]以及聲學[13]等領域的應用.如Hou 等[14]研究了凹角六邊形蜂窩的能量吸收特性,與傳統蜂窩相比,凹角六邊形蜂窩具有更高的吸能效率.而在低速準靜態壓縮荷載作用下,蜂窩結構的胞壁會發生屈曲破壞從而喪失穩定性.

已有研究表明,蜂窩結構在壓縮載荷作用下發生屈曲的特性,并不總是作為一種結構缺陷的存在.蜂窩結構的微觀不穩定模式作為一種調制周期性結構宏觀力學性能的方法已經被廣泛采用,比如手性、波傳播、聲子特性以及光學特性等.Bertoldi 等[15]揭示了多胞結構由于彈性失穩而產生負Poisson 比效應的機理.Yang等[16]利用多胞結構的屈曲提出了一種屈曲致動器用于制作軟體機器人.由此可知,對蜂窩結構的屈曲分析尤其是凹角蜂窩這種拉脹蜂窩的屈曲分析是十分必要的.Jiménez 和Triantafyllidis[17]研究了矩形和正六邊形蜂窩在軸向壓縮和橫向剪切組合荷載作用下的屈曲,發現屈曲模式高度依賴于載荷類型.Combescure 等[18]利用群論方法研究了有限應變圓形蜂窩的變形模式及其穩定性.梁觀坡等[19]基于數值模擬和理論分析研究了含周期性橢圓孔二維結構的屈曲行為,研究發現改變孔的幾何參數,橢圓孔結構的屈曲模態隨之發生轉換.Peng等[20]將變分漸近法擴展到凹角蜂窩夾層板的面外屈曲研究,得出在相同條件下,凹角蜂窩夾層板的屈曲載荷大于傳統的蜂窩夾層板.Gavazzoni 等[21]基于材料和幾何非線性以及制造缺陷之間的耦合作用,研究了軟多胞結構局部不穩定性的循環響應.然而,卻少見針對負Poisson 比凹角蜂窩結構面內屈曲模態力學性能的研究報道.

本文系統研究了凹角蜂窩的面內屈曲特性,在第1 節中,建立了凹角蜂窩的有限元模型,對凹角蜂窩在單軸壓縮下的屈曲行為進行了數值模擬.以此為基礎,在第2 節中建立了凹角蜂窩對應于不同屈曲模態的屈曲強度解析公式,揭示凹角蜂窩不同屈曲模態的產生機理.隨后,在第3 節中通過實驗,對以上數值模擬和理論分析結果進行了驗證.

1 有限元分析

利用有限元軟件ABAQUS 對凹角蜂窩在宏觀應力下的屈曲進行模擬.建立的有限元模型在水平方向上設置5 個胞元,在豎直方向上設置7 個胞元,將胞壁材料設置為理想的彈塑性模型,質量密度為1230 kg/m3,彈性模量為35 MPa,Poisson 比為0.45.有限元模型是三維殼結構,為了研究凹角蜂窩面內屈曲性能,限制模型的面外位移和旋轉.構造的三維凹角蜂窩模型胞壁的面內長細比均大于1/10,受力情況類似于板殼結構,實體單元在模擬板殼結構時容易形成剪力自鎖從而影響計算結果,所以本文在進行凹角蜂窩面內屈曲的有限元模擬時采用殼單元.模擬沿著平行于水平胞壁和垂直于水平胞壁方向的準靜態壓縮時,在模型底部的每一個節點上設置固定約束限制模型底部的線位移和轉角位移.在模型頂部的每一個節點上設置滑動約束限制模型頂部的水平位移和轉角位移,如圖1所示.為了描述方便,將平行于水平胞壁和垂直于水平胞壁的方向分別定義為x方向和y方向.用X=σx/(E(t/L)3)和Y=σy/(E(t/L)3)分別表示沿著x方向和y方向的應力系數.X=0,Y=1表示沿著y方向的單軸壓縮,X=1,Y=0表示沿著x方向的單軸壓縮.

圖1 凹角蜂窩的有限元模型和邊界條件:(a)沿著x 方向準靜態壓縮的邊界條件;(b)沿著y 方向準靜態壓縮的邊界條件Fig.1 The finite element model and boundary conditions for the re-entrant honeycomb:(a)boundary conditions for the quasi-static compression along the x direction;(b)boundary conditions for the quasi-static compression along the y direction

對凹角蜂窩有限元模型劃分網格時,網格類型設置為S4R,該單元類型是通用的殼單元類型,并且這種單元類型既適用于薄殼結構,也適用于厚殼結構,比較適合用來計算凹角蜂窩這種胞壁長細比不同的結構.為了確定網格的全局尺寸,進行了收斂性分析.收斂性分析選用的有限元模型的幾何設置、材料屬性、邊界條件和加載條件與前文所述的有限元模型保持一致.根據圖2所示的收斂性計算結果,網格尺寸設置為5 mm,4 mm,3 mm 時,計算結果誤差較大;網格尺寸為2 mm 時,曲線開始趨于平緩.為確保有限元模擬結果最終收斂,我們計算了網格尺寸為1 mm,0.5 mm,0.3 mm,0.2 mm,0.1 mm 的情況,發現網格尺寸為0.5 mm,0.3 mm,0.2 mm,0.1 mm 時,計算結果基本一致,考慮到有限元模擬研究的簡潔高效,將網格尺寸設置為0.5 mm.

圖2 收斂性計算結果Fig.2 Convergence calculation results

圖3 列出了凹角蜂窩在單軸壓縮荷載作用下的前十階屈曲模態.當凹角蜂窩承受沿著y方向的壓縮載荷時(X=0,Y=1),前五階屈曲模態的單胞變形是一致的,而整體屈曲模態表現為加載方向上的波峰數逐階遞增,這與單桿壓縮的高階屈曲模態類似[22].從第六階開始出現不同的單胞變形,這是因為所選結構的尺寸達到所能承受的波峰數的極限從而導致多階屈曲模態偶連.當凹角蜂窩承受沿著x方向的壓縮載荷時(X=1,Y=0),前八階屈曲模態的單胞變形保持一致,而整體屈曲模態沒有觀察到波峰數遞增的現象.從第九階開始出現不同的單胞變形,同樣考慮是多階屈曲模態偶連導致的.因此分別提取第一階模態進行理論分析,將沿著x方 向準靜態壓縮得到的屈曲模態命名為屈曲模態Ⅰ,沿著y方向準靜態壓縮得到的屈曲模態命名為屈曲模態Ⅱ,如圖4所示.

圖3 凹角蜂窩在單軸荷載作用下的前十階屈曲模態:(a)沿著y 方向單軸壓縮;(b)沿著x 方向單軸壓縮Fig.3 The 1st 10 buckling modes of the re-entrant honeycomb under uniaxial loading:a)uniaxial compression along the y direction;(b)uniaxial compression along the x direction

圖4 凹角蜂窩的第一階模態:(a)沿著x 方向準靜態壓縮得到屈曲模態Ⅰ;(b)沿著y 方向準靜態壓縮得到屈曲模態ⅡFig.4 The 1st buckling modes of the re-entrant honeycomb:(a)buckling mode Ⅰ obtained under quasi-static compression along the x direction;(b)buckling mode Ⅱ obtained under quasi-static compression along the y direction

為了研究尺寸效應對凹角蜂窩的屈曲模態和屈曲強度的影響,設置了三種不同胞元數量的凹角蜂窩模型:3 × 4 構型、5 × 7 構型、7 × 10 構型,第一個數字代表水平方向的胞元數量,第二個數字代表豎直方向的胞元數量.這三種構型的面內長寬比均等于1,這樣就排除了面內長寬比對凹角蜂窩屈曲行為的影響.結果表明:三種構型的一階屈曲模態是相同的;隨著胞元數量的增加,單胞變形模態轉變的現象會得到延緩,如圖3和圖5所示;凹角蜂窩的胞元數量對屈曲強度的影響可以忽略不計,如圖6所示,因此后續研究采用5 × 7 構型.

圖5 不同單胞數量的凹角蜂窩的屈曲模態:(a)3×4 構型沿著y 方向單軸壓縮;(b)3×4 構型沿著x 方向單軸壓縮;(c)7×10 構型沿著y 方向單軸壓縮;(d)7×10 構型沿著x 方向單軸壓縮Fig.5 Buckling modes of the re-entrant honeycomb with different numbers of unit cells:(a)the 3×4 configuration under uniaxial compression in the y direction;(b)the 3×4 configuration under uniaxial compression in the x direction;(c)the 7×10 configuration under uniaxial compression in the y direction;(d)the 7×10 configuration under uniaxial compression in the x direction

圖6 不同單胞數量的凹角蜂窩的屈曲強度Fig.6 Buckling strengths of the re-entrant honeycomb with different numbers of unit cells

2 理論分析

通過有限元模擬凹角蜂窩在不同加載方式下的屈曲行為,分別獲得凹角蜂窩沿著x方向和y方向的單軸壓縮屈曲變形模式,并從中提取具有代表性的單胞變形模式為理論分析模型(圖4).依據 Timoshenko 和Gere 提出的梁柱理論[23],構建承受軸力作用的桿件桿端彎矩和桿端轉角位移的關系.即梁柱結構在軸向力P作用下,當其外荷載為一對順時針桿端彎矩Ma和Mb時(見圖7,下標a和b分別對應桿件的兩端結點),梁柱的桿端轉角θa和 θb可以使用梁柱方程計算得到:

圖7 梁柱的桿端轉角和桿端彎矩Fig.7 End rotations and end bending moments of the beam-column

利用梁柱理論公式對凹角蜂窩單胞的每一個胞壁的桿端彎矩和轉角位移建立相應的梁柱方程,對于具有桿端線位移的胞壁還需建立彎矩平衡方程.將這些方程用統一的矩陣公式表示,根據屈曲臨界條件,即相應位移不全為零,可得方程組的系數矩陣的行列式為零,解行列式得凹角蜂窩屈曲穩定方程.

2.1 屈曲模態Ⅰ

凹角蜂窩整體與胞元的應力狀態和幾何結構如圖8所示,水平桿的長度為La,斜桿的長度為Lb,水平桿和斜桿的夾角為θ,面外厚度為c.對于任意的應力 σx和 σy,可等效為水平桿和斜桿上的集中力F,P和W,分別為

圖8 一般宏觀應力狀態下的凹角蜂窩結構:(a)凹角蜂窩整體的幾何結構和應力狀態;(b)基礎單胞的幾何結構和應力狀態Fig.8 The re-entrant honeycomb structure under a general macroscopic stress state:(a)the geometric structure and the stress state of the re-entrant honeycomb;(b)the geometric structure and the stress state of the unit cell

由式(2)可得OA,OB和OC桿的軸力分別為Fa,Fb和Fc:

這里軸力的符號規定為壓力為正,拉力為負.從有限元模擬的屈曲模態圖4(a)中提取凹角蜂窩的胞元屈曲模態如圖9所示,虛線表示屈曲前的位置,實線表示屈曲后的位置.由于OA桿的屈曲變形是中心對稱的,所以OA桿兩端的彎矩大小和方向均相等,而OB,OC桿的屈曲變形是關于過桿中點的法線對稱的,所以OB,OC桿兩端的彎矩大小相等方向相反.OA桿兩端有相對線位移,而OB,OC桿沒有.

圖9 凹角蜂窩屈曲模態Ⅰ的單胞模態和OA,OC 桿的內力示意圖Fig.9 The unit cell modes for buckling mode Ⅰ of the re-entrant honeycomb and the internal forces of rods OA and OC

將OA,OB和OC三根桿的梁柱理論關系表達式以及OA桿的彎矩平衡方程表示為式(5)中的矩陣形式:

其中前三行分別表示OA,OB和OC桿的梁-柱理論關系式,第四行表示O點的力矩平衡,第五行表示OA桿關于O點的力矩平衡.式(5)中方程的非齊次項不包含桿端彎矩、桿端轉角和軸力,所包含的剪力只會引起靜態撓度而非屈曲,不影響系數矩陣的結果.

由圖9 可知,當發生屈曲時,α和β不全為0,故式(5)中的方程組存在非零解,系數行列式等于零.

化簡后可得

式 (8)是關于qa,qb和qc的關系式,表示凹角蜂窩在宏觀應力狀態下產生屈曲模態Ⅰ的應力條件.

2.2 屈曲模態Ⅱ

從有限元模擬的屈曲模態圖4(b)中提取凹角蜂窩的胞元屈曲模態如圖10所示,虛線表示屈曲前的位置,實線表示屈曲后的位置.由于每根桿的屈曲變形是中心對稱的,所以OA桿兩端的彎矩大小和方向均相等,OB和OC桿同理.OB,OC桿兩端有相對線位移,而OA桿沒有,這與屈曲模態Ⅰ不同.

將OA,OB和OC三根桿的梁柱理論關系表達式以及OB和OC桿的彎矩平衡方程表示為式(9)中的矩陣形式:

其中前三行分別表示OA,OB和OC桿的梁-柱理論關系式,第四行表示O點的力矩平衡,第五行表示OB桿關于O點的力矩平衡,第六行表示OC桿關于O點的力矩平衡.

由圖10 可知,當發生屈曲時,α和β不全為0,故式(9)中的方程組存在非零解,系數行列式等于零.

圖10 凹角蜂窩屈曲模態Ⅱ的單胞模態和OA,OB 桿的內力示意圖Fig.10 The unit cell modes for buckling mode Ⅱ of the re-entrant honeycomb and the internal forces of rods OA and OC

化簡后可得

式(12)是關于qa,qb和qc的關系式,表示凹角蜂窩在宏觀應力狀態下產生屈曲模態Ⅱ的應力條件.

2.3 理論結果分析

根據式(8)和式(12)的關系繪制曲線,如圖11所示.圖11 給出了凹角蜂窩在宏觀應力狀態下屈曲模態Ⅰ和屈曲模態Ⅱ對應的應力關系,水平軸和豎直軸分別代表無量綱化的宏觀應力,分別用X=σx/(E(t/L)3)和Y=σy/(E(t/L)3)表示.點線代表屈曲模態Ⅰ,點劃線代表屈曲模態Ⅱ.

圖11 凹角蜂窩屈曲模態Ⅰ和Ⅱ的失效界面Fig.11 Failure surfaces for buckling modes Ⅰ and Ⅱ of the re-entrant honeycomb

兩條曲線的交點出現在第一象限,在交點右邊的區域,同樣的Y值,屈曲模態Ⅰ的X值小于屈曲模態Ⅱ的X值.在這塊區域包含的任意應力狀態下,凹角蜂窩的屈曲模態為屈曲模態Ⅰ.在交點左邊的區域,同樣的X值屈曲模態Ⅱ的Y值小于屈曲模態Ⅰ的Y值.在這塊區域包含的任意應力狀態下,凹角蜂窩的屈曲模態為屈曲模態Ⅱ.也即表明,當沿著x方 向的應力為主應力時,凹角蜂窩的屈曲變形呈現為屈曲模態Ⅰ;當沿著y方向的應力為主應力時,凹角蜂窩的屈曲變形呈現為屈曲模態Ⅱ.凹角蜂窩在雙軸應力狀態下的優先屈曲模式對雙軸應力的比例關系表現出敏感性.同時我們還發現,點劃線延伸至第三象限,這表明凹角蜂窩在雙軸拉伸的應力狀態下依然會發生屈曲變形,而傳統正六邊形蜂窩僅在受壓的狀態下才會發生屈曲變形,這可歸結為凹角蜂窩不同于一般蜂窩的凹角幾何特征所引起的拉脹效應.根據圖11 可得凹角蜂窩在x方向單軸壓縮下的屈曲強度為σx/E=0.076(t/L)3(圖11 中三角形標注),在y方向單軸壓縮下的屈曲強度為σy/E=0.158(t/L)3(圖11中菱形標注),分別對應著在單軸壓縮下屈曲模態Ⅰ和屈曲模態Ⅱ的屈曲強度.這個結果表明屈曲模態Ⅰ的屈曲強度低于屈曲模態Ⅱ的屈曲強度,因此屈曲模態Ⅰ可認定為整體屈曲模態,屈曲模態Ⅱ可認定為局部屈曲模態.,

3 實驗

3.1 實驗設置

采用增材制造技術制備實驗樣本,3D 打印機型號為RAISED3D E2,打印材料采用PolyFlex TPU95A.樣本的長寬高為129.8 mm×129.2mm×18.3mm,單胞水平桿長度為18.3mm,斜桿長度為10.2mm,水平桿和斜桿的夾角壁厚為1 mm,打印層厚為0.025mm.實驗樣本如圖12所示.使用萬能試驗機分別對實驗樣本沿著x方向和y方向進行準靜態壓縮實驗.行程位移設置為20 mm,壓縮速度設置為1 mm/min,溫度設置為23 ℃,可認定為常溫條件下的準靜態壓縮.

圖12 凹角蜂窩沿著x方向和y 方向單軸壓縮的實驗結果和數值結果對比Fig.12 Comparison of experimental and numerical results of re-entrant honeycombs under uniaxial compression along the x and y direction

3.2 實驗結果和理論研究、數值研究的對比

圖12 展示了凹角蜂窩結構在準靜態壓縮條件下的變形情況.提取實驗結果中凹角蜂窩完全屈曲變形的模態圖與第1 節中有限元模擬的屈曲模態進行對比,發現凹角蜂窩結構的整體屈曲變形模態和單胞屈曲變形模態在實驗和有限元模擬中保持高度的一致性,如圖12所示.實驗結果表明凹角蜂窩沿著x方向單軸壓縮時,屈曲變形模態表現為屈曲模態Ⅰ,而沿著y方向單軸壓縮時,屈曲變形模態表現為屈曲模態Ⅱ,這與理論預測的結果一致.

從圖13 中可以看到,實驗樣本分別承受沿著x方向和y方向的單軸壓縮,在應變到0.039 時,沿著x方向壓縮的應力-應變曲線由近似線性的彈性階段過渡為平臺階段,即桿件開始發生屈曲變形;而沿著y方向壓縮的應力-應變曲線的彈性階段在應變到0.044 時才結束.根據應力-應變曲線,可以觀察到在發生屈曲時,沿著y方向壓縮的曲線會出現短暫的下降,隨后呈現平臺關系.且沿著y方向壓縮的平臺應力明顯高于沿著x方向壓縮的平臺應力,也就是說凹角蜂窩沿著y方向的屈曲強度高于沿著x方向的屈曲強度,這與圖11 的結果一致.將應力-應變曲線從彈性階段結束過渡為穩定平臺階段的這部分應力的平均值視為屈曲應力,并與數值和理論研究得到的屈曲強度進行對比,如表1所示.結果表明理論計算的屈曲強度與實驗和數值研究得到的屈曲強度存在一定的誤差,理論結果略低于實驗和數值研究的結果.這是因為理論研究是以凹角蜂窩單胞承受理想的應力狀態為前提,即圖8所示的一般應力狀態,并將這種一般應力狀態等效為單胞所受的集中荷載(即圖8(b)中所示的P,W,F),以此構建的矩陣公式(5)和(9).在實驗和數值研究中凹角蜂窩所承受的荷載為均勻荷載,這就使得實驗和數值研究得到的屈曲應力值與理論分析得到的應力值存在誤差.同時,還可以看到實驗和數值研究得到的屈曲應力值存在微小誤差,這是因為數值研究中增加了限制凹角蜂窩面外位移的約束而導致的.

圖13 凹角蜂窩沿著x方向和y方向單軸壓縮的應力-應變曲線Fig.13 Stress-strain curves of the re-entrant honeycomb in uniaxial compression along the x and y directions

表1 實驗、數值和理論研究得到的凹角蜂窩屈曲強度對比Table 1 Comparison of buckling strengths of re-entrant honeycombs obtained from experimental,numerical and theoretical studies

4 結論

通過有限元仿真分析發現了凹角蜂窩在宏觀應力狀態下的兩種屈曲模態.為了研究這兩種屈曲模態的屈曲強度以及產生的機理,本文采用梁柱理論對其進行了理論分析.根據梁-柱方程和平衡關系建立包含桿端彎矩和桿端轉角的方程組,利用屈曲臨界條件得穩定方程,進而得到屈曲強度的解析表達式.采用增材制造技術打印凹角蜂窩試件,進而對其屈曲性能進行實驗驗證.將理論結果與實驗和有限元結果進行對比分析.結果表明,通過改變加載條件,分別獲得凹角蜂窩的兩種屈曲模態,即屈曲模態Ⅰ和屈曲模態Ⅱ.當雙軸應力中x方向上的應力占據主導地位時,即σx>σy,凹角蜂窩優先產生屈曲模態Ⅰ.當雙軸應力中y方向上的應力占據主導地位時,即σx<σy,凹角蜂窩優先產生屈曲模態Ⅱ.凹角蜂窩具有不同于一般蜂窩的凹角幾何特征,從而產生了負Poisson 比效應;并且這種負Poisson 比效應使得凹角蜂窩在雙軸拉伸的應力狀態下依然會發生屈曲變形,而傳統正六邊形蜂窩卻不會出現這種現象.本文提出的凹角蜂窩兩種屈曲模態的穩定方程式(8)和(12),能夠有效預測凹角蜂窩在任意應力狀態下的屈曲模式,對凹角蜂窩因失穩而破壞以及利用凹角蜂窩的屈曲實現特殊力學性能的研究具有一定的指導意義.

猜你喜歡
模態方向有限元
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
磨削淬硬殘余應力的有限元分析
位置與方向
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 午夜福利网址| 久久黄色毛片| 黄色片中文字幕| 免费无码一区二区| 色综合综合网| 欧美在线导航| 欧美激情二区三区| 欧美日韩在线亚洲国产人| 亚洲午夜18| 国产18在线播放| 久久亚洲中文字幕精品一区| 91无码网站| 在线观看欧美精品二区| 亚洲综合在线最大成人| 色偷偷一区| 国产欧美视频在线| 亚洲伊人久久精品影院| 午夜福利免费视频| 美女一级毛片无遮挡内谢| 日韩福利视频导航| 精品欧美一区二区三区在线| 亚洲人成影视在线观看| 91区国产福利在线观看午夜| 手机在线国产精品| 午夜久久影院| 人妖无码第一页| 人妻一区二区三区无码精品一区| 真实国产乱子伦视频| 欧美黄网站免费观看| 国产乱人伦偷精品视频AAA| 九九免费观看全部免费视频| 亚洲中文字幕23页在线| 美美女高清毛片视频免费观看| 成人毛片免费在线观看| 国产在线91在线电影| 98超碰在线观看| av性天堂网| 亚洲伦理一区二区| 国产亚洲精品97在线观看| 国产一区二区免费播放| 国产伦精品一区二区三区视频优播| 在线亚洲小视频| www.91中文字幕| 国产女人18毛片水真多1| 97在线碰| 青草国产在线视频| 久草性视频| 国产成人精品视频一区二区电影| 国产欧美日本在线观看| 亚洲AV色香蕉一区二区| 老司机午夜精品网站在线观看| 亚洲人成高清| 久久婷婷色综合老司机| 亚洲男人天堂网址| 毛片网站在线看| 国产欧美网站| 欧美日韩高清在线| 伊人蕉久影院| 99re在线观看视频| 国产男女XX00免费观看| 老司机精品一区在线视频| 91小视频在线| 高清无码一本到东京热| 午夜国产理论| 欧美成人午夜影院| 中文成人无码国产亚洲| 亚洲视频影院| 亚洲性视频网站| 日韩欧美在线观看| 亚洲国产理论片在线播放| 久久青草免费91观看| 黄色网在线| 亚洲日韩Av中文字幕无码| 青青操视频在线| 日本一本正道综合久久dvd| 538精品在线观看| 亚洲大尺度在线| 国产精品极品美女自在线看免费一区二区 | 国产不卡一级毛片视频| 亚洲天堂在线免费| 日韩中文无码av超清| 激情综合五月网|