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

爆炸載荷下固支矩形板的大撓度塑性動力響應(yīng)

2013-06-07 02:53:46劉敬喜
中國艦船研究 2013年1期
關(guān)鍵詞:有限元效應(yīng)變形

顏 豐,劉敬喜

華中科技大學(xué)船舶與海洋工程學(xué)院,湖北武漢 430074

爆炸載荷下固支矩形板的大撓度塑性動力響應(yīng)

顏 豐,劉敬喜

華中科技大學(xué)船舶與海洋工程學(xué)院,湖北武漢 430074

從矩形板的小撓度運動方程出發(fā),通過引入膜力因子,給出四邊固支矩形板在大撓度變形情況下的運動方程,分析矩形板大撓度塑性動力響應(yīng),并根據(jù)運動方程導(dǎo)出在矩形脈沖載荷作用下四邊固支矩形板的運動微分方程,求解矩形板的最大殘余變形計算式。同時,通過假設(shè)的應(yīng)變率效應(yīng)系數(shù)選取方法,解決大撓度加載情況下材料屈服應(yīng)力的增加問題。使用有限元仿真手段驗證了帶有移行鉸線的變形機(jī)構(gòu),對已有的實驗樣本和補充的有限元模型進(jìn)行計算,并將計算出的理論結(jié)果與已有實驗結(jié)果和有限元結(jié)果進(jìn)行了比對,吻合較好。

矩形板;塑性動力響應(yīng);爆炸載荷;矩形脈沖;膜力因子;應(yīng)變率效應(yīng)

0 引 言

爆炸載荷作用下矩形板結(jié)構(gòu)的塑性動力響應(yīng)研究很早便受到理論界與工程界的廣泛關(guān)注。對于矩形板這種二維的非軸對稱結(jié)構(gòu),因其主應(yīng)力方向事先不知道,因而使得有關(guān)矩形板的理論分析相當(dāng)復(fù)雜和困難,尤其在大撓度情形下,既有材料非線性,又有幾何非線性。盡管如此,近年來仍有不少學(xué)者在研究有關(guān)爆炸載荷下矩形板結(jié)構(gòu)的大撓度塑性變形問題。例如,朱錫等[1]對爆炸載荷下的固支方板進(jìn)行了理論分析和試驗分析,導(dǎo)出了固支方板在爆炸載荷作用下的應(yīng)變場,分析了方板的破裂形式,并給出了破裂臨界壓力值。Chung和 Nurick[2-3]對爆炸載荷下方板、單加筋、雙加筋、十字加筋和雙十字筋方板的塑性動力響應(yīng)進(jìn)行了一系列實驗研究,獲得了相當(dāng)數(shù)量的爆炸載荷下加筋板大撓度變形和撕裂的實驗結(jié)果,并對實驗樣本進(jìn)行了有限元仿真。Park和Cho[4]對爆炸載荷下矩形板和加筋矩形板的最終變形提出了新的經(jīng)驗公式,采用MSC/DYTRAN軟件進(jìn)行了有限元仿真,并與數(shù)值結(jié)果和實驗結(jié)果進(jìn)行了比對。章媛和茍瑞君[5]根據(jù)塑性變形理論,推導(dǎo)出了矩形板在水下爆炸載荷作用下中心處撓度的解析解,利用LS-DYNA有限元軟件對矩形板的水下爆炸沖擊過程進(jìn)行了數(shù)值仿真,并進(jìn)行了相關(guān)的水下爆炸試驗。何建等[6]運用能量守恒原理分析了四邊固支矩形鋼板的剛塑性大撓度變形的動力響應(yīng),推導(dǎo)了最大塑性變形計算公式,并與LS-DYNA的有限元仿真結(jié)果和前人的試驗結(jié)果進(jìn)行了比對,獲得了令人滿意的結(jié)果。

早在1992年,余同希和陳發(fā)良[7]就曾考慮過大撓度響應(yīng)的移行塑性鉸相,并在角動量守恒方程中同時考慮了彎矩和大撓度誘導(dǎo)的膜力產(chǎn)生的合力矩的作用,進(jìn)行了較為完整的理論分析,但只給出了瞬動載荷作用下的算例,而沒有給出矩形脈沖載荷作用下的理論分析結(jié)果。

