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

圓弧形蜂窩夾芯板在低速沖擊下的動力響應研究

2024-03-19 07:19:40陽,付
振動與沖擊 2024年5期
關鍵詞:有限元理論模型

余 陽,付 濤

(昆明理工大學 機電工程學院,昆明 650500)

蜂窩夾芯板具有質量輕、強度高、剛度大、能量吸收及減振等良好的性能,廣泛應用于航空航天和船舶等沖擊防護領域。根據結構的功能、使用壽命、可用性和成本的不同,蜂窩夾芯板的結構設計可以千變萬化。憑借這些優勢,蜂窩夾芯板受到研究者的青睞。在工程應用領域,蜂窩夾芯板容易受到鳥類或石子等物體的撞擊,給蜂窩夾芯板造成損傷,導致蜂窩夾芯板的剛度、穩定性和承載能力降低,進而影響結構的安全性[1-2]。因此,能夠準確地預測低速沖擊下蜂窩夾芯板動力響應的模型是必要的。在低速沖擊響應研究中,理論模型從彈道沖擊響應到空氣爆炸、水下爆炸沖擊響應,使用現有的理論和數值模型進行蜂窩夾芯板低速沖擊分析的理論推導,采用質量-彈簧單元、能量平衡原理和哈密爾頓原理等理論對模型進行求解,隨著理論的不斷豐富,蜂窩夾芯板的低速沖擊模型越來越精確[3-4]。

在過去幾十年里,有限元分析方法已經發展成為一種強大的技術,能夠模擬各種幾何形狀和不同沖擊載荷條件下的沖擊響應[5],將試驗結果和有限元仿真結果進行對比分析的研究方法受到研究者的青睞。Liu等[6]利用ANSYS基于時間步長積分的有限差分算法計算了沖擊位移、應力和應變隨時間的變化。姚鵬等[7]采用LS-DYNA軟件對夾層結構 (sandwich plate system,SPS)的沖擊吸能特性進行了模擬,并通過試驗數據驗證有限元仿真結果的可靠性。辛亞軍等[8]通過試驗和有限元模型對鋁蜂窩夾芯板的面外剪切行為和力學性能進行了研究。楊姝等[9]采用勢能駐值原理獲得集中力與局部變形關系,通過試驗、有限元及理論分析相結合方法對金屬面夾芯板抗沖擊性能進行研究。Zhang等[10]建立了考慮微觀結構的三維有限元模型,結合試驗和數值方法,推導了球面壓頭對夾芯板低速沖擊的壓痕響應和能量平衡解析模型。

在低速沖擊下夾芯板動力響應的數學模型研究中,Zhu等[11]基于能量法建立了數值分析模型,將分析模型結果與試驗結果進行對比,驗證模型的有效性。能量平衡模型能夠預測沖擊過程中夾芯板結構受到的最大沖擊力,但是這種方法不能提供沖擊力隨時間變化的趨勢[12-13],這是研究夾芯板結構損傷的重要信息。而利用最小勢能原理[14],通過總勢能的極小化,可以得到壓痕載荷與橫向撓度和變形長度的關系。在沖擊過程中,當夾芯板的總勢能達到其最小勢能時,凹陷深度最大[15]。隨著理論的不斷豐富,Andrews等[16]考慮動力學效應,將夾芯板建模為單自由度質量-彈簧系統,并與準靜態試驗結果進行了比較。金其多等[17]通過結合四階龍格庫塔法拓展了兩步攝動-伽遼金法,以此研究了含黏彈性夾芯的功能梯度石墨烯增強復合材料(functionally graded graphene reinforced composite,FG-GRC) 后屈曲梁在低速沖擊下的動力響應特性。朱秀杰等[18]基于精確板理論提供了一種分析復合材料格柵/波紋夾芯板屈曲特性的解析模型。Shu等[19]基于Mindlin正交各向異性板理論,得到指定邊界條件下波紋夾芯板在彎曲荷載作用下的 位移閉式解。Mishra等[20]利用模態求解技術,確定了沖擊點處的接觸力、板與沖擊器的側向位移和速度以及板內部的應力狀態。運動方程適用于特殊正交各向異性層合板的小撓度彈性響應。Sun等[21]利用哈密頓原理,提出了一種研究帶金屬面板的蜂窩夾芯板低速沖擊響應的解析建模方法,將哈密頓原理從彈性體擴展到具有塑性的面板。Malekzadeh等[22]基于改進的高階夾層板理論,采用彈簧-質量-阻尼器-減震器(spring-mass-damper-dashpot,SMDD)和彈簧-質量-阻尼器(spring-mass-damper,SMD)兩種新型三自由度系統對沖擊器與面板之間的相互作用進行了建模。

