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

非線性系統(tǒng)帶集中質(zhì)量懸臂梁易損件跌落沖擊特性

2016-01-15 03:35:27郝蒙,陳安軍
振動與沖擊 2015年15期
關(guān)鍵詞:有限元

非線性系統(tǒng)帶集中質(zhì)量懸臂梁易損件跌落沖擊特性

郝蒙1,陳安軍1,2

(1.江南大學(xué)包裝工程系,江蘇無錫2141222; 2.國家輕工業(yè)包裝制品質(zhì)量監(jiān)督檢測中心,江蘇無錫214122)

摘要:以三次非線性緩沖包裝系統(tǒng)為研究對象,建立系統(tǒng)跌落沖擊動力學(xué)方程;綜合龍格-庫塔和有限元法設(shè)計動力學(xué)方程的求解算法;探討了集中質(zhì)量和主體振動頻率對易損件響應(yīng)的影響。數(shù)值分析表明:系統(tǒng)跌落沖擊過程,易損件最大位移和加速度響應(yīng)均位于懸臂梁自由端,易損件與主體連接部內(nèi)應(yīng)力最大;易損件相對主體質(zhì)量較小時,系統(tǒng)耦合作用對易損件響應(yīng)影響不明顯;隨集中質(zhì)量增加或主體振動頻率向易損件第一階固有頻率接近時,易損件內(nèi)應(yīng)力響應(yīng)、自由端相對位移響應(yīng)幅值顯著增大。在含彈性易損件產(chǎn)品系統(tǒng)包裝設(shè)計中,連接部的內(nèi)應(yīng)力和易損件相對位移是關(guān)注的重要參數(shù)。研究結(jié)論可為帶集中質(zhì)量懸臂梁產(chǎn)品緩沖包裝設(shè)計提供理論依據(jù)。

關(guān)鍵詞:集中質(zhì)量;懸臂梁;易損件;有限元;三次非線性;跌落沖擊

中圖分類號:TB485.3;O322文獻標志碼:A

基金項目:國家自然科學(xué)基金青年基金資助項目(11002084);上海市教委創(chuàng)新基金資助項目(12YZ092, 12YZ074) ;滬江基金(D14005);上海理工大學(xué)研究生創(chuàng)新基金項目

收稿日期:2014-03-26修改稿收到日期:2014-07-23

Dropping shock characteristics of a cubic nonlinear system with a cantilever beam type elastic critical component attaching a concentrated tip mass

HAOMeng1,CHENAn-jun1,2(1. Department of Packaging Engineering, Jiangnan University, Wuxi 214122, China; 2. China National Light Industry Packaging Quality Supervision & Inspection Center, Wuxi 214122, China)

Abstract:In order to investigate the dropping shock characteristics of a cubic nonlinear system with a cantilever beam critical component attaching a concentrated tip mass, the dynamic model of the whole system was built. Runge-Kutta method and the finite element method were used to solve the system’s dynamic equations, and the effects of the system parameters, such as, concentrated tip mass and natural frequencies of the main part on the responses of the critical component were studied. The case study showed that the maximum dropping displacement and acceleration responses of the critical component occur at its free end, while the maximum internal stress appears at the joint point between the critical component and the main part; the effects of the coupling action of the system on the responses of the critical component can be ignored when its mass is much less than that of the main component; with increase in concentrated tip mass and/or the natural frequency of the main part close to the 1st natural frequencu of the critical component, the free end of the critical component’s response amplitudes increase obviously; in the packaging design of systems with elastic critical components, the internal stress at the joint point between the critical component and the main part and the relative displacement of the critical component are be parameters to evaluate if products are damaged. The results provided a theoretical foundation for the package cushioning design of this type systems.

Key words:concentrated tip mass; cantilever beam; critical component; finite element; cubic nonlinear; dropping impact