本文參照文獻(xiàn)[7]提出的膜力因子法,將對固支矩形板的塑性動力響應(yīng)進(jìn)行研究,推導(dǎo)同時考慮彎矩和膜力聯(lián)合作用下的固支矩形板運動模式及運動方程,給出矩形脈沖載荷作用下的計算結(jié)果,并將與有限元仿真和Nurick的試驗結(jié)果進(jìn)行對比,以得出一些有意義的結(jié)論。

1 矩形板小撓度變形分析

假設(shè)一寬為a,長為b,厚度為H的矩形板,其周邊為固支,單位面積質(zhì)量為μ,其上承受均布動力載荷P(t)的作用,不計自身板重,如圖1所示。板的單位長度塑性極限彎矩為

式中,σs為屈服應(yīng)力。

圖1 四邊固支矩形板Fig.1 Fully clamped rectangular plate

1.1 靜力極限狀態(tài)下變形機(jī)構(gòu)的選取

首先討論四邊固支矩形板處于靜力極限狀態(tài)時的破壞機(jī)構(gòu),圖2和圖3所示分別為兩種可能的變形機(jī)構(gòu)形式。

圖2 變形機(jī)構(gòu)IFig.2 Deformationmechanism,pattern I

圖3 變形機(jī)構(gòu)IIFig.3 Deformationmechanism,pattern II

變形機(jī)構(gòu)I的靜力極限方程為

變形機(jī)構(gòu)II的靜力極限方程為

將方程組(3)中的ξp消去并聯(lián)立方程,得

1.2 小撓度變形的運動控制方程

假設(shè)在某一瞬間均布載荷P(0)=P*>P0,此時,板開始發(fā)生變形。一般而言,板變形應(yīng)優(yōu)先考慮變形機(jī)構(gòu)II,即在板中部有一個發(fā)生平動運動的塑性區(qū),在它四周的四塊剛性板則分別繞各自的邊界作剛體轉(zhuǎn)動。

矩形板小撓度變形的變形機(jī)構(gòu)II的運動方程組為

式中,τ為無因次時間

?1為無因次轉(zhuǎn)角

p(τ)為無因次均布載荷

2 固支矩形板大撓度變形分析

2.1 引入膜力因子建立大撓度變形運動控制方程

當(dāng)固支矩形板變形的撓度達(dá)到甚至超過板厚,同時平板邊緣在平面方向上為固定而無法移動時,就需要計入大撓度誘導(dǎo)的膜力效應(yīng),因為膜力將會通過耗散變形能量增強(qiáng)結(jié)構(gòu)剛性而對變形產(chǎn)生影響。

引入隨矩形板變形而增大的膜力因子 f1和f2

式中,w0為中央平臺區(qū)的位移;η為無因次化后的中央平臺區(qū)位移。

則考慮膜力因素后,矩形板大撓度變形的變形機(jī)構(gòu)II的運動方程組為:

在本文討論的爆炸載荷下,矩形板先以變形機(jī)構(gòu)II的方式運動,并發(fā)生鉸線移行,當(dāng)鉸線移行至δ=時,轉(zhuǎn)化為機(jī)構(gòu)I的運動方式。

2.2 矩形脈沖加載下矩形板的動態(tài)微分方程

可將板的變形分解為若干階段。

2.2.1 第1階段

0≤τ≤τT,τT為脈沖加載時間,p(τ)=p* ,I(τ)=p*τ,運動微分方程為

由此,ξ和δ的初值ξ0和δ0的近似解便可由方程組(15)用數(shù)值方法求得。再將 ξ0,δ0和式(13)聯(lián)立,使用龍格—庫塔積分獲得ξ和δ的值。

此外,矩形板區(qū)域①繞邊界軸的轉(zhuǎn)角為

2.2.2 第2階段

ξ和δ的初值由第1階段計算獲得,和式(18)聯(lián)立使用龍格—庫塔積分獲得ξ和δ的值。

矩形板區(qū)域①繞邊界軸的轉(zhuǎn)角及轉(zhuǎn)角速度計算同式(16)和式(17)。

2.2.3 第3階段

τ2<τ≤τ3,其中 τ=τ3時,ξ=0 。運動微分方程組為

其中,η的初值為

ξ和?1的初值由第2階段獲得,和式(20)聯(lián)立使用龍格—庫塔積分獲得ξ,φ1和η的值。

矩形板區(qū)域①繞邊界軸的轉(zhuǎn)角計算同式(16)。

2.2.4 第4階段

τ3<τ≤ τf,ξ=0 即 ξ=ξf,常數(shù)。式(20)可簡化為

