• <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国产精品免费福利视频| 久久97久久精品| 欧美精品人与动牲交sv欧美| 老女人水多毛片| 天美传媒精品一区二区| 欧美精品一区二区大全| 少妇精品久久久久久久| 天堂8中文在线网| 久久99一区二区三区| 免费黄网站久久成人精品| 99久国产av精品国产电影| 成人漫画全彩无遮挡| 久久久久精品久久久久真实原创| 国产毛片在线视频| 一级av片app| 亚洲精品国产av成人精品| 桃花免费在线播放| 中文字幕人妻丝袜制服| 午夜视频国产福利| 午夜激情福利司机影院| 丰满迷人的少妇在线观看| 久久国产乱子免费精品| 熟女电影av网| 久久婷婷青草| 欧美xxxx性猛交bbbb| 国产精品免费大片| 久久久久久久久久久久大奶| 亚洲精品一区蜜桃| 亚洲精品自拍成人| 国产在线男女| 国产精品一二三区在线看| 国产高清有码在线观看视频| 亚洲国产av新网站| 免费观看av网站的网址| 亚洲国产毛片av蜜桃av| av不卡在线播放| 中国三级夫妇交换| 99视频精品全部免费 在线| 人妻 亚洲 视频| 曰老女人黄片| 国产精品秋霞免费鲁丝片| 欧美激情国产日韩精品一区| 久久久久精品性色| 精品久久久久久电影网| 我的老师免费观看完整版| 免费少妇av软件| av卡一久久| 免费观看a级毛片全部| 男女无遮挡免费网站观看| 我要看黄色一级片免费的| 97超视频在线观看视频| 国产日韩欧美在线精品| 自线自在国产av| 久久久久国产精品人妻一区二区| 九色成人免费人妻av| 亚洲精品aⅴ在线观看| 午夜免费观看性视频| 91精品国产九色| 亚洲欧美清纯卡通| 黄色毛片三级朝国网站 | 纯流量卡能插随身wifi吗| 女性被躁到高潮视频| 在线观看美女被高潮喷水网站| 亚洲精品色激情综合| 日产精品乱码卡一卡2卡三| 美女脱内裤让男人舔精品视频| 久久毛片免费看一区二区三区| 伊人久久国产一区二区| av女优亚洲男人天堂| 国产av国产精品国产| 国产亚洲午夜精品一区二区久久| 日韩精品有码人妻一区| 日日啪夜夜撸| 成人亚洲欧美一区二区av| 91精品一卡2卡3卡4卡| tube8黄色片| 99热这里只有是精品在线观看| 99久久精品国产国产毛片| 久久ye,这里只有精品| 国产精品国产三级国产av玫瑰| 91aial.com中文字幕在线观看| 国产女主播在线喷水免费视频网站| 伊人久久精品亚洲午夜| 久久久精品免费免费高清| 一个人免费看片子| 精品99又大又爽又粗少妇毛片| 国产黄片视频在线免费观看| 亚洲欧美一区二区三区黑人 | av专区在线播放| 久久韩国三级中文字幕| 成人影院久久| 国产精品成人在线| 制服丝袜香蕉在线| 日韩电影二区| 新久久久久国产一级毛片| 亚洲av综合色区一区| 麻豆成人午夜福利视频| 日日摸夜夜添夜夜爱| 大码成人一级视频| 蜜臀久久99精品久久宅男| 精品熟女少妇av免费看| 十八禁网站网址无遮挡 | 亚洲激情五月婷婷啪啪| 国产精品99久久久久久久久| 伊人久久国产一区二区| 久久国内精品自在自线图片| 国产成人午夜福利电影在线观看| a级毛片免费高清观看在线播放| 欧美日韩视频精品一区| av天堂中文字幕网| 免费少妇av软件| 中文字幕制服av| 国产黄色视频一区二区在线观看| 成人无遮挡网站| 91久久精品国产一区二区三区| 七月丁香在线播放| 少妇被粗大猛烈的视频| 国产精品无大码| 一区二区三区精品91| 亚洲第一区二区三区不卡| 九九在线视频观看精品| 国产在线男女| 日韩一区二区视频免费看| 久久久久久久亚洲中文字幕| 涩涩av久久男人的天堂| 中文字幕亚洲精品专区| 国产日韩欧美视频二区| 五月天丁香电影| 插阴视频在线观看视频| 麻豆精品久久久久久蜜桃| 国产精品久久久久久av不卡| 多毛熟女@视频| 高清av免费在线| 爱豆传媒免费全集在线观看| 成人国产av品久久久| 亚洲美女黄色视频免费看| 岛国毛片在线播放| 国产精品国产av在线观看| tube8黄色片| 黑人巨大精品欧美一区二区蜜桃 | 欧美变态另类bdsm刘玥| 高清在线视频一区二区三区| 国产精品国产三级国产专区5o| 久久 成人 亚洲| av一本久久久久| 久久人人爽人人爽人人片va| 久久久久久久久大av| 亚洲欧美精品专区久久| 亚洲精品日韩在线中文字幕| 亚洲欧美一区二区三区黑人 | 国产av国产精品国产| 成人黄色视频免费在线看| 91午夜精品亚洲一区二区三区| 日本与韩国留学比较| 最新的欧美精品一区二区| 男女啪啪激烈高潮av片| 久久99一区二区三区| 人体艺术视频欧美日本| 黑丝袜美女国产一区| 欧美日韩综合久久久久久| 欧美激情极品国产一区二区三区 | 亚洲天堂av无毛| 中文欧美无线码| 91aial.com中文字幕在线观看| 国产精品福利在线免费观看| 男的添女的下面高潮视频| 欧美精品人与动牲交sv欧美| 久久久久久伊人网av| 丰满乱子伦码专区| 激情五月婷婷亚洲| 欧美国产精品一级二级三级 | 一级av片app| 另类精品久久| 久久精品国产亚洲网站| 久久久久国产精品人妻一区二区| 秋霞在线观看毛片| 97精品久久久久久久久久精品| 日韩亚洲欧美综合| 在现免费观看毛片| 99热6这里只有精品| 丰满饥渴人妻一区二区三| 嘟嘟电影网在线观看| 中文字幕人妻熟人妻熟丝袜美| 国产亚洲av片在线观看秒播厂| 久久久久久久亚洲中文字幕| 久久国产精品男人的天堂亚洲 | 久久午夜综合久久蜜桃| 有码 亚洲区| 人人澡人人妻人| av卡一久久| 大香蕉97超碰在线| a级毛色黄片| 亚洲美女搞黄在线观看| 中文字幕免费在线视频6| 高清在线视频一区二区三区| 男人添女人高潮全过程视频| 日韩视频在线欧美| a级片在线免费高清观看视频| 国产伦理片在线播放av一区| 久久亚洲国产成人精品v| 99re6热这里在线精品视频| 女人久久www免费人成看片| 国产成人精品久久久久久| 亚洲一区二区三区欧美精品| 自拍偷自拍亚洲精品老妇| 国产一区二区在线观看日韩| 最近最新中文字幕免费大全7| 免费久久久久久久精品成人欧美视频 | 国产白丝娇喘喷水9色精品| 永久网站在线| 久久99一区二区三区| 十八禁网站网址无遮挡 | 精品国产一区二区三区久久久樱花| 人妻一区二区av| 男女无遮挡免费网站观看| 免费看日本二区| 国产视频内射| 免费看日本二区| 九九爱精品视频在线观看| 精品亚洲成a人片在线观看| 成人美女网站在线观看视频| 国产av一区二区精品久久| kizo精华| 女性生殖器流出的白浆| 国产男女超爽视频在线观看| 日本-黄色视频高清免费观看| 男的添女的下面高潮视频| 日产精品乱码卡一卡2卡三| 少妇被粗大猛烈的视频| 18禁裸乳无遮挡动漫免费视频| 一级片'在线观看视频| 狂野欧美激情性xxxx在线观看| 国产综合精华液| 久久99热这里只频精品6学生| 日韩欧美精品免费久久| 国产淫片久久久久久久久| 国产精品熟女久久久久浪| 99热这里只有是精品在线观看| 久久久a久久爽久久v久久| 国产精品人妻久久久久久| 亚洲性久久影院| 老司机影院毛片| av网站免费在线观看视频| 在线观看免费日韩欧美大片 | 黄色欧美视频在线观看| 亚洲精品一二三| 三级经典国产精品| 日韩人妻高清精品专区| 观看av在线不卡| 亚洲av成人精品一二三区| 高清av免费在线| 欧美性感艳星| 美女福利国产在线| 99久久精品国产国产毛片| 精品一区二区三区视频在线| 另类精品久久| 成人综合一区亚洲| 久久婷婷青草| 久久女婷五月综合色啪小说| 亚洲精品视频女| 亚洲天堂av无毛| 各种免费的搞黄视频| 国产高清国产精品国产三级| 嫩草影院入口| a级毛色黄片| av在线老鸭窝| 国产成人精品无人区| 久久午夜综合久久蜜桃| 在线 av 中文字幕| 超碰97精品在线观看| 黄色欧美视频在线观看| 大香蕉久久网| 汤姆久久久久久久影院中文字幕| 国产成人91sexporn| 亚洲精品国产色婷婷电影| 精品少妇内射三级| 韩国高清视频一区二区三区| 国产色婷婷99| av福利片在线观看| av播播在线观看一区| 亚洲av福利一区| 在线观看美女被高潮喷水网站| 国产老妇伦熟女老妇高清| 香蕉精品网在线| 国产午夜精品一二区理论片| 伦理电影大哥的女人| 日产精品乱码卡一卡2卡三| 国产精品国产三级国产av玫瑰| 亚洲,一卡二卡三卡| 下体分泌物呈黄色| 国产精品人妻久久久影院| 午夜免费鲁丝| 黄色视频在线播放观看不卡| 国产在线一区二区三区精| 97超碰精品成人国产| 国产有黄有色有爽视频| 只有这里有精品99| 精品久久久久久久久av| 亚洲va在线va天堂va国产| 亚洲国产欧美日韩在线播放 | 日本av免费视频播放| av免费在线看不卡| 毛片一级片免费看久久久久| 久久精品国产亚洲网站| 精品一区二区三卡| 色网站视频免费| 久久精品国产自在天天线| 日本猛色少妇xxxxx猛交久久| 亚洲国产精品一区三区| tube8黄色片| 成人毛片a级毛片在线播放| 制服丝袜香蕉在线| 观看美女的网站| 天堂8中文在线网| 国产成人a∨麻豆精品| 国产中年淑女户外野战色| 日日啪夜夜爽| 中国美白少妇内射xxxbb| 永久免费av网站大全| 看非洲黑人一级黄片| 一边亲一边摸免费视频| 免费观看性生交大片5| 久久久久精品性色| 狠狠精品人妻久久久久久综合| 亚洲丝袜综合中文字幕| 国产av码专区亚洲av| 精品久久久久久久久亚洲| 狂野欧美激情性xxxx在线观看| 亚洲国产色片| 国产色爽女视频免费观看| 三级经典国产精品| 2022亚洲国产成人精品| 有码 亚洲区| 亚洲精品国产av成人精品| 亚洲av综合色区一区| 男女无遮挡免费网站观看| 国产av码专区亚洲av| 国内揄拍国产精品人妻在线| 99九九线精品视频在线观看视频| 婷婷色综合www| av线在线观看网站| 看免费成人av毛片| 久久久久久久久久久久大奶| 日韩三级伦理在线观看| 免费观看无遮挡的男女| 超碰97精品在线观看| 人人妻人人爽人人添夜夜欢视频 | 草草在线视频免费看| 久久青草综合色| 国产精品伦人一区二区| 日韩成人伦理影院| 汤姆久久久久久久影院中文字幕| 日本猛色少妇xxxxx猛交久久| 国产精品熟女久久久久浪| 中文在线观看免费www的网站| 黄色日韩在线| 亚洲欧洲国产日韩| av在线播放精品| 嘟嘟电影网在线观看| 在线精品无人区一区二区三| 亚洲真实伦在线观看| 亚洲精品国产成人久久av| 狂野欧美激情性bbbbbb| 免费人成在线观看视频色| 国产免费福利视频在线观看| 爱豆传媒免费全集在线观看| 日韩欧美精品免费久久| 国产深夜福利视频在线观看| 久久久久久久久久久免费av| 最黄视频免费看| 一区二区三区免费毛片| 老熟女久久久| 最近中文字幕2019免费版| 久久久久视频综合| 欧美 亚洲 国产 日韩一| 亚洲av成人精品一区久久| 亚洲第一区二区三区不卡| 黑人高潮一二区| 一区二区三区乱码不卡18| 十八禁高潮呻吟视频 | 成人18禁高潮啪啪吃奶动态图 | 美女xxoo啪啪120秒动态图| 亚洲精品456在线播放app| 免费大片黄手机在线观看| 亚洲欧美成人综合另类久久久| 91久久精品国产一区二区成人| 十八禁高潮呻吟视频 | 久久精品久久精品一区二区三区| 最近的中文字幕免费完整| 成人午夜精彩视频在线观看| 日韩av不卡免费在线播放| 亚洲精品国产色婷婷电影| 日韩 亚洲 欧美在线| 蜜桃在线观看..| 国产精品福利在线免费观看| 综合色丁香网| 自线自在国产av| 最近中文字幕2019免费版| 亚洲av.av天堂| 亚洲四区av| 黄色欧美视频在线观看| 22中文网久久字幕| 精品一区在线观看国产| 国产视频内射| 麻豆成人av视频| 搡女人真爽免费视频火全软件| 午夜福利视频精品| 涩涩av久久男人的天堂| 男女免费视频国产| 26uuu在线亚洲综合色| 久久6这里有精品| av线在线观看网站| 欧美最新免费一区二区三区| 国产成人免费观看mmmm| 热99国产精品久久久久久7| 欧美97在线视频| 成人毛片a级毛片在线播放| 日本猛色少妇xxxxx猛交久久| 99热6这里只有精品| 黄色配什么色好看| 桃花免费在线播放| 久久精品夜色国产| 伦精品一区二区三区| 又大又黄又爽视频免费| 91精品国产国语对白视频| av在线老鸭窝| 精品人妻偷拍中文字幕| 免费大片18禁| 高清在线视频一区二区三区| 亚洲欧美一区二区三区黑人 | 欧美bdsm另类| 亚洲va在线va天堂va国产| 日本欧美视频一区| 免费高清在线观看视频在线观看| 黄色怎么调成土黄色| 国产欧美日韩一区二区三区在线 | 国产在线视频一区二区| 日韩制服骚丝袜av| 欧美 亚洲 国产 日韩一| 成年女人在线观看亚洲视频| 国产精品福利在线免费观看| 青青草视频在线视频观看| 老熟女久久久| 在线看a的网站| 9色porny在线观看| 免费看日本二区| av天堂中文字幕网| 国产男人的电影天堂91| 国产精品一二三区在线看| 国产精品一区二区在线观看99| av又黄又爽大尺度在线免费看| 秋霞在线观看毛片| 成人免费观看视频高清| 亚洲av福利一区| 久久人人爽人人片av| 亚洲在久久综合| 久久国产乱子免费精品| 九色成人免费人妻av| 99热这里只有是精品在线观看| 亚洲成色77777| 国产精品免费大片| 日韩欧美精品免费久久| 99热这里只有精品一区| 久久人人爽人人爽人人片va| 一区二区三区四区激情视频| 欧美精品高潮呻吟av久久| 亚洲精品视频女| 久久 成人 亚洲| 99久久精品一区二区三区| 国产欧美日韩综合在线一区二区 | 亚洲精品乱码久久久v下载方式| 我的老师免费观看完整版| 久久久久久久精品精品| 一级爰片在线观看| 中文欧美无线码| 2021少妇久久久久久久久久久| 日韩成人伦理影院| 丰满人妻一区二区三区视频av| 99热这里只有是精品在线观看| 老司机影院毛片| 成人18禁高潮啪啪吃奶动态图 | 啦啦啦视频在线资源免费观看| 在线观看一区二区三区激情| 亚洲天堂av无毛| 日韩亚洲欧美综合| 午夜福利,免费看| 亚州av有码| 亚洲精品国产av蜜桃| 欧美精品一区二区大全| 久久韩国三级中文字幕| 少妇的逼水好多| av免费在线看不卡| 男人狂女人下面高潮的视频| 日本-黄色视频高清免费观看| 亚洲精品久久午夜乱码| 日韩电影二区| 亚洲怡红院男人天堂| 麻豆乱淫一区二区| 高清午夜精品一区二区三区| 日韩成人av中文字幕在线观看| 亚洲国产毛片av蜜桃av| 国产熟女午夜一区二区三区 | 一本大道久久a久久精品| 亚洲久久久国产精品| 啦啦啦在线观看免费高清www| 成年av动漫网址| 久久久午夜欧美精品| 亚洲国产最新在线播放| 看十八女毛片水多多多| 激情五月婷婷亚洲| 久久精品国产亚洲av天美| av在线观看视频网站免费| 大香蕉97超碰在线| 看十八女毛片水多多多| 国产亚洲最大av| 老女人水多毛片| 日韩中字成人| 国产日韩欧美视频二区| 亚洲欧美一区二区三区黑人 | 多毛熟女@视频| 国产日韩欧美亚洲二区| 欧美bdsm另类| a级毛色黄片| 成人毛片a级毛片在线播放| 国模一区二区三区四区视频| 97超碰精品成人国产| 又粗又硬又长又爽又黄的视频| 久久久久久久久久久丰满| 深夜a级毛片| 99热6这里只有精品| 午夜福利在线观看免费完整高清在| 91久久精品国产一区二区成人| 免费看不卡的av| 777米奇影视久久| 国产深夜福利视频在线观看| 国产真实伦视频高清在线观看| av女优亚洲男人天堂| 肉色欧美久久久久久久蜜桃| 一二三四中文在线观看免费高清| 国产欧美日韩综合在线一区二区 | 免费看日本二区| 亚洲精品视频女| 久久国内精品自在自线图片| h日本视频在线播放| 青春草亚洲视频在线观看| 亚洲国产毛片av蜜桃av| 亚洲美女视频黄频| 国产亚洲一区二区精品| 男的添女的下面高潮视频| 国产精品久久久久成人av| 黑丝袜美女国产一区| 日本av手机在线免费观看| 日韩欧美一区视频在线观看 | 特大巨黑吊av在线直播| 22中文网久久字幕| 一本大道久久a久久精品| 老司机影院毛片| 久久韩国三级中文字幕| 国产伦在线观看视频一区| 亚洲欧洲精品一区二区精品久久久 | 国产亚洲一区二区精品| 大陆偷拍与自拍| 亚洲婷婷狠狠爱综合网| 人妻制服诱惑在线中文字幕| 在线观看美女被高潮喷水网站| 久久亚洲国产成人精品v| 国产男女内射视频| 精品亚洲成a人片在线观看| 哪个播放器可以免费观看大片| 韩国高清视频一区二区三区| 日韩人妻高清精品专区| 老司机亚洲免费影院| av天堂中文字幕网| 亚洲高清免费不卡视频| 欧美激情极品国产一区二区三区 | 久久久精品免费免费高清| 2018国产大陆天天弄谢| 欧美一级a爱片免费观看看| 各种免费的搞黄视频| av在线老鸭窝| 日韩三级伦理在线观看| av在线老鸭窝| 日韩中文字幕视频在线看片| av在线app专区| 十分钟在线观看高清视频www | 国产在线视频一区二区| 人人妻人人添人人爽欧美一区卜| 国产伦在线观看视频一区| 一级二级三级毛片免费看| 国产毛片在线视频| 久久人人爽av亚洲精品天堂| 国产伦在线观看视频一区| 综合色丁香网| 国产一区二区三区综合在线观看 | 熟女av电影| 美女xxoo啪啪120秒动态图| 日日摸夜夜添夜夜添av毛片| 一级二级三级毛片免费看| 高清视频免费观看一区二区| 少妇裸体淫交视频免费看高清| 99国产精品免费福利视频| 最近中文字幕高清免费大全6| 亚洲欧美日韩东京热| 91精品国产国语对白视频| 女性生殖器流出的白浆| 热re99久久国产66热| 人体艺术视频欧美日本| 免费观看的影片在线观看| 永久网站在线| 精品人妻熟女av久视频| 精品人妻偷拍中文字幕| 国产中年淑女户外野战色| 久久99蜜桃精品久久| 久热久热在线精品观看|