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

    平穩(wěn)/非平穩(wěn)激勵(lì)下中厚圓柱殼隨機(jī)振動(dòng)響應(yīng)的基準(zhǔn)解1)

    2022-04-07 06:56:42霍慧陳國(guó)海王文培楊迪雄
    力學(xué)學(xué)報(bào) 2022年3期
    關(guān)鍵詞:方根振型圓柱

    霍慧 陳國(guó)海 王文培 楊迪雄

    (大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧大連 116024)

    引言

    作為工程實(shí)際中常見的結(jié)構(gòu)形式,圓柱殼得益于其優(yōu)良的結(jié)構(gòu)特性和力學(xué)性能,被廣泛應(yīng)用于潛艇、火箭、導(dǎo)彈、飛機(jī)、儲(chǔ)液罐等工程結(jié)構(gòu)中[1-5].其中,考慮橫向剪切變形和轉(zhuǎn)動(dòng)慣量影響的中厚圓柱殼結(jié)構(gòu)不可避免地會(huì)承受各類隨機(jī)激勵(lì)作用,如隨機(jī)地震動(dòng)、風(fēng)載、噪聲激勵(lì)等[6-7].因此,中厚殼隨機(jī)振動(dòng)分析和不確定性傳播研究對(duì)其設(shè)計(jì)和安全服役具有重要的理論意義及應(yīng)用價(jià)值.

    自Love 建立殼體基本方程以來,已有許多學(xué)者和工程師對(duì)中厚圓柱殼結(jié)構(gòu)的自由振動(dòng)特性進(jìn)行了分析[8-9].Leissa[10]對(duì)殼體自由振動(dòng)理論進(jìn)行了深入廣泛的論述.然而,目前僅有簡(jiǎn)支邊圓柱殼能獲得精確解析的自振頻率和振型函數(shù).為此,各種數(shù)值解法被相繼提出,如動(dòng)力剛度法[11]、傅里葉級(jí)數(shù)法[12]、有限元法[13-14]、波函數(shù)法[15]等.

    實(shí)際工程中,殼體結(jié)構(gòu)經(jīng)常承受環(huán)境動(dòng)載荷作用.最近幾年,確定性動(dòng)載荷作用下圓柱殼的動(dòng)力響應(yīng)研究已經(jīng)得到了廣泛關(guān)注,載荷形式包括黏性流體[16]、移動(dòng)載荷[17]、溫度或濕度場(chǎng)[18]、磁場(chǎng)作用[19]等.由于服役環(huán)境的復(fù)雜性,動(dòng)力載荷普遍具有隨機(jī)性,通常采用隨機(jī)過程對(duì)其進(jìn)行表征.基于狀態(tài)協(xié)方差分配法,To和Chen[20]給出了非平穩(wěn)隨機(jī)激勵(lì)作用下壓電殼體元件的最優(yōu)隨機(jī)振動(dòng)控制方案.針對(duì)邊界隨機(jī)激勵(lì)作用下的壓電空心圓柱厚殼,Ying等[21]基于伽遼金法計(jì)算了隨機(jī)響應(yīng)的統(tǒng)計(jì)矩.Esmailzadeh和Lakis[22]采用有限元方法求解了湍流邊界層激勵(lì)作用下的開口薄殼振動(dòng)響應(yīng).基于FPK (Fokker-Planck-Kolmogorov)方程,Asnafi[23]獲得了隨機(jī)激勵(lì)作用下圓柱淺殼的概率密度函數(shù).Li 等[24]利用辛對(duì)偶方法研究了層合圓柱殼在軸向壓縮和湍流邊界層激勵(lì)下的隨機(jī)振動(dòng)響應(yīng).

    此外,Crandall和Elishkoff[25]以圓柱殼為研究對(duì)象,揭示了振型互相關(guān)項(xiàng)對(duì)隨機(jī)振動(dòng)響應(yīng)具有不可忽視的作用.振型互相關(guān)項(xiàng)的重要作用也在球殼[26]、多跨梁[27]中被深入考察.對(duì)于離散多自由度系統(tǒng),林家浩等[28]提出了高效的虛擬激勵(lì)法,該方法精確計(jì)入振型互相關(guān)項(xiàng).Chen 等[29]將虛擬激勵(lì)法推廣到連續(xù)體,獲得了薄殼平穩(wěn)隨機(jī)振動(dòng)響應(yīng)的解析精確解.相比于薄殼,中厚殼在實(shí)際工程中的應(yīng)用更加廣泛,還會(huì)受到非平穩(wěn)隨機(jī)激勵(lì)作用.因此,本文針對(duì)中厚圓柱殼,開展平穩(wěn)、非平穩(wěn)隨機(jī)振動(dòng)分析,以獲得其隨機(jī)響應(yīng)的精確基準(zhǔn)解.

    為了保證結(jié)果的準(zhǔn)確性,本文首先通過自由振動(dòng)分析,獲得了簡(jiǎn)支中厚圓柱殼的精確頻率和解析的振型函數(shù).考慮點(diǎn)、線、面隨機(jī)激勵(lì),構(gòu)造平穩(wěn)和非平穩(wěn)虛擬激勵(lì),結(jié)合振型疊加法,將平穩(wěn)和非平穩(wěn)隨機(jī)振動(dòng)分別轉(zhuǎn)化為簡(jiǎn)諧振動(dòng)和時(shí)程分析,推導(dǎo)了中厚圓柱殼隨機(jī)響應(yīng)功率譜的解析表達(dá)式.將功率譜在頻域上積分即可得到均方響應(yīng)及均方根.此外,在解析解的推導(dǎo)過程中還涉及對(duì)振型函數(shù)空間域積分和時(shí)間域的Duhamel 積分.基于符號(hào)運(yùn)算的解析法在中厚殼隨機(jī)振動(dòng)響應(yīng)解析求解中無法實(shí)現(xiàn)多空間點(diǎn)同時(shí)輸出,且計(jì)算難度和效率隨參振振型的增加而顯著增加.為此,進(jìn)一步提出了離散解析法.通過對(duì)空間域先解析積分后離散,頻域和時(shí)域數(shù)值積分,將解析運(yùn)算轉(zhuǎn)化為矩陣運(yùn)算,不僅顯著提高了計(jì)算效率,而且能批量獲得殼體隨機(jī)響應(yīng)的分布,便于研究參數(shù)變化對(duì)圓柱殼體隨機(jī)響應(yīng)全局影響.此外,基于精確基準(zhǔn)解,深入討論了隨機(jī)激勵(lì)類型及結(jié)構(gòu)參數(shù)對(duì)中厚殼振動(dòng)響應(yīng)的影響.

    1 中厚圓柱殼自由振動(dòng)精確解

    對(duì)于殼結(jié)構(gòu),隨著厚徑比(厚度h與曲率半徑R之比)的增大,剪切變形和轉(zhuǎn)動(dòng)慣量對(duì)力學(xué)行為的影響越來越大,Kirchhoff 薄殼理論不再適用.針對(duì)中等厚度圓柱殼(如厚徑比≥1/36~ 1/20),分析應(yīng)采用考慮橫向剪切變形和轉(zhuǎn)動(dòng)慣量影響的中厚殼理論.本文分別討論了如圖1(a)所示的封閉殼(φT=2π)和圖1(b)所示的開口角度φT的開口殼結(jié)構(gòu).R,L和h分別為殼的半徑、長(zhǎng)度和厚度,坐標(biāo)系(x,θ,z),x,θ,z分別代表軸向、環(huán)向和徑向.引入5 個(gè)廣義獨(dú)立位移描述殼體的中面變形:軸向位移u1、環(huán)向位移u2、徑向位移u3、軸向轉(zhuǎn)角u4及環(huán)向轉(zhuǎn)角u5.考慮橫向剪切變形和轉(zhuǎn)動(dòng)慣量的中厚圓柱殼隨機(jī)偏微分方程可以描述為[10]

    圖1 中厚圓柱殼幾何模型及載荷情況Fig.1 Geometric model and load cases for moderately thick cylindrical shell

    其中c為等效阻尼系數(shù),ρ為體密度,i=1,2,3 分別代表軸向、環(huán)向和徑向,j=4,5 分別代表軸向轉(zhuǎn)角和環(huán)向轉(zhuǎn)角方向,qi為圖1 所示施加在第i個(gè)方向上的隨機(jī)激勵(lì),L1~L5為微分算子

    其中,E和υ分別為材料彈性模量和泊松比,κ為剪切修正系數(shù),本文取為5/6[10].薄膜內(nèi)力Nx,Nθ,Nxθ和Nθx,彎矩Mx和Mθ,扭矩Mxθ和Mθx具體形式分別寫為[30]

    式中,C=Eh/(1-υ2)表示拉伸剛度,D=Eh3/[12(1-υ2)]為彎曲剛度.

    為了得到圓柱殼的自振頻率和振型,需進(jìn)行中厚圓柱殼的自由振動(dòng)分析.去掉式(1)中的阻尼力項(xiàng)和外力項(xiàng),即可得到中厚圓柱殼自由振動(dòng)的微分方程

    對(duì)于如圖1 (a)所示的簡(jiǎn)支封閉圓柱殼,其自由振動(dòng)精確解可設(shè)為分離變量形式[10]

    其中,Us,mn(x,θ)為與第s個(gè)廣義位移相應(yīng)的第mn階振型函數(shù),s=1,2,3,4,5 分別代表軸向、環(huán)向、徑向、軸向轉(zhuǎn)角及環(huán)向轉(zhuǎn)角方向,m和n分別為軸向和環(huán)向的半波數(shù),αm=mπ/L,A,B,F,M和N為待求常數(shù).將式(7)代入式(6)可以得到

    此外,由式(6)和式(7)可以得到矩陣形式的對(duì)稱齊次方程組

    式(9)中的各個(gè)系數(shù)分別為

    式中,Ω=ωR為待求的無量綱頻率,λ=mπR/L,k=h2/(12R2),Iκ=κ(1-υ)/2.為得到式(9)的非平凡解,令其系數(shù)矩陣行列式為0,則可以獲得關(guān)于Ω2的五次方的頻率方程

    其中,a1-a5為代數(shù)方程的系數(shù).求解該方程從而獲得殼結(jié)構(gòu)的精確頻率ω.將求得的自振頻率回代到式(9),就可以最終確定式(7)中振型函數(shù)的系數(shù)A,B,F,M和N之間的比值關(guān)系.

    類似地,對(duì)于如圖1 (b)所示的四邊簡(jiǎn)支的開口形式中厚圓柱殼,其自由振動(dòng)封閉形式解可以表達(dá)為如下形式

    其中βn=nπ/φT.采用與封閉圓柱殼類似的計(jì)算步驟,將式(12)代入式(6),得到開口中厚圓柱殼自由振動(dòng)的頻率方程,同樣可以獲得該結(jié)構(gòu)的精確頻率及解析振型函數(shù).

    2 中厚圓柱殼隨機(jī)振動(dòng)響應(yīng)的解析法

    由振型疊加法可知,軸向、環(huán)向及徑向上的線彈性中厚圓柱殼位移可以表達(dá)為

    其中ηmn(t)為第mn階振型的正則坐標(biāo).將式(13)代入式(1)可以得到

    對(duì)于i=1,2,3,分別在式(14)兩端乘上第lk階振型Ui,lk,然后將i=1,2,3 的3 個(gè)方程在圓柱殼曲面內(nèi)進(jìn)行積分,綜合考慮式(8)及振型正交化,最終得到一系列解耦的單自由度系統(tǒng)

    其中 ζmn=c/(2ρhωmn) 為第mn階的中厚殼模態(tài)阻尼比,

    針對(duì)隨機(jī)激勵(lì)作用下中厚殼的振動(dòng)分析,基于功率譜的傳統(tǒng)頻域分析方法采用維納-辛欽關(guān)系獲得響應(yīng)的功率譜密度.其中完全二次項(xiàng)組合(CQC)法能夠計(jì)及所有參振振型耦合項(xiàng),但向量乘法運(yùn)算量過大;而忽略振型耦合項(xiàng)的平方和開平方(SRSS)法計(jì)算誤差較大,且僅適用于阻尼比很小、參振頻率稀疏分布的均質(zhì)系統(tǒng).因此,本文利用虛擬激勵(lì)法和連續(xù)體結(jié)構(gòu)解析振型函數(shù),在保證結(jié)果準(zhǔn)確性的同時(shí)只需要計(jì)算一次向量乘法,計(jì)算效率得到顯著提高.

    基于虛擬激勵(lì)法[28],構(gòu)造虛擬激勵(lì)

    對(duì)于平穩(wěn)隨機(jī)激勵(lì)qi(x,θ,t)=Γi(x,θ)X(t),其中Γi(x,θ)為作用在第i個(gè)方向上的平穩(wěn)隨機(jī)過程X(t)的空間分布函數(shù),式(16)成為

    其中SXX(ω) 為X(t)的功率譜密度(PSD)函數(shù).

    將虛擬激勵(lì)式(17)代入式(15)可以得到

    其中

    因此,平穩(wěn)隨機(jī)振動(dòng)分析轉(zhuǎn)化為確定性的簡(jiǎn)諧振動(dòng)分析,從而得到虛擬正則坐標(biāo)

    除平穩(wěn)隨機(jī)激勵(lì)外,本文還考察統(tǒng)計(jì)信息隨時(shí)間變化的非平穩(wěn)隨機(jī)激勵(lì).為描述非平穩(wěn)隨機(jī)激勵(lì),Priestley 定義了具有明確物理意義的演變功率譜[31].若平穩(wěn)隨機(jī)過程X(t)表示為:非平穩(wěn)隨機(jī)過程Y(t)可進(jìn)一步描述為

    其中,f(ω,t) 為確定性慢變的非均勻調(diào)制函數(shù).Y(t)的方差可表達(dá)為

    進(jìn)而得到Y(jié)(t)的演變功率譜密度

    然而,通常難以計(jì)算式(23)所示的非平穩(wěn)激勵(lì)下結(jié)構(gòu)隨機(jī)振動(dòng)響應(yīng)的演變功率譜密度.當(dāng)激勵(lì)功率譜密度函數(shù)中的頻率特性不隨時(shí)間變化時(shí),將非均勻調(diào)制函數(shù)f(ω,t)簡(jiǎn)化近似為均勻調(diào)制函數(shù)f(t),此時(shí)的均勻調(diào)制非平穩(wěn)隨機(jī)過程可以寫為

    對(duì)于施加在中厚殼第i個(gè)方向的均勻調(diào)制非平穩(wěn)激勵(lì)qi(x,θ,t)=Γi(x,θ)Y(t)=Γi(x,θ)f(t)X(t),基于虛擬激勵(lì)法的非平穩(wěn)虛擬激勵(lì)表示為

    將式(25)描述的非平穩(wěn)虛擬激勵(lì)代入解耦的單自由度系統(tǒng)式(15),可得到

    由Duhamel 時(shí)域積分,得到虛擬正則坐標(biāo)

    其中hmn(t-τ) 為單位脈沖響應(yīng)函數(shù),即

    基于平穩(wěn)及非平穩(wěn)隨機(jī)激勵(lì)作用下得到的虛擬正則坐標(biāo)如式(20)和式(27)所示,進(jìn)而中厚圓柱殼第i個(gè)方向的虛擬位移表示為

    進(jìn)一步,可以基于虛擬激勵(lì)法得到各個(gè)方向位移、速度及加速度響應(yīng)功率譜密度函數(shù)如下

    其中,上標(biāo)“*”表示復(fù)共軛.

    此外,在虛擬位移響應(yīng)的基礎(chǔ)上,根據(jù)殼結(jié)構(gòu)的幾何關(guān)系及材料本構(gòu)關(guān)系,可以由式(4)和式(5)獲得虛擬薄膜內(nèi)力、彎矩及扭矩

    同樣地,基于虛擬激勵(lì)法,將虛擬內(nèi)力響應(yīng)的復(fù)共軛乘上自己本身,可以解析獲得中厚圓柱殼薄膜內(nèi)力及彎曲內(nèi)力等待求響應(yīng)的功率譜密度函數(shù).

    此節(jié)解析法推導(dǎo)過程沒有引入任何近似,解析地獲得了平穩(wěn)及非平穩(wěn)隨機(jī)激勵(lì)作用下中厚圓柱殼的各類響應(yīng)功率譜密度函數(shù)精確解.

    3 中厚圓柱殼隨機(jī)振動(dòng)響應(yīng)的離散解析法

    在上一節(jié)的解析推導(dǎo)中,盡管采用了虛擬激勵(lì)法,解析解難以滿足求解隨機(jī)振動(dòng)響應(yīng)多點(diǎn)輸出的需求,且計(jì)算難度和效率隨參振振型的增加而增長(zhǎng).為了在保證解析法空間域上精確性的前提下提高計(jì)算效率,本節(jié)基于空間域解析積分和微分運(yùn)算后離散化策略,提出了中厚圓柱殼隨機(jī)振動(dòng)響應(yīng)分析的離散解析法.

    針對(duì)解析法推導(dǎo)過程中涉及到的振型函數(shù)在空間域的積分和微分運(yùn)算(見式(15)、式(19)、式(31)和式(32)),離散解析法同樣采取了解析求解,即

    由于引入了圓柱殼的封閉振型函數(shù)Ui,mn,如式(7)和式(12)所示,上述積分和微分均可獲得解析解.相比功率譜分析有限元法,采用空間域積分和微分解析運(yùn)算的離散解析法能保證在空間域的精確性,具有數(shù)值方法無法比擬的優(yōu)勢(shì).

    此時(shí),圓柱殼的隨機(jī)偏微分方程便轉(zhuǎn)化為求解式(18)和式(26)所示的常微分方程.在離散化系統(tǒng)中,虛擬激勵(lì)法的核心是在響應(yīng)功率譜密度計(jì)算中將向量乘法減少到一次,進(jìn)而提高計(jì)算效率.而在解析的功率譜密度公式(30)中,只涉及符號(hào)運(yùn)算.

    為此,將中厚殼結(jié)構(gòu)空間域進(jìn)行離散,各點(diǎn)位置坐標(biāo)寫為矩陣形式:r=[r1r2···rNG],rα=(xα,θα)T,α=1,2,···,NG,NG為空間離散點(diǎn)總數(shù),則式(13)描述的中厚殼位移公式可寫為向量形式

    其中Ui(r)為第i個(gè)方向上的振型矩陣,矩陣維度為NG×NF,NF表示參振頻率階數(shù),i=1,2,3,η 代表NF維的正則坐標(biāo)列向量.

    作用在中厚殼上的虛擬激勵(lì)可以離散化為

    因此,對(duì)于平穩(wěn)隨機(jī)振動(dòng),第i個(gè)方向上中厚殼虛擬位移改寫為向量形式

    其中,P(ω)為與式(19)對(duì)應(yīng)的廣義激勵(lì)的振幅矩陣,H(ω)為頻率響應(yīng)函數(shù)矩陣.由于離散解析法中針對(duì)振型函數(shù)積分的解析求解(見式(33)中的I1和I2),可知P(ω)在幾何空間域內(nèi)不涉及離散化操作.將多個(gè)節(jié)點(diǎn)坐標(biāo)代入式(36)中振型Ui(r)中,使得離散解析法實(shí)現(xiàn)了多空間點(diǎn)虛擬響應(yīng)的同時(shí)輸出,且虛擬位移響應(yīng)精度與空間離散點(diǎn)多寡無關(guān).

    對(duì)于非平穩(wěn)隨機(jī)振動(dòng),由于Duhamel 積分解析求解通常耗時(shí)較長(zhǎng),所以采用高精度的精細(xì)積分法取代Duhamel 積分進(jìn)行時(shí)間歷程分析可提高計(jì)算效率.將解耦得到的單自由度系統(tǒng)式(26)改寫為狀態(tài)空間形式

    式(37)求解過程中需要對(duì)時(shí)域離散化處理,對(duì)任意t∈[tk,tk + 1],采用精細(xì)積分法對(duì)式(37)進(jìn)行精確高效計(jì)算[32],其中tk=kΔt,k=0,1,2,···,Nt,Δt=T/(Nt),T為激勵(lì)持續(xù)時(shí)間,Nt為離散時(shí)間點(diǎn)數(shù).從而最終得到中厚圓柱殼在非平穩(wěn)隨機(jī)激勵(lì)作用下的虛擬位移響應(yīng)向量注意到,針對(duì)非平穩(wěn)隨機(jī)振動(dòng)響應(yīng)分析的離散解析法同樣保證了在空間域的精確性.

    基于虛擬激勵(lì)法,任一時(shí)刻N(yùn)G個(gè)空間點(diǎn)第i個(gè)方向的位移響應(yīng)功率譜密度矩陣為

    需要說明的是,式(39)需要在一系列離散頻率點(diǎn)ωk上進(jìn)行計(jì)算,其中ωk=ωL+kΔω,Δω=(ωU-ωL)/(Nω),k=0,1,···,Nω,ωU和ωL為截?cái)囝l率上下限,Nω為離散頻率點(diǎn)數(shù).

    響應(yīng)功率譜密度是關(guān)于頻率的單變量光滑函數(shù),通過頻率帶寬離散化操作,對(duì)響應(yīng)功率譜密度簡(jiǎn)單數(shù)值積分即可獲得響應(yīng)的均方根

    同樣地,基于獲得的虛擬位移響應(yīng),亦可容易獲得其他待求響應(yīng)的功率譜密度矩陣及響應(yīng)均方根.需要強(qiáng)調(diào)的是,在基于有限元計(jì)算圓柱殼內(nèi)力功率譜密度時(shí),內(nèi)力的計(jì)算需要對(duì)位移形函數(shù)求導(dǎo),使得連續(xù)性減低,導(dǎo)致內(nèi)力的計(jì)算精度比位移低,出現(xiàn)如線性三角形單元中常應(yīng)力情形.而離散解析法采用位移關(guān)于空間坐標(biāo)的偏導(dǎo)數(shù)先解析求導(dǎo),后將導(dǎo)函數(shù)離散化的方式,獲得各類內(nèi)力分布的精確解,避免了空間離散化對(duì)結(jié)果精度的影響,所得的內(nèi)力功率譜和均方根與位移是同等精確的.

    因此,本節(jié)提出的離散解析法是對(duì)解析法的離散化求解,具有極高精度,可作為基準(zhǔn)解驗(yàn)證中厚圓柱殼隨機(jī)振動(dòng)分析的其他數(shù)值方法.

    4 數(shù)值算例

    4.1 中厚圓柱殼自由振動(dòng)精確解

    采用文獻(xiàn)[12]中的封閉中厚圓柱殼結(jié)構(gòu),彈性模量E=210 GPa,密度ρ=7800 kg/m3,長(zhǎng)度L=0.502 m,半徑R=0.063 5 m,厚度h=0.003 26 m,泊松比υ=0.28.表1 展示了前5 階精確頻率.文獻(xiàn)[12]及有限元(FEM)所得自振頻率也列于表1.有限元結(jié)果基于ABAQUS 軟件,結(jié)構(gòu)模型分別劃分為500,1000,2000,4000,6000 及8000 個(gè)S4 R 單元(考慮縮減積分的4 節(jié)點(diǎn)彎曲薄殼或厚殼單元) 或S8 R 單元(考慮縮減積分的8 節(jié)點(diǎn)雙重彎曲厚殼單元).可見有限元計(jì)算精度受限于幾何網(wǎng)格疏密程度和單元類型選擇等因素,存在一定誤差,而文獻(xiàn)級(jí)數(shù)解[12]與本文所得精確頻率更加接近.

    表1 封閉中厚殼的前5 階自振頻率(Hz)Table 1 The first 5 natural frequencies (Hz) of closed moderately thick cylindrical shells

    除封閉中厚殼外,考慮一個(gè)開口角度為φT=π 的中厚圓柱殼,彈性模量E=27.466 GPa,密度ρ=7850 kg/m3,長(zhǎng)度L=14.4 m,半徑R=3.6 m,厚度h=0.36 m,泊松比υ=0.2.圖2 分別給出了開口中厚殼的第1 階、第2 階、第3 階及第5 階自振頻率和振型,有限元結(jié)果由ABAQUS 軟件對(duì)結(jié)構(gòu)剖分為4000 個(gè)S8 R 單元得到.有限元結(jié)果與本文精確分析獲得的自由振動(dòng)自振頻率及振型吻合良好.

    圖2 開口中厚殼的自振頻率(Hz)及振型(精確解和有限元解)Fig.2 Natural frequencies (Hz) and modal functions of open moderately thick shells (exact solutions and FEM results)

    4.2 隨機(jī)點(diǎn)激勵(lì)作用下中厚圓柱殼振動(dòng)響應(yīng)

    本節(jié)考慮如下封閉中厚圓柱殼結(jié)構(gòu):楊氏模量E=2.96 GPa、密度ρ=0.733 kg/m3、長(zhǎng)度L=14.4 m、半徑R=3.6 m、泊松比υ=0.25、厚度h=0.36 m、阻尼比ζ=0.05.分別在兩個(gè)空間點(diǎn) case 1:(L/2,0);case 2:(L/4,π/4)對(duì)中厚殼上施加徑向隨機(jī)點(diǎn)激勵(lì).該隨機(jī)激勵(lì)為限帶白噪聲,其頻率范圍為[20,2000] Hz,功率譜密度為1 N2/Hz.

    考慮這兩種隨機(jī)點(diǎn)激勵(lì),利用解析法(analytical method,AM)、離散解析法(discrete analytical method,DAM)及ABAQUS 軟件計(jì)算,表2 分別給出了(L/2,0)處隨機(jī)響應(yīng)的均方根與各方法相應(yīng)的CPU 運(yùn)行時(shí)間,其中ABAQUS 軟件結(jié)果由基于傳統(tǒng)功率譜方法的隨機(jī)振動(dòng)分析模塊給出,中厚殼結(jié)構(gòu)被劃分為4000 (40 × 100)個(gè)S8 R 單元.為便于比較計(jì)算效率,基于離散解析法的中厚殼結(jié)構(gòu)同樣劃分為4000 (40 × 100) 個(gè)網(wǎng)格,空間離散點(diǎn)總數(shù)NG為4141 (41 × 101).三類方法均采用激勵(lì)帶寬內(nèi)的全部自振模態(tài)參與隨機(jī)振動(dòng)分析.對(duì)本小節(jié)所考慮的中厚殼結(jié)構(gòu),頻帶[20,2000] Hz 內(nèi)存在14 階自振模態(tài),因此式(34)中離散解析法的正則坐標(biāo)η維度NF取為14.

    由表2可知,兩類隨機(jī)點(diǎn)激勵(lì)作用下,離散解析法所得位移、加速度及彎矩Mx響應(yīng)及均方根與解析法完全吻合,而ABAQUS 軟件結(jié)果與解析法相比存在一定誤差.這驗(yàn)證了時(shí)域、頻域數(shù)值積分的精度是非常高的.由于解析法單次只能計(jì)算一個(gè)空間點(diǎn)響應(yīng),為了公平比較三類方法的計(jì)算效率,將能夠多點(diǎn)運(yùn)算的離散解析解及ABAQUS 軟件結(jié)果的CPU 運(yùn)行時(shí)間均攤至各節(jié)點(diǎn).可以看到,離散解析解法在保證計(jì)算精度的前提下,計(jì)算效率得到了極大提升,展示了該方法的計(jì)算優(yōu)越性.值得注意的是,基于空間域離散建模的ABAQUS 軟件結(jié)果計(jì)算精度依賴于結(jié)構(gòu)的離散形式,而離散解析法基于精確振型函數(shù)的解析積分或微分后再將節(jié)點(diǎn)坐標(biāo)代入式(39)中,其空間離散點(diǎn)多寡并不影響離散解析法的計(jì)算精度.

    表2 徑向隨機(jī)點(diǎn)激勵(lì)作用下封閉中厚殼(L/2,0)處響應(yīng)均方根Table 2 Response RMSs of (L/2,0) for the closed moderately thick shells under random point excitation

    圖3 給出了兩類隨機(jī)點(diǎn)激勵(lì)作用下封閉中厚殼徑向位移響應(yīng)功率譜密度曲線.除位移外,中厚圓柱殼的內(nèi)力響應(yīng)同樣值得關(guān)注.隨機(jī)點(diǎn)激勵(lì)(case 2)下相應(yīng)的結(jié)構(gòu)薄膜內(nèi)力Nx和彎矩Mx的響應(yīng)功率譜在圖4 中給出.結(jié)果表明,離散解析法和解析法得到的響應(yīng)功率譜密度在[20,2000] Hz 內(nèi)完全吻合,ABAQUS 軟件結(jié)果在低頻范圍內(nèi)與解析解吻合良好,但是隨著激勵(lì)頻率增加,ABAQUS 軟件結(jié)果呈現(xiàn)了一定的偏差,且ABAQUS 軟件內(nèi)力計(jì)算結(jié)果的精度低于位移計(jì)算結(jié)果.由于解析法結(jié)果采用了精確振型,保證了沒有錯(cuò)根和漏根.因此,這種偏差顯然是由ABAQUS 軟件計(jì)算高階自振頻率的誤差引起的.

    圖3 隨機(jī)點(diǎn)激勵(lì)下封閉中厚殼(L/2,0)處徑向位移響應(yīng)功率譜密度曲線Fig.3 Deflection response PSD curves of (L/2,0) for the closed moderately thick shells under random point excitation

    圖4 隨機(jī)點(diǎn)激勵(lì)(case 2)下封閉中厚殼(L/2,0)處薄膜內(nèi)力Nx和彎矩Mx 功率譜密度曲線Fig.4 Response PSD curves of membrane force Nx and bending moment Mx of (L/2,0) for the closed moderately thick shells under random point excitation (case 2)

    4.3 隨機(jī)環(huán)形線激勵(lì)作用下中厚圓柱殼振動(dòng)響應(yīng)

    本節(jié)對(duì)于上節(jié)中討論的半徑R=3.6 m 的封閉中厚圓柱殼,分別取厚度h為0.1 m,0.18 m,0.36 m,也即將厚徑比(h/R)分別取為1/36,1/20,1/10.考慮在封閉中厚殼的環(huán)線 (L/2,0~2π)上施加徑向隨機(jī)線激勵(lì),該激勵(lì)同樣是限帶寬高斯白噪聲,頻帶為[20,2000] Hz,功率譜密度為1 N2/Hz,針對(duì)圓柱殼的三種厚徑比情形,采用激勵(lì)帶寬內(nèi)的全部振型參與隨機(jī)振動(dòng)分析(NF分別取為46,27,14).針對(duì)上述三種厚徑比,表3 分別給出了中厚殼在(L/2,0)處的徑向位移、速度和加速度均方根.結(jié)果表明,在徑向隨機(jī)線激勵(lì)下,隨著中厚殼厚徑比的增加,結(jié)構(gòu)的抗彎剛度增大,從而導(dǎo)致徑向位移、速度及加速度響應(yīng)均方根減小.

    表3 隨機(jī)環(huán)形線激勵(lì)作用下封閉中厚殼(L/2,0)處徑向響應(yīng)均方根Table 3 Radial response RMSs of (L/2,0) for the closed moderately thick shells under random annular line excitation

    而且,也討論了徑向隨機(jī)線激勵(lì)下中厚殼結(jié)構(gòu)的厚徑比對(duì)彎矩均方根的影響.圖5 分別給出了厚徑比為1/36,1/20,1/10 下的中厚殼的彎矩Mx和Mθ的功率譜密度函數(shù).由式(5)可知,彎曲剛度D的表達(dá)式中包含厚度的三次方項(xiàng).在虛擬激勵(lì)法中,彎矩響應(yīng)Mx和Mθ的功率譜密度可由相應(yīng)的虛擬彎矩與其自共軛相乘,即與D2呈正比.因此,隨著中厚殼厚徑比的增加,彎矩響應(yīng)Mx和Mθ的功率譜密度隨之增大.

    圖5 隨機(jī)環(huán)形線激勵(lì)下中厚殼(L/2,0)處彎矩響應(yīng)功率譜密度曲線Fig.5 Bending moment response PSD curves of (L/2,0) for the closed moderately thick shells under random annular line excitation

    4.4 分布隨機(jī)激勵(lì)下中厚圓柱殼振動(dòng)響應(yīng)

    針對(duì)4.1 小節(jié)討論的開口中厚圓柱殼結(jié)構(gòu),在其表面施加徑向分布的隨機(jī)激勵(lì).考慮激勵(lì)具有如下兩類空間分布形式case 1:p1(x,θ)=x;case 2:p2(x,θ)=5θ.隨機(jī)激勵(lì)為[0,100] Hz 內(nèi)的限帶寬白噪聲,其功率譜密度Spp=1 Pa2/Hz,采用激勵(lì)帶寬內(nèi)的全部振型參與隨機(jī)振動(dòng)分析(NF=32).阻尼比ζ=0.05.圖6 給出了這兩類分布激勵(lì)在開口中厚圓柱殼上的具體空間分布形式.

    圖6 兩類隨機(jī)分布載荷的空間分布形式Fig.6 Spatial distribution forms of two random distributed loads

    針對(duì)這兩類隨機(jī)分布載荷,圖7~圖9 分別給出了開口中厚殼徑向位移、速度及加速度響應(yīng)的均方根空間分布.對(duì)于分布形式case 1 的隨機(jī)載荷,其關(guān)于θ=π/2 的軸對(duì)稱,激勵(lì)峰值為x=L=14.4 m.由圖7 (a)、圖8 (a)及圖9 (a)可以看到,位移、速度及加速度響應(yīng)均方根峰值逐漸向激勵(lì)峰值方向移動(dòng),這些均方根關(guān)于θ=π/2 的軸是對(duì)稱的,并且位移、速度響應(yīng)均方根關(guān)于x=L/2 的軸呈現(xiàn)一定的對(duì)稱性,而加速度響應(yīng)均方根關(guān)于x=L/2 的軸沒有對(duì)稱性.對(duì)于隨機(jī)載荷case 2,其關(guān)于x=L/2 的軸對(duì)稱,激勵(lì)峰值為θ=φT=π.由圖7 (b)、圖8 (b)及圖9 (b)可知,位移、速度及加速度響應(yīng)均方根峰值逐漸向θ最大值位置移動(dòng),這些均方根關(guān)于x=L/2 的軸也是對(duì)稱,并且位移、速度響應(yīng)均方根關(guān)于θ=π/2 的軸呈現(xiàn)一定的對(duì)稱性,而加速度響應(yīng)均方根關(guān)于θ=π/2 的軸沒有對(duì)稱性.其原因在于,與其他空間位置相比,激勵(lì)峰值作用位置往往能夠激起結(jié)構(gòu)更高階振型.由式(30)可知,高階頻率對(duì)加速度響應(yīng)影響最大,其次是速度響應(yīng).由此可知,由于高階參振振型的影響,激勵(lì)峰值作用位置附近的加速度響應(yīng)均方根出現(xiàn)較大值,導(dǎo)致圖9 (a)、圖9 (b)所示的加速度響應(yīng)均方根分別關(guān)于x=L/2,θ=π/2 的軸呈現(xiàn)非對(duì)稱分布.分布隨機(jī)激勵(lì)的空間分布形式對(duì)中厚圓柱殼的加速度響應(yīng)均方根分布影響較大.由于隨著響應(yīng)均方根的增大,結(jié)構(gòu)響應(yīng)的變異性隨之增大.因此在實(shí)際工程中,為了保證中厚殼結(jié)構(gòu)在隨機(jī)環(huán)境激勵(lì)作用下的安全性,除徑向位移及速度響應(yīng)外,應(yīng)尤其關(guān)注結(jié)構(gòu)的加速度響應(yīng)分布.

    圖7 兩類分布載荷作用下開口圓柱殼的位移響應(yīng)均方根分布Fig.7 RMS distributions of deflection for open shell under two distributed random loads

    圖8 兩類分布載荷作用下開口殼的速度響應(yīng)均方根分布Fig.8 RMS distributions of velocity for open shell under two distributed random loads

    圖9 兩類分布載荷作用下開口殼的加速度響應(yīng)均方根分布Fig.9 RMS distributions of acceleration for open shell under two distributed random loads

    4.5 均勻調(diào)制非平穩(wěn)激勵(lì)下中厚圓柱殼振動(dòng)響應(yīng)

    除平穩(wěn)隨機(jī)激勵(lì)外,均勻調(diào)制非平穩(wěn)激勵(lì)下中厚圓柱殼的隨機(jī)振動(dòng)響應(yīng)同樣值得重視.對(duì)于如下開口中厚圓柱殼(φT=π,E=27.466 GPa,ρ=7850 kg/m3,L=14.4 m,R=3.6 m,υ=0.2),考慮在開口殼上施加均勻調(diào)制非平穩(wěn)的隨機(jī)徑向均布載荷q3(x,θ,t)=f(t)X(t).隨機(jī)激勵(lì)X(t)假設(shè)為[0,100] Hz 內(nèi)的限帶寬高斯白噪聲激勵(lì),激勵(lì)PSDSxx(ω)=1 Pa2/Hz,采用激勵(lì)帶寬內(nèi)的全部振型參與隨機(jī)振動(dòng)分析.

    采用如圖10 所示的時(shí)間調(diào)制函數(shù)f(t)[33]

    圖10 時(shí)間調(diào)制函數(shù)f(t)Fig.10 Time-modulated function f(t)

    其中,A為時(shí)間調(diào)制函數(shù)幅值,?(t)為Heaviside 函數(shù).本文中,取A=105,激勵(lì)持續(xù)時(shí)間T=9.6 s.可知,t=4.8 s 時(shí),時(shí)間調(diào)制函數(shù)f(t)取最大值.

    表4 給出了t=4.8 s 時(shí)開口殼中心點(diǎn)(L/2,π/2)處的徑向位移、速度及加速度響應(yīng)均方根,其中結(jié)構(gòu)厚徑比分別采用1/20和1/10,對(duì)應(yīng)的參振頻率階數(shù)NF分別為65和32.為了驗(yàn)證結(jié)果,基于10 000次并行蒙特卡洛模擬(MCS)結(jié)果也被給出.利用譜表達(dá)方法,隨機(jī)激勵(lì)描述為

    表4 非平穩(wěn)隨機(jī)均布激勵(lì)下中厚殼中心點(diǎn)響應(yīng)均方根 (t=4.8 s)Table 4 Stochastic response RMSs of the center of moderately thick shell under nonstationary random distributed excitations(t=4.8 s)

    其中,Sqq(x,θ,t,ωk)為非平穩(wěn)隨機(jī)載荷在第k個(gè)頻率點(diǎn)的功率譜密度函數(shù),φk是[0,2π]內(nèi)獨(dú)立均勻分布的隨機(jī)相位角.對(duì)于頻帶[0,100] Hz 的離散頻率點(diǎn)數(shù)Nω,離散解析法和蒙特卡洛模擬均將其取為500 以保證公平比較計(jì)算效率.此外,若假設(shè)0~9.8 s 內(nèi)f(t)均等于A,則該均勻調(diào)制非平穩(wěn)激勵(lì)轉(zhuǎn)變?yōu)槠椒€(wěn)隨機(jī)激勵(lì).表4 同時(shí)給出了中厚殼在該平穩(wěn)均布激勵(lì)下的解析解.可以看到,非平穩(wěn)隨機(jī)振動(dòng)的離散解析解法在采用精細(xì)積分法進(jìn)行時(shí)程分析時(shí),取了較大的時(shí)間步長(zhǎng)(Δt=0.2 s),得到的計(jì)算結(jié)果與平穩(wěn)隨機(jī)響應(yīng)的解析解吻合良好,且計(jì)算效率很高.而蒙特卡洛模擬即使時(shí)間步長(zhǎng)已經(jīng)足夠小,其計(jì)算誤差仍大于離散解析解,這主要?dú)w因于蒙特卡洛的隨機(jī)收斂性.此外,從耗費(fèi)的CPU 計(jì)算時(shí)間,也可看出,即使蒙特卡洛模擬采用了并行算法,計(jì)算效率仍遠(yuǎn)低于離散解析法.

    非平穩(wěn)均布激勵(lì)下,圖11 分別繪制了厚徑比為1/20和1/10 的兩種情況下位移響應(yīng)的均方根曲線.可以看到,兩種方法的計(jì)算結(jié)果吻合良好,但蒙特卡洛模擬曲線存在一定數(shù)值振蕩.對(duì)于均勻調(diào)制非平穩(wěn)激勵(lì)下的中厚殼隨機(jī)振動(dòng)響應(yīng),其徑向位移、速度及加速度響應(yīng)均方根隨時(shí)間的變化趨勢(shì)與圖10所示的時(shí)間調(diào)制函數(shù)趨勢(shì)一致.

    圖11 非平穩(wěn)隨機(jī)激勵(lì)下中厚殼中心點(diǎn)位移響應(yīng)均方根Fig.11 Stochastic deflection responses RMSs of the center of moderately thick shell under nonstationary random excitations

    圖11 非平穩(wěn)隨機(jī)激勵(lì)下中厚殼中心點(diǎn)位移響應(yīng)均方根(續(xù))Fig.11 Stochastic deflection responses RMSs of the center of moderately thick shell under nonstationary random excitations(continued)

    為了討論開口中厚圓柱殼的厚徑比對(duì)位移響應(yīng)功率譜密度的影響,圖12 分別繪制了厚徑比h/R=1/36,1/20,1/15 下不同時(shí)刻(t=4.8 s,t=7 s)的中厚殼中心點(diǎn)處的位移響應(yīng)功率譜密度曲線.由圖12可知,徑向位移響應(yīng)的功率譜密度曲線的第一個(gè)峰值出現(xiàn)位置從約7.4 Hz 后移到約14.2 Hz.這是由于中厚圓柱殼的基頻隨著中厚殼厚徑比的增加而增大.這表明,在線性隨機(jī)振動(dòng)分析中,結(jié)構(gòu)的自振頻率精確與否將直接影響響應(yīng)功率譜密度峰值出現(xiàn)的位置.此外,在徑向激勵(lì)作用下,中厚殼徑向位移響應(yīng)功率譜密度的幅值整體減小.因此,隨著結(jié)構(gòu)厚徑比的增加,徑向位移的均方根逐漸減小,也即響應(yīng)的變異性減小,圓柱殼的安全性得到提高.

    圖12 非平穩(wěn)隨機(jī)激勵(lì)下中厚殼中心點(diǎn)位移響應(yīng)功率譜密度Fig.12 Stochastic deflection responses PSD curves of the center of moderately thick shell under nonstationary random excitations

    圖12 非平穩(wěn)隨機(jī)激勵(lì)下中厚殼中心點(diǎn)位移響應(yīng)功率譜密度(續(xù))Fig.12 Stochastic deflection responses PSD curves of the center of moderately thick shell under nonstationary random excitations(continued)

    5 結(jié)論

    本文考察了考慮轉(zhuǎn)動(dòng)慣量和剪切變形的中厚圓柱殼的隨機(jī)振動(dòng),獲得了平穩(wěn)、非平穩(wěn)隨機(jī)振動(dòng)響應(yīng)的基準(zhǔn)解.通過引入中厚圓柱殼的精確自振頻率和封閉振型函數(shù),基于虛擬激勵(lì)法,解析推導(dǎo)了中厚圓柱殼在平穩(wěn)、非平穩(wěn)隨機(jī)激勵(lì)作用下的響應(yīng)功率譜密度函數(shù).為了充分發(fā)揮虛擬激勵(lì)法在矩陣運(yùn)算中的高效性優(yōu)勢(shì),將空間域積分解析求解,頻域和時(shí)域積分?jǐn)?shù)值求解,提出了離散解析法高效計(jì)算了中厚殼在各類隨機(jī)激勵(lì)下響應(yīng)均方根的基準(zhǔn)解.通過與解析解、ABAQUS 軟件結(jié)果及蒙特卡洛模擬結(jié)果的比較,展示了離散解析法的準(zhǔn)確性和高效性.通過算例分析,得到以下主要結(jié)論:

    (1) 徑向隨機(jī)線激勵(lì)下,隨著圓柱殼厚度的增加,結(jié)構(gòu)徑向抗壓強(qiáng)度增大,位移、速度及加速度響應(yīng)均方根等比例降低;結(jié)構(gòu)彎曲剛度增大,導(dǎo)致彎矩響應(yīng)均方根增加.

    (2) 分布隨機(jī)激勵(lì)的空間分布形式對(duì)中厚圓柱殼的加速度響應(yīng)均方根分布影響較大.為保證結(jié)構(gòu)安全性,中厚殼結(jié)構(gòu)的加速度響應(yīng)也應(yīng)受到關(guān)注.

    (3) 非平穩(wěn)隨機(jī)激勵(lì)下,中厚圓柱殼位移響應(yīng)均方根隨時(shí)間變化與時(shí)間調(diào)制函數(shù)相一致.此外,隨著圓柱殼厚度的增加,位移響應(yīng)功率譜的幅值整體呈減小趨勢(shì),相應(yīng)的均方根(即變異性)減小,因此結(jié)構(gòu)安全性得到提高.

    猜你喜歡
    方根振型圓柱
    方根拓展探究
    關(guān)于模態(tài)綜合法的注記
    縱向激勵(lì)下大跨鋼桁拱橋高階振型效應(yīng)分析
    圓柱的體積計(jì)算
    “圓柱與圓錐”復(fù)習(xí)指導(dǎo)
    塔腿加過渡段輸電塔動(dòng)力特性分析
    均方根嵌入式容積粒子PHD 多目標(biāo)跟蹤方法
    揭開心算方根之謎
    削法不同 體積有異
    結(jié)構(gòu)振型幾何辨識(shí)及應(yīng)用研究
    山西建筑(2015年14期)2015-06-05 09:37:07
    熟女少妇亚洲综合色aaa.| 日本熟妇午夜| 每晚都被弄得嗷嗷叫到高潮| 丰满的人妻完整版| 精品免费久久久久久久清纯| 国产国语露脸激情在线看| 国产又黄又爽又无遮挡在线| 婷婷精品国产亚洲av| 色在线成人网| 日本a在线网址| 久久这里只有精品19| 欧美一级a爱片免费观看看 | 老司机在亚洲福利影院| 精品久久久久久,| 免费高清在线观看日韩| 1024手机看黄色片| 级片在线观看| 国产视频一区二区在线看| 日韩精品免费视频一区二区三区| 夜夜爽天天搞| 午夜久久久久精精品| 手机成人av网站| 美女免费视频网站| videosex国产| 久久久精品国产亚洲av高清涩受| 亚洲国产日韩欧美精品在线观看 | 久久精品国产99精品国产亚洲性色| 国产成人精品久久二区二区免费| 亚洲精品一区av在线观看| 亚洲无线在线观看| 欧美日韩一级在线毛片| 淫妇啪啪啪对白视频| 女同久久另类99精品国产91| 88av欧美| 欧美日韩一级在线毛片| 国产99久久九九免费精品| 国产午夜福利久久久久久| 国产欧美日韩一区二区三| 男女视频在线观看网站免费 | 欧美绝顶高潮抽搐喷水| 婷婷精品国产亚洲av在线| 成人亚洲精品一区在线观看| 国产97色在线日韩免费| 欧美乱码精品一区二区三区| 国产精品久久久av美女十八| 狠狠狠狠99中文字幕| 国产单亲对白刺激| 国产伦在线观看视频一区| 免费看a级黄色片| 久久久精品欧美日韩精品| 无人区码免费观看不卡| 国产激情偷乱视频一区二区| av中文乱码字幕在线| 日韩视频一区二区在线观看| 可以免费在线观看a视频的电影网站| 亚洲一区二区三区不卡视频| 美女午夜性视频免费| 大型av网站在线播放| 成人三级做爰电影| 亚洲最大成人中文| 国产精品久久久av美女十八| 中亚洲国语对白在线视频| 亚洲国产精品久久男人天堂| 久久精品91蜜桃| 2021天堂中文幕一二区在线观 | 成人18禁在线播放| 手机成人av网站| 日韩欧美一区视频在线观看| 18美女黄网站色大片免费观看| 欧美国产日韩亚洲一区| 国产黄片美女视频| 久久精品aⅴ一区二区三区四区| 亚洲专区国产一区二区| 一进一出好大好爽视频| 亚洲精品久久国产高清桃花| 国产精品一区二区免费欧美| 日本免费一区二区三区高清不卡| 亚洲九九香蕉| 国产精品久久电影中文字幕| 精品一区二区三区视频在线观看免费| 黄色视频不卡| 女人被狂操c到高潮| 欧美绝顶高潮抽搐喷水| 中文字幕人成人乱码亚洲影| av视频在线观看入口| 美女扒开内裤让男人捅视频| 一边摸一边做爽爽视频免费| 免费高清视频大片| 久99久视频精品免费| 久久国产亚洲av麻豆专区| 亚洲 国产 在线| 听说在线观看完整版免费高清| 国产视频内射| 国产亚洲精品av在线| 亚洲美女黄片视频| 性色av乱码一区二区三区2| 精品日产1卡2卡| 美女国产高潮福利片在线看| 国产精品,欧美在线| 给我免费播放毛片高清在线观看| 高潮久久久久久久久久久不卡| 久久久久久人人人人人| 亚洲无线在线观看| 亚洲成a人片在线一区二区| 国产伦在线观看视频一区| 亚洲国产精品sss在线观看| 国产精品久久久久久人妻精品电影| √禁漫天堂资源中文www| 久久国产精品影院| 久久精品夜夜夜夜夜久久蜜豆 | 18禁黄网站禁片午夜丰满| 亚洲国产欧美网| 最好的美女福利视频网| 亚洲av第一区精品v没综合| 午夜久久久在线观看| 久热爱精品视频在线9| 88av欧美| 99热6这里只有精品| 热re99久久国产66热| 伊人久久大香线蕉亚洲五| 99在线人妻在线中文字幕| 99久久国产精品久久久| 国产成+人综合+亚洲专区| 99re在线观看精品视频| 最近最新免费中文字幕在线| 国产激情偷乱视频一区二区| 一区福利在线观看| 18禁裸乳无遮挡免费网站照片 | 国产蜜桃级精品一区二区三区| a在线观看视频网站| 啦啦啦观看免费观看视频高清| 久久久久亚洲av毛片大全| 性色av乱码一区二区三区2| svipshipincom国产片| 成人特级黄色片久久久久久久| 长腿黑丝高跟| 99久久综合精品五月天人人| 国产一级毛片七仙女欲春2 | 亚洲中文字幕一区二区三区有码在线看 | 久久草成人影院| 色婷婷久久久亚洲欧美| 丰满的人妻完整版| 久久久国产成人免费| 在线国产一区二区在线| 巨乳人妻的诱惑在线观看| 视频在线观看一区二区三区| 日韩大尺度精品在线看网址| 午夜福利高清视频| 香蕉丝袜av| svipshipincom国产片| 国产成人欧美| 欧美日韩亚洲综合一区二区三区_| 老司机福利观看| 国产精品亚洲一级av第二区| 俄罗斯特黄特色一大片| 午夜成年电影在线免费观看| 日韩欧美 国产精品| 露出奶头的视频| 可以免费在线观看a视频的电影网站| 少妇 在线观看| 亚洲熟妇中文字幕五十中出| 精品午夜福利视频在线观看一区| 久久香蕉精品热| 欧美激情极品国产一区二区三区| 法律面前人人平等表现在哪些方面| 日本黄色视频三级网站网址| 日韩视频一区二区在线观看| 国产精品一区二区免费欧美| 国产精品久久久久久亚洲av鲁大| 亚洲av熟女| 男女那种视频在线观看| av中文乱码字幕在线| 一级片免费观看大全| 50天的宝宝边吃奶边哭怎么回事| 精华霜和精华液先用哪个| www.www免费av| 一区二区三区国产精品乱码| 婷婷六月久久综合丁香| 久久九九热精品免费| 精品卡一卡二卡四卡免费| 午夜福利欧美成人| 午夜精品在线福利| 国产熟女xx| 国产高清激情床上av| 制服诱惑二区| 色老头精品视频在线观看| 岛国视频午夜一区免费看| 亚洲欧美激情综合另类| 老司机在亚洲福利影院| 88av欧美| 午夜久久久在线观看| 一区二区三区高清视频在线| 男女视频在线观看网站免费 | 欧美日韩黄片免| 日韩视频一区二区在线观看| 亚洲第一电影网av| 国产爱豆传媒在线观看 | 在线av久久热| 又紧又爽又黄一区二区| 19禁男女啪啪无遮挡网站| 国产单亲对白刺激| 国产精品1区2区在线观看.| 亚洲熟妇熟女久久| 国产日本99.免费观看| 久久久久国产精品人妻aⅴ院| 黄色女人牲交| 欧美黄色淫秽网站| 中文亚洲av片在线观看爽| 免费高清在线观看日韩| 级片在线观看| 午夜福利欧美成人| 免费在线观看黄色视频的| 色老头精品视频在线观看| 日本 欧美在线| 一级毛片精品| 国产精品免费一区二区三区在线| 欧美一级毛片孕妇| 老汉色∧v一级毛片| 亚洲九九香蕉| 成熟少妇高潮喷水视频| 又紧又爽又黄一区二区| 国产精品久久久人人做人人爽| 亚洲欧美一区二区三区黑人| 成年免费大片在线观看| 久久精品亚洲精品国产色婷小说| 99久久久亚洲精品蜜臀av| 天堂动漫精品| 成人18禁在线播放| 国产精品1区2区在线观看.| 午夜免费观看网址| 国产成人精品无人区| www.熟女人妻精品国产| 97碰自拍视频| 国产欧美日韩一区二区精品| 又黄又粗又硬又大视频| www.999成人在线观看| 精品日产1卡2卡| 999久久久国产精品视频| a级毛片在线看网站| 久久久久久久久中文| 久久99热这里只有精品18| 国产极品粉嫩免费观看在线| 老司机在亚洲福利影院| 免费高清在线观看日韩| 欧美久久黑人一区二区| 男人舔女人的私密视频| 18美女黄网站色大片免费观看| 91九色精品人成在线观看| 91麻豆精品激情在线观看国产| 99热6这里只有精品| 黄网站色视频无遮挡免费观看| 性色av乱码一区二区三区2| 99精品在免费线老司机午夜| 久久久久久久久久黄片| 特大巨黑吊av在线直播 | 久久久精品国产亚洲av高清涩受| 国产成人欧美| 久久性视频一级片| 757午夜福利合集在线观看| 国产一区二区激情短视频| 欧美国产日韩亚洲一区| 精品久久久久久,| 一级作爱视频免费观看| 国产精品免费视频内射| 叶爱在线成人免费视频播放| 看黄色毛片网站| 99热只有精品国产| 国产亚洲欧美精品永久| 欧美又色又爽又黄视频| 91国产中文字幕| 欧美国产精品va在线观看不卡| 美女高潮喷水抽搐中文字幕| 在线观看日韩欧美| 国产精品乱码一区二三区的特点| 国产精品av久久久久免费| 黄色a级毛片大全视频| 亚洲国产看品久久| 国产亚洲欧美98| 身体一侧抽搐| 午夜福利视频1000在线观看| av免费在线观看网站| 丝袜美腿诱惑在线| 国产精品电影一区二区三区| 男女床上黄色一级片免费看| 欧美大码av| 搡老岳熟女国产| 精品无人区乱码1区二区| 亚洲电影在线观看av| 久久亚洲真实| 白带黄色成豆腐渣| 久久中文字幕人妻熟女| 黄色毛片三级朝国网站| 亚洲,欧美精品.| 国产精品久久电影中文字幕| 波多野结衣av一区二区av| 亚洲成人久久爱视频| 成人特级黄色片久久久久久久| 色播亚洲综合网| 女人爽到高潮嗷嗷叫在线视频| 久久精品影院6| 天堂影院成人在线观看| 国产精品久久久av美女十八| 禁无遮挡网站| 91九色精品人成在线观看| 婷婷六月久久综合丁香| 久久香蕉激情| 国产熟女xx| 国产主播在线观看一区二区| tocl精华| 久久久久免费精品人妻一区二区 | 最近最新免费中文字幕在线| 亚洲欧美日韩高清在线视频| 亚洲自拍偷在线| 日韩高清综合在线| 真人做人爱边吃奶动态| 国产精品亚洲一级av第二区| 欧美乱码精品一区二区三区| 2021天堂中文幕一二区在线观 | 欧美日本视频| 在线观看一区二区三区| 久久香蕉激情| 又黄又粗又硬又大视频| 国产真人三级小视频在线观看| 男女视频在线观看网站免费 | 国产亚洲精品一区二区www| 午夜a级毛片| bbb黄色大片| 婷婷精品国产亚洲av| 淫秽高清视频在线观看| 成人三级做爰电影| 黄色 视频免费看| 满18在线观看网站| 日日爽夜夜爽网站| 一区福利在线观看| 欧美成狂野欧美在线观看| 精品卡一卡二卡四卡免费| 美女午夜性视频免费| 欧美成狂野欧美在线观看| 成人国产综合亚洲| 女生性感内裤真人,穿戴方法视频| 国产区一区二久久| 亚洲精品在线美女| 亚洲精品美女久久av网站| 国产麻豆成人av免费视频| 亚洲国产精品合色在线| 琪琪午夜伦伦电影理论片6080| 亚洲第一电影网av| 麻豆久久精品国产亚洲av| 91av网站免费观看| 九色国产91popny在线| 97碰自拍视频| 日韩大码丰满熟妇| 国产99久久九九免费精品| 亚洲第一av免费看| 国产真实乱freesex| 又黄又爽又免费观看的视频| 波多野结衣高清作品| 亚洲一区高清亚洲精品| 身体一侧抽搐| 18禁美女被吸乳视频| 亚洲国产中文字幕在线视频| 亚洲av五月六月丁香网| 成人免费观看视频高清| 亚洲精品国产一区二区精华液| 色播亚洲综合网| 日本一本二区三区精品| 午夜影院日韩av| 久久九九热精品免费| 免费在线观看视频国产中文字幕亚洲| 宅男免费午夜| 亚洲一区高清亚洲精品| 国产午夜福利久久久久久| 大型av网站在线播放| 亚洲av成人一区二区三| 久久天堂一区二区三区四区| 成人三级做爰电影| 国产三级在线视频| 久久香蕉激情| 午夜两性在线视频| 黄频高清免费视频| 波多野结衣巨乳人妻| 国产精品,欧美在线| 首页视频小说图片口味搜索| 国产亚洲av高清不卡| 制服丝袜大香蕉在线| 日本在线视频免费播放| 一区二区三区高清视频在线| av视频在线观看入口| 亚洲av成人一区二区三| 日本精品一区二区三区蜜桃| 欧美一级毛片孕妇| 久久久久久人人人人人| 热99re8久久精品国产| 一级毛片高清免费大全| 亚洲激情在线av| 首页视频小说图片口味搜索| 麻豆久久精品国产亚洲av| 免费搜索国产男女视频| 夜夜躁狠狠躁天天躁| 国产精品久久久av美女十八| 久99久视频精品免费| 亚洲成人免费电影在线观看| 少妇裸体淫交视频免费看高清 | 久久婷婷人人爽人人干人人爱| 人人妻人人澡人人看| 欧美日韩瑟瑟在线播放| 亚洲,欧美精品.| 给我免费播放毛片高清在线观看| 久久久久国产精品人妻aⅴ院| 一个人免费在线观看的高清视频| 淫妇啪啪啪对白视频| 国产在线观看jvid| 香蕉国产在线看| 777久久人妻少妇嫩草av网站| 国产精品精品国产色婷婷| 在线看三级毛片| 日本黄色视频三级网站网址| 国产成人精品久久二区二区91| 丁香六月欧美| 免费电影在线观看免费观看| 婷婷丁香在线五月| 亚洲精品一区av在线观看| 1024手机看黄色片| 老司机在亚洲福利影院| 亚洲av电影在线进入| 看免费av毛片| 男女做爰动态图高潮gif福利片| 夜夜夜夜夜久久久久| 精品电影一区二区在线| 给我免费播放毛片高清在线观看| 亚洲精品一卡2卡三卡4卡5卡| 满18在线观看网站| 亚洲av中文字字幕乱码综合 | www国产在线视频色| 伊人久久大香线蕉亚洲五| 中国美女看黄片| 国产成+人综合+亚洲专区| 1024香蕉在线观看| 99热这里只有精品一区 | 色精品久久人妻99蜜桃| 午夜日韩欧美国产| 精品乱码久久久久久99久播| 国产精品久久电影中文字幕| 久久精品人妻少妇| 亚洲成av人片免费观看| 亚洲欧美精品综合一区二区三区| 亚洲av熟女| 搡老熟女国产l中国老女人| 亚洲av电影在线进入| 日日摸夜夜添夜夜添小说| 1024香蕉在线观看| 法律面前人人平等表现在哪些方面| 国产真人三级小视频在线观看| 少妇的丰满在线观看| 国产亚洲av嫩草精品影院| 久99久视频精品免费| 美国免费a级毛片| 视频区欧美日本亚洲| 欧洲精品卡2卡3卡4卡5卡区| 国产亚洲精品第一综合不卡| 无人区码免费观看不卡| 丝袜在线中文字幕| 亚洲精品久久国产高清桃花| 亚洲av五月六月丁香网| 亚洲欧美一区二区三区黑人| 亚洲国产日韩欧美精品在线观看 | 99久久无色码亚洲精品果冻| 韩国av一区二区三区四区| 免费高清在线观看日韩| 国产1区2区3区精品| 欧美黑人欧美精品刺激| 最新美女视频免费是黄的| 亚洲第一av免费看| 十分钟在线观看高清视频www| 此物有八面人人有两片| 可以在线观看的亚洲视频| 不卡一级毛片| 欧美成狂野欧美在线观看| 欧美日韩黄片免| 青草久久国产| 黄网站色视频无遮挡免费观看| 欧美成人一区二区免费高清观看 | 亚洲精品久久国产高清桃花| 99久久99久久久精品蜜桃| 久久中文看片网| 亚洲欧洲精品一区二区精品久久久| 久久精品aⅴ一区二区三区四区| 人妻久久中文字幕网| 国产在线观看jvid| 久久久久久亚洲精品国产蜜桃av| 日韩免费av在线播放| 91国产中文字幕| 国产精品亚洲一级av第二区| 俺也久久电影网| 亚洲av五月六月丁香网| 最近最新中文字幕大全免费视频| 又黄又粗又硬又大视频| 少妇裸体淫交视频免费看高清 | 操出白浆在线播放| 日本免费一区二区三区高清不卡| 久久九九热精品免费| 啦啦啦 在线观看视频| 淫妇啪啪啪对白视频| 国产私拍福利视频在线观看| www.熟女人妻精品国产| 国产免费av片在线观看野外av| 日本在线视频免费播放| 色综合亚洲欧美另类图片| √禁漫天堂资源中文www| 久久精品夜夜夜夜夜久久蜜豆 | a级毛片在线看网站| 两性午夜刺激爽爽歪歪视频在线观看 | 中文亚洲av片在线观看爽| a级毛片在线看网站| 国产精品一区二区精品视频观看| 国产成人精品无人区| 激情在线观看视频在线高清| 欧美激情高清一区二区三区| 欧美乱码精品一区二区三区| 亚洲性夜色夜夜综合| 成人一区二区视频在线观看| 宅男免费午夜| 久久久久精品国产欧美久久久| 在线观看一区二区三区| 精品久久久久久成人av| 国产av一区二区精品久久| 熟女少妇亚洲综合色aaa.| 亚洲片人在线观看| 人人妻人人澡人人看| 成人国产一区最新在线观看| 亚洲成人久久爱视频| 欧美+亚洲+日韩+国产| 可以在线观看的亚洲视频| 久久欧美精品欧美久久欧美| 人人妻人人澡人人看| 日本免费一区二区三区高清不卡| 无人区码免费观看不卡| 日本免费一区二区三区高清不卡| aaaaa片日本免费| 亚洲欧美日韩无卡精品| e午夜精品久久久久久久| 久久欧美精品欧美久久欧美| 久99久视频精品免费| 韩国精品一区二区三区| 夜夜夜夜夜久久久久| 亚洲欧美日韩无卡精品| 精品久久久久久成人av| 国产精品,欧美在线| 男女那种视频在线观看| 岛国视频午夜一区免费看| 国产成人av激情在线播放| 十八禁网站免费在线| 国产精品电影一区二区三区| 久久精品影院6| 色播在线永久视频| 深夜精品福利| 国内精品久久久久久久电影| 亚洲电影在线观看av| 国产成人精品无人区| 真人做人爱边吃奶动态| 亚洲av美国av| 波多野结衣巨乳人妻| 国产一区二区三区在线臀色熟女| 亚洲自偷自拍图片 自拍| 最近最新中文字幕大全免费视频| 手机成人av网站| 男女视频在线观看网站免费 | 国产精品永久免费网站| 免费在线观看完整版高清| 亚洲 欧美一区二区三区| av有码第一页| 搡老熟女国产l中国老女人| 中文字幕另类日韩欧美亚洲嫩草| 天天一区二区日本电影三级| cao死你这个sao货| 日本成人三级电影网站| 成年人黄色毛片网站| 99久久久亚洲精品蜜臀av| 一进一出抽搐动态| 亚洲自拍偷在线| 日韩欧美在线二视频| 波多野结衣高清作品| 国产亚洲精品综合一区在线观看 | 日本a在线网址| 亚洲国产精品合色在线| a级毛片在线看网站| 夜夜看夜夜爽夜夜摸| 日本一区二区免费在线视频| 亚洲国产毛片av蜜桃av| 国产成人啪精品午夜网站| 一进一出好大好爽视频| 国产在线精品亚洲第一网站| 亚洲精品粉嫩美女一区| 中文资源天堂在线| 搡老妇女老女人老熟妇| 又黄又爽又免费观看的视频| 亚洲国产精品成人综合色| 亚洲五月色婷婷综合| 1024手机看黄色片| 97碰自拍视频| 人妻久久中文字幕网| 午夜福利欧美成人| 最新美女视频免费是黄的| 啦啦啦 在线观看视频| 精品卡一卡二卡四卡免费| 欧美成人性av电影在线观看| 变态另类成人亚洲欧美熟女| 真人一进一出gif抽搐免费| 无遮挡黄片免费观看| 午夜久久久久精精品| 亚洲无线在线观看| 手机成人av网站| 久99久视频精品免费| 日韩精品免费视频一区二区三区| 欧美最黄视频在线播放免费|