王彥斌等[23]對具有負泊松比特性的圓弧形蜂窩結構進行了力學分析研究,推導出圓弧形蜂窩胞元的等效參數。通過查閱相關文獻發現,對于蜂窩夾芯板在低速沖擊下的動力響應理論分析模型研究較少,而且圓弧形蜂窩夾芯板在低速沖擊下的動力響應還未被研究。因此為了能更好的研究圓弧形蜂窩夾芯板在低速沖擊中的動力響應,本文基于哈密頓原理和一階剪切變形理論來推導蜂窩夾芯板的運動方程,建立兩自由度的質量-彈簧模型,來研究沖擊器與蜂窩夾芯板之間的接觸力,采用Navier法和Duhamel積分對蜂窩夾芯板的運動方程進行解析求解。將模型計算結果與ABAQUS有限元仿真計算結果進行對比分析,驗證理論模型的準確性和有效性,同時討論不同的蜂窩胞元參數對蜂窩夾芯板動力響應的影響,可供其他研究者參考。

1 蜂窩夾芯板的低速沖擊理論模型

在建立圓弧形蜂窩夾芯板的低速沖擊理論模型之前,本文先假設蜂窩夾芯板是由上下蒙皮層和中間夾芯層組成,其中夾芯層是由負泊松比特性的圓弧形蜂窩胞元經過線性陣列而組成。在蜂窩夾芯板的中間平面建立笛卡爾坐標系(x,y,z),沿著x、y和z軸方向分別為蜂窩夾芯板的長度、寬度和高度,蜂窩夾芯板的長度、寬度和高度分別記為Lx、Ly和h。蜂窩夾芯板受到初始速度為vs,半徑為Rs的球形沖擊器的沖擊載荷,低速沖擊下的圓弧形蜂窩夾芯板的模型如圖1所示。為了更好地計算分析,假設頂部蒙皮層的高度ht和底部蒙皮層的高度hb相等,夾芯層的高度為hc。蜂窩夾芯板的圓弧形蜂窩胞元結構如圖2所示。圖2中:d為胞元的壁厚;Rc為胞元的半徑;θ為胞元的角度。引用王彥斌等所推導的圓弧形蜂窩結構的等效參數為

圖1 低速沖擊下蜂窩夾芯板模型

圖2 蜂窩胞元結構

(1)

式中:a=2(2Rc+d)sinθ,b=2(2Rc+d)cosθ,S1=8π(2Rcd+d2)θ/360°+4dRc,S2=2ab;E、G和μ分別為基體材料的彈性模量、剪切模量和泊松比。

1.1 蜂窩夾芯板的運動方程

為了獲得蜂窩夾芯板的位移場,采用一階剪切變形理論,位移場方程為

(2)

式中:u,v,w分別為中間平面相對于x軸,y軸和z軸平面內的橫向位移;φx和φy為旋轉位移分量;t為時間。

與位移場有關的任意位置的橫向應變-位移關系為

(3)

式中:εx、εy為正應變;γxy、γxz、γyz為剪切應變。

蜂窩夾芯板的線性本構關系

(4)

(5)

(6)

式中,標量符號的上角標t、c和b分別為頂部蒙皮層、圓弧形夾芯層和底部蒙皮層所對應的參數值。

其剛度系數為

(7)

根據哈密頓原理,蜂窩夾芯板一階剪切變形理論的運動方程為

(8)

式中:q為沖擊載荷;I0、I1和I2為質量慣性矩;Nx、Ny和Nxy為軸向力;Qx和Qy為剪切力;Mx、My和Mxy為力矩。

軸向力、力矩和剪切力分別為

(9)

(10)

1.2 接觸力的求解方法

為了能夠計算出球形沖擊器與蜂窩夾芯板之間的接觸力,兩自由度質量-彈簧模型如圖3所示。球形沖擊器的質量為m1,蜂窩夾芯板的質量為m2,F(t)為球形沖擊器在沖擊下的接觸力,δ(t)為接觸后蜂窩夾芯板的接觸位移。接觸位移由沖擊器和蜂窩夾芯板之間的位移差確定,如圖4所示。

圖3 兩自由度彈簧質量模型

圖4 球形沖擊器沖擊下的接觸力和接觸位移示意圖