研究包裝系統(tǒng)的動力學(xué)模型及跌落條件下系統(tǒng)的動力學(xué)特性是產(chǎn)品緩沖包裝設(shè)計的理論基礎(chǔ)。Mindlin[1]首次提出易損件的概念,建立兩種包裝系統(tǒng)的動力學(xué)模型(單自由度和二自由度集中質(zhì)量-彈簧系統(tǒng)模型),研究線性和非線性系統(tǒng)跌落沖擊動力學(xué)特性。此后,考慮易損件的二自由度集中質(zhì)量-彈簧模型的動力學(xué)特性研究得到不斷豐富和發(fā)展[2-5]。孫勇[6]、王軍[7]和高德[8]分別對多自由度包裝系統(tǒng)振動和沖擊動力學(xué)特性進行了分析,Suhir[9]和Wong[10-12]針對帶彈性板易損件的產(chǎn)品包裝系統(tǒng),建立動力學(xué)模型,分析彈性板的跌落沖擊特性,結(jié)果表明對于帶彈性易損件的產(chǎn)品,由于易損件為無限自由度,傳統(tǒng)集中質(zhì)量-彈簧模型應(yīng)用有失準確性;奚德昌等[13-14]針對具有簡支梁彈性部件工藝品包裝系統(tǒng),建立剛體-彈性體模型,探究了緩沖材料參數(shù)對易損件和系統(tǒng)響應(yīng)的影響;高德等[15-17]分別以桿式、簡支梁和懸臂梁3種不同類型易損件產(chǎn)品與發(fā)泡聚乙烯緩沖材料組成的耦合系統(tǒng)為研究對象,建立了系統(tǒng)的動力學(xué)方程,利用顯式有限差分法研究了跌落沖擊過程中易損件位移響應(yīng)、加速度響應(yīng)和應(yīng)力分布的特點。對包含集中質(zhì)量懸臂梁彈性易損件產(chǎn)品包裝系統(tǒng)的研究現(xiàn)未見報道,且此類易損件在機電產(chǎn)品中較為常見,其包裝防護理論分析有待進一步探討。

本文以三次非線性包裝系統(tǒng)為研究對象,產(chǎn)品易損件為自由端有集中質(zhì)量的懸臂梁,建立系統(tǒng)跌落沖擊下的動力學(xué)模型;龍格-庫塔法與有限元法結(jié)合,設(shè)計系統(tǒng)動力學(xué)方程的數(shù)值分析算法。探討易損件位移、加速度響應(yīng)以及內(nèi)應(yīng)力變化規(guī)律,并分析其影響因素。研究方法和結(jié)論可為此類產(chǎn)品緩沖包裝設(shè)計提供理論依據(jù)。

1系統(tǒng)動力學(xué)模型

圖1 跌落條件下考慮易 損件包裝系統(tǒng)力學(xué)模型 Fig.1 The dropping model of packaging system with critical component

圖2 帶集中質(zhì)量懸臂梁易損件有限單元劃分 Fig.2 The finite element division of cantilever beam component with concentrated tip mass

(1a)

(1b)

式中:E、I、ρ和A分別為易損件彈性模量、橫截面慣性矩、材料密度和截面面積。

易損件動能T和應(yīng)變能U為

(2)

以u、w分別表示易損件節(jié)點絕對位移和相對位移向量,即

單元節(jié)點位移向量與總節(jié)點位移向量關(guān)系為

ui=Ciu,wi=Ciw

(3)

式中:Ci為4×(2n+2)的矩陣,記為第i個單元位移提取矩陣,具體為

將式(3)和Ci代入式(2),易損件動能和應(yīng)變能分別為

(4)

式中:M和K均為2(n+1)×2(n+1)的矩陣,分別表示易損件質(zhì)量矩陣和剛度矩陣,即

(5)

考慮自由端的集中質(zhì)量影響,將m2添加到總質(zhì)量矩陣的第(2n+1)×(2n+1)元素上。

由拉格朗日方程建立易損件的動力學(xué)方程

(6)

式中:w=u-y1(t)v,v=[1010…1010]T。

系統(tǒng)跌落過程中,可將易損件在主體位移激勵下的響應(yīng)分析轉(zhuǎn)化為求解易損件在慣性力作用下的相對位移響應(yīng)問題。將u=w+y1(t)v代入方程(6)得

(7)

考慮易損件對主體的影響,跌落沖擊條件下主體的動力學(xué)方程為

(8)

若不考慮易損件對主體的影響,初始條件不變,主體動力學(xué)方程為

(9)

同時,由變分迭代法[19]可得的主體頻率近似表達式

(10)

2動力學(xué)方程求解