其中,η的初值為

?1的初值由第3階段計算獲得,將?1和η的初值與式(22)聯(lián)立,使用龍格—庫塔積分獲得?1和η的值。

3 計算實例與結(jié)果分析

3.1 矩形板變形的有限元仿真

采用有限元分析軟件ABAQUS/Explicit v6.8-1進(jìn)行數(shù)值模擬。模型參照文獻(xiàn)[2]中的實驗樣本建立,其面板的長和寬均取原文的實驗樣本數(shù)據(jù),而外加均布載荷則為

式中,T為載荷加載時間。本文亦給出了14.5μs的估計值。

自建模型,其面板厚度均取1.6 mm,面積約16 000mm2。建模時,在各矩形板和加筋板的面板外加一外框,再如圖4所示將邊界約束施加于外框上。自建模型的外加均布載荷為矩形脈沖,幅值150 MPa,持續(xù)時間14.5 μs。采用彈塑性材料模型,楊氏模量E=210GPa,材料密度 ρ=7 850 kg/m3,靜態(tài)屈服應(yīng)力為242MPa;考慮材料應(yīng)變率效應(yīng)的影響,D=40.4 s-1,q=5[8]。

圖4 有限元模型邊界約束和加載示意圖Fig.4 Boundary conditions and loadsof finite elementmodel

為了驗證前文所提出的固支矩形板的兩種變形機(jī)構(gòu)及其轉(zhuǎn)化,取長寬比為1.6∶1的有限元模型U02作為參照,有限元變形仿真結(jié)果如圖5和圖6所示。

圖5 模型U02有限元仿真變形總體示意圖Fig.5 Deformation of finite elementmodel U02

圖6 模型U02有限元仿真中點截面變形圖Fig.6 Section deformation of finite elementmodel U02

在以上矩形板變形圖中,從變形起始至111 μs左右,中點截面圖的板截面呈平頂狀突起,這說明在固支矩形板的變形過程中出現(xiàn)了中央平臺區(qū)。而由變形總體示意圖則可以看出,其中央平臺區(qū)大致呈矩形,并且其大小是隨時間的后移逐漸變小,這說明鉸線移行在有限元仿真中也是存在的。

在111μs和165μs之間的某個時刻,中點截面圖中的平頂轉(zhuǎn)變?yōu)榻萍忭敚从吵鲎冃螜C(jī)構(gòu)II的中央平臺區(qū)消失,轉(zhuǎn)換為對應(yīng)尖頂?shù)闹醒脬q線,這說明在有限元仿真中同樣存在變形機(jī)構(gòu)II向變形機(jī)構(gòu)I的轉(zhuǎn)換。

3.2 應(yīng)變率效應(yīng)系數(shù)的選取

當(dāng)對矩形板進(jìn)行實際動力加載時,其材料的屈服應(yīng)力將會隨著應(yīng)變率的增大而增大,稱之為應(yīng)變率效應(yīng)。若應(yīng)變率效應(yīng)較為顯著,則此時式(1)中的屈服應(yīng)力應(yīng)采用動態(tài)屈服應(yīng)力,其表達(dá)式為

式中,σs為材料的靜態(tài)屈服應(yīng)力;α為應(yīng)變率效應(yīng)系數(shù),表征動態(tài)屈服應(yīng)力較靜態(tài)屈服應(yīng)力的放大倍數(shù)。

關(guān)于應(yīng)變率效應(yīng)系數(shù)的取值,Cloete等[9]曾針對周邊固支中央簡支的圓板,將具體實驗結(jié)果與塑性功—動能守恒條件相結(jié)合,推導(dǎo)出其參考值為2.8。本文的理論計算所依據(jù)的變形機(jī)構(gòu)因含有平臺區(qū)及移行鉸線,而文獻(xiàn)[9]在固支圓板的理論分析中并未考慮平臺區(qū)和移行鉸線,所以不宜直接采用上述參考值。此外,由于移行鉸線的存在,使得對矩形板在不同時刻的整體應(yīng)變推導(dǎo)較為困難。

針對應(yīng)變率效應(yīng)系數(shù)推導(dǎo)的難點,本文擬從相反方向著手解決,采用1.0~3.0范圍內(nèi)的各應(yīng)變率效應(yīng)系數(shù)同時進(jìn)行理論計算,并將獲得的不同計算值與實驗數(shù)據(jù)或有限元結(jié)果進(jìn)行比對,從而選擇出合適的應(yīng)變率效應(yīng)系數(shù)。

