• <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
    亚洲熟女精品中文字幕| 青草久久国产| 国产亚洲一区二区精品| 视频区图区小说| 久久人妻熟女aⅴ| 国产精品二区激情视频| 我要看黄色一级片免费的| 亚洲五月婷婷丁香| a级片在线免费高清观看视频| 又黄又粗又硬又大视频| 老司机午夜福利在线观看视频 | 人人妻人人澡人人看| 大香蕉久久网| 飞空精品影院首页| 精品久久久精品久久久| av天堂久久9| 日本一区二区免费在线视频| 欧美在线黄色| 国产日韩欧美视频二区| 中国国产av一级| 在线观看免费午夜福利视频| 国产精品久久久久久人妻精品电影 | 女性生殖器流出的白浆| 国产精品香港三级国产av潘金莲| 纯流量卡能插随身wifi吗| 婷婷成人精品国产| 国产精品一区二区精品视频观看| 成人三级做爰电影| 母亲3免费完整高清在线观看| 精品国产乱码久久久久久男人| 女性生殖器流出的白浆| 热99久久久久精品小说推荐| 欧美国产精品va在线观看不卡| 制服人妻中文乱码| 亚洲国产欧美在线一区| 免费在线观看完整版高清| 丝袜在线中文字幕| 国产日韩欧美视频二区| 免费久久久久久久精品成人欧美视频| 国产欧美亚洲国产| 国产又爽黄色视频| 啦啦啦中文免费视频观看日本| 精品久久蜜臀av无| av又黄又爽大尺度在线免费看| 天天添夜夜摸| 精品国产一区二区三区四区第35| 亚洲专区国产一区二区| 真人做人爱边吃奶动态| 国产成人欧美在线观看 | 亚洲精品久久午夜乱码| 免费av中文字幕在线| tocl精华| 中国美女看黄片| 老鸭窝网址在线观看| 交换朋友夫妻互换小说| 国产免费福利视频在线观看| 大型av网站在线播放| 自线自在国产av| 丁香六月欧美| 高潮久久久久久久久久久不卡| 久久久久精品国产欧美久久久 | 亚洲天堂av无毛| 人成视频在线观看免费观看| 国产在线免费精品| 在线 av 中文字幕| 成人18禁高潮啪啪吃奶动态图| 国产男人的电影天堂91| 国产av又大| 蜜桃在线观看..| 老熟女久久久| 搡老熟女国产l中国老女人| 50天的宝宝边吃奶边哭怎么回事| 91大片在线观看| 精品一区二区三区av网在线观看 | 9热在线视频观看99| 亚洲欧洲精品一区二区精品久久久| 精品国内亚洲2022精品成人 | 日本av手机在线免费观看| 精品久久蜜臀av无| 久久精品久久久久久噜噜老黄| av网站免费在线观看视频| 人人妻人人澡人人爽人人夜夜| 80岁老熟妇乱子伦牲交| 午夜免费成人在线视频| 老熟女久久久| 久久久久久人人人人人| 国产高清视频在线播放一区 | 国产精品九九99| 人人妻人人澡人人爽人人夜夜| 国产一区二区三区综合在线观看| 岛国毛片在线播放| 日韩中文字幕欧美一区二区| 精品国产国语对白av| 两个人看的免费小视频| 青青草视频在线视频观看| 日日夜夜操网爽| 欧美日韩视频精品一区| 中文字幕另类日韩欧美亚洲嫩草| 亚洲av美国av| 91麻豆精品激情在线观看国产 | 老司机亚洲免费影院| 热99re8久久精品国产| 色精品久久人妻99蜜桃| 久久久国产精品麻豆| 精品一区二区三卡| 成年动漫av网址| 搡老岳熟女国产| 亚洲精品成人av观看孕妇| 亚洲专区字幕在线| 亚洲精品久久午夜乱码| 搡老乐熟女国产| 欧美日韩亚洲国产一区二区在线观看 | 91麻豆精品激情在线观看国产 | 天天躁夜夜躁狠狠躁躁| 妹子高潮喷水视频| 一边摸一边做爽爽视频免费| 手机成人av网站| 咕卡用的链子| 少妇猛男粗大的猛烈进出视频| 黑人欧美特级aaaaaa片| 免费av中文字幕在线| 精品欧美一区二区三区在线| 日本vs欧美在线观看视频| 少妇人妻久久综合中文| 精品熟女少妇八av免费久了| 99久久国产精品久久久| e午夜精品久久久久久久| 国产97色在线日韩免费| 色婷婷av一区二区三区视频| 亚洲国产欧美在线一区| 青青草视频在线视频观看| 少妇人妻久久综合中文| av超薄肉色丝袜交足视频| 久久久久国内视频| 制服诱惑二区| 免费看十八禁软件| 久久人妻福利社区极品人妻图片| 亚洲精品美女久久av网站| 久久毛片免费看一区二区三区| 人妻 亚洲 视频| 国产日韩欧美亚洲二区| av欧美777| 亚洲精品一二三| 国产精品久久久人人做人人爽| 久久国产精品男人的天堂亚洲| 免费少妇av软件| 99久久人妻综合| 女人被躁到高潮嗷嗷叫费观| 男女午夜视频在线观看| 成年动漫av网址| 在线观看一区二区三区激情| 日韩中文字幕欧美一区二区| 亚洲五月色婷婷综合| 日韩 亚洲 欧美在线| 国产精品99久久99久久久不卡| 国产精品熟女久久久久浪| 搡老熟女国产l中国老女人| 亚洲av片天天在线观看| 亚洲自偷自拍图片 自拍| 午夜91福利影院| 久久人人爽av亚洲精品天堂| av天堂久久9| 巨乳人妻的诱惑在线观看| 午夜两性在线视频| 建设人人有责人人尽责人人享有的| 狂野欧美激情性bbbbbb| 97精品久久久久久久久久精品| a 毛片基地| 黄色视频在线播放观看不卡| 久久久久网色| 曰老女人黄片| 亚洲欧美日韩另类电影网站| 欧美日韩亚洲高清精品| 成年av动漫网址| 91成人精品电影| 日韩欧美一区二区三区在线观看 | 女人久久www免费人成看片| 国产成人精品久久二区二区91| 90打野战视频偷拍视频| 国产精品九九99| 一二三四社区在线视频社区8| 欧美日韩亚洲综合一区二区三区_| 色精品久久人妻99蜜桃| 日韩大码丰满熟妇| 麻豆av在线久日| 欧美亚洲 丝袜 人妻 在线| 桃红色精品国产亚洲av| 老司机午夜十八禁免费视频| 亚洲熟女精品中文字幕| 久久久久国内视频| 精品久久久久久电影网| 欧美另类亚洲清纯唯美| 免费在线观看黄色视频的| 香蕉国产在线看| 成人黄色视频免费在线看| 国产成人精品无人区| 国产精品成人在线| 日本vs欧美在线观看视频| 亚洲精品国产一区二区精华液| 青草久久国产| 中文欧美无线码| 国产欧美日韩综合在线一区二区| 亚洲国产av影院在线观看| 国产区一区二久久| 桃花免费在线播放| 亚洲avbb在线观看| 伦理电影免费视频| 啦啦啦中文免费视频观看日本| 少妇猛男粗大的猛烈进出视频| 亚洲国产欧美网| 久久精品亚洲熟妇少妇任你| 手机成人av网站| 丰满饥渴人妻一区二区三| 国产成人欧美| 亚洲精品一区蜜桃| 国产精品久久久av美女十八| 80岁老熟妇乱子伦牲交| 人妻一区二区av| 国产亚洲欧美在线一区二区| 国产有黄有色有爽视频| 母亲3免费完整高清在线观看| 桃红色精品国产亚洲av| 99热国产这里只有精品6| av片东京热男人的天堂| 各种免费的搞黄视频| 国产成人啪精品午夜网站| 久久国产精品人妻蜜桃| 啦啦啦在线免费观看视频4| 亚洲成人手机| 手机成人av网站| 午夜免费成人在线视频| 国产一卡二卡三卡精品| 午夜福利在线免费观看网站| 一级,二级,三级黄色视频| 国产精品秋霞免费鲁丝片| 久久青草综合色| 成年人午夜在线观看视频| 99国产极品粉嫩在线观看| 19禁男女啪啪无遮挡网站| 亚洲精品美女久久久久99蜜臀| 十八禁网站网址无遮挡| 久久精品国产亚洲av高清一级| 欧美日韩一级在线毛片| 久久影院123| 少妇猛男粗大的猛烈进出视频| 91九色精品人成在线观看| 黑人巨大精品欧美一区二区mp4| 飞空精品影院首页| 欧美午夜高清在线| 国产免费av片在线观看野外av| h视频一区二区三区| 在线看a的网站| 日韩有码中文字幕| 亚洲精品国产色婷婷电影| 亚洲精品av麻豆狂野| 真人做人爱边吃奶动态| 一边摸一边做爽爽视频免费| 日韩欧美国产一区二区入口| 搡老岳熟女国产| 亚洲av日韩精品久久久久久密| 国产成人影院久久av| 自拍欧美九色日韩亚洲蝌蚪91| 黄色a级毛片大全视频| 99国产精品一区二区三区| 国产老妇伦熟女老妇高清| 亚洲欧美精品综合一区二区三区| 日韩欧美国产一区二区入口| 午夜免费观看性视频| 女性生殖器流出的白浆| 中亚洲国语对白在线视频| 最黄视频免费看| 日本猛色少妇xxxxx猛交久久| 久久ye,这里只有精品| 欧美人与性动交α欧美软件| 午夜精品国产一区二区电影| 十八禁网站网址无遮挡| 黄片大片在线免费观看| √禁漫天堂资源中文www| 岛国在线观看网站| 欧美亚洲日本最大视频资源| videos熟女内射| 国产欧美日韩一区二区三 | 母亲3免费完整高清在线观看| 国产麻豆69| 午夜成年电影在线免费观看| 9色porny在线观看| 不卡av一区二区三区| 一区二区三区激情视频| 91成人精品电影| 精品人妻在线不人妻| 国内毛片毛片毛片毛片毛片| 国产欧美日韩综合在线一区二区| 国产成人av激情在线播放| 亚洲av成人一区二区三| 欧美日韩av久久| 少妇精品久久久久久久| 成人亚洲精品一区在线观看| 国产成人影院久久av| 免费久久久久久久精品成人欧美视频| av电影中文网址| 丝袜喷水一区| 成在线人永久免费视频| 亚洲美女黄色视频免费看| 两人在一起打扑克的视频| 看免费av毛片| 黄色视频在线播放观看不卡| 亚洲少妇的诱惑av| 99国产精品免费福利视频| 国产免费福利视频在线观看| 视频区图区小说| 精品一区二区三区av网在线观看 | 美女主播在线视频| 久久久精品94久久精品| 又黄又粗又硬又大视频| 少妇粗大呻吟视频| 热99久久久久精品小说推荐| 欧美激情高清一区二区三区| av福利片在线| 丰满迷人的少妇在线观看| 老司机影院成人| 亚洲精品一区蜜桃| 一二三四社区在线视频社区8| av网站免费在线观看视频| 欧美黄色淫秽网站| 男男h啪啪无遮挡| 国产成人欧美| 亚洲成人国产一区在线观看| 久久99一区二区三区| 少妇人妻久久综合中文| 脱女人内裤的视频| 建设人人有责人人尽责人人享有的| 国产一区有黄有色的免费视频| 黄色视频不卡| 久久性视频一级片| 啦啦啦中文免费视频观看日本| 老汉色∧v一级毛片| 99re6热这里在线精品视频| 亚洲熟女精品中文字幕| 丁香六月欧美| 日韩欧美国产一区二区入口| 叶爱在线成人免费视频播放| 老司机午夜福利在线观看视频 | 18禁黄网站禁片午夜丰满| 午夜激情久久久久久久| 国产成人a∨麻豆精品| 建设人人有责人人尽责人人享有的| 亚洲第一青青草原| av网站免费在线观看视频| 男女国产视频网站| 97精品久久久久久久久久精品| 精品人妻在线不人妻| 国产精品av久久久久免费| 两人在一起打扑克的视频| 啦啦啦在线免费观看视频4| 国产亚洲av高清不卡| 亚洲 国产 在线| 精品久久久久久久毛片微露脸 | 久久久久国产一级毛片高清牌| 真人做人爱边吃奶动态| 欧美精品亚洲一区二区| 一区二区日韩欧美中文字幕| 最近最新免费中文字幕在线| 亚洲国产欧美在线一区| 不卡av一区二区三区| 久久亚洲精品不卡| a在线观看视频网站| 日本av手机在线免费观看| 欧美日韩亚洲高清精品| 搡老岳熟女国产| 精品免费久久久久久久清纯 | 一边摸一边做爽爽视频免费| 80岁老熟妇乱子伦牲交| 母亲3免费完整高清在线观看| 大码成人一级视频| 天堂俺去俺来也www色官网| 国产亚洲欧美在线一区二区| 精品欧美一区二区三区在线| 夜夜骑夜夜射夜夜干| 涩涩av久久男人的天堂| 叶爱在线成人免费视频播放| av片东京热男人的天堂| 成人三级做爰电影| 一边摸一边抽搐一进一出视频| 搡老岳熟女国产| 亚洲av日韩在线播放| 中文精品一卡2卡3卡4更新| 热99国产精品久久久久久7| 日本欧美视频一区| 中文欧美无线码| 国产免费视频播放在线视频| 99国产极品粉嫩在线观看| 少妇裸体淫交视频免费看高清 | 高清欧美精品videossex| 日韩欧美一区二区三区在线观看 | 日韩制服骚丝袜av| 一区二区三区乱码不卡18| 久久精品熟女亚洲av麻豆精品| 国产黄色免费在线视频| 丁香六月欧美| 丁香六月天网| 99国产精品一区二区蜜桃av | 精品亚洲成国产av| 肉色欧美久久久久久久蜜桃| 法律面前人人平等表现在哪些方面 | a级毛片黄视频| 国产伦理片在线播放av一区| 国产麻豆69| 老司机在亚洲福利影院| 人妻 亚洲 视频| 国产精品偷伦视频观看了| 精品高清国产在线一区| 午夜福利视频在线观看免费| 亚洲av电影在线进入| 成年女人毛片免费观看观看9 | 男女午夜视频在线观看| 精品国产一区二区三区久久久樱花| 亚洲免费av在线视频| 一区二区三区乱码不卡18| 国产真人三级小视频在线观看| 最近最新免费中文字幕在线| 国产日韩一区二区三区精品不卡| 亚洲精品国产区一区二| 90打野战视频偷拍视频| 少妇的丰满在线观看| 欧美精品亚洲一区二区| 人妻一区二区av| 黄色片一级片一级黄色片| 亚洲专区字幕在线| 日韩欧美一区视频在线观看| 亚洲人成电影观看| 国精品久久久久久国模美| 午夜免费成人在线视频| 高清黄色对白视频在线免费看| 亚洲avbb在线观看| 久久久精品区二区三区| 女性生殖器流出的白浆| 久久久久久亚洲精品国产蜜桃av| 在线 av 中文字幕| 国产伦理片在线播放av一区| 女警被强在线播放| 亚洲va日本ⅴa欧美va伊人久久 | 日韩视频一区二区在线观看| 自拍欧美九色日韩亚洲蝌蚪91| 久久国产精品影院| 亚洲精品国产精品久久久不卡| 亚洲国产av新网站| 9热在线视频观看99| 日本精品一区二区三区蜜桃| 免费少妇av软件| 亚洲性夜色夜夜综合| 亚洲精品国产色婷婷电影| 亚洲伊人色综图| 国产黄色免费在线视频| netflix在线观看网站| 欧美黄色片欧美黄色片| 久久这里只有精品19| 国产精品久久久久久精品电影小说| 人妻一区二区av| 狂野欧美激情性xxxx| 蜜桃国产av成人99| 最近最新中文字幕大全免费视频| 一级毛片女人18水好多| avwww免费| 精品国产乱码久久久久久小说| 王馨瑶露胸无遮挡在线观看| 免费黄频网站在线观看国产| 嫁个100分男人电影在线观看| 男男h啪啪无遮挡| 午夜福利一区二区在线看| 制服诱惑二区| 老司机影院成人| 日韩三级视频一区二区三区| www.精华液| 亚洲av成人不卡在线观看播放网 | 久久人人爽av亚洲精品天堂| 久久久久久人人人人人| 久久中文字幕一级| √禁漫天堂资源中文www| 少妇粗大呻吟视频| 亚洲精华国产精华精| 国产成人精品久久二区二区免费| 免费少妇av软件| 国产av又大| 丝袜脚勾引网站| 久久久国产精品麻豆| 国产人伦9x9x在线观看| 日本91视频免费播放| 又大又爽又粗| 国产成人精品久久二区二区免费| 巨乳人妻的诱惑在线观看| 欧美日韩视频精品一区| 一边摸一边做爽爽视频免费| 久久久久久久久久久久大奶| 亚洲国产精品成人久久小说| 俄罗斯特黄特色一大片| 高清在线国产一区| 99国产极品粉嫩在线观看| 黄片大片在线免费观看| 欧美人与性动交α欧美软件| 成人18禁高潮啪啪吃奶动态图| 又紧又爽又黄一区二区| 91成年电影在线观看| 中国国产av一级| av在线老鸭窝| 国产熟女午夜一区二区三区| 欧美日韩av久久| 人人妻人人澡人人爽人人夜夜| 多毛熟女@视频| 嫩草影视91久久| 超色免费av| 91九色精品人成在线观看| 久久精品国产a三级三级三级| 亚洲欧美精品自产自拍| 我的亚洲天堂| 日韩精品免费视频一区二区三区| 亚洲熟女精品中文字幕| 亚洲精品久久午夜乱码| 国产野战对白在线观看| 亚洲精品日韩在线中文字幕| 搡老乐熟女国产| 乱人伦中国视频| 嫁个100分男人电影在线观看| 国产成+人综合+亚洲专区| 捣出白浆h1v1| 欧美中文综合在线视频| 极品少妇高潮喷水抽搐| 成人黄色视频免费在线看| 777久久人妻少妇嫩草av网站| 国产精品一区二区在线观看99| 国产高清videossex| 国产一区有黄有色的免费视频| 成人黄色视频免费在线看| 777久久人妻少妇嫩草av网站| 国产精品一区二区在线观看99| 青春草亚洲视频在线观看| 日韩视频在线欧美| 国产麻豆69| 国产成人欧美| a在线观看视频网站| 亚洲一区中文字幕在线| 最近中文字幕2019免费版| 精品国内亚洲2022精品成人 | 久久人妻熟女aⅴ| 丝袜美足系列| 亚洲精品美女久久av网站| 50天的宝宝边吃奶边哭怎么回事| 日韩欧美一区视频在线观看| 久久久久国产一级毛片高清牌| 亚洲精品一区蜜桃| 久久精品久久久久久噜噜老黄| 国产无遮挡羞羞视频在线观看| 久久天躁狠狠躁夜夜2o2o| 美女国产高潮福利片在线看| 免费观看av网站的网址| 一本大道久久a久久精品| 99国产综合亚洲精品| 亚洲成人手机| 老汉色av国产亚洲站长工具| 久久久久久久久免费视频了| 1024香蕉在线观看| 欧美少妇被猛烈插入视频| 一区福利在线观看| 国产精品免费大片| 久久狼人影院| 高清av免费在线| 欧美在线黄色| 国产精品亚洲av一区麻豆| av在线app专区| 国产一区二区在线观看av| 日韩制服丝袜自拍偷拍| 一进一出抽搐动态| 天天躁夜夜躁狠狠躁躁| 欧美日韩一级在线毛片| 老熟妇仑乱视频hdxx| 午夜精品久久久久久毛片777| 动漫黄色视频在线观看| 1024视频免费在线观看| 日韩 亚洲 欧美在线| 久久久久久久久久久久大奶| 亚洲一区中文字幕在线| 亚洲av欧美aⅴ国产| 水蜜桃什么品种好| 老鸭窝网址在线观看| 欧美日韩av久久| 欧美在线一区亚洲| 久久久久网色| av网站免费在线观看视频| 精品福利观看| 电影成人av| 狂野欧美激情性bbbbbb| 午夜福利在线观看吧| 国产精品一区二区精品视频观看| 69av精品久久久久久 | 久久久久久久精品精品| 男女床上黄色一级片免费看| 久久久久久久久免费视频了| 一本一本久久a久久精品综合妖精| 成人影院久久| 亚洲精品一二三| 国产日韩一区二区三区精品不卡| 久久毛片免费看一区二区三区| 日本猛色少妇xxxxx猛交久久| 日本欧美视频一区| 久久久国产一区二区| 亚洲精品久久成人aⅴ小说| 国产成人啪精品午夜网站| 久热这里只有精品99| 久久久久精品国产欧美久久久 | 俄罗斯特黄特色一大片| 人妻人人澡人人爽人人| 日韩欧美免费精品| 欧美国产精品va在线观看不卡| 日本a在线网址| 亚洲国产精品成人久久小说| 亚洲国产欧美一区二区综合|