對二階常微分方程,龍格-庫塔數(shù)值分析方法具有很好的精度[20]。方程(9)可直接采用龍格-庫塔法求解,由于方程(8)包含易損件對主體的影響,需有限差分離散后再使用龍格-庫塔法求解。有限差分離散方法如下。

(11)

(12)

式中:B=EIτ2/(ρAl4)。進一步

(13)

(14)

應(yīng)用Newmark法[23]求解;對應(yīng)廣義特征值方程可表示為

K*φ=ω2M*φ

(15)

采用子空間迭代法[23]求解。

(16)

式中:J(n×2n)為橫向相對位移提取矩陣,即

易損件橫向絕對位移向量可表示為

u2=w2+y1(:,t)

(17)

易損件內(nèi)應(yīng)力為

(18)

式中:r為易損件橫截面半徑。

3數(shù)值計算結(jié)果及其分析

某小型機電產(chǎn)品,主體質(zhì)量m1=10 kg;懸臂梁(45號鋼)橫截面半徑r=0.003 m,長L=0.1 m,材料密度ρ=7 850 kg/m3,彈性模量E=200 GPa,彈性極限σe=180 MPa,材料安全系數(shù)取1.2,許用內(nèi)應(yīng)力為150 MPa;集中質(zhì)量m2=0.02 kg。緩沖材料初始彈性常數(shù)k0=100 N/cm,非線性常數(shù)e=72 N/cm3,系統(tǒng)跌落高度H=0.6 m。

圖3 易損件跌落沖擊響應(yīng) Fig.3 The dropping shock response of critical component

表1 易損件跌落沖擊響應(yīng)最值

考慮易損件與主體耦合作用,易損件位移、加速度響應(yīng)和懸臂梁內(nèi)應(yīng)力分布如圖3所示,相對位移、絕對加速度及內(nèi)應(yīng)力最值如表1所示。易損件絕對位移響應(yīng)和絕對加速度響應(yīng)沿懸臂梁長度方向變化不明顯;相對位移響應(yīng)和絕對加速度響應(yīng)最值出現(xiàn)在自由端,分別為0.41 mm和-560.02 m/s2;內(nèi)應(yīng)力最大值81.60 MPa出現(xiàn)在與主體連接部。傳統(tǒng)脆值理論以加速度是否超過許用值作為判斷產(chǎn)品失效標準,而對于含彈性易損件產(chǎn)品,根據(jù)彈性材料失效準則,當最大應(yīng)力超過其彈性極限,材料將發(fā)生塑性變形,產(chǎn)品失效,同時考慮到產(chǎn)品內(nèi)部構(gòu)件間隙,避免碰撞,應(yīng)綜合考慮最大應(yīng)力是否超過許用值和易損件相對位移是否超過許用值作為判定該產(chǎn)品失效的準則。

圖4 不同集中質(zhì)量易損件跌落沖擊響應(yīng) Fig.4 The dropping shock response of critical component with different concentrated tip mass

考慮易損件和主體耦合作用,集中質(zhì)量取不同值(0.01、0.02、0.03和0.04 kg)時易損件自由端相對位移、相對加速度和根部內(nèi)應(yīng)力變化如圖4。隨著集中質(zhì)量增加,易損件響應(yīng)顯著增加;集中質(zhì)量為0.04 kg時,內(nèi)應(yīng)力及相對位移最大值分別為138.44 MPa和0.73 mm,分別是集中質(zhì)量為0.01 kg時的2.53倍和2.7倍。最大值出現(xiàn)時間和響應(yīng)周期基本不變。

考慮耦合作用,易損件響應(yīng)最值及出現(xiàn)時間如表2所示。集中質(zhì)量為0.02 kg時,易損件總質(zhì)量(0.042 3 kg)遠小于主體質(zhì)量(10 kg),耦合作用對易損件的響應(yīng)影響不明顯;當集中質(zhì)量為0.08 kg時,相對位移及絕對加速度最大值誤差達到5%以上,耦合作用影響不能忽略。

表2 考慮和不考慮耦合作用易損件跌落沖擊響應(yīng)最值

表3 主體振動頻率

圖5 不同主體頻率易損件跌落沖擊響應(yīng) Fig.5 The dropping shock response of critical component with different main component frequency