文獻(xiàn)[9]在進(jìn)行應(yīng)變率效應(yīng)系數(shù)的驗證時,使用了 Bodner和 Symonds[10],以及 Nurick[11]在實驗中所獲得的達(dá)到變形峰值的響應(yīng)時間。本文相應(yīng)地也對計算樣本建立了有限元模型,而后將通過有限元仿真獲得的達(dá)到變形峰值的響應(yīng)時間作為主要參考值,再將其和各應(yīng)變率效應(yīng)系數(shù)所對應(yīng)的板變形峰值響應(yīng)時間計算值進(jìn)行比對,以選出適合該計算模型的應(yīng)變率效應(yīng)系數(shù)。此外,由于在有限元仿真中將板材料定義為彈塑性模型而未設(shè)置其阻尼系數(shù),使得矩形板模型在變形后期發(fā)生了反復(fù)的小幅振動,因此,此處取其第1次達(dá)到變形峰值的響應(yīng)時間。同時,由于本文的理論計算所采用的板材料假設(shè)為剛塑性模型,故其變形的總響應(yīng)時間即為平板達(dá)到變形峰值的響應(yīng)時間。

以文獻(xiàn)[2]中的試驗樣本S01為例,以有限元計算獲得的板變形峰值時間為參考值,另選取應(yīng)變率效應(yīng)系數(shù)1.0~3.0進(jìn)行計算,獲得了對應(yīng)的板變形響應(yīng)時間理論計算值—系數(shù)關(guān)系曲線如圖7所示。

圖7 模型S01板變形響應(yīng)時間的理論計算與有限元結(jié)果Fig.7 Theoreticaland finite element results of deformation time ofmodel S01

由圖7可以看出,當(dāng)理論計算采用的應(yīng)變率系數(shù)為2.3時,平板變形峰值響應(yīng)時間的理論計算值與有限元結(jié)果較為近似。所以,在計算板中點位移時,選取2.3作為樣本S01的應(yīng)變率效應(yīng)系數(shù)。

3.3 與試驗結(jié)果的比對及分析

按照以上提出的應(yīng)變率效應(yīng)系數(shù)選取方法計算文獻(xiàn)[2]的方板樣本,并與實驗值進(jìn)行比對,結(jié)果如表1所示。

表1 固支方板樣本理論計算與實驗結(jié)果對比Tab.1 Resu ltsofm id-point disp lacem ent of theoretical calcu lation and experim ents

由表1的結(jié)果對比可以看出,在選取適當(dāng)應(yīng)變率效應(yīng)系數(shù)的前提下,方板大撓度變形的中點位移的理論計算值與實驗值吻合良好。

3.4 有限元補充模型計算結(jié)果

文獻(xiàn)[7]運用膜力因子法對Jones在1971年的試驗[12]中獲得的、無因次變形范圍在0~5內(nèi)的數(shù)據(jù)進(jìn)行了理論計算與試驗數(shù)據(jù)比對,而前文所計算的文獻(xiàn)[2]中的試驗樣本均為長寬比為1∶1的方板,為進(jìn)一步考察固支矩形板在大撓度變形下的響應(yīng),本文利用有限元仿真分析另給出了長寬比在1.3~4范圍且無因次變形均大于5的5個模型作為補充算例,所得計算結(jié)果如表2所示。

表2 固支矩形板模型理論計算與有限元仿真結(jié)果對比Tab.2 Resu lts ofm id-poin t disp lacem en t of theoretical calcu lation and finite elem ent sim ulation

由表2的結(jié)果對比可以看出,在選取適當(dāng)應(yīng)變率效應(yīng)系數(shù)的前提下,即使在模型長寬比達(dá)4∶1的情況下,固支矩形板變形的中點位移的理論計算值與有限元計算值仍然符合較好。

4 結(jié) 論

在矩形板的塑性動力響應(yīng)過程中,當(dāng)其撓度達(dá)到或超過厚度量級時,由于膜力的作用,使得僅考慮彎矩作用的小撓度理論不再適用,此時,必須引入膜力作用的影響來研究矩形板的大撓度塑性變形響應(yīng)。