接觸位移δ(t)為

δ(t)=w1(t)-w2(t)

(11)

式中,w1(t)和w2(t)為兩個質量塊在撞擊后時間的位移響應。

根據參考文獻[24],沖擊器在沖擊期間的非線性赫茲沖擊力可推導為

(12)

式中,λ為彈性恢復常數,取1.5。

(13)

式中,Es、μs和Et、μt分別為球形沖擊器和蜂窩夾芯板頂部蒙皮層的楊氏模量和泊松比。

使用等效剛度K1解析求解,接觸力為

F(t)=δ(t)K1

(14)

引用參考文獻[25]等效接觸剛度為

(15)

式中,Γ(·)為Gamma函數。

δm最大接觸變形為

(16)

應用牛頓第二運動定律,兩自由度質量-彈簧系統的平衡方程為

(17)

平衡方程的初始條件為

(18)

由式(11)、(14)、(17)和(18)可以推導出球形沖擊器與蜂窩夾芯板之間的接觸力為

(19)

其中

1.3 橫向位移的求解方法

本文將蜂窩夾芯板視為四邊簡支邊的矩形板,在這種情況下,其邊界條件為

v=w=φy=Nx=Mx=0 當x=0或x=Lx

(20)

u=w=φx=Ny=My=0 當y=0或y=Ly

(21)

四邊簡支夾芯板的求解可以通過Navier法獲得,將位移假定為如下雙三角級數

(22)

式中:α=mπ/Lx,β=nπ/Ly;Umn(t)、Vmn(t)、Wmn(t)、Xmn(t)和Ymn(t)為位移振幅分量;m和n為模數。

沖擊載荷q可表示為

(23)

其中,系數Qmn(t)為

(24)

將式(22)、(23)代入式(8)中,并應用Galerkin方法,忽略平面內和旋轉慣性,可推導為

(25)

根據式(25),可得到一個線性二階微分方程為

(26)

式中,ωmn為蜂窩夾芯板的固有頻率。

使用Duhamel積分,將式(24)代入式(26),可推導出Wmn(t)為

(27)

使用式(19)、(22)和(27),可以得到蜂窩夾芯板的橫向位移方程

(28)

式中,e1和e2為代替系數。

1.4 理論模型驗證

在本節中,將理論模型計算結果與ABAQUS有限元仿真結果或已發表文獻的結果進行對比分析,來驗證當前理論模型可靠性和有效性。首先將理論模型與有限元仿真模型所計算的蜂窩夾芯板中心最大橫向位移結果進行對比驗證,蜂窩夾芯板的幾何尺寸為:蜂窩夾芯板的長度與寬度分別為Lx=300 mm和Ly=300 mm,蜂窩夾芯板的總高度h=12 mm,夾芯層高度hc=10 mm,蜂窩胞元壁厚d=2 mm,蜂窩胞元角度θ=45°,蜂窩胞元半徑Rc=6 mm。材料參數為:蒙皮層和芯層為相同材料,材料的彈性模量E=69 GPa,泊松比μ=0.3,密度ρ=2 700 kg/m3。球形沖擊器的半徑為10 mm,其彈性模量E=200 GPa,泊松比μ=0.3,密度ρ=7 871.8 kg/m3。

建立ABAQUS有限元仿真模型,首先在有限元仿真軟件中對球形沖擊器和圓弧形蜂窩夾芯板進行3D建模,并將材料屬性賦予模型。其次對沖擊器和蜂窩夾芯板兩個部件進行裝配,將球形沖擊器參考點P與蜂窩夾芯板頂部蒙皮層中心點相接觸,其接觸點坐標為(Lx/2,Ly/2,h),同時將參考點P與球形沖擊器進行剛體約束,在分析開始時將參考點調整到質心,然后對參考點施加速度載荷,通過改變球形沖擊器的速度參數就可以達到不同速度的沖擊效果。為了防止低速沖擊過程中發生干預穿透問題,采用通用接觸,面面接觸在法向設置為硬接觸,在切向設置為罰函數接觸,其摩擦因數為0.2。然后創建動力顯示分析步,對蜂窩夾芯板施加四邊簡支載荷邊界條件,分別在x=0和x=Lx邊,限制y方向與z方向的位移和繞y軸的旋轉;在y=0和y=Ly邊,限制x方向與z方向的位移和繞x軸的旋轉。采用十結點二次四面體單元(C3D10)對圓弧形蜂窩夾芯板進行網格劃分,其邊界網格尺寸為5 mm,為了更好的分析蜂窩夾芯板橫向位移,對蜂窩夾芯板的中心20 mm×20 mm區域網格尺寸設置為2.5 mm。球形沖擊器同樣采用十結點二次四面體單元(C3D10),其網格尺寸為1 mm。有限元軟件建立的圓弧形蜂窩夾芯板低速沖擊仿真模型網格劃分如圖5所示。