初始條件和產(chǎn)品主體質(zhì)量不變,緩沖材料參數(shù)e取不同的值,相應(yīng)主體的響應(yīng)頻率通過式(10)求得如表3所示;集中質(zhì)量m2=0.02 kg時,由式(14)模態(tài)分析可求得懸臂梁易損件第一階頻率為195.63 Hz。圖5為主體運動頻率取不同值時,易損件自由端相對位移、相對加速度和根部內(nèi)應(yīng)力變化。當主體振動頻率為26.91時,內(nèi)應(yīng)力最大值為185.46 MPa已超過許用應(yīng)力150 MPa,緩沖材料需重新設(shè)計;且自由端相對位移最大值為0.94 mm,為主體振動頻率14.32 Hz時的2.29倍。主體振動頻率向易損件一階頻率靠近時,易損件自由端相對位移、相對加速度響應(yīng)以及根部內(nèi)應(yīng)力幅值均顯著增大,易損件響應(yīng)對主體振動頻率變化敏感。

4結(jié)論

以易損件為帶集中質(zhì)量的懸臂梁產(chǎn)品包裝系統(tǒng)為研究對象,建立了多自由度剛體-彈性體動力學(xué)模型和運動方程;結(jié)合龍格-庫塔法和有限元法提出解決此類問題的有效數(shù)值算法,主要結(jié)論如下:

(1)帶集中質(zhì)量懸臂梁易損件的包裝系統(tǒng),易損件自由端位移和加速度響應(yīng)最大,內(nèi)應(yīng)力最小;而與主體連接部內(nèi)應(yīng)力最大。

(2)集中質(zhì)量對易損件位移、加速度響應(yīng)及內(nèi)應(yīng)力峰值影響顯著。當集中質(zhì)量為0.04 kg時,內(nèi)應(yīng)力及相對位移最大值分別為138.44 MPa和0.73 mm,是集中質(zhì)量為0.01 kg時的2.53倍和2.7倍。

(3)考慮動力學(xué)方程耦合作用具有實際意義。集中質(zhì)量為0.02 kg時,耦合作用對易損件的響應(yīng)影響不明顯;當集中質(zhì)量為0.08 kg時,相對位移及絕對加速度最大值誤差達到5%以上,耦合作用影響不能忽略。

(4)易損件響應(yīng)峰值與周期對主體振動頻率敏感。當主體頻率向易損件的第一階固有頻率靠近時,易損件相對位移和內(nèi)應(yīng)力峰值增加尤為顯著,峰值出現(xiàn)時間及響應(yīng)周期同樣變化明顯。

對含彈性易損件產(chǎn)品包裝設(shè)計中,需根據(jù)彈性材料失效準則,將彈性易損件的最大內(nèi)應(yīng)力是否超過許用值,作為判斷該類產(chǎn)品是否失效的標準之一,同時考慮相對位移是否超過允許值,作為產(chǎn)品是否失效的判定標準,使得產(chǎn)品包裝系統(tǒng)設(shè)計及評價更具有針對性。緩沖包裝設(shè)計中,緩沖材料的選擇應(yīng)使主體的振動頻率遠離易損件的固有頻率。

參考文獻

[1]Mindlin R D. Dynamics of packaging cushioning[J]. Bell System Technical Journal, 1945, 25:353-445.

[2]王志偉, 胡長鷹. 洗衣機包裝系統(tǒng)在矩形波沖擊下的特性研究[J]. 包裝工程, 1997, 18(2,3):14-16.

WANG Zhi-wei, HU Chang-ying. The Shock characteristics of washing machine under the action of rectangular pulse[J]. Packaging Engineering, 1997, 18(2,3):14-16.

[3]Wang Zhi-wei, Hu Chang-ying. Shock spectra and damage boundary curves for nonlinear package cushioning system[J]. Packaging Technology and Science, 1999, 12(5): 207-217.

[4]王軍,王志偉.半正弦脈沖激勵下考慮易損件的正切型包裝系統(tǒng)沖擊特性研究[J].振動與沖擊,2008,27(1):167-168.

WANG Jun, WANG Zhi-wei.3-Dimensional shock response spectra characterizing shock response of a tangent packaging system with critical component[J]. Journal of Vibration and Shock, 2008, 27(1):167-168.