本文從矩形板的小撓度變形出發(fā),將膜力因子引入小撓度變形運動方程,獲得了大撓度情況下的變形運動方程。推導(dǎo)了在加載矩形脈沖情況下矩形板的運動微分方程。與文獻(xiàn)[7]相比,本文是直接使用在推導(dǎo)時所用的矩形脈沖載荷,而非再另外推導(dǎo)相應(yīng)的瞬動載荷進(jìn)行加載計算,而且所用的實驗樣本和有限元模型也具有更大的無因次位移和多種長寬比。

本文通過理論分析、有限元仿真和實驗對比,得到了以下結(jié)論:

1)在理論分析中所采用的帶有中央平臺區(qū)和移行鉸線,并且發(fā)生轉(zhuǎn)換的變形機(jī)構(gòu)是能夠通過有限元仿真手段進(jìn)行驗證的。

2)本文所假設(shè)的通過有限元仿真方法獲得板變形峰值對應(yīng)的響應(yīng)時間,進(jìn)而反推得到理論計算可用的應(yīng)變率效應(yīng)系數(shù)的方法是可行的。

3)在選取適當(dāng)?shù)膽?yīng)變率效應(yīng)系數(shù)前提下,膜力因子法能對無因次位移達(dá)到5以上的大撓度變形和矩形板長寬比在1~4范圍內(nèi)的實例進(jìn)行良好的中點最終位移預(yù)報。

[1] 朱錫,馮剛,張振華.爆炸載荷作用下固支方板的應(yīng)變場及破壞分析[J].船舶力學(xué),2005,9(2):83-89.

ZHU Xi,F(xiàn)ENG Gang,ZHANG Zhenhua.Strain field and damage analysis of clamped square plates subjected to explosive loading[J].Journal of Ship Mechanics,2005,9(2):83-89.

[2] CHUNG K Y S,NURICK G N.Experimental and numerical studies on the response of quadrangular stiffened plates.PartⅠ :subjected to uniform blast load[J].International Journal of Impact Engineering,2005,31(1):55-83.

[3]BONORCHISD,NURICK G N.The analysis and simulation ofwelded stiffener plates subjected to localised blast loading[J].International Journal of Impact Engineering,2010,37(3):260-273.

[4]PARK BW,CHO SR.Simple design formulae for predicting the residual damage of unstiffened and stiffened plates under explosion loadings[J].International Journalof Impact Engineering,2006,32(10):1721-1736.

[5]章媛,茍瑞君.矩形板在水下爆炸沖擊載荷作用下的動力響應(yīng)分析[J].四川兵工學(xué)報,2006,27(4):40-42.

[6]何建,肖玉鳳,陳振勇,等.空爆載荷作用下固支矩形鋼板的塑性極限變形[J].哈爾濱工業(yè)大學(xué)學(xué)報,2007,39(2):310-313,329.

HE Jian,XIAO Yufeng,CHEN Zhenyong,et al.Plastic limited deformation analysis of the clamped rectangular steel plate subjected to air non-contactexplosions[J].Journal of Harbin Institute of Technology,2007,39(2):310-313,329.

[7] YU Tongxi,CHEN Faliang.The large deflection dynamic plastic response of rectangular plates[J].International Journal of Impact Engineering,1992,12(4):605-616.

[8] SYMONDSP S,JONESN.Impulsive loading of fully clamped beams with finite plastic deflections and strain-rate sensitivity[J].International Journal of Mechanical Sciences,1972,14(1):49-69.

[9]CLOETE T J,NURICK G N,PALMER R N.The deformation and shear failure of peripherally clamped centrally supported blast loaded circular plates[J].International Journal of Impact Engineering,2005,31(1/4):92-117.

[10]BODNER SR,SYMONDSPS.Experimentson viscoplastic response of circular plates to impulsive loading[J].Journal of the Mechanics and Physics of Solids,1979,27(2):91-113.

[11] NURICK G N.A new techniquemeasure the deflection-time history of a structure subjected to high strain rates[J].International Journal of Impact Engineering,1985,3(1):17-26.

[12] JONESN.A theoretical study of the dynamic plastic behaviour of beams and plates with finite deflections[J].International Journal of Solids and Structures,1971,7(8):1007-1029.

The Large Deflection Dynam ic Plastic Response of Rectangu lar Plates Sub jected to Blast Load

YAN Feng,LIU Jingxi
Schoolof Naval Architecture and Ocean Engineering,Huazhong University of Science and Technology,Wuhan 430074,China