圖5 有限元仿真模型網格劃分

采用相對收斂性檢查的方法來驗證當前網格的收斂性,逐步調整中心網格尺寸,每次以近似2∶1的比例調整網格尺寸。在10 m/s的沖擊速度下,中心網格尺寸分別設計為5 mm、2.5 mm、1.2 mm和0.6 mm,蜂窩夾芯板中心最大橫向位移分別為0.209 mm、0.305 mm、0.307 mm和0.307 mm,有限元數值模擬計算的時間分別為90 s、90 s、141 s和360 s。使用MATLAB軟件中Smoothing Spline函數對離散數據進行曲線擬合,蜂窩夾芯板的中心網格尺寸對中心最大橫向位移和計算時間的影響如圖6所示。由圖6可知,隨著網格尺寸的逐漸減小,夾芯板中心最大橫向位移也逐漸趨于一個穩定值0.307 mm,其5%的誤差下橫向位移為0.262 mm,網格尺寸為2.5 mm時可以滿足網格收斂性的要求,同時此網格尺寸下,有限元數值模擬計算出結果所花費的時間也相對較少。通過上述分析可以得出當中心網格尺寸為2.5 mm時,滿足網格收斂性的要求,且有限元數值模擬計算效率較高。

圖6 中心網格尺寸對最大橫向位移和計算時間的影響

在不同的沖擊載荷下,采用上述有限元仿真模型與本文理論模型所計算的中心最大橫向位移如表1所示。由表1可知,有限元仿真結果與本文理論模型計算結果最大的相對誤差為4.9%,驗證了本文方法的正確性。

表1 圓弧形夾芯板中心最大橫向位移

最后利用本文所建立的理論模型求解方法來計算球形沖擊器和矩形板的接觸力,并與Yang等[25]和Wu等[26]計算的接觸力進行對比研究。矩形板結構的長寬高分別為200 mm、200 mm和8 mm,球形沖擊器半徑為10 mm,沖擊器與矩形板的材料密度ρ=7 971.8 kg/m3,彈性模量E=200 GPa,泊松比μ=0.3,沖擊器以1 m/s的速度沖擊矩形板。本文研究預測的接觸力分布與文獻[25]和文獻[26]得到的接觸力分布對比如圖7所示,從圖中可以看,本文的分析模型預測的最大接觸力高于文獻[26]的試驗結果,低于文獻[25]的試驗結果,矩形板與球形沖擊器之間的接觸力總體變化趨勢相同,理論模型與參考文獻的接觸力最大誤差為8%,驗證了當前理論模型的有效性。

圖7 理論模型與參考文獻的接觸力隨時間變化的對比

2 蜂窩夾芯板的參數分析

在本章中,本文通過所建立的理論模型來研究低速沖擊下蜂窩胞元的不同幾何參數對圓弧形蜂窩夾芯板動力響應的影響。假設圓弧形蜂窩夾芯板材料的密度ρ=2 700 kg/m3,彈性模量E=69 GPa和泊松比μ=0.3,蜂窩夾芯板的長度和寬度為Lx=Ly=100 mm,高度hc=10 mm,頂部蒙皮層和底部蒙皮層高度相等且為1 mm。蜂窩胞元角度θ=45°,蜂窩胞元半徑Rc=6 mm,蜂窩胞元壁厚d=2 mm。球形沖擊器的半徑為10 mm,沖擊速度為8 m/s,沖擊器材料的密度ρ=7 971.8 kg/m3,彈性模量E=200 GPa,泊松比μ=0.3。

2.1 蜂窩胞元半徑對蜂窩夾芯板動力響應的影響