[5]段寧寧,陳安軍.矩形脈沖激勵下斜支承系統(tǒng)易損件的沖擊特性研究[J].包裝工程,2012,34(7):21-24.

DUAN Ning-ning, CHEN An-jun. Shock characteristics of vulnerable components of tilted support system under rectangular pulse excitation[J]. Packaging Engineering, 2012, 34(7):21-24.

[6]孫勇, 熊彪, 王石. 多自由度緩沖包裝系統(tǒng)的平均能量設(shè)計法[J]. 包裝工程, 1997, 18(1):17-20.

SUN Yong, XIONG Biao, WANG Shi. The averaging energy design method used for a multi-freedom-degree cushioning packaging system[J]. Packaging Engineering, 1997, 18(1):17-20.

[7]王軍, 王志偉. 多層堆碼包裝系統(tǒng)沖擊動力學(xué)特性研究(I): 沖擊譜[J]. 振動與沖擊, 2008, 27(8):106-107.

WANG Jun, WANG Zhi-wei. Combined shock spectrum of linear stacking packaging system[J]. Journal of Vibration and Shock, 2008, 27(8):106-107.

[8]高德, 盧富德. 聚乙烯緩沖材料多自由度跌落包裝系統(tǒng)優(yōu)化設(shè)計[J]. 振動與沖擊, 2012, 31(3):69-72.

GAO De, LU Fu-de.Optimization design of MDOF package cushioning system made of polyethylene[J]. Journal of Vibration and Shock, 2012, 31(3):69-72.

[9]Suhir E, Burke R. Dynamic response of a rectangular plate to a shock load, with application to portable electronic products[J]. Components, Packaging and Manufacturing Technology, Part B: Advanced Packaging, IEEE Transactions on, 1994,17(3):449-460.

[10]Wong E H. Dynamics of board-level drop impact[J]. Journal of Electronic Packaging, 2005, 127(3): 200-207.

[11]Wong E H, Mai Y W, Seah S K. Board level drop impact-fundamental and parametric analysis[J]. Journal of Electronic Packaging, 2005, 127(4):496-502.

[12]Wong E H, Mai Y W. New insights into board level drop impact[J]. Microelectronics Reliability, 2006, 46(5):930-938.

[13]奚德昌, 彭南陵. 具有簡支梁彈性部件包裝物品的防震[J]. 包裝工程, 1989, 10(1):42-48.

XI De-chang, PENG Nan-ling. Shock proof of the packaging for elastic parts with simple beam[J]. Packaging Engineering, 1989, 10(1):42-48.

[14]彭南陵, 奚德昌. 具有簡支梁彈性部件包裝物品的緩沖[J]. 南昌大學(xué)學(xué)報(工程技術(shù)版), 1994, 16(1):43-48.

PENG Na-ling, XI De-chang. Cushioning of packaged articles with an elastic element of simply supported beam[J]. Journal of Nanchang University, 1994, 16(1):43-48.

[15]盧富德,陶偉明,高德.具有簡支梁式易損部件的產(chǎn)品包裝系統(tǒng)跌落沖擊研究[J].振動與沖擊,2012,31(5):79-81.

LU Fu-de, TAO Wei-ming, GAO De.Drop impact analysis on item packaging system with beam type elastic critical component[J]. Journal of Vibration and Shock, 2012, 31(15):79-81.

[16]高德,盧富德.基于桿式彈性易損部件的非線性系統(tǒng)跌落沖擊研究[J].振動與沖擊,2012,31(15):7-49,59.

GAO De, LU Fu-de.Drop impact analysis of packaging system with bar type elastic critical components[J]. Journal of Vibration and Shock, 2012, 31(15):47-49.

[17]Gao D, Lu F D, Chen S J. Drop impact analysis of cushioning system with an elastic critical component of cantilever beam type[J]. Mathematical Problems in Engineering,2013, Article(379068):1-5.

[18]Rao S S. The finite element method in engineering[M]. Butterworth-heinemann, 2005.

[19]陳安軍. 非線性包裝系統(tǒng)跌落沖擊問題變分迭代法[J]. 振動與沖擊, 2013, 32(18):105-107.

CHEN An-jun. Variational iteration method for droppingshock problem of a cubic non-linear packaging system[J]. Journal of Vibration and Shock, 2013, 32(18):105-107.

