亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于各向同性本構(gòu)關(guān)系薄層單元的螺栓連接參數(shù)識別

        2014-05-25 00:34:03吳邵慶史勤豐費慶國
        振動與沖擊 2014年22期
        關(guān)鍵詞:本構(gòu)薄層螺栓

        姜 東,吳邵慶,史勤豐,費慶國

        (1.東南大學(xué)工程力學(xué)系,南京 210096;2.江蘇省工程力學(xué)分析重點實驗室,南京 210096)

        基于各向同性本構(gòu)關(guān)系薄層單元的螺栓連接參數(shù)識別

        姜 東1,2,吳邵慶1,2,史勤豐1,2,費慶國1,2

        (1.東南大學(xué)工程力學(xué)系,南京 210096;2.江蘇省工程力學(xué)分析重點實驗室,南京 210096)

        針對螺栓連接建模與參數(shù)識別問題開展研究?;诟飨蛲员緲?gòu)關(guān)系薄層單元理論,提出螺栓連接結(jié)構(gòu)接觸面力學(xué)性能識別方法。對單個螺栓搭接及多個螺栓搭接兩種結(jié)構(gòu)分別進行有限元建模,忽略螺栓質(zhì)量、螺孔影響,搭接界面采用各向同性本構(gòu)關(guān)系的薄層單元模擬。根據(jù)試驗?zāi)B(tài)參數(shù)構(gòu)造優(yōu)化問題,識別搭接界面薄層單元材料參數(shù)。結(jié)果表明,兩種螺栓連接結(jié)構(gòu)前四階彎曲模態(tài)頻率計算精度較高,薄層單元能準(zhǔn)確反映接觸界面力學(xué)性能。該方法適用于單個螺栓及螺栓較密集工況動力學(xué)精確模擬。

        螺栓連接;各向同性;薄層單元;參數(shù)識別

        機械結(jié)構(gòu)中的螺栓等連接件會嚴(yán)重影響結(jié)構(gòu)局部剛度及整體動力學(xué)性能。連接處動力學(xué)參數(shù)對建立準(zhǔn)確的動力學(xué)模型起決定性作用,是準(zhǔn)確計算結(jié)構(gòu)動態(tài)特性與響應(yīng)的前提。

        螺栓連接結(jié)構(gòu)有限元模擬方法主要分兩類,即非線性接觸算法與界面單元法[1-2]。前者需不斷通過接觸算法判別接觸狀態(tài),迭代計算結(jié)構(gòu)響應(yīng),計算量較大;而后者可在保證計算精度前提下降低計算量。彈簧單元、一般單元(Generic element)、零厚度單元(Zero thickness element)及薄層單元(Thin layer element)為較常用接觸界面力學(xué)性能參數(shù)化方法。Ren等[3]僅考慮接觸界面平動自由度,轉(zhuǎn)動自由度通過平動自由度實現(xiàn)。Ahmadian等[4]采用一般單元方法模擬螺栓連接搭接情況,通過試驗頻響函數(shù)識別單元參數(shù);Mayer等[5]利用零厚度單元及薄層單元提出能有效模擬接觸面方法。薄層單元理論最早來自模擬巖石接觸力學(xué)分析[6-7]。已有研究將薄層單元成功用于焊接、螺栓連接等形式的機械連接建模[8-9]。金峰等[10]用薄層單元分析重力壩體的抗震穩(wěn)定性;黃耀英等[11]分析基于橫觀各向同性體本構(gòu)方程簡化獲得薄層單元本構(gòu)方程。楊紅平等[12]基于分形幾何及接觸力學(xué)理論表征粗糙表面微凸體參數(shù),提出機械結(jié)合面法向接觸剛度計算模型。田紅亮等[13]采用虛擬各向同性材料模擬固定連接界面,基于赫茲接觸理論與分形理論推導(dǎo)等效虛擬材料參數(shù),有限元分析頻率與試驗值誤差小于9%。黃開放等[14]基于虛擬材料方法研究預(yù)緊力變化的連接結(jié)構(gòu)動力學(xué)仿真,結(jié)果與試驗值誤差在±6%以內(nèi)。

        本文針對螺栓連接建模與參數(shù)識別問題在線性范圍內(nèi)開展研究?;诒訂卧碚?,對單、多螺栓搭接兩種結(jié)構(gòu)進行有限元建模與參數(shù)識別。接觸界面采用基于各向同性本構(gòu)關(guān)系的薄層單元模擬。進而提出基于薄層單元的螺栓連接結(jié)構(gòu)接觸面力學(xué)性能模擬方法。

        1 薄層單元基本理論

        對于固定接觸界面的螺栓連接結(jié)構(gòu),結(jié)合面始終保持接觸,接觸狀態(tài)為粘附或者僅發(fā)生微小的局部相對位移(微觀滑移和微觀拍擊)。在該前提下,接觸面之間的作用力僅存在線性的法向接觸力與切向摩擦力。薄層單元能夠較準(zhǔn)確的反映接觸面力學(xué)性能。

        薄層單元最早由Desai等[6]提出并用于巖土結(jié)構(gòu)接觸分析,在相鄰接觸體間定義一層單元等效模擬連接界面接觸力學(xué)特征。考慮尺寸l1×l2×d薄層單元,據(jù)虛位移原理得虛功方程為

        圖1 薄層單元等參變換Fig.1 Isoparametric transformation of thin-layer element

        對薄層單元,設(shè)厚度d遠小于另兩方向特征尺寸l1,l2。研究表明[6],單元面內(nèi)應(yīng)變分量(εx,εy,γxy)、應(yīng)力分量(σx,σy,τxy)會被忽略。用單元形函數(shù)[11]分析時則有?Ni/?z遠大于?Ni/?x及?Ni/?y(Ni為單元形函數(shù)),可認為?Ni/?x=?Ni/?y≈0,從而得出應(yīng)變分量εx=εy=γxy≈0。因此薄層單元在高斯點的應(yīng)變分量有三個不為零,應(yīng)變分量可簡化為ε=[εzγyzγzx]T。若接觸面法向{e}n與兩切向{e}t分別定義為薄層單元局部坐標(biāo)系z,x,y方向,據(jù)以上分析,設(shè)連接界面法向、切向接觸性能相互獨立,兩切向接觸性能具有一致性,表征界面接觸性能的薄層單元本構(gòu)方程為

        此時法向彈性常數(shù)與切向彈性常數(shù)為非獨立的。

        2 薄層單元參數(shù)確定

        2.1 單元厚度

        薄層單元厚度選擇對計算結(jié)果影響較大。厚度過大則單元會有6個應(yīng)變分量,難以準(zhǔn)確體現(xiàn)接觸界面力學(xué)特征,不能用于模擬接觸情況;厚度過小則會導(dǎo)致雅可比矩陣行列式值趨向于零,致矩陣病態(tài)、求逆困難、無法計算位移-應(yīng)變關(guān)系。對薄層單元建模單元厚度選擇定義比例系數(shù)[1,6-7]為

        文獻[6]認為Ratio取值為10~100時能獲得較準(zhǔn)確結(jié)果。

        2.2 接觸剛度

        接觸剛度的確定為基于薄層單元螺栓連接結(jié)構(gòu)有限元模擬關(guān)鍵。接觸面保持線性粘合狀態(tài)時彈性矩陣中G,E均為常數(shù)[1],需試驗獲得。通過測試接觸界面應(yīng)力與位移間關(guān)系得到法向、切向接觸剛度[6]為

        式中:kn為法向接觸剛度;kτ為切向接觸剛度;d為薄層厚度;ur,vr為薄層法向及切向位移;σ,τ為薄層法向及切向應(yīng)力。

        該方法雖基于混凝土接觸試驗,但對機械連接結(jié)構(gòu)界面接觸問題亦具有一定參考意義。通過螺栓連接結(jié)構(gòu)切向受力試驗獲得接觸面等效剪切模量與切向剛度之關(guān)系[15]為

        式中:A為實際接觸面積;k為連接切向剛度,與接觸界面面壓及摩擦性能等有關(guān),其值為試驗所得連接處切向力與切向位移關(guān)系曲線斜率。

        以上兩種界面接觸剛度測試方法可獲得近似薄層單元彈性參數(shù),將材料本構(gòu)關(guān)系集成于薄層單元中能對連接結(jié)構(gòu)進行準(zhǔn)確有限元動力學(xué)分析。采用該方法對不同方向接觸剛度需分別進行試驗;對試驗裝置、試驗方法要求較高。本文采用動態(tài)試驗數(shù)據(jù)識別螺栓連接結(jié)構(gòu)結(jié)合面薄層單元參數(shù)。

        3 算例研究

        分別對單個螺栓及多個螺栓搭接結(jié)構(gòu)開展研究。結(jié)構(gòu)Ⅰ由一個螺栓與兩塊鋁板搭接;結(jié)構(gòu)Ⅱ由兩塊鋁板通過4個螺栓搭接。螺栓型號均為M10。結(jié)構(gòu)幾何尺寸及搭接見圖2、圖3,搭接板材料參數(shù)見表1。

        圖2 螺栓搭接結(jié)構(gòu)Ⅰ(單位:mm)Fig.2 Bolted joint structureⅠ

        圖3 螺栓搭接結(jié)構(gòu)Ⅱ(單位:mm)Fig.3 Bolted joint structureⅡ

        表1 搭接板材料參數(shù)Tab.1 Material properties of the lapp ing p late

        3.1 模態(tài)試驗結(jié)果

        對兩種螺栓搭接模型分別進行模態(tài)試驗。螺栓預(yù)緊力對連接結(jié)構(gòu)動態(tài)特性有較大影響,預(yù)緊力大到一定程度時結(jié)構(gòu)固有頻率不會隨預(yù)緊力的增加而變化。用錘擊法以螺栓搭接結(jié)構(gòu)Ⅱ為例進行模態(tài)試驗說明,測試與測點布置見圖4、圖5。采用彈簧繩懸掛方法模擬自由-自由邊界條件,傳感器布置于結(jié)構(gòu)端部以避開振型節(jié)點;為獲得較光滑振型曲線共布置13個測點。螺栓搭接結(jié)構(gòu)Ⅰ的試驗條件與此類似。通過試驗得螺栓搭接結(jié)構(gòu)Ⅰ,Ⅱ前四階試驗?zāi)B(tài)數(shù)據(jù)見圖6。

        圖4 螺栓搭接結(jié)構(gòu)Ⅱ測試情況Fig.4 The experiment condition of bolted structureⅡ

        圖5 測點布置Fig.5 Test point arrangement

        圖6 螺栓搭接結(jié)構(gòu)Ⅰ,Ⅱ前四階試驗?zāi)B(tài)數(shù)據(jù)Fig.6 The first four experimentalmode shapes and modal frequencies of bolted structureⅠandⅡ

        3.2 薄層單元參數(shù)識別

        本文對模型Ⅰ,Ⅱ建立有限元模型,忽略螺栓質(zhì)量及螺孔影響,搭接板用實體單元模擬,界面用各向同性薄層單元模擬。將薄層單元材料參數(shù)識別轉(zhuǎn)化為優(yōu)化問題。待識別參數(shù)為接觸面薄層單元材料彈性模量E與剪切模量G。建立目標(biāo)函數(shù)為前四階彎曲振型計算、試驗結(jié)果殘差加權(quán)平方和最小,即

        式中:J(p)為目標(biāo)函數(shù)定義在結(jié)構(gòu)待識別參數(shù)合理取值范圍p1≤p≤p2內(nèi),試驗與計算模態(tài)參數(shù)殘差加權(quán)平方和取極小值;ε為模態(tài)參數(shù)殘差;zm,za(p)分別為試驗與計算的模態(tài)參數(shù);W為反映各模態(tài)參數(shù)殘差相對權(quán)重的對角陣。

        設(shè)定待修正參數(shù)初值,用靈敏度分析法迭代求解式(14)。第j個迭代步問題可描述為

        式中:Sj=W1/2?zj/?pj為模態(tài)參數(shù)對待修正參數(shù)的加權(quán)靈敏度矩陣。迭代直到待識別參數(shù)p收斂。計算模態(tài)參數(shù)精度滿足要求時可得準(zhǔn)確的接觸面材料參數(shù)。

        3.2.1 螺栓搭接結(jié)構(gòu)Ⅰ

        螺栓搭接結(jié)構(gòu)Ⅰ結(jié)合面薄層單元示意圖見圖7。由于靠近螺栓部位接觸剛度高于遠離螺栓部位,將接觸面薄層單元等分為三部分,采用兩種不同的各向同性本構(gòu)關(guān)系模擬及參數(shù)識別??拷菟ǎ▓D中深色)區(qū)域待識別參數(shù)為彈性模量E1與剪切模量G1;遠離螺栓(圖中淺色)區(qū)域待識別參數(shù)為彈性模量E2及剪切模量G2。薄層比例系數(shù)Ratio=10,材料密度為0。

        圖7 螺栓搭接結(jié)構(gòu)Ⅰ結(jié)合面薄層單元示意圖Fig.7 Thin-layer element in the contact surface of bolted structureⅠ

        圖8為螺栓搭接結(jié)構(gòu)Ⅰ結(jié)合面薄層單元材料參數(shù)迭代收斂曲線,表2為薄層單元材料參數(shù)初始值與識別值,表3為參數(shù)識別前后計算與試驗?zāi)B(tài)參數(shù)誤差比較。由識別結(jié)果看出,固有頻率最大誤差不超過1.0%,識別精度較高。表2中楊氏模量E大于剪切模量G,能反映接觸界面法向剛度大于切向接觸剛度特征;靠近螺栓區(qū)域結(jié)合面材料參數(shù)識別結(jié)果均遠大于遠離螺栓區(qū)域,能較好反映實際情況。因此基于各向同性材料的薄層單元能準(zhǔn)確反映結(jié)合面切向、法向接觸性能,計算模態(tài)參數(shù)精度高。

        表2 結(jié)合面薄層單元材料參數(shù)初始值和識別值Tab.2 The initial and the identified material parameters of the thin-layer element

        表3 參數(shù)識別前后計算與試驗?zāi)B(tài)參數(shù)誤差比較Tab.3 Com parison betw een the com putational and experimental modal data before and after parameter identification

        3.2.2 螺栓搭接結(jié)構(gòu)Ⅱ

        螺栓搭接結(jié)構(gòu)Ⅱ由兩塊鋁板通過4個螺栓搭接而成。接觸界面螺栓分布更密集,與結(jié)構(gòu)Ⅰ相比結(jié)構(gòu)Ⅱ的接觸情況更復(fù)雜。由于螺栓預(yù)緊力產(chǎn)生的接觸面面壓分布非均勻,且4個螺栓距離較近,導(dǎo)致接觸剛度分布不均勻。結(jié)構(gòu)Ⅱ在接觸界面定義兩種材料,在約兩倍螺桿直徑的矩形區(qū)域內(nèi)定義各向同性材料1,待識別參數(shù)為彈性模量E1及剪切模量G1,其它區(qū)域定義各向同性材料2,識別參數(shù)為彈性模量E2及剪切模量G2。參數(shù)設(shè)置方法與結(jié)構(gòu)Ⅰ相同。

        圖9為螺栓搭接結(jié)構(gòu)Ⅱ結(jié)合面等效材料參數(shù)迭代收斂曲線,表4為結(jié)合面薄層單元材料參數(shù)初始值與識別值,表5為參數(shù)識別前后計算與試驗?zāi)B(tài)參數(shù)誤差比較。由識別結(jié)果看出,計算模態(tài)頻率最大誤差不超過2.5%,識別精度較高;薄層單元材料彈性參數(shù)在螺栓附近明顯高于其它區(qū)域,且法向彈性模量大于剪切模量,與實際搭接結(jié)構(gòu)接觸面接觸性能一致。因此基于各向同性材料的薄層單元能較好模擬多個螺栓搭接結(jié)構(gòu)的連接性能。

        圖8 螺栓搭接結(jié)構(gòu)Ⅰ結(jié)合面薄層單元材料參數(shù)迭代收斂曲線Fig.8 Material parameter convergence of the thin-layer element in the contact surface of bolted structureⅠ

        圖9 螺栓搭接結(jié)構(gòu)Ⅱ結(jié)合面薄層單元材料參數(shù)迭代收斂曲線Fig.9 Parameter convergence of the thin-layer element in the contact surface of bolted structureⅡ

        表4 結(jié)合面薄層單元材料參數(shù)初始值和識別值Tab.4 The initial and the identified material parameters of the thin-layer element

        表5 參數(shù)識別前后計算與試驗?zāi)B(tài)參數(shù)誤差比較Tab.5 Com parison betw een the com putational and experimental modal data before and after param eter identification

        4 結(jié) 論

        (1)本文在各向同性本構(gòu)關(guān)系薄層單元理論基礎(chǔ)上提出螺栓連接結(jié)構(gòu)接觸面力學(xué)性能識別方法。通過固定結(jié)合面單、多個螺栓搭接兩種結(jié)構(gòu)對方法進行驗證。

        (2)采用線性各向同性本構(gòu)關(guān)系薄層單元模擬接觸面,所提方法建模與參數(shù)識別效率較高,識別后參數(shù)能準(zhǔn)確反映螺栓搭接結(jié)構(gòu)結(jié)合面接觸性能。近螺栓區(qū)域接觸面材料參數(shù)識別結(jié)果大于遠離螺栓區(qū)域,能較好反映接觸界面剛度分布不均勻性。

        (3)該方法能準(zhǔn)確描述螺栓搭接結(jié)構(gòu)接觸面力學(xué)性能,適用于單、多個螺栓搭接等多種工況有限元精確模擬。

        [1]Bograd S,Reuss P,Schmidt A,et al.Modeling the dynamics of mechanical joints[J].Mechanical Systems and Signal Processing,2011,25:2801-2826.

        [2]Mackerle J.Finite element analysis of fastening and joining:a bibliography(1990-2002)[J].International Journal of Pressure Vessels and Piping,2003,80(4):253-271.

        [3]Ren Y,Beards C F.Identification of effective linear joints using coupling and joint identification techniques[J].Journal of Vibration and Acoustics,1998,120(2):331-338.

        [4]Ahmadian H,Jalali H.Generic element formulation for modelling bolted lap joints[J].Mechanical Systems and Signal Processing,2007,21:2318-2334.

        [5]Mayer M,Gaul L.Segment-to-segment contact elements for modeling joint interfaces in finite element analysis[J].Mechanical Systems and Signal Processing,2007,21:724-734.

        [6]Desai C S,Zaman MM,Lightner J G,et al.Thin-layer element for interfaces and joints[J].International Journal for Numerical and Analytical Methods in Geomechanics,1984,8(1):19-43.

        [7]Sharma K,Desai C.Analysis and implementation of thinlayer element for interfaces and joints[J].Journal of Engineering Mechanics,1992,118(12):2442-2462.

        [8]Ahmadian H,Jalali H.Identification of bolted lap joints parameters in assembled structures[J].Mechanical Systems and Signal Processing,2007,21:1041-1050.

        [9]Jalali H,Hedayati A,Ahmadian H.Modeling mechanical interfaces experiencingmicro-slip/slap[J].Inverse Problems in Science and Engineering,2011,19(6):751-764.

        [10]金峰,邵偉,張立翔,等.模擬軟弱夾層動力特性的薄層單元及其工程應(yīng)用[J].工程力學(xué),2002,19(2):36-40.

        JIN Feng,SHAOWei,ZHANG Li-xiang,et al.A thin-layer element for simulation of static and dynamic characteristics of soft interlayer and its application[J].Engineering Mechanics,2002,19(2):36-40.

        [11]黃耀英,吳中如,王德信.薄層單元基本假設(shè)和簡化探討[J].力學(xué)與實踐,2008,130(2):49-52.

        HUANG Yao-ying,WU Zhong-ru,WANG De-xin.Discuss on fundamental assumption and simplification of thin-layer element[J].Mechanics in Engineering,2008,130(2):49-52.

        [12]楊紅平,傅衛(wèi)平,王雯,等.基于分形幾何與接觸力學(xué)理論的結(jié)合面法向接觸剛度計算模型[J].機械工程學(xué)報,2013,49(1):102-107.

        YANG Hong-ping,F(xiàn)U Wei-ping,WANG Wen,et al.Calculation model of the normal contact stiffness of joints based on the fractal geometry and contact theory[J].Journal of Mechanical Engineering,2013,49(1):102-107.

        [13]Tian Hong-liang,Li Bin,Liu Hong-qi,et al.A new method of virtual material hypothesis-based dynamic modeling on fixed joint interface in machine tools[J].International Journal of Machine Tools&Manufacture,2011,51:239-249.

        [14]黃開放,金建新.基于虛擬材料方法的螺栓預(yù)緊力模擬的研究[J].機械設(shè)計與制造,2012(8):148-150.

        HUANG Kai-fang,JIN Jian-xin.Research on bolt preload simulation based on virtual material method[J].Machine Design&Manufacture,2012(8):148-150.

        [15]Schmidt A,Bograd S,Gaul L.Measurement of join patch properties and their integration into finite-element calculations of assembled structures[J].Shock and Vibration,2012,19(5):1125-1133.

        Parameter identification of bolted-joint based on themodelwith thin-layer elements with isotropic constitutive relationship

        JIANG Dong1,2,WU Shao-qing1,2,SHIQin-feng1,2,F(xiàn)EIQing-guo1,2
        (1.Department of Engineering Mechanics,Southeast University,Nanjing 210096,China;2.Jiangsu Key Laboratory of Engineering Mechanics,Nanjing 210096,China)

        The finite elementmodeling and parameter identification of bolted jointswere concerned.According to the basic theory of thin-layer element with isotropic constitutive relationship,a parameter identification method for recognizing the mechanical characteristics of contact surface in bolted structures was proposed.The finite element modelingmethods for single andmultiple bolted structureswere investigated respectively,ignoring the hole and themass of the bolt,and modeling the contact surface by using thin-layer elements.Experimentalmodal datawere used for identifying the constitutive parameters of thin-layer elementwith an optimization procedure.Applying the identified parameters in the finite element model,the maximum error between the computational and experimental modal frequencies was reduced reasonably.It is shown that thin layer elementwith identified parameters can be used for accurately simulating the normal and tangential stiffness of contact surface.The proposed approach is available for precisely simulating the single and multiple bolted structures.

        bolted-joint connection;isotropic;thin-layer element;parameter identification

        TU318

        :A

        10.13465/j.cnki.jvs.2014.22.007

        國家自然科學(xué)基金(10902024);教育部新世紀(jì)優(yōu)秀人才支持計劃(NCET-11-0086);江蘇省自然科學(xué)基金(BK2010397);江蘇高校優(yōu)勢學(xué)科建設(shè)工程資助項目(1105007001)

        2013-07-08 修改稿收到日期:2013-08-16

        姜東男,博士生,1985年12月生

        費慶國男,教授,博士生導(dǎo)師,1977年生

        猜你喜歡
        本構(gòu)薄層螺栓
        M16吊耳螺栓斷裂失效分析
        預(yù)緊力衰減對摩擦型高強螺栓群承載力的影響
        四川建筑(2020年1期)2020-07-21 07:26:08
        離心SC柱混凝土本構(gòu)模型比較研究
        螺栓緊固雜談
        鋸齒形結(jié)構(gòu)面剪切流變及非線性本構(gòu)模型分析
        維藥芹菜根的薄層鑒別
        一種新型超固結(jié)土三維本構(gòu)模型
        SiN_x:H膜沉積壓強與擴散薄層電阻的匹配性研究
        參芪苓口服液的薄層色譜鑒別
        芪參清幽膠囊的薄層鑒別研究
        在线精品国产亚洲av蜜桃| 在线av野外国语对白| 蜜桃av区一区二区三| 自拍成人免费在线视频| 欧美变态另类刺激| 无码人妻精品一区二区三区66| 欧美在线观看www| 九九精品国产亚洲av日韩| 亚洲精品美女久久777777| 无码熟熟妇丰满人妻啪啪| 国产盗摄XXXX视频XXXX| 男女视频网站在线观看| 无码国内精品久久人妻| 国产内射合集颜射| 人妻丰满熟妇av一区二区| 亚洲av少妇高潮喷水在线| 国产av无码专区亚洲av毛网站 | 高潮抽搐潮喷毛片在线播放| 日产精品久久久久久久蜜臀| 久久99久久99精品观看| 中文av字幕一区二区三区| 亚洲日韩av一区二区三区中文 | 亚洲精品无码专区| 日产精品久久久久久久| 高清av一区二区三区在线| 亚洲精品国产精品乱码视色| 粉嫩虎白女毛片人体| 成人无码区免费AⅤ片WWW| 在线不卡av一区二区| 亚洲成aⅴ人片久青草影院| 国产成人精品日本亚洲18| 日本一区二区三区四区在线看| 美女用丝袜脚玩我下面| 中国农村妇女hdxxxx| 巨臀中文字幕一区二区| 男男做h嗯啊高潮涩涩| 99久久精品国产一区二区| 中文亚洲日韩欧美| 精品蜜臀国产av一区二区| 国产亚洲精品美女久久久m| 欧美黑人乱大交|