本文假設蜂窩胞元半徑分別為5 mm、6 mm和7 mm,在上述其它幾何參數不變的情況下,不同蜂窩胞元半徑對蜂窩夾芯板動力響應的影響如圖8所示。由圖8(a)可知,隨著蜂窩胞元半徑的增加,蜂窩夾芯板的中心最大橫向位移也在不斷增加,表明在相同的沖擊速度下,蜂窩夾芯板的抗沖擊性能反而減小。由圖8(b)可知,在相同的沖擊速度下,隨著蜂窩胞元半徑的增大,沖擊器與蜂窩夾芯板之間的接觸力也在增大。在蜂窩胞元其它參數確定的情況下,隨著蜂窩胞元半徑的增加,蜂窩胞元呈現出擴張的趨勢,進而導致蜂窩胞元的抗沖擊性能減弱。在相同的沖擊載荷下,當蜂窩胞元半徑從5 mm增加至7 mm時,蜂窩夾芯板的中心最大橫向位移從0.211 mm增加至0.296 mm,夾芯板的抗沖擊性能減少40.28%。沖擊器與蜂窩夾芯板之間的接觸力從12 161.2 N增加至12 389.3 N,接觸力增加1.88%。

(a) 橫向位移隨時間變化曲線

2.2 蜂窩胞元角度對蜂窩夾芯板動力響應的影響

本文假設蜂窩胞元的角度分別為30°、45°和60°,在其它幾何參數不變的情況下,不同蜂窩胞元角度對蜂窩夾芯板動力響應的影響如圖9所示。

(a) 橫向位移隨時間變化曲線

由圖9可知,在相同的沖擊載荷下,圓弧形蜂窩夾芯板的中心最大橫向位移隨著蜂窩胞元角度的增大而增大。在相同的蜂窩胞元參數下,蜂窩胞元角度的減小可以提升胞元結構的緊密性,使蜂窩胞元呈現出一種收縮的現象,進而提升圓弧形蜂窩夾芯板的抗沖擊特性。當蜂窩胞元角度從30°增加至60°時,蜂窩夾芯板的中心最大橫向位移從0.165 mm增加至0.303 mm,蜂窩夾芯板的抗沖擊性能減少83.64%;蜂窩夾芯板與沖擊器之間的接觸力從12 184.5 N增加至12 545.5 N,接觸力增加2.96%。

2.3 蜂窩胞元壁厚對蜂窩夾芯板動力響應的影響

本文假設蜂窩胞元的壁厚分別為1 mm、2 mm和3 mm,在其它幾何參數保持不變的情況下,不同蜂窩胞元壁厚對蜂窩夾芯動力響應的影響如圖10所示。由圖10可知,在相同的沖擊載荷下,圓弧形蜂窩夾芯板的中心最大橫向位移隨著蜂窩胞元壁厚的增大而減小,當蜂窩胞元壁厚從1 mm增加至3 mm時,蜂窩夾芯板中心最大橫向位移0.368 mm減少至0.149 mm,夾芯板的抗沖擊性能提升59.51%;蜂窩夾芯板與沖擊器之間的接觸力從12 524.8 N減少至11 758.8 N,接觸力減小6.12%。在蜂窩胞元其它參數確定情況下,蜂窩胞元的壁厚增加可以提升夾芯層與蒙皮層之間的接觸面積,進而提升圓弧形蜂窩夾芯板的抗沖擊特性。

(a) 橫向位移隨時間變化曲線

2.4 材料密度對蜂窩夾芯板動力響應的影響

假設圓弧形蜂窩夾芯板的長度和寬度為Lx=Ly=100 mm,夾芯層高度hc=10 mm,頂部和底部蒙皮層高度相等且為1 mm。蜂窩胞元角度θ=45°,胞元半徑Rc=6 mm,胞元壁厚d=2 mm。本文選用常見的軋制鋁、灰鑄鐵和合金鋼作為蜂窩夾芯板的材料,三種材料的參數如表2所示。球形沖擊器材料的密度ρ=7 971.8 kg/m3,彈性模量E=200 GPa,泊松比μ=0.3,沖擊器的半徑為10 mm,沖擊速度為8 m/s。

表2 蜂窩夾芯板的材料參數設置

在上述蜂窩夾芯板的材料參數和幾何參數設計下,材料密度對蜂窩夾芯板低速沖擊下動力響應的影響如圖11所示。由圖11(a)可知,隨著蜂窩夾芯板材料密度的增大,蜂窩夾芯板中心最大橫向位移在減小,軋制鋁、灰鑄鐵和合金鋼三種蜂窩材料的中心最大橫向位移分別為0.261 mm、0.166 mm和0.11 mm,可以得出隨著蜂窩夾芯板密度的增大,蜂窩夾芯板的抗沖擊特性明顯增強。由圖11(b)可知,隨著材料密度的增大,蜂窩夾芯板與球形沖擊器之間的接觸力在明顯增大。