[20]陳安軍. 矩形脈沖激勵下斜支承彈簧系統(tǒng)沖擊特性的研究[J]. 振動與沖擊, 2010, 29(10):225-227.

CHEN An-jun. Shock characteristics of a tilted support spring system under action of a rectangular pulse[J]. Journal of Vibration and Shock, 2010, 29(10):225-227.

[21]杜紹洪. 梁自由橫振動方程的有限差分方法[J]. 重慶交通大學(xué)學(xué)報:自然科學(xué)版, 2012, 31(1):6-10.

DU Shao-hong. Finite difference method of free transverse vibration equations of beams[J]. Journal of Chongqing Jiaotong University:Natural Science, 2012, 31(1):6-10.

[22]丁科, 殷水平. 有限單元法[M]. 2版.北京:北京大學(xué)出版社, 2012.

[23]王勖成, 邵敏. 有限單元法基本原理和數(shù)值方法[M]. 2版.北京:清華大學(xué)出版社, 1991.

第一作者胡育佳男,副教授,1979年生

猜你喜歡
有限元
基于擴展有限元的疲勞裂紋擴展分析
非線性感應(yīng)加熱問題的全離散有限元方法
TDDH型停車器制動過程有限元分析
新型有機玻璃在站臺門的應(yīng)用及有限元分析
基于I-DEAS的履帶起重機主機有限元計算
基于有限元模型對踝模擬扭傷機制的探討
10MN快鍛液壓機有限元分析
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 午夜丁香婷婷| 久久久久久久蜜桃| 91亚洲视频下载| 国产一区二区三区免费| 国产精品污视频| 91久久偷偷做嫩草影院| 无码专区在线观看| 99视频在线看| 性网站在线观看| 国产丝袜无码精品| 狂欢视频在线观看不卡| 久热99这里只有精品视频6| 国产导航在线| 日本日韩欧美| 亚洲综合久久成人AV| 欧类av怡春院| 伊人婷婷色香五月综合缴缴情| 久久天天躁狠狠躁夜夜2020一| 免费人成网站在线高清| 综合色区亚洲熟妇在线| 亚洲日韩在线满18点击进入| 国产美女丝袜高潮| 久久精品只有这里有| 国产亚洲高清在线精品99| 性欧美精品xxxx| 午夜福利视频一区| 欧美国产日韩在线| 永久成人无码激情视频免费| 国产一区二区三区精品久久呦| 国产AV无码专区亚洲A∨毛片| 3D动漫精品啪啪一区二区下载| 日韩小视频在线观看| 高清欧美性猛交XXXX黑人猛交| 最新亚洲人成无码网站欣赏网| 狠狠干欧美| 国产色婷婷| 国产精品第一区在线观看| 国产美女久久久久不卡| а∨天堂一区中文字幕| 中文精品久久久久国产网址 | 国产精品不卡永久免费| 亚洲成aⅴ人在线观看| 日本三级黄在线观看| 欧美日韩专区| 亚洲综合精品香蕉久久网| 伊人久久大香线蕉影院| 亚洲无码视频一区二区三区| 亚洲二区视频| 日韩亚洲综合在线| 日韩高清中文字幕| 免费国产黄线在线观看| 欧美综合区自拍亚洲综合绿色| 国产成人亚洲精品蜜芽影院| 亚洲一级毛片| 久久久久久久久亚洲精品| 国产第一页亚洲| 国国产a国产片免费麻豆| 波多野结衣AV无码久久一区| 国产视频只有无码精品| 日韩精品视频久久| 亚洲av日韩av制服丝袜| 国产在线自乱拍播放| 日本高清成本人视频一区| 手机在线国产精品| 日韩精品专区免费无码aⅴ| 十八禁美女裸体网站| 亚洲香蕉在线| 1024你懂的国产精品| 国产亚洲精品资源在线26u| 免费观看国产小粉嫩喷水| 亚洲香蕉在线| 黄色网站不卡无码| 欧美日韩国产在线人| 精品国产Av电影无码久久久| 色老二精品视频在线观看| 欧美精品伊人久久| 国产91丝袜在线播放动漫 | 欧美人与动牲交a欧美精品| 全部免费毛片免费播放| 亚洲国产第一区二区香蕉| 精品综合久久久久久97超人| 亚洲综合色婷婷|