By introducing the membrane force factors,the motion equations of rectangular plates in the case of small deflection aremodified to fit in the case of large deflection,and the corresponding dynamic plastic response is analyzed.Meanwhile,the motion differential equations of fully-clamped rectangular plates subjected to blast loads are derived,with the loads being rectangular pulses.Therefore,the formulas ofmaximum permanent deformation of rectangular plates can be solved.In addition,this paper investigates the problem of yield stress increasing under high strain rate by employing an assumed selectionmethod for the strain rate effect coefficient.Finally,the proposed motion patterns are verified with a finite elementmethod,where the existing experimental samplesand the new finite elementmodels are both calculated using the described theoreticalmethod.The results are then compared with those obtained from previous experiments aswell as finite element simulations,and good agreement between the two can be clearly observed.

rectangular plate;dynamic plastic response;blast loading;rectangular pulse;membrane force factor;strain rate effect

U661.41

A

1673-3185(2013)01-47-07

10.3969/j.issn.1673-3185.2013.01.008

http://www.cnki.net/kcms/detail/42.1755.TJ.20130116.1428.007.htm l

2012-05-08 網(wǎng)絡(luò)出版時間:2013-01-16 14:28

顏 豐(1987-),男,碩士。研究方向:沖擊動力學(xué)。E-mail:545225243@qq.com

劉敬喜(1975-),男,博士,副教授。研究方向:船舶結(jié)構(gòu)振動與噪聲控制。E-mail:liu_jing_xi@m(xù)ail.hust.edu.cn

劉敬喜。

[責(zé)任編輯:盧圣芳]

猜你喜歡
有限元效應(yīng)變形
鈾對大型溞的急性毒性效應(yīng)
懶馬效應(yīng)
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
“我”的變形計
例談拼圖與整式變形
應(yīng)變效應(yīng)及其應(yīng)用
會變形的餅
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 91探花国产综合在线精品| 免费视频在线2021入口| 久久熟女AV| 啪啪免费视频一区二区| 国产一区自拍视频| 国产精品亚洲一区二区三区在线观看| 大陆精大陆国产国语精品1024| 欧美日本不卡| 免费无码一区二区| 亚洲精品桃花岛av在线| 欧美日韩精品一区二区在线线| 亚洲中文字幕av无码区| 国产91精品调教在线播放| 亚洲午夜天堂| 国产农村1级毛片| 亚洲欧美成aⅴ人在线观看| 91在线日韩在线播放| 在线观看av永久| 97视频免费在线观看| 中文字幕 日韩 欧美| 波多野结衣无码视频在线观看| 色综合a怡红院怡红院首页| 综合色在线| 国产精品亚洲日韩AⅤ在线观看| 久久特级毛片| 中文字幕色站| 欧美日韩中文字幕在线| 国产成人福利在线| 亚洲欧洲一区二区三区| 在线a视频免费观看| 91精品福利自产拍在线观看| 综合色区亚洲熟妇在线| 美女无遮挡被啪啪到高潮免费| 国产精品一老牛影视频| 精品国产免费观看一区| 亚洲一级色| www.youjizz.com久久| 国产成人一区免费观看| 国产成人调教在线视频| 亚洲va在线∨a天堂va欧美va| 亚洲色图欧美视频| 亚洲国产成人精品一二区| 久久精品国产999大香线焦| 久久大香伊蕉在人线观看热2| 午夜老司机永久免费看片| 欧美一级黄色影院| 久久久久久久蜜桃| 久久综合九色综合97婷婷| 久久亚洲天堂| 久久永久精品免费视频| 91福利片| 国产微拍一区二区三区四区| 九色在线观看视频| 91激情视频| 国产精品爆乳99久久| 99久久精品国产精品亚洲| 日韩亚洲综合在线| 国产成人亚洲无码淙合青草| 99在线观看国产| 精品国产网| 亚洲天堂免费观看| 国产好痛疼轻点好爽的视频| 亚洲综合久久成人AV| 欧美三级视频网站| 欧美国产在线看| 88国产经典欧美一区二区三区| 美女扒开下面流白浆在线试听| 久久国产精品电影| 亚洲免费人成影院| 9999在线视频| 欧美成人a∨视频免费观看 | 精品亚洲国产成人AV| 日韩小视频在线观看| 91在线无码精品秘九色APP| 四虎综合网| 无码高清专区| 亚洲第一色视频| 91国内在线观看| 日韩国产高清无码| 免费看久久精品99| 午夜福利在线观看成人| 欧美精品另类|