• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Z型折疊板內(nèi)共振下非線性振動(dòng)特性研究

    2018-06-14 14:54郭翔鷹張楊張偉
    振動(dòng)工程學(xué)報(bào) 2018年2期
    關(guān)鍵詞:數(shù)值模擬

    郭翔鷹 張楊 張偉

    摘要:Z型折疊板是一類復(fù)雜多體結(jié)構(gòu),在工程實(shí)際中應(yīng)用廣泛。以大型的空間可展開(kāi)結(jié)構(gòu)為實(shí)際工程背景,考慮了Z型折疊板折疊過(guò)程中各部件的運(yùn)動(dòng)耦合關(guān)系,基于Reddy經(jīng)典板理論和von-Karman幾何非線性應(yīng)變位移關(guān)系,利用Hamilton原理建立了Z型折疊板的非線性動(dòng)力學(xué)控制方程。利用有限元分析得到了Z型折疊板結(jié)構(gòu)的模態(tài)函數(shù),運(yùn)用二階Galerkin截?cái)?,得到了Z型折疊板二自由度的非線性常微分方程。考慮系統(tǒng)主參數(shù)共振-1∶2內(nèi)共振的情況,通過(guò)攝動(dòng)分析得到系統(tǒng)四維直角坐標(biāo)形式的平均方程,最后利用數(shù)值模擬方法研究了橫向激勵(lì)對(duì)Z型折疊板非線性動(dòng)力學(xué)特性的影響。

    關(guān)鍵詞: 非線性動(dòng)力學(xué); Z型折疊板; 內(nèi)共振; 數(shù)值模擬

    中圖分類號(hào): O322文獻(xiàn)標(biāo)志碼: A文章編號(hào): 1004-4523(2018)02-0183-15

    DOI:10.16385/j.cnki.issn.1004-4523.2018.02.001

    引言

    可展開(kāi)(折疊)結(jié)構(gòu)具有悠久的研究歷史和工程應(yīng)用。折疊板結(jié)構(gòu)是可展結(jié)構(gòu)體系中最常用的一種,近年來(lái)在航空航天和建筑領(lǐng)域中得到了廣泛的應(yīng)用。在折疊狀態(tài)下,結(jié)構(gòu)體積較小,可用于運(yùn)輸或存儲(chǔ);在外力作用或系統(tǒng)內(nèi)部驅(qū)動(dòng)下,結(jié)構(gòu)逐步展開(kāi),最終達(dá)到完全展開(kāi)的工作狀態(tài),然后鎖定為穩(wěn)定狀態(tài)[1-2]。其中,Z型折疊板結(jié)構(gòu)是此類折疊結(jié)構(gòu)中最為經(jīng)典的一種,在航空領(lǐng)域主要被應(yīng)用于可變體飛行器機(jī)翼結(jié)構(gòu),而其中運(yùn)動(dòng)過(guò)程中由于變形引起的結(jié)構(gòu)非線性動(dòng)態(tài)特性問(wèn)題是結(jié)構(gòu)平穩(wěn)運(yùn)行的關(guān)鍵。

    在Z型折疊板結(jié)構(gòu)振動(dòng)特性的研究方面,國(guó)內(nèi)外學(xué)者進(jìn)行了大量的研究,如:Lee等[3-4]分別利用有限元方法和高階板理論研究了復(fù)合材料折疊板的振動(dòng)特性。Pal等[5-8]分別使用有限元方法,高階層合板理論研究了復(fù)合材料層合板折疊結(jié)構(gòu)的自由振動(dòng)特性。Topal等[9] 利用一階板理論建立了具有對(duì)稱折疊角度的復(fù)合材料折疊板的動(dòng)力學(xué)模型,通過(guò)MFD方法進(jìn)行頻率優(yōu)化,數(shù)值分析得出層合板的長(zhǎng)厚度比、夾角、板的長(zhǎng)度以及邊界條件會(huì)影響到其頻率特性。Jian等[10]建立了有限元模型,用一種條狀網(wǎng)格劃分單元,模擬了折疊板的靜態(tài)特征和動(dòng)力學(xué)特性分析。Liu等[11]對(duì)折疊板多體系統(tǒng)進(jìn)行了動(dòng)力學(xué)建模,并用數(shù)值方法研究了其動(dòng)態(tài)特性。張偉等[12]建立了變角度Z型梁的動(dòng)力學(xué)方程,計(jì)算了結(jié)構(gòu)的固有頻率和振動(dòng)模態(tài),并通過(guò)數(shù)值仿真和實(shí)驗(yàn)驗(yàn)證。此外,Zhang等[13-14]研究了復(fù)合材料正交鋪設(shè)層合板結(jié)構(gòu)在外激勵(lì)作用下的非線性振動(dòng)問(wèn)題,分析了系統(tǒng)不同參數(shù)下的振動(dòng)響應(yīng)。浙江大學(xué)的趙孟良和關(guān)富玲等[15-16]研究了空間可展結(jié)構(gòu)在外載荷作用下的運(yùn)動(dòng)特性。

    然而,Z型折疊結(jié)構(gòu)屬于多體結(jié)構(gòu),在外激勵(lì)作用下將產(chǎn)生大幅的非線性耦合振動(dòng)響應(yīng),這部分的研究成果目前還很有限。本文根據(jù)實(shí)際工程背景,研究了一類Z型折疊板結(jié)構(gòu)在一定頻率的外激勵(lì)作用下發(fā)生1∶2內(nèi)共振情況下的非線性動(dòng)力學(xué)特性。

    1動(dòng)力學(xué)模型和方程〖2〗1.1動(dòng)力學(xué)模型以大型的空間可展開(kāi)結(jié)構(gòu)為實(shí)際工程背景,建立如圖1所示的力學(xué)模型,所研究的Z型折疊板采用碳纖維復(fù)合材料層合板,折疊結(jié)構(gòu)由三部分組成,靠近固定端的部分稱為內(nèi)板,用Ω1表示,完成翻轉(zhuǎn)動(dòng)作的稱為中間板,用Ω2表示,中間板外側(cè)部分為外板,用Ω3表示。圖中O1x1y1,O2x2y2,O3x3y3分別為內(nèi)板、中間板、外板的局部坐標(biāo)系,全局坐標(biāo)系Oxy與局部坐標(biāo)系O1x1y1重合,各個(gè)板表面受到的橫向簡(jiǎn)諧激勵(lì)表示為P1,P2和P3,如圖1所示。

    模型中,各個(gè)部分之間通過(guò)剛性鉸鏈相連接,通過(guò)內(nèi)部的作動(dòng)器和機(jī)械結(jié)構(gòu)進(jìn)行折疊與展開(kāi)動(dòng)作,機(jī)構(gòu)在變形時(shí),內(nèi)板與固定端相連接,中間板在一定的角度范圍內(nèi)轉(zhuǎn)動(dòng),外板始終與內(nèi)板保持平行。

    本文研究過(guò)程中,滿足以下條件:

    1)Z型折疊板3塊板是等厚度,等寬度的;

    2)Z型折疊板3塊板在外激勵(lì)作用下不會(huì)產(chǎn)生材料本身的破壞。

    圖1Z型折疊板的動(dòng)力學(xué)模型

    Fig.1Mechanical model of the Z-type folding plates第2期郭翔鷹,等: Z型折疊板內(nèi)共振下非線性振動(dòng)特性研究振 動(dòng) 工 程 學(xué) 報(bào)第31卷1.2動(dòng)力學(xué)方程的建立

    本文選取折疊板材料為正交鋪設(shè)的碳纖維復(fù)合材料層合板,不考慮剪切效應(yīng),因此(σz )z=+ /-h2=0,(τxz)z=+/-h2=0,(τyz)z=+/-h2=0,其中h為層合板結(jié)構(gòu)厚度。

    材料的本構(gòu)關(guān)系如下

    σxx

    σyy

    τyz

    τzx

    τxy=1112000

    2122000

    004400

    000550

    000066εxx

    εyy

    γyz

    γzx

    γxy(1)

    式中ij為轉(zhuǎn)換彈性常數(shù),定義為ij=T-1QijT-1T(2)并且T-1=cos2θsin2θ0

    sin2θ〖〗cos2θ0

    00cos2θ-sin2θ(3)正交鋪設(shè)復(fù)合材料層合板的彈性模量Qij(i=1,2,4,5,6; j=1,2,4,5,6)可以表示為如下形式:Q11=Q22=E1-ν2,Q12=Q21=νE1-ν2,

    Q44=Q55=Q66=E2(1+ν)(4)內(nèi)力和彎矩的關(guān)系可表示為:Nxxi

    Nyyi

    Nzzi=∑Nk=1∫zk+1zkσxxi

    σyyi

    σzzidz (5)

    Mxxi

    Myyi

    Mzzi=∑Nk=1∫zk+1zkσxxi

    σyyi

    σzzizdz(6)本文所研究的Z型可折疊結(jié)構(gòu)在橫向外載荷的作用下會(huì)產(chǎn)生較大的變形,因此這里用von-Karman大變形理論進(jìn)行分析,各個(gè)板的非線性應(yīng)變表達(dá)式為εx

    εy

    γxy=ε0+zε1(7)式中ε0=u0ixi+12w0ixi2

    v0iyi+12w0iyi2

    v0ixi+u0iyi+w0ixiw0iyi(8)

    ε1=-ω20ix2i

    -ω20iy2i

    -2ω20ixiyi (9)式中下標(biāo)i表示第i個(gè)板(i=1,2,3)。

    根據(jù)正交鋪設(shè)碳纖維增強(qiáng)復(fù)合材料的結(jié)構(gòu)特點(diǎn),可得到內(nèi)力與應(yīng)變關(guān)系為N

    M=AB

    BDε1

    ε2 (10)式中剛度矩陣Aij, Bij, Dij, 用如下形式表示Aij,Bij,Dij=∫h2-h2ij1,z,z2dz(11)根據(jù)經(jīng)典的層合板理論[17],3塊板上任意一點(diǎn)在局部坐標(biāo)系中的位移分別表示為以下形式:ui=u0i(xi,yi,zi,t)-ziw0ixi(12a)

    vi=v0i(xi,yi,zi,t)-ziw0iyi(12b)

    wi=w0i(xi,yi,zi,t),(i=1,2,3)(12c)式中u0i,v0i,w0i是在板Ωi中性層上任意一點(diǎn)在X, Y, Z方向的位移。

    根據(jù)本文模型中所建立的坐標(biāo)轉(zhuǎn)換關(guān)系,Z型折疊板結(jié)構(gòu)在全局坐標(biāo)系下的位移場(chǎng)可描述為:r1=R1u1

    v1

    w1 (13a)

    r2 = r1 L1 cosw1′(L1 )

    v2

    w2 + R2 u2

    v2

    w2 (13b)

    r3 = r2 L2 cosw2′(L2 )

    v2

    w2 (L2 ) + u3

    v3

    w3 (13c)式中Li為各個(gè)板的長(zhǎng)度, Ri 為坐標(biāo)轉(zhuǎn)換矩陣,形式如下R1=1

    1

    1 (14a)

    R2=cosθ0-sinθ

    0〖〗10

    sinθ0cosθ (14b)將上述表達(dá)式(1)~(14)代入Hamilton原理中,整理展開(kāi)后得到用廣義位移表示的Z型折疊板結(jié)構(gòu)的非線性動(dòng)力學(xué)方程如下。

    內(nèi)板的非線性動(dòng)力學(xué)方程為:

    A11(2u01x21+w01x12w01x21)+A12(2u01x21+

    w01x12w01x21)+A66(v01x1+u01y1+w01x1w01y1)=

    I02u01t2-I13w01t2x1(15a)

    A12(2v01y21+w01y12w01y21)+A22(2v01y21+

    w01y12w01y21)+A66(v01x1+u01y1+w01x1w01y1)=

    I02v01t2-I1(3w01〖〗t2y1)(15b)

    A11[2u01x21w01x1+u01x12w01x21+3〖〗2w01x122w01x21]+

    A12[2u01x21w01x1+u01x12w01x21+32w01x122w01x21]+

    A12[2v01y21w01y1+v01y12w01y21+32w01y122w01y21]+

    A22[2v01y21w01y1+v01y12w01y21+32w01y122w01y21]-

    (D11+D12)4w01x41-(D22+D12)4w01y41-

    8D664ω01x21y21+A66[2v01x12ω01x1y1+2u01y12ω01x1y1+

    4w01x1w01y12ω01x1y1+v201x1y1w01x1+2u01x1y1w01x1+ 2v01x21w01y1+2u01x1y1w01y1+2w01x21w01y12+

    2w01y21w01x12]+P1=I02w01t2+I1(3u01t2x1+

    3v01t2y1)-I2(3w01t2x1+3w01t2y1)(15c)

    中間板的非線性動(dòng)力學(xué)方程為:

    A11(2u02x22+w02x22w02x22)+A12(2u02x22+

    w02x22w02x22)+A66(v02x2+u02y2+w02x2w02y2)=

    I0(sinkθ+coskθ)2u02t2+I0(sinkθ+

    coskθ)(2w01(L1)tx1+)u02-I1(sinkθ+

    coskθ)3w01t2x2-I1(sinkθ+coskθ)A12(2v02〖〗y(tǒng)22+

    w02y22w02y22)+A22(2v02y22+w02y22w02y22)+A66(v02x2+

    u02y2+w02x2w02y2)=I02v02t2-I1(3w02t2y2)(15d)

    A12(2v02y22+w02y22w02y22)+A22(2v02y22+

    w02y22w02y22)+A66(v02〖〗x2+u02y2+w02x2w02y2)=

    I02v02t2-I1(3w02t2y2) (15e)

    A11[2u02x22w02x2+u02x22w02x22+w02x222w02x22]+

    A12[2u02x22w02x2+u02〖〗x22w02x22+

    32w02x222w02x22]+A12[2v02y22w02y2+

    v02y22w02y22+32w02y222w02y22]+A22[2v02y22w02y2+

    v02y22w02y22+32w02y222w02y22]-(D11+

    D12)4w02x42-(D22+D12)4w02y42-8D664ω02x22y22+

    A66[2v02x22ω02x2y2+2u02y22ω02x2y2+

    4w02x2w02y22ω02x2y2+v202x2y2w02〖〗x2+

    2u02x2y2w02x2+2v02x22w02y2+2u02x2y2w02y2+

    2w02x22w02y22+2w02y22w02x22]+P2=

    -I02sinkθ2w02t2+I02coskθ(2w01(L1)tx1+

    )w02+I1(sinkθ+coskθ)3u02t2x2+I1(sinkθ+

    coskθ)(2w01(L1)tx1+)u02x2+I1+3w02t2y2-

    I2(sinkθ+coskθ)4w02t2x22+I2(sinkθ+

    coskθ)(2w01(L1)tx1+)w202x22-I24w02t2y22 (15f)

    外板的非線性動(dòng)力學(xué)方程為:

    A11(2u03x23+w03x32w03x23)+A12(2u03x23+

    w03〖〗x32w03x23)+A66(v03x3+u03y3+w03〖〗x3w03y3)=

    I02u03t2-I13w03t2x3+(-θ)(15g)

    A12(2v03y23+w03y32w03〖〗y(tǒng)23)+A22(2v03y23+

    w03y32w03y23)+A66(v03x3+u03y3+w03x3w03y3)=

    I02v03t2-I1(3w03t2y3)(15h)

    A11[2u03x23w03x3+u03x32w03x23+32w03x322w03x23]+

    A12[2u03x23w03x3+u03x32w03x23+32w03x322w03x23]+

    A12[2v03y23w03y3+v03y32w03y23+32w03y322w03y23]+

    A22[2v03y23w03y3+v03y32w03y23+32w03y322w03y23]-

    (D11+D12)4w03x43-(D22+D12)4w03y43-D664ω03x23y23+

    A66[2v03x32ω03x3y3+2u03y32ω03x3y3+

    4w03x3w03y32ω03x3y3+v203x3y3w03x3+2u03x3y3w03x3+

    2v03x23w03y3+2u03x3y3w03y3+2w03x23w03y32+

    2w03y23w03x32]+P3=I02w03t2+I13u03t2x3+

    I13u03t2y3-I24w03t2x23-I24w03t2y23(15i)

    結(jié)構(gòu)整體的邊界條件滿足如下等式:u1(0)=0,u1 x1 = L1 = u2 (0) = L1 (16a)

    u2 x2 = L2 = u3 (0) = L1 + L2 cosθ(16b)

    u3 x3 = L3 = L1 + L2 cosθ + L3 (16c)

    w1(0)=0,w1(L1)=w2(0)=0(16d)

    w3 (0) = w2 (L2 ) = L2 sinθ(16e)

    w1′(0)=0,w2″(0)=0,w3(0)=0(16f)

    w01 (0)y1 x1 = 0 = 0 (16g)

    2w2(0)y22x1 = L1 = 2w3(0)y23x2 = L2 = 0(16h)

    w01 (b)y1 x1 = 0 = 0(16i)

    2w2(b)y22x1 = L1 = 2w3(b)y23x2 = L2 = 0(16j)

    2w03 tx3 x3 = L3 = L1 + L2 cosθ + L3 (16k)

    (A12+A22)[v0iyi+12(w0iyi)2]y=0,b=0 (16l)

    2w02 tx2 x2 = L2 = L1 + L2 cosθ (16m)2有限元模態(tài)分析

    上述所建立的動(dòng)力學(xué)控制方程是偏微分方程,利用數(shù)學(xué)方法直接求解極為困難,因此本文將采用Galerkin離散將偏微分方程轉(zhuǎn)換到常微分方程后進(jìn)行求解。對(duì)于單一結(jié)構(gòu)模型,其模態(tài)函數(shù)通常是可以直接根據(jù)經(jīng)驗(yàn)公式假設(shè)的,但Z型折疊板結(jié)構(gòu)為多體結(jié)構(gòu),很難確定折疊板在外激振力作用下結(jié)構(gòu)的振動(dòng)模態(tài)。因此首先通過(guò)ANSYS有限元方法,對(duì)Z型折疊板進(jìn)行模態(tài)分析和諧響應(yīng)分析,得到Z型折疊板結(jié)構(gòu)的固有頻率和模態(tài)。通過(guò)研究結(jié)構(gòu)模態(tài)振型,確定系統(tǒng)的模態(tài)函數(shù)形式。最后,將無(wú)量綱形式的Z型折疊板結(jié)構(gòu)的非線性動(dòng)力學(xué)控制方程通過(guò)Galerkin方法進(jìn)行二階離散,得到可求解的常微分方程組。

    2.1有限元模型

    航空領(lǐng)域?qū)嶋H應(yīng)用中,Z型折疊板結(jié)構(gòu)在折疊展開(kāi)運(yùn)動(dòng)過(guò)程中的角速度很小,且折疊角度在0~150°的范圍內(nèi),因此本文將板折疊過(guò)程分解,看作是不同折疊角度板結(jié)構(gòu)的準(zhǔn)靜態(tài)慢變過(guò)程,選取幾個(gè)特定折疊角度來(lái)建立Z型碳纖維復(fù)合材料層合板結(jié)構(gòu)的力學(xué)模型。

    本章研究中取折疊角度為60°,90°和120°作為典型參數(shù)值,進(jìn)行對(duì)比分析和討論。

    有限元模型統(tǒng)一采用四邊形板單元,如圖2所示;設(shè)置彈性模量為5×105 MPa,泊松比為0.3。Z型折疊板的有限元模型尺寸參數(shù)如表1所示。

    圖2模型網(wǎng)格劃分示意圖

    Fig.2Meshing of the Z-type folding plate wing of the three angles

    表1Z型折疊板有限元模型的尺寸參數(shù)

    Tab.1Geometric parameters of the finite element model for the Z-type folding wing

    板長(zhǎng)度/m厚度/m寬度/m內(nèi)板20.012.1中間板10.012.1外板40.012.12.2模態(tài)分析

    通過(guò)對(duì)建立的有限元模型施加初始條件,進(jìn)行模態(tài)分析,得到Z型折疊板結(jié)構(gòu)不同折疊角度下前5階固有頻率如表2所示。

    表2不同折疊角度前5階固有頻率(單位:Hz)

    Tab.2The first five order natural frequencies of different folding angles (Unit:Hz)

    階數(shù)60°90°120°10.560180.617540.6944322.711402.358002.6580033.225903.319903.6326043.550903.731604.2836057.367306.220905.96140

    下面列出Z型折疊板結(jié)構(gòu)在60°,90°,120°的折疊角度下的模態(tài)振型圖,如圖3~5所示。

    由有限元模態(tài)分析結(jié)果可以看出,不同折疊角度下的Z型折疊板結(jié)構(gòu)前5階模態(tài)可表示為: 第1圖3折疊角度為60°時(shí)Z型折疊板前5階模態(tài)振型圖

    Fig.3The first five order mode shapes of the Z-type folding angle 60°圖4折疊角度為90°時(shí)Z型折疊板前5階模態(tài)振型圖

    Fig.4The first five order mode shapes of the Z-type folding angle 90°圖5折疊角度為120°時(shí)Z型折疊板前5階模態(tài)振型圖

    Fig.5The first five order mode shapes of the Z-type folding angle 120°

    階為彎曲振動(dòng),第2階為扭轉(zhuǎn)振動(dòng),第3階為彎曲振動(dòng),第4階為彎扭耦合振動(dòng),第5階為彎曲振動(dòng)。通過(guò)以上分析可以發(fā)現(xiàn),Z型折疊板結(jié)構(gòu)的前5階振動(dòng)模態(tài)的形式與懸臂板結(jié)構(gòu)前5階振動(dòng)模態(tài)的形式相似。

    3Galerkin離散

    由于結(jié)構(gòu)在共振情況下會(huì)發(fā)生劇烈的大幅振動(dòng),容易產(chǎn)生失穩(wěn)及結(jié)構(gòu)整體破壞等。因此,對(duì)于結(jié)構(gòu)在共振情況下動(dòng)力學(xué)響應(yīng)的研究是十分必要的。

    在有限元分析結(jié)果中可以發(fā)現(xiàn),Z型折疊板結(jié)構(gòu)在折疊角度為60°時(shí),第4階固有頻率幾乎是第5階固有頻率1/2,存在1∶2內(nèi)共振的情況,因此,本文將對(duì)Z型折疊板結(jié)構(gòu)在主參數(shù)共振-1∶2內(nèi)共振的情況下的動(dòng)力學(xué)特性進(jìn)行深入分析。

    根據(jù)上述的分析結(jié)果,可選取懸臂板結(jié)構(gòu)的振動(dòng)模態(tài)函數(shù)形式作為Z型折疊板在折疊角度為60°時(shí)的模態(tài)函數(shù),利用Galerkin方法對(duì)方程(15)進(jìn)行二階截?cái)?,得到系統(tǒng)常微分形式的非線性動(dòng)力學(xué)方程。在滿足位移邊界條件的前提下,選取3個(gè)方向的基本函數(shù)分別如下:

    Ui(x,y,t) = u1 (t)sinx/2πLi cosπyb +

    u2 (t)sin3x/2πLi cos2πyb(17a)

    Vi (x,y,t) = v1 (t)cosx/2πLi sinπyb +

    v2 (t)cos3x/2πLi sin2πyb(17b)

    Wi(x,y,t)=w1(t)Xi1(x)Yi1(y)+

    w2(t)Xi2(x)Yi2(y) (17c)

    根據(jù)4階常微分方程的通解,可將方程(17c)中的Xij(x)和Yi1(y)取為Xij (x) = Ai1 coshki x-Ai2 coski x-

    Ai3 βi1 sinhki x + Ai4 sinki x,

    Yij (x) = Bi1 coshki x-Bi2 coski x-

    Bi3 βi1 sinhki x + Bi4 sinki x(18)式中Xij(x)是沿x軸方向的固支-自由梁函數(shù),Yij(y)為y軸方向的自由-自由梁函數(shù),ki1和ki2為特征方程的根,并有如下關(guān)系coski1coshki1+1=0

    coski2coshki2-1=0(19)將系統(tǒng)實(shí)際參數(shù)L1=2 m,L2=1 m,L3=4 m代入式(17)可求得X,Y方向的模態(tài)函數(shù),將方程(18),(19)代入邊界條件求得模態(tài)參數(shù)Aij和Bij,得到Xij(x)和Yij(y)的函數(shù)表達(dá)式,如下:

    X11 =-4.9×10 - 5coshk11 x +

    2.1×10 - 5cosk11 x - 6.6×10 - 5β11 ·

    sinhk11 x-sink11 x (20a)

    X12 = 0.25coshk12 x + 27.3cosk12 x -

    0.25β12 sinhk12 x - sink12 x(20b)

    X21 = -1.4×10 - 6coshk21 x +

    1.3×10 - 6cosk21 x - 2.5×10 - 7β21 ·

    sinhk21 x-sink21 x(20c)

    X22 = 1.2×10 - 4coshk22 x-1.4cosk22 x -

    0.12β22 sinhk22 x-sink22 x(20d)

    X31 = -1.5×10 - 4coshk31 x +

    3.4×10 - 4cosk31 x - 3×10 - 4β31 ·

    sinhk31 x-sink31 x(20e)

    X32 =0.25coshk32 x + 27.3cosk32 x -

    0.25β32 sinhk32 x - sink32 x (20f)

    同時(shí),將外激勵(lì)也進(jìn)行離散,表示為

    Pi=Fi1sin3πxlisinπyb+Fi2sinπxlisin3πyb (21)

    用二階Galerkin方法離散方程(15),并將離散后的面內(nèi)位移u0,v0用橫向位移w0表示,整理可得到Z型折疊板橫向振動(dòng)的常微分運(yùn)動(dòng)控制方程:

    內(nèi)板方程:

    11+μ111+21w11+α112w12+α113w311+α114w312+ α115w12w211+α116w11w212=f11cos(ξt)(22a)

    12+μ212+α121w11+22w12+α123w311+α123w312+α125w12w211+α126w11w212=f12cos(ξt)(22b)

    中間板方程:

    21+μ321+α21121w21+α21222w21+α21321w22+α21422w22+21w21+α216w22+α217w21u21+

    α218w21u22+α219w22u21+α2110w22u22+

    (α2111/θ)w321+(α2112/θ)w322+α2113(sinθ+

    cosθ)w22w221+α2114(sinθ+cosθ)w21w222=

    f21cos(ξt) (22c)

    22+μ422+α22121w21+α22222w21+α22321w22+α22422w22+α225w21+22w22+α227w21u21+

    α228w21u22+α229w22u21+α2210w22u22+w321+(α2212/θ)w322+α2213(sinθ+cosθ)w22w221+

    α2214(sinθ+cosθ)w21w222=f22cos(ξt) (22d)

    外板方程:

    31+μ531+21w31+(α311+F1cosΩt)ω31+

    α312w32+α313w331+α314w332+α315w32w231+

    α316w31w232=f31cos(ξt)(22e)

    32+μ632+22w32+(α322+F2cosΩt)ω32+

    α321w31+α323w331+α324w332+α325w32w231+

    α326w31w232=f32cos(ξt)(22f)

    式中折疊角θ為3塊板的運(yùn)動(dòng)方程的連接參數(shù),體現(xiàn)了3塊板之間的耦合運(yùn)動(dòng)關(guān)系。

    4攝動(dòng)分析

    對(duì)于較復(fù)雜的非線性常微分方程,很難求出其精確解,需要用近似解析的方法求其漸近解來(lái)替代精確解。攝動(dòng)分析是近似解析的一種方法,包括直接攝動(dòng)法、多尺度法以及KBM法等。

    因此,本章基于系統(tǒng)主參數(shù)共振-1∶2內(nèi)共振的共振關(guān)系,使用多尺度方法進(jìn)行攝動(dòng)分析。系統(tǒng)共振關(guān)系表示如下:21=14Ω2+εσ1, 22=Ω2+εσ2(23)式中1,2是相應(yīng)線性系統(tǒng)的第1階、第2階固有頻率。σ1,σ2為系統(tǒng)的調(diào)諧參數(shù),為了方便分析,令Ω=1。

    經(jīng)過(guò)計(jì)算得到系統(tǒng)的直角坐標(biāo)下平均方程為:

    內(nèi)板平均方程:

    11=-12u1x11-σ11x12-3α113x211x12- 3α113x312-2α116x12x213+x214(24a)

    12=-12u1x12+σ11x11+3α113x11x212+x311+ 2α116x11x213+x214(24b)

    13=-12u2x13-12σ12x14-32α124x213x14+x314-α125x14x211+x212(24c)

    14=-12u2x14+12σ12x13+32α124x214x13+x313+α125x13x211+x212-14f12(24d)

    中間板平均方程:

    21=-12u3x21-σ11x22+14α212x21x24-x22x23- 3α214x22x221+x222-2α217x21x223+x224(25a)

    22=-12u3x22+σ11x21-14α212x21x23-x22x24+ 3α214x21x221+x222+2α217x21x223+x224(25b)

    23=-12u4x23-32α225x24x223+x224-

    α226x24x221+x222-14σ2x4(25c)

    24=-12u4x24+32α225x23x223+x224+

    α226x23x221+x222+14σ2x3-14f22 (25d)

    外板方程:

    31=-12u5x31-σ1x32-3α313x231x32+x332-

    2α316x32x233+x234+12f31x32(26a)

    32=-12u5x32+σ1x31+3α313x31x232+x331+

    2α316x31x233+x234+12f31x32(26b)

    33=-12u6x33-12σ2x34-32α324x233x34+x334-α325x34x231+x232 (26c)

    34=-12u6x34+12σ2x33+32α324x234x33+x333+α325x33x231+x232-12f32(26d)

    5數(shù)值模擬

    根據(jù)數(shù)值分析結(jié)果發(fā)現(xiàn),Z型折疊板內(nèi)板的振動(dòng)幅值很小且多為周期性顫振[18],考慮到本文篇幅的限制,不再詳述,主要討論Z型折疊機(jī)板在橫向激勵(lì)作用下中間板和外板的非線性振動(dòng)特性。

    5.1幅頻響應(yīng)特性分析

    通過(guò)數(shù)值求解系統(tǒng)的四維平均方程,利用matlab軟件繪制3塊板的幅頻響應(yīng)曲線,選取外激勵(lì)幅值和系統(tǒng)阻尼系數(shù)為控制參數(shù),研究參數(shù)對(duì)系統(tǒng)幅頻特性的影響。

    首先,根據(jù)實(shí)際參數(shù)的取值范圍,經(jīng)過(guò)無(wú)量綱處理后,選取參數(shù)為μ1=0.25,μ3=0.22,μ4=0.14, σ3=-0.015,σ4=-0.014,α214=0.31,α217=0.0002,α225=-0.108,α226=3.5,μ5=0.37,μ6=0.56,σ5=1.77,σ6=1.88,α313=2.49,α316=-3.27,α324=7.16,α325=6.81,α313=2.49。

    將2塊板的初始條件均設(shè)為x10=1.44,x20=1.55,x30=1.35,x40=-1.799。令外激勵(lì)幅值fi (i=2,3)的值分別為50和100,研究結(jié)構(gòu)幅頻響應(yīng)曲線的變化。圖中藍(lán)色曲線和紅色曲線分別表示結(jié)構(gòu)第4階模態(tài)和第5階模態(tài)的幅頻響應(yīng)曲線,橫坐標(biāo)為調(diào)頻參數(shù)σi (i=2,3),縱坐標(biāo)為振動(dòng)幅值ai (i=1,2), a1, a2 分別表示第4階和第5階的振動(dòng)幅值。

    由圖6和7可知,隨著外激勵(lì)幅值的增加,幅頻

    圖6不同外激勵(lì)幅值下的中間板Ω2幅頻響應(yīng)曲線

    Fig.6The frequency-response curves of the middle plate Ω2 to the external excitation amplitude f2

    圖7不同外激勵(lì)幅值下的外板Ω3幅頻響應(yīng)曲線

    Fig.7The frequency-response curves of the outer plate Ω3 to the external excitation amplitude f3

    響應(yīng)曲線的形態(tài)發(fā)生了不規(guī)律變化。隨著調(diào)頻參數(shù)的增加,系統(tǒng)振幅除中間板第5階模態(tài)以外都呈現(xiàn)減小的趨勢(shì),并且會(huì)出現(xiàn)多值現(xiàn)象和跳躍現(xiàn)象。

    為研究系統(tǒng)阻尼系數(shù)對(duì)幅頻響應(yīng)曲線的影響,分別令阻尼系數(shù)μi (i=3,5)的值為0.3和0.8。繪制如下幅頻響應(yīng)曲線。

    觀察圖8和9發(fā)現(xiàn),隨著阻尼的增加,系統(tǒng)的振動(dòng)幅值降低,幅頻響應(yīng)曲線形態(tài)發(fā)生變化,且對(duì)第4階模態(tài)的頻響特性影響較大。隨著調(diào)頻參數(shù)σ的增加,系統(tǒng)會(huì)出現(xiàn)多值和跳躍的現(xiàn)象。

    圖8不同阻尼下的中間板Ω2幅頻響應(yīng)曲線

    Fig.8The frequency-response curves of the middle plate Ω2 to the damping coefficient μ3

    圖9不同阻尼下的外板Ω3幅頻響應(yīng)曲線

    Fig.9The frequency-response curves of the outer plate Ω3 to the damping coefficient μ5

    5.2非線性振動(dòng)響應(yīng)分析

    為了研究Z型折疊板系統(tǒng)在主參數(shù)共振-1∶2內(nèi)共振的情況下的非線性振動(dòng)響應(yīng)特性,選取外激勵(lì)幅值fi為控制參數(shù),研究橫向激勵(lì)幅值對(duì)系統(tǒng)產(chǎn)生周期運(yùn)動(dòng)和混沌運(yùn)動(dòng)的影響。

    固定上述結(jié)構(gòu)參數(shù),改變外激勵(lì)的幅值,運(yùn)用四階龍格庫(kù)塔法對(duì)系統(tǒng)運(yùn)動(dòng)方程進(jìn)行數(shù)值求解,得到中間板和外板的混沌分叉圖,如圖10和11所示。

    圖10中間板Ω2隨外激勵(lì)幅值變化的分叉圖

    Fig.10The bifurcation diagram of the middle plate Ω2 with the transverse excitation f2

    圖11外板隨外激勵(lì)幅值變化的分叉圖

    Fig.11The bifurcation diagram of the outer plate Ω3 with the transverse excitation f3

    由圖中可以看出,在外激勵(lì)幅值增大的過(guò)程中,結(jié)構(gòu)會(huì)出現(xiàn)周期運(yùn)動(dòng)-混沌運(yùn)動(dòng)-周期運(yùn)動(dòng)-混沌運(yùn)動(dòng)的變化,說(shuō)明共振情況下,外激勵(lì)幅值的變化對(duì)于系統(tǒng)的運(yùn)動(dòng)穩(wěn)定性具有重要的影響。

    為了更好地描述上述分叉圖的特性,首先對(duì)中間板隨外激勵(lì)幅值變化的振動(dòng)特性給出了具體的分析:圖12,13和14分別給出了中間板在外激勵(lì)幅值為3,10和35時(shí)的波形圖、三維相圖和Poincaré截面,此時(shí)中間板的振動(dòng)是從單倍周期進(jìn)入短暫的混沌運(yùn)動(dòng)之后又變?yōu)楸吨芷谶\(yùn)動(dòng)。繼續(xù)增大外激勵(lì)幅值到,中間板的振動(dòng)形式趨于復(fù)雜,逐漸變?yōu)槿鐖D15所示的概周期運(yùn)動(dòng),最后進(jìn)入圖16所示的混沌運(yùn)動(dòng),且不會(huì)伴隨外激勵(lì)幅值繼續(xù)增大而改變運(yùn)動(dòng)形態(tài)。

    圖12f2=3時(shí)中間板的周期運(yùn)動(dòng)

    Fig.12The periodic motion of the middle plate when f2=3圖13f2=10時(shí)中間板的混沌運(yùn)動(dòng)

    Fig.13The chaotic motion of the middle plate when f2=10

    外板隨外激勵(lì)幅值變化的振動(dòng)特性如圖17~20所示。當(dāng)外激勵(lì)幅值小于10時(shí),外板呈現(xiàn)不規(guī)律的混沌運(yùn)動(dòng),隨著外激勵(lì)幅值的增大,外板的振動(dòng)形式趨于平穩(wěn),呈現(xiàn)周期運(yùn)動(dòng)形式;當(dāng)外激勵(lì)幅值增大至25~38之間時(shí),再次進(jìn)入混沌運(yùn)動(dòng),繼續(xù)增大外激勵(lì)幅值會(huì)發(fā)現(xiàn)外板再次變?yōu)橐?guī)律的周期運(yùn)動(dòng)。圖17~20分別是外激勵(lì)幅值為3, 15, 35和45時(shí)外板的波形圖、三維相圖和Poincaré截面。圖14f2=35時(shí)中間板的4倍周期運(yùn)動(dòng)

    Fig.14The period-4 motion of the middle plate when f2=35

    圖15f2=65時(shí)中間板的概周期運(yùn)動(dòng)

    Fig.15The quasi-period motion of the middle plate when f2=65圖16f2=80時(shí)中間板的混沌運(yùn)動(dòng)

    Fig.16The chaotic motion of the middle plate when f2=80

    圖17f3=3時(shí)外板的混沌運(yùn)動(dòng)

    Fig.17The chaotic motion of the outer plate when f3=3圖18f3=15時(shí)外板的周期運(yùn)動(dòng)

    Fig.18The periodic motion of the outer plate when f3=15

    圖19f3=35時(shí)外板的混沌運(yùn)動(dòng)

    Fig.19The chaotic motion of the outer plate when f3=35圖20f3=45時(shí)外板的周期運(yùn)動(dòng)

    Fig.20The periodic motion of the outer plate when f3=45此外,中間板第5階模態(tài)的振動(dòng)幅值比第4階模態(tài)的振動(dòng)幅值大,外板第4,5階模態(tài)的振動(dòng)幅值基本相差不大,這是由于系統(tǒng)在與第5階固有頻率對(duì)應(yīng)的外激勵(lì)作用下,在這兩階模態(tài)存在1∶2的關(guān)系時(shí)發(fā)生了耦合的內(nèi)共振現(xiàn)象。

    6結(jié)論

    本文利用Hamilton原理建立了在外激勵(lì)作用下Z型折疊板的幾何非線性動(dòng)力學(xué)方程,并對(duì)系統(tǒng)在主參數(shù)共振-1∶2內(nèi)共振情況下的非線性動(dòng)力學(xué)行為進(jìn)行了攝動(dòng)分析,得到系統(tǒng)4自由度的平均方程。利用數(shù)值方法分析了系統(tǒng)的幅頻響應(yīng)特性和混沌分叉特性。

    數(shù)值結(jié)果表明,不同的外激勵(lì)幅值和阻尼系數(shù)會(huì)對(duì)系統(tǒng)的頻響特性產(chǎn)生一定的影響,且隨著調(diào)頻參數(shù)σ的增大,對(duì)應(yīng)的系統(tǒng)振動(dòng)幅值會(huì)出現(xiàn)多值和跳躍的現(xiàn)象。

    選取一定的參數(shù)和初始條件,通過(guò)數(shù)值模擬發(fā)現(xiàn),在主參數(shù)共振-1∶2內(nèi)共振的共振關(guān)系下,當(dāng)外激勵(lì)的頻率與系統(tǒng)第5階固有頻率相同時(shí),只改變外激勵(lì)幅值時(shí),中間板會(huì)出現(xiàn)單倍周期-混沌-概周期-混沌運(yùn)動(dòng),外板會(huì)出現(xiàn)混沌-單倍周期-概周期-混沌運(yùn)動(dòng),由此可見(jiàn)外激勵(lì)的改變會(huì)對(duì)系統(tǒng)的非線性動(dòng)力學(xué)特性產(chǎn)生顯著的影響,且系統(tǒng)的第4階模態(tài)對(duì)應(yīng)的幅值也會(huì)產(chǎn)生明顯的變化,說(shuō)明此非線性系統(tǒng)的不同模態(tài)振動(dòng)之間存在復(fù)雜的耦合關(guān)系。

    因此,在研究Z型折疊板這一類結(jié)構(gòu)的非線性動(dòng)力學(xué)行為時(shí),不應(yīng)該只考慮單一的模態(tài)振動(dòng),還應(yīng)考慮多階模態(tài)之間的相互作用,以便更好地利用或控制其運(yùn)動(dòng)形式,為實(shí)際工程提供重要的理論依據(jù)。

    參考文獻(xiàn):

    [1]韓運(yùn)龍. 折疊板殼結(jié)構(gòu)的設(shè)計(jì)與分析[D]. 南京:東南大學(xué), 2011:02.

    Han Yunlong. Design and analysis of foldable plates structures[D]. Nanjing: Southeast University, 2011:02.

    [2]陳務(wù)軍, 關(guān)富玲, 陳向陽(yáng). 可折疊航天結(jié)構(gòu)展開(kāi)動(dòng)力學(xué)分析[J]. 計(jì)算力學(xué)學(xué)報(bào), 1999,16:4.

    Chen Wujun, Guan Fulin, Chen Xingyang. Dynamic analysis for deployment process of foldable aerospace structures[J]. Chinese Journal of Computational Mechanics, 1999, 16:4.

    [3]Lee S Y, Wooh S C. Finite element vibration analysis of composite box structures using the high order plate theory[J]. Journal of Sound and Vibration, 2004,277(4-5): 801—814.

    [4]Lee S Y, Wooh S C, Yhim S S. Dynamic behavior of folded composite plates analyzed by the third order plate theory[J]. International Journal of Solids and Structures, 2004,41(7):1879 —1892.

    [5]Pal S, Gu H, Niyogi A. Application of folded plate formulation in analyzing stiffened laminated composite and sandwich folded plate vibration[J]. Journal of Reinforced Plastics and Composites, 2008,27(7):693 —710.

    [6]Haldar S, Sheikh A H. Free vibration analysis of isotropic and composite folded plates using a shear flexible element[J]. Finite Elements in Analysis and Design, 2005,42(3):208—226.

    [7]Hernández E, Hervella-Nieto L. Finite element approximation of free vibration of folded plates[J]. Computer Methods in Applied Mechanics and Engineering, 2009,198(15-16):1360—1367.

    [8]Peng L X, Kitipornchai S, Liew K M. Free vibration analysis of folded plate structures by the FSDT mesh-free method[J]. Computational Mechanics, 2006,39(6):799—814.

    [9]Topal U, Uzman . Frequency optimization of laminated folded composite plates[J]. Materials & Design, 2009,30(3):494 —501.

    [10]Jiang R J, Au F T K. A general finite strip for the static and dynamic analyses of folded plates[J]. Thin-Walled Structures, 2011,49(10):1288—1294.

    [11]Liu C, Tian Q, Hu H. Dynamics of a large scale rigid-flexible multibody system composed of composite laminated plates[J]. Multibody System Dynamics, 2011,26(3):283—305.

    [12]Zhang W, Hu W H, Cao D X, et al. Vibration frequencies and modes of a Z-shaped beam with variable folding angles[J]. Journal of Vibration and Acoustics, 2016,138(4):041004.

    [13]Guo X Y, Zhang W, Yao M H. Nonlinear dynamics of angle-ply composite laminated thin plate with third-order shear deformation[J]. Science China Technological Sciences, 2010,53(3):612—622.

    [14]Zhang W, Guo X Y, Lai S K. Research on periodic and chaotic oscillations of composite laminated plates with one-to-one internal resonance[J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2009,10(11-12):1567—1583.

    [15]趙孟良, 關(guān)富玲, 吳開(kāi)成. 空間可展板殼結(jié)構(gòu)的展開(kāi)分析[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版), 2006,40(11):1837 —1841.

    Zhao Mengliang, Guang Fulin, Wu Kaicheng. Deployment analysis of deployable space panel and shell structure[J]. Journal of Zhejiang University(Engineering Science), 2006,40(11):1837—1841.

    [16]關(guān)富玲, 張惠峰, 韓克良. 二維可展板殼結(jié)構(gòu)展開(kāi)過(guò)程分析[J]. 工程設(shè)計(jì)學(xué)報(bào), 2008,15(5):351—356.

    Guan Fuling, Zhang Huifeng, Han Keliang. Deployment analysis of two-dimensional deployable panel and shell structures[J]. Journal of Engineering Design, 2008,15(5):351—356.

    [17]Reddy J N. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis[M]. CRC Press, 2004.

    [18]樸金麗.Z 型折疊機(jī)翼的非線性動(dòng)力學(xué)研究[D].北京:北京工業(yè)大學(xué),2016:05.

    Piao Jinli. Nonlinear vibrations for the Z-type folding wings of the morphing aircraft[D]. Beijing: Beijing University of Technology, 2016:05.

    Nonlinear vibration characteristics of Z-type folding plates

    with internal resonance

    GUO Xiang-ying, ZHANG Yang, ZHANG Wei

    (Beijing Key Laboratory of Nonlinear Vibrations and Strength of Mechanical Structures,

    Beijing University of Technology, Beijing 100124, China)

    Abstract: The Z-type folding plate is a kind of complex multi-body structure, which is widely applied to many engineering fields. Based on the classical laminated plate theory and the von Karman type equation, the nonlinear dynamic equations of the Z-type folding plates are obtained by using the Hamilton′s principle. The mode functions of the Z-type folding plates are analyzed with the ANSYS. Then, the Galerk in procedure is used to obtain the normal differential governing equations of the nonlinear system. The case of primary parametric resonance 1∶2 inner resonance is considered. Based on the averaged equation obtained with the method of multiple scales, the numerical simulation is performed to indicate the nonlinear dynamical characteristics of the system.

    Keywords: nonlinear dynamics; Z-type folding plates; inner resonance; numerical simulation

    猜你喜歡
    數(shù)值模擬
    基于AMI的雙色注射成型模擬分析
    錐齒輪精密冷擺輾成形在“材料成型數(shù)值模擬”課程教學(xué)中的應(yīng)用
    西南地區(qū)氣象資料測(cè)試、預(yù)處理和加工研究報(bào)告
    張家灣煤礦巷道無(wú)支護(hù)條件下位移的數(shù)值模擬
    張家灣煤礦開(kāi)切眼錨桿支護(hù)參數(shù)確定的數(shù)值模擬
    跨音速飛行中機(jī)翼水汽凝結(jié)的數(shù)值模擬研究
    雙螺桿膨脹機(jī)的流場(chǎng)數(shù)值模擬研究
    一種基于液壓緩沖的減震管卡設(shè)計(jì)與性能分析
    蒸汽發(fā)生器一次側(cè)流阻數(shù)值模擬研究
    亚洲激情五月婷婷啪啪| 99久久无色码亚洲精品果冻| 尾随美女入室| 99久国产av精品国产电影| 国产男人的电影天堂91| 亚洲成人av在线免费| 久久人妻av系列| 日本一二三区视频观看| 久久久成人免费电影| 国产高清不卡午夜福利| 成人综合一区亚洲| 午夜福利视频1000在线观看| 亚州av有码| 大话2 男鬼变身卡| 麻豆成人av视频| 女的被弄到高潮叫床怎么办| 国产免费视频播放在线视频 | 色综合色国产| 中文字幕亚洲精品专区| 欧美成人一区二区免费高清观看| 国产成人91sexporn| 在线观看av片永久免费下载| 欧美日本视频| 免费黄网站久久成人精品| 久久国内精品自在自线图片| 国国产精品蜜臀av免费| 精华霜和精华液先用哪个| 国产高清不卡午夜福利| 青春草视频在线免费观看| 亚洲国产精品成人久久小说| 我要看日韩黄色一级片| 91av网一区二区| 国产精品久久久久久久电影| 久久这里有精品视频免费| 国产精品一二三区在线看| 我要搜黄色片| 亚洲av成人精品一二三区| 人妻夜夜爽99麻豆av| 免费看光身美女| 波多野结衣巨乳人妻| 久热久热在线精品观看| 一区二区三区免费毛片| 国产成人福利小说| 国内揄拍国产精品人妻在线| 国产熟女欧美一区二区| 可以在线观看毛片的网站| 国产精品伦人一区二区| 亚洲va在线va天堂va国产| 欧美性猛交黑人性爽| 久久亚洲精品不卡| 永久免费av网站大全| 欧美日韩国产亚洲二区| 99视频精品全部免费 在线| 久久鲁丝午夜福利片| 久久久精品大字幕| 久久久国产成人精品二区| 中文乱码字字幕精品一区二区三区 | 日韩一区二区三区影片| 亚洲欧美日韩高清专用| 一个人免费在线观看电影| 国产精品,欧美在线| 国产精品乱码一区二三区的特点| 最近手机中文字幕大全| 亚洲真实伦在线观看| 99久久成人亚洲精品观看| 色哟哟·www| 精品久久久噜噜| 日韩成人av中文字幕在线观看| 国产免费一级a男人的天堂| 99久久人妻综合| 天美传媒精品一区二区| 欧美精品一区二区大全| 欧美区成人在线视频| 九草在线视频观看| 亚洲欧美成人精品一区二区| 免费人成在线观看视频色| 精品99又大又爽又粗少妇毛片| 免费在线观看成人毛片| 国产亚洲av嫩草精品影院| 亚洲av熟女| 18禁在线播放成人免费| 日本wwww免费看| 国产精品无大码| 嫩草影院入口| 亚洲乱码一区二区免费版| 麻豆精品久久久久久蜜桃| 天堂√8在线中文| 国产av在哪里看| 日本免费在线观看一区| 日本熟妇午夜| 日本一本二区三区精品| 99热这里只有精品一区| 国产精品嫩草影院av在线观看| 校园人妻丝袜中文字幕| 99久久精品热视频| 亚洲国产日韩欧美精品在线观看| 精品99又大又爽又粗少妇毛片| 一本久久精品| 免费大片18禁| 亚洲av免费在线观看| 久久久午夜欧美精品| 三级经典国产精品| 国产 一区 欧美 日韩| 国产精品国产三级国产专区5o | 精品久久久噜噜| a级毛色黄片| av国产久精品久网站免费入址| 日韩三级伦理在线观看| 国产真实伦视频高清在线观看| 小说图片视频综合网站| 欧美97在线视频| 久久久久性生活片| 国产精品熟女久久久久浪| 国产高清视频在线观看网站| 日本午夜av视频| 波多野结衣高清无吗| 日本-黄色视频高清免费观看| 精品少妇黑人巨大在线播放 | 永久免费av网站大全| 色网站视频免费| 熟妇人妻久久中文字幕3abv| 国产av一区在线观看免费| 麻豆成人午夜福利视频| 国语对白做爰xxxⅹ性视频网站| 久久久精品欧美日韩精品| 亚州av有码| 日韩中字成人| 久久欧美精品欧美久久欧美| 日韩,欧美,国产一区二区三区 | 我的老师免费观看完整版| 久久久久国产网址| 久久久久久久久中文| 岛国毛片在线播放| 天堂av国产一区二区熟女人妻| 又爽又黄无遮挡网站| 插阴视频在线观看视频| 我要看日韩黄色一级片| 亚洲内射少妇av| 国产成年人精品一区二区| 九九久久精品国产亚洲av麻豆| 久久精品人妻少妇| 身体一侧抽搐| 插逼视频在线观看| 天天一区二区日本电影三级| 日韩成人伦理影院| 欧美潮喷喷水| 观看免费一级毛片| 国产日韩欧美在线精品| 欧美一级a爱片免费观看看| 久久久久性生活片| 在现免费观看毛片| 国产精华一区二区三区| 欧美xxxx黑人xx丫x性爽| 久久精品久久精品一区二区三区| 亚洲av电影在线观看一区二区三区 | 久久久成人免费电影| 午夜a级毛片| 2021天堂中文幕一二区在线观| 在线观看av片永久免费下载| 边亲边吃奶的免费视频| 人妻夜夜爽99麻豆av| 亚洲精品aⅴ在线观看| 中文天堂在线官网| 在线观看av片永久免费下载| 国产高潮美女av| av在线蜜桃| 汤姆久久久久久久影院中文字幕 | 久久99热6这里只有精品| 成年免费大片在线观看| 国产v大片淫在线免费观看| 嘟嘟电影网在线观看| 99热6这里只有精品| 日韩亚洲欧美综合| 国产片特级美女逼逼视频| 亚洲国产精品成人久久小说| 简卡轻食公司| 国产精品国产三级国产专区5o | 亚洲精品自拍成人| 日本五十路高清| 日产精品乱码卡一卡2卡三| 久久99精品国语久久久| 国产三级中文精品| 久久久久久九九精品二区国产| 熟女电影av网| 少妇高潮的动态图| 成人高潮视频无遮挡免费网站| 成人午夜精彩视频在线观看| 联通29元200g的流量卡| 99热精品在线国产| 国产精品电影一区二区三区| 真实男女啪啪啪动态图| 欧美成人精品欧美一级黄| 中文字幕制服av| 亚洲精品456在线播放app| 又爽又黄无遮挡网站| 欧美成人免费av一区二区三区| 国产淫片久久久久久久久| 免费观看的影片在线观看| 亚洲精品国产av成人精品| 亚洲国产精品成人久久小说| 精品一区二区三区视频在线| 欧美丝袜亚洲另类| 国产精品蜜桃在线观看| 国产日韩欧美在线精品| 亚洲av成人精品一区久久| 精品久久久噜噜| 看片在线看免费视频| 在线播放无遮挡| 搞女人的毛片| 国产成人一区二区在线| 麻豆成人av视频| 国产白丝娇喘喷水9色精品| 免费无遮挡裸体视频| 一个人看视频在线观看www免费| 一个人看的www免费观看视频| 国内精品一区二区在线观看| 国语自产精品视频在线第100页| 色综合色国产| 亚洲电影在线观看av| 三级国产精品欧美在线观看| 2021天堂中文幕一二区在线观| 亚洲欧美成人综合另类久久久 | 欧美xxxx黑人xx丫x性爽| 18禁在线无遮挡免费观看视频| 一本一本综合久久| 久久精品久久久久久久性| 蜜桃亚洲精品一区二区三区| 97超碰精品成人国产| 色网站视频免费| 大话2 男鬼变身卡| 亚洲精品,欧美精品| 国产精品国产三级专区第一集| 亚洲欧洲国产日韩| 日韩三级伦理在线观看| 听说在线观看完整版免费高清| 99久久无色码亚洲精品果冻| 春色校园在线视频观看| 精品人妻熟女av久视频| 免费在线观看成人毛片| 18禁动态无遮挡网站| 久久久欧美国产精品| www.av在线官网国产| 国产精品一区二区性色av| 91aial.com中文字幕在线观看| 免费av毛片视频| 欧美一区二区亚洲| 欧美xxxx黑人xx丫x性爽| 视频中文字幕在线观看| 国产单亲对白刺激| 亚洲性久久影院| 亚洲欧美日韩无卡精品| 免费大片18禁| 亚洲欧美精品综合久久99| 国产精品国产高清国产av| 国产不卡一卡二| 99久久无色码亚洲精品果冻| 特大巨黑吊av在线直播| 岛国毛片在线播放| 久久这里只有精品中国| 日韩高清综合在线| 国产综合懂色| 亚洲欧美精品自产自拍| 国产亚洲91精品色在线| 1024手机看黄色片| 如何舔出高潮| 国产极品天堂在线| 亚洲欧美中文字幕日韩二区| 亚洲中文字幕一区二区三区有码在线看| 欧美zozozo另类| 国产精品人妻久久久久久| 亚洲精品乱久久久久久| 级片在线观看| 九九久久精品国产亚洲av麻豆| 欧美3d第一页| 亚洲精品久久久久久婷婷小说 | 国产精品一区二区三区四区免费观看| 一级毛片久久久久久久久女| 日韩视频在线欧美| 亚洲自偷自拍三级| 真实男女啪啪啪动态图| 高清日韩中文字幕在线| 久久久久久久国产电影| 亚洲欧美日韩无卡精品| 三级国产精品片| 永久免费av网站大全| 国产精品久久久久久av不卡| 少妇高潮的动态图| 岛国在线免费视频观看| 国产成人午夜福利电影在线观看| 高清av免费在线| 国产在视频线在精品| 日韩欧美精品v在线| 成人av在线播放网站| 色网站视频免费| 成人午夜高清在线视频| 国产精品乱码一区二三区的特点| 欧美性猛交黑人性爽| 久久久久久久亚洲中文字幕| 日韩欧美三级三区| 国产 一区 欧美 日韩| 男人舔女人下体高潮全视频| 99热这里只有是精品在线观看| 中文天堂在线官网| 国产真实乱freesex| 免费观看在线日韩| 毛片女人毛片| 久久精品91蜜桃| av视频在线观看入口| 亚洲成av人片在线播放无| 久久人人爽人人爽人人片va| 中文字幕久久专区| 国产成年人精品一区二区| 久久久久久久久久久丰满| 精品久久久久久久久亚洲| 麻豆精品久久久久久蜜桃| 日本免费一区二区三区高清不卡| 亚洲成人久久爱视频| 身体一侧抽搐| 男人的好看免费观看在线视频| 色综合色国产| 一个人观看的视频www高清免费观看| 熟妇人妻久久中文字幕3abv| 偷拍熟女少妇极品色| 小蜜桃在线观看免费完整版高清| 久久久精品94久久精品| 国产淫语在线视频| 变态另类丝袜制服| 久久这里有精品视频免费| 亚洲18禁久久av| 中文在线观看免费www的网站| 少妇被粗大猛烈的视频| 女的被弄到高潮叫床怎么办| eeuss影院久久| 免费看av在线观看网站| 国产精品福利在线免费观看| 超碰97精品在线观看| 国产极品天堂在线| 久久精品久久久久久久性| 欧美又色又爽又黄视频| 国产中年淑女户外野战色| 欧美一区二区精品小视频在线| 国产精品嫩草影院av在线观看| 日本免费在线观看一区| 国产精品熟女久久久久浪| 在线天堂最新版资源| 国模一区二区三区四区视频| 精品国产一区二区三区久久久樱花 | 亚洲国产精品成人久久小说| 男的添女的下面高潮视频| 国产老妇伦熟女老妇高清| 在线天堂最新版资源| 亚洲国产最新在线播放| 久久久久久久久久久免费av| 18禁在线无遮挡免费观看视频| 看十八女毛片水多多多| 成人午夜高清在线视频| 亚洲国产精品合色在线| av卡一久久| 亚洲色图av天堂| 69人妻影院| 国产精品乱码一区二三区的特点| 一二三四中文在线观看免费高清| 黄色日韩在线| 性色avwww在线观看| 国产老妇女一区| av国产久精品久网站免费入址| 亚洲成色77777| 嫩草影院入口| 中文亚洲av片在线观看爽| 国语对白做爰xxxⅹ性视频网站| 精品不卡国产一区二区三区| 插阴视频在线观看视频| 大话2 男鬼变身卡| 亚洲国产色片| 麻豆久久精品国产亚洲av| 一个人观看的视频www高清免费观看| 亚洲精华国产精华液的使用体验| 免费大片18禁| 一个人看视频在线观看www免费| 天美传媒精品一区二区| 黄色一级大片看看| 好男人视频免费观看在线| 欧美不卡视频在线免费观看| 欧美xxxx黑人xx丫x性爽| 国产高清不卡午夜福利| 亚洲最大成人手机在线| 国产老妇女一区| videossex国产| 九九在线视频观看精品| 老司机影院成人| 边亲边吃奶的免费视频| 免费在线观看成人毛片| 日日啪夜夜撸| 狂野欧美白嫩少妇大欣赏| 22中文网久久字幕| 桃色一区二区三区在线观看| 免费观看a级毛片全部| 老司机福利观看| 久久久久精品久久久久真实原创| 中文字幕人妻熟人妻熟丝袜美| 精品国产三级普通话版| 高清在线视频一区二区三区 | 一级av片app| 五月伊人婷婷丁香| 亚洲精品一区蜜桃| 国产亚洲91精品色在线| 国产成人午夜福利电影在线观看| 日本免费a在线| 欧美激情久久久久久爽电影| 能在线免费看毛片的网站| 夫妻性生交免费视频一级片| 一本一本综合久久| 午夜日本视频在线| 国产单亲对白刺激| 日韩制服骚丝袜av| 日本一本二区三区精品| 久久久久久国产a免费观看| kizo精华| av视频在线观看入口| 亚洲经典国产精华液单| 亚洲四区av| 91久久精品电影网| 啦啦啦啦在线视频资源| 亚洲欧洲国产日韩| 国产色爽女视频免费观看| 国产av不卡久久| 赤兔流量卡办理| 我要看日韩黄色一级片| 国产精品一区二区三区四区免费观看| 插阴视频在线观看视频| 久久久久久九九精品二区国产| 男人和女人高潮做爰伦理| 国产女主播在线喷水免费视频网站 | 人妻系列 视频| 91av网一区二区| 午夜精品国产一区二区电影 | 在线a可以看的网站| 久久国内精品自在自线图片| 91午夜精品亚洲一区二区三区| 成年免费大片在线观看| 亚洲欧美成人综合另类久久久 | 久久久精品欧美日韩精品| 精品少妇黑人巨大在线播放 | 91久久精品电影网| 欧美3d第一页| 国内揄拍国产精品人妻在线| 国产成人a区在线观看| 久久精品影院6| 日韩在线高清观看一区二区三区| 69av精品久久久久久| av福利片在线观看| 好男人在线观看高清免费视频| 99视频精品全部免费 在线| 老司机影院成人| 神马国产精品三级电影在线观看| 国产v大片淫在线免费观看| 精品久久久久久电影网 | 国产精品日韩av在线免费观看| 久久久久免费精品人妻一区二区| av卡一久久| 国产成人精品久久久久久| 久久精品夜色国产| 久久精品国产自在天天线| 中文在线观看免费www的网站| 天堂√8在线中文| 国产亚洲精品av在线| 69av精品久久久久久| 欧美日韩一区二区视频在线观看视频在线 | 日韩人妻高清精品专区| 国产视频内射| 久久精品夜夜夜夜夜久久蜜豆| 亚洲天堂国产精品一区在线| 丰满乱子伦码专区| 伦理电影大哥的女人| 毛片女人毛片| 欧美日本亚洲视频在线播放| 国产在线一区二区三区精 | 欧美丝袜亚洲另类| 欧美日韩精品成人综合77777| 午夜久久久久精精品| 久久久久久国产a免费观看| 久久久久免费精品人妻一区二区| 午夜福利在线观看吧| 91久久精品电影网| 欧美成人免费av一区二区三区| 亚洲一级一片aⅴ在线观看| av专区在线播放| 免费大片18禁| 国产熟女欧美一区二区| 中文字幕精品亚洲无线码一区| 久久精品国产自在天天线| 国产精品麻豆人妻色哟哟久久 | 少妇熟女欧美另类| 亚洲美女视频黄频| 青春草视频在线免费观看| 我要搜黄色片| 夜夜看夜夜爽夜夜摸| 少妇熟女欧美另类| 久久99热这里只有精品18| 秋霞伦理黄片| 一二三四中文在线观看免费高清| 我的老师免费观看完整版| 啦啦啦啦在线视频资源| 国产私拍福利视频在线观看| 在线播放无遮挡| 美女内射精品一级片tv| 一区二区三区高清视频在线| 欧美成人精品欧美一级黄| 亚州av有码| 亚洲无线观看免费| 亚洲精品国产av成人精品| 日日啪夜夜撸| 免费看av在线观看网站| 久久久久久久久久黄片| av在线亚洲专区| 日韩高清综合在线| 只有这里有精品99| 乱码一卡2卡4卡精品| 国产男人的电影天堂91| 两性午夜刺激爽爽歪歪视频在线观看| 大香蕉97超碰在线| 亚洲国产精品sss在线观看| 真实男女啪啪啪动态图| 久久综合国产亚洲精品| 亚洲国产高清在线一区二区三| a级毛色黄片| 七月丁香在线播放| 亚洲无线观看免费| 99久久精品国产国产毛片| 午夜福利在线在线| 老师上课跳d突然被开到最大视频| 国产成人精品婷婷| 国产探花极品一区二区| 亚洲欧美日韩高清专用| 日韩中字成人| 久久精品熟女亚洲av麻豆精品 | 大香蕉97超碰在线| 国产 一区精品| 国产成人精品一,二区| 日本免费a在线| 亚洲中文字幕一区二区三区有码在线看| 免费看av在线观看网站| 亚洲国产欧洲综合997久久,| 国内精品一区二区在线观看| 成人性生交大片免费视频hd| 看非洲黑人一级黄片| 日韩视频在线欧美| 欧美一区二区亚洲| 国产精品久久久久久久电影| 超碰97精品在线观看| 亚洲自偷自拍三级| 人妻制服诱惑在线中文字幕| 人妻少妇偷人精品九色| 亚洲美女视频黄频| 亚洲精华国产精华液的使用体验| videossex国产| 国产免费一级a男人的天堂| 嘟嘟电影网在线观看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 国产精华一区二区三区| 美女大奶头视频| 成人美女网站在线观看视频| 最近中文字幕高清免费大全6| 人妻夜夜爽99麻豆av| 午夜精品国产一区二区电影 | 波多野结衣高清无吗| 国产在视频线在精品| 欧美一区二区亚洲| 亚洲av日韩在线播放| 成人亚洲精品av一区二区| 国产精品日韩av在线免费观看| 欧美最新免费一区二区三区| 国产精品一二三区在线看| 国产欧美另类精品又又久久亚洲欧美| 亚洲最大成人中文| 性色avwww在线观看| 国产一区二区三区av在线| 日韩视频在线欧美| 亚洲国产精品成人综合色| 22中文网久久字幕| 亚洲精品乱码久久久v下载方式| 噜噜噜噜噜久久久久久91| 日本免费在线观看一区| 18禁动态无遮挡网站| 久久久久久久午夜电影| 九九在线视频观看精品| 久久精品国产99精品国产亚洲性色| 高清av免费在线| 欧美最新免费一区二区三区| 成年免费大片在线观看| 长腿黑丝高跟| 国产精品.久久久| 免费黄色在线免费观看| 久久久久久久久久黄片| 黄色一级大片看看| 少妇熟女欧美另类| 精品午夜福利在线看| 天天一区二区日本电影三级| 国产黄色视频一区二区在线观看 | 国产精品野战在线观看| 一边摸一边抽搐一进一小说| 精品人妻视频免费看| 国产淫语在线视频| 亚洲国产最新在线播放| av免费观看日本| 国产午夜精品一二区理论片| 免费观看的影片在线观看| 最新中文字幕久久久久| 亚洲国产日韩欧美精品在线观看| 老司机影院成人| 亚洲伊人久久精品综合 | 国产亚洲午夜精品一区二区久久 | 亚洲国产精品合色在线| 精品无人区乱码1区二区| 国产高清视频在线观看网站|