(a) 橫向位移隨時間變化曲線

3 結 論

(1) 在6種不同沖擊速度下理論模型與ABAQUS有限元仿真軟件計算出的圓弧形蜂窩夾芯板中心最大橫向位移的最大誤差為4.9%;在沖擊速度為1 m/s時,理論模型與參考文獻[25-26]求解出的最接觸力的最大誤差為8%,驗證了當前理論模型的有效性。

(2) 通過理論模型計算出,隨著球形沖擊器的沖擊速度遞增,圓弧形蜂窩夾芯板中心最大橫向位移也呈現出遞增的規律。當沖擊速度從4 m/s提升至14 m/s時,夾芯板中心最大橫向位移從0.107 mm增加至0.452 mm。

(3) 圓弧形蜂窩夾芯板的抗沖擊特性隨著蜂窩胞元半徑和角度的增大而減小,隨著蜂窩胞元壁厚的增大而增大。當蜂窩胞元半徑從5 mm增加至7 mm時,蜂窩夾芯板的抗沖擊特性減少40.32%;當蜂窩胞元角度從30°增加至60°時;蜂窩夾芯板的抗沖擊特性減少80.82%;當蜂窩胞元壁厚從1 mm增加至3 mm時,蜂窩夾芯板的抗沖擊特性提升59.57%。

(4) 材料密度對圓弧形蜂窩夾芯板的抗沖擊特性也有顯著的影響,隨著蜂窩夾芯板密度的增加,蜂窩夾芯板的抗沖擊特性也會增強。

猜你喜歡
有限元理論模型
一半模型
堅持理論創新
當代陜西(2022年5期)2022-04-19 12:10:18
神秘的混沌理論
理論創新 引領百年
相關于撓理論的Baer模
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 国模私拍一区二区三区| 九色在线观看视频| 亚洲女同欧美在线| 超碰色了色| 亚洲国产精品日韩av专区| 成年A级毛片| 欧美一级大片在线观看| 欧美成人二区| 久久精品人人做人人爽97| 91探花在线观看国产最新| 国产精女同一区二区三区久| 日韩成人在线一区二区| 欧美精品亚洲二区| 国产精品视频观看裸模| 成年片色大黄全免费网站久久| 国产精品亚洲一区二区三区z | 国产精品短篇二区| 亚洲精品视频免费看| 性喷潮久久久久久久久| 亚洲色图欧美| 男女男免费视频网站国产| 无码电影在线观看| 超碰91免费人妻| 91在线高清视频| 国产一二视频| 亚洲视频欧美不卡| 欧美高清国产| 亚洲综合中文字幕国产精品欧美| 91无码人妻精品一区二区蜜桃| 国产精品13页| 美女国产在线| 亚洲欧美不卡视频| 欧美日韩在线亚洲国产人| 99免费视频观看| 99在线视频免费观看| 欧美午夜精品| 国产精品综合色区在线观看| 在线另类稀缺国产呦| 精品久久香蕉国产线看观看gif| 色爽网免费视频| 久久人搡人人玩人妻精品| 97超级碰碰碰碰精品| 欧美亚洲一二三区 | 色综合久久久久8天国| 91精品在线视频观看| 青草免费在线观看| 呦女精品网站| 色网站免费在线观看| 婷婷99视频精品全部在线观看| 日韩欧美在线观看| 日韩国产黄色网站| 97se综合| 欧美黄色a| 久久综合丝袜长腿丝袜| 亚洲人成色77777在线观看| 国国产a国产片免费麻豆| 国产成人高清精品免费软件| 亚洲日韩精品无码专区| www.99精品视频在线播放| 久久精品免费国产大片| 台湾AV国片精品女同性| 无码区日韩专区免费系列| 国产在线观看一区二区三区| 91无码人妻精品一区二区蜜桃| 国产高清色视频免费看的网址| 国产国产人成免费视频77777 | 九色91在线视频| 高清欧美性猛交XXXX黑人猛交 | 国产亚洲成AⅤ人片在线观看| 伊人国产无码高清视频| 九九九国产| 欧美第二区| 国产成人艳妇AA视频在线| 中文字幕在线观| 狠狠做深爱婷婷久久一区| 91精品在线视频观看| 国产伦片中文免费观看| 中文国产成人久久精品小说| 亚洲成人在线网| 欧美翘臀一区二区三区| 亚洲欧洲日韩久久狠狠爱| 91年精品国产福利线观看久久|