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

    熱梯度環(huán)境下梁高頻振動的能量流模型

    2022-07-04 07:20:06劉知輝牛軍川賈睿昊
    航空學(xué)報 2022年5期
    關(guān)鍵詞:熱應(yīng)力軸向密度

    劉知輝,牛軍川,2,*,賈睿昊

    1. 山東大學(xué) 機(jī)械工程學(xué)院,濟(jì)南 250061

    2. 山東大學(xué)高效潔凈機(jī)械制造教育部重點實驗室,濟(jì)南 250061

    航空領(lǐng)域中高速飛行器常工作于復(fù)雜、嚴(yán)苛的極端熱環(huán)境中,并承受來自動力源、氣動湍流、噪聲等多種形式的外界激勵。在這些外界振源的激勵下,結(jié)構(gòu)產(chǎn)生的振動會影響系統(tǒng)的平穩(wěn)運行,同時也對精密儀器儀表、電子設(shè)備等造成危害。此外,結(jié)構(gòu)振動引起的聲輻射還會影響設(shè)備中操作人員的舒適性。因此在系統(tǒng)設(shè)計階段針對結(jié)構(gòu)動力學(xué)特性進(jìn)行預(yù)測分析,從而實現(xiàn)結(jié)構(gòu)聲振特性的優(yōu)化設(shè)計具有重要的理論與工程價值。

    在高速飛行器等航空航天設(shè)備的聲振特性分析中,往往需考慮極端熱環(huán)境對結(jié)構(gòu)動態(tài)特性的影響。極端熱環(huán)境一方面可顯著改變結(jié)構(gòu)的材料屬性,從而影響結(jié)構(gòu)的剛度等關(guān)鍵動力學(xué)參數(shù);另一方面,處于極端熱環(huán)境中的結(jié)構(gòu)熱膨脹被約束時,結(jié)構(gòu)內(nèi)還將產(chǎn)生熱應(yīng)力。楊雄偉等以X43A超高速飛行器為例,研究了熱效應(yīng)對飛行器結(jié)構(gòu)各部分材料的物性及聲振耦合特性的影響規(guī)律。胡君逸等對熱環(huán)境下正交各向異性板的動力學(xué)特性進(jìn)行了理論和試驗研究,結(jié)果表明熱應(yīng)力和熱變形會顯著影響結(jié)構(gòu)的動力學(xué)特性。

    在航空航天領(lǐng)域,梁單元被廣泛用于高速飛行器中機(jī)翼、壁板、航空葉片等結(jié)構(gòu)的動力學(xué)建模分析。在這些應(yīng)用中,結(jié)構(gòu)常處于極端熱環(huán)境。Gu等建立了Timoshenko梁在熱環(huán)境中的動力學(xué)理論分析模型。Avsec和Oblak對自由和簡支邊界條件下的梁在熱環(huán)境下的振動進(jìn)行了研究。Chen等在考慮參數(shù)不確定性的情況下,使用基于間接Trefftz方法的波函數(shù)法(Wave Based Method,WBM)分析了均勻熱環(huán)境下梁的中低頻動力學(xué)響應(yīng)。在高速飛行器運行時,氣動湍流、外界噪聲等激勵往往頻率較高,且飛行器中結(jié)構(gòu)一般較為輕薄,易表現(xiàn)為短波高頻振動特性。上述以位移為基本變量的確定性分析方法由于計算量及參數(shù)敏感性等原因難以適用于高頻動力學(xué)分析。

    高頻分析時一般采用以聲振能量為基本變量的統(tǒng)計性分析方法。為考慮高速飛行器在極端溫度場中的高頻響應(yīng),陳強(qiáng)等建立了熱環(huán)境下的統(tǒng)計能量分析(Statistical Energy Analysis,SEA)模型,并就L型折板的長寬比對SEA模型關(guān)鍵參數(shù)的影響規(guī)律進(jìn)行了研究。胡婉璐等研究了軸向力對SEA模型中耦合梁耦合損耗因子的影響規(guī)律。

    近年來,以能量密度為基本變量的能量流方法在結(jié)構(gòu)高頻分析領(lǐng)域得到了廣泛應(yīng)用。相比于SEA,該方法可得到振動能量的空間分布信息,從而可提供更為詳細(xì)的高頻動力學(xué)預(yù)測結(jié)果。代文強(qiáng)等使用能量有限元方法對高速列車車內(nèi)的噪聲進(jìn)行了理論預(yù)測。王懷志等使用能量有限元法對衛(wèi)星的整流罩和儀器艙的中頻聲振響應(yīng)進(jìn)行了分析。陳兆林等使用能量有限元法和虛擬模態(tài)綜合法(EFEM-VMSS)分析了航天器分離等沖擊狀態(tài)下的瞬態(tài)響應(yīng)。滕曉艷等通過采用大阻尼假設(shè)建立了分析自由阻尼梁高頻響應(yīng)的解析模型。為實現(xiàn)對熱環(huán)境中梁的高頻分析,Zhang等考慮軸向熱應(yīng)力和材料屬性的變化,建立了均勻熱環(huán)境下梁的能量流模型。Wang等基于板的平面波疊加假設(shè)研究了板在熱環(huán)境中具有各向不同薄膜力時的能量流模型。

    飛行器在超音速飛行時,氣動加熱等效應(yīng)往往在壁板等結(jié)構(gòu)內(nèi)外表面產(chǎn)生較大的溫差,壁板等結(jié)構(gòu)的內(nèi)外表面之間也因熱傳導(dǎo)而存在較大的溫度梯度。此時結(jié)構(gòu)內(nèi)的材料參數(shù)以及熱應(yīng)力等將具有明顯的空間變化。而在上述熱環(huán)境高頻分析模型中,結(jié)構(gòu)被認(rèn)為處于均勻的溫度場中,此時建立的高頻分析模型也因此無法準(zhǔn)確反映溫度場的空間變化對結(jié)構(gòu)高頻振動的影響。

    二維壁板是飛行器結(jié)構(gòu)振動問題研究中廣泛采用的模型。當(dāng)考慮此類結(jié)構(gòu)的柱形彎曲(Cylindrical Bending)而忽略展向彎曲(Spanwise Bending)時,可采用梁單元進(jìn)行建模分析。為分析飛行器中此類結(jié)構(gòu)因氣動加熱等效應(yīng)而處于非均勻熱環(huán)境時的高頻動力學(xué)特性,建立熱梯度環(huán)境下梁的能量流解析模型。通過求解穩(wěn)態(tài)熱傳導(dǎo)方程得到梁內(nèi)溫度的空間分布,考慮溫度對材料屬性的改變及熱應(yīng)力,引入物理中性層以消除軸向與橫向彎曲變形之間的耦合。使用哈密頓原理導(dǎo)出梁在非均勻熱梯度環(huán)境下的運動控制方程,得到波動的頻散關(guān)系?;谖⒃w的功率平衡關(guān)系得到梁的能量流解析模型,并通過與基準(zhǔn)解的對比驗證其正確性。

    1 熱梯度環(huán)境下梁的模型

    1.1 梁中的溫度場

    圖1為處于熱梯度環(huán)境中的梁,其長度為,截面寬度為,截面高度為。梁的上下表面溫度分別為和。由于外界溫度在梁的上下表面存在差異,熱傳導(dǎo)將使梁在厚度方向上的溫度呈梯度分布。外界溫度在梁軸向上的變化率一般較小,因此不考慮軸向的溫度變化。

    圖1 處于熱梯度環(huán)境中的梁Fig.1 Beam exposed in thermal gradient environment

    首先對梁內(nèi)的溫度場進(jìn)行建模。以梁的幾何中性層為參考,坐標(biāo)為處的溫度滿足一維穩(wěn)態(tài)熱傳導(dǎo)方程

    (1)

    式中:()為坐標(biāo)處的材料熱傳導(dǎo)系數(shù),對于均質(zhì)材料,()為常數(shù);為溫度,K。式(1)的邊界條件可表示為

    (2)

    式中:為參考溫度,取為室溫300 K,此時材料不產(chǎn)生熱膨脹;Δ和Δ分別為梁的上下表面處溫度與參考溫度的差值。將式(2)代入式(1),對于均勻材料可求解得到在厚度方向的溫度為

    (3)

    根據(jù)Touloukian提出的高階多項式擬合材料屬性-溫度擬合關(guān)系,材料參數(shù)與溫度的關(guān)系可表示為

    =(+1+++)

    (4)

    式中:表示密度、彈性模量、泊松比、熱膨脹系數(shù)等材料物性參數(shù);(=-1,0,1,2,3)為材料參數(shù)的溫度系數(shù),需根據(jù)具體材料確定。在厚度方向的熱傳導(dǎo)導(dǎo)致溫度在厚度方向連續(xù)變化,由式(4) 可知,此時梁的材料參數(shù)在厚度方向也將呈現(xiàn)連續(xù)性變化。

    1.2 運動控制方程

    厚度方向存在的熱梯度使梁的材料屬性在空間中具有非均勻的特點。同時,若軸向變形被約束,熱膨脹還將產(chǎn)生軸向應(yīng)力。這兩方面因素使熱梯度環(huán)境下梁的運動控制方程異于常溫下梁的運動控制方程。

    已有的一些研究表明,基于等效單層理論(Equivalent Single Layer,ESL)建立非均勻材料結(jié)構(gòu)的位移場時,由于在厚度方向彈性模量等力學(xué)參數(shù)的空間變化,結(jié)構(gòu)橫向彎曲變形與軸向伸展之間將存在耦合。一些學(xué)者為此引入了物理中性層的概念以消除這一耦合。圖2中熱梯度環(huán)境下的梁由于在厚度方向的彈性模量也具有非均勻性,因此亦可使用物理中性層的概念建立熱梯度環(huán)境中梁的運動控制方程。

    圖2 梁在熱梯度環(huán)境中的變形Fig.2 Deformation of beam in thermal gradient environment

    當(dāng)梁的軸向跨度遠(yuǎn)大于截面高度時,一般可忽略其截面轉(zhuǎn)動慣量及剪切變形的影響而采用歐拉-伯努利梁對其進(jìn)行分析。在飛行器中,薄壁類壁板等結(jié)構(gòu)一般具有這一特性?;跉W拉-伯努利梁理論的位移場可表示為

    (5)

    式中:為梁的軸向位移;為彎曲變形對應(yīng)的橫向位移;為物理中性層與幾何中性層的距離,表示為

    (6)

    式中:()為幾何中性層坐標(biāo)下處的彈性模量,可通過式(3)和式(4)確定。一般情況下高溫將減小材料的彈性模量,因此對于圖2所示的熱梯度環(huán)境下的梁,物理中性層相對于幾何中性層將向溫度較低的方向偏移?;谑?5)給定的位移場,梁的彈性變形勢能可表示為

    (7)

    當(dāng)梁處于超靜定約束狀態(tài)、熱膨脹完全被約束時,由熱梯度環(huán)境產(chǎn)生的軸向熱應(yīng)力可表示為

    ()=-()()Δ()

    (8)

    式中:()為梁的熱膨脹系數(shù);Δ()為溫度差。

    對于小變形情況,變形后梁的微弧段長度d與未變形狀態(tài)長度d之間的差值可表示為

    (9)

    從而軸向上初始熱應(yīng)力對應(yīng)的勢能密度可表示為

    (10)

    由于熱應(yīng)力在結(jié)構(gòu)超靜定狀態(tài)下產(chǎn)生壓縮軸向力,式(10)中的勢能密度將是負(fù)的,這實質(zhì)上將減小梁的固有頻率。當(dāng)熱應(yīng)力使一階固有頻率為0時,梁將失穩(wěn)發(fā)生熱屈曲,此時梁將產(chǎn)生非線性變形。能量流分析是基于線性小變形假設(shè)的,因此不討論梁在熱梯度環(huán)境中處于后屈曲狀態(tài)時的能量流模型。

    基于等效單層理論,熱梯度環(huán)境下梁的動能可表示為

    (11)

    式中:()為梁的材料密度。

    由外力對處于熱梯度環(huán)境中的梁所做的功可表示為

    (12)

    式中:(,)為外力;為時間變量。

    當(dāng)外力為集中力形式時,可采用δ函數(shù)進(jìn)行描述。由式(7)、式(10)~式(12)可得熱梯度環(huán)境下梁彎曲變形的能量泛函。對能量泛函使用哈密頓原理并收集變分項δ的系數(shù),可得熱梯度環(huán)境下梁的彎曲振動控制方程:

    (13)

    式中:為在熱梯度環(huán)境下梁的截面有效彎曲剛度;為處于超靜定狀態(tài)的梁在熱環(huán)境中產(chǎn)生的有效軸向力;為等效截面密度。它們可分別表示為

    (14)

    (15)

    (16)

    當(dāng)梁處于均勻的溫度場時,式(13)將轉(zhuǎn)換為Zhang等使用的梁運動控制方程。進(jìn)一步地,當(dāng)不考慮梁的軸向熱應(yīng)力時,式(13)將退化為均勻梁的運動控制方程。

    式(15)中由熱應(yīng)力產(chǎn)生的軸向力是在默認(rèn)熱膨脹被完全約束時得到的。在實際工程中,薄壁結(jié)構(gòu)的屈服極限往往較小,為避免過大的溫度應(yīng)力引起結(jié)構(gòu)的熱屈曲,在實際工程中細(xì)長薄壁類結(jié)構(gòu)的兩端可能留有伸縮縫以削弱對熱膨脹的約束。如圖3所示,可使用端點的軸向彈簧模擬軸向伸縮縫的效應(yīng)。當(dāng)該軸向彈簧剛度系數(shù)為0時,將不存在針對軸向熱變形的約束,此時梁中也將不存在熱應(yīng)力,溫度將僅改變梁的材料屬性。當(dāng)彈簧剛度系數(shù)遠(yuǎn)大于軸向剛度時,熱膨脹可視為被完全約束從而產(chǎn)生熱應(yīng)力。其他剛度值可類似設(shè)定。

    圖3 軸向未完全約束的梁Fig.3 Beam without full axial constraint

    對于圖3中的梁,根據(jù)彈簧與梁的軸向力相等可得

    (17)

    式中:為軸向彈簧的剛度系數(shù);Δ為彈簧的伸縮量,Δ一般遠(yuǎn)小于梁的尺寸,因此不計這一微小變形對梁幾何尺寸的影響。通過在式(17)中求解Δ便可得到梁的軸向力。將代入式(13) 便可得此時對應(yīng)的運動控制方程。

    2 熱梯度環(huán)境下梁的能量流模型

    2.1 波動頻散關(guān)系

    在結(jié)構(gòu)的能量流研究中需建立能量密度控制方程,而能量密度控制方程基于振動的波動解。因此需要首先求解式(13)波動形式的解。

    為考慮在高頻時結(jié)構(gòu)阻尼的影響,可在勢能表達(dá)式中引入復(fù)剛度,此時導(dǎo)出的式(13)將被改寫為

    (18)

    式(18)具有時間簡諧和空間簡諧形式的波動解可表示為

    (,)=ejej

    (19)

    式中:為彈性波幅值;為波數(shù),表征彈性波的空間振蕩特性;為圓頻率,表征彈性波的時間振蕩特性。將式(19)代入運動控制方程式(18)對應(yīng)的齊次方程,可得

    (20)

    式(20)建立了波數(shù)與圓頻率之間波動頻散關(guān)系。通過求解式(20)便可得到熱梯度環(huán)境下彈性波時間振蕩與空間振蕩之間的關(guān)系。數(shù)學(xué)形式上式(20)為關(guān)于波數(shù)的四次多項式方程,共有4個根:

    (21)

    式中:為遠(yuǎn)場波波數(shù)的實部;為近場波波數(shù)的虛部。

    式(21)中的根可根據(jù)多項式求根公式得到:

    (22)

    (23)

    將式(21)中的前兩項代入式(19)后,得到以波數(shù)為頻率、在空間振蕩的遠(yuǎn)場波,式(21)中的后兩項對應(yīng)于指數(shù)衰減的近場波。能量流中一般不考慮近場波的影響,從而式(19)可表示為

    (,)=(ej+e-j)ej

    (24)

    式中:為第個波的幅值。

    由于使用復(fù)剛度模擬結(jié)構(gòu)遲滯阻尼效應(yīng),式(24) 中的波數(shù)也將變?yōu)閺?fù)數(shù),也即=+j?;谛∽枘峒僭O(shè)(?1),對式(22)以j為未知量、以0為原點進(jìn)行一階泰勒展開,可得

    (25)

    常見金屬材料的阻尼系數(shù)一般是小量,此時?1,因此式(25)中與同階的項可忽略,從而可得

    (26)

    (27)

    式中:和分別為復(fù)波數(shù)的實部和虛部。

    由式(26)可得彈性波的波長

    (28)

    式(26)表明在小阻尼假設(shè)下導(dǎo)出的表征波動空間振蕩特性的波數(shù)實部部分與無阻尼的波數(shù)一致,也即

    (29)

    為得到表征遠(yuǎn)場彈性波能量傳遞速度的群速度,對式(26)兩邊應(yīng)用隱函數(shù)求導(dǎo)法則,可得

    (30)

    式中:為彎曲波的群速度。

    類似的,由式(29)還可得

    (31)

    比較式(27)和式(30),表征波動隨空間衰減的項可表示為

    (32)

    2.2 能量密度控制方程

    時間周期平均后,熱梯度環(huán)境下梁的動能密度可表示為

    (33)

    將式(24)代入式(33),可得動能密度為

    (34)

    同理,由梁的軸向熱應(yīng)力以及彈性變形對應(yīng)的勢能密度表達(dá)式為

    (35)

    將式(24)代入式(35),忽略阻尼系數(shù)的高階項,可得勢能密度為

    (36)

    在式(34)和式(36)中,均包含指數(shù)衰減項和三角函數(shù)項。在物理意義上,前者對應(yīng)于遠(yuǎn)場波的衰減,而后者對應(yīng)于正負(fù)向遠(yuǎn)場波的干涉。在能量流分析中,基本變量為局部空間平均后的能量密度,因此式(34)和式(36)中的三角函數(shù)項在經(jīng)過波長局部空間平均后可被消除。注意到式(26) 給出的頻散關(guān)系,時間與空間平均后總的能量密度可表示為

    (37)

    圖4為熱梯度環(huán)境下梁的受力示意圖。振動功率流將同時由彎矩、剪力和熱應(yīng)力傳遞。時間周期平均意義下的振動能量流可表示為

    圖4 熱梯度環(huán)境下梁的受力示意圖Fig.4 Schematic diagram of forces acting on beam in thermal gradient environment

    (38)

    將式(24)給出的波動解代入式(38),進(jìn)一步經(jīng)局部空間平均以消除行波之間干涉項引起的功率流,可得時間與局部空間平均后的功率流表達(dá)式:

    (39)

    比較式(37)與式(39),可得

    (40)

    將式(31)和式(32)代入式(40),可得

    (41)

    在能量流分析中,能量密度控制方程可根據(jù)微元體控制體積內(nèi)的能量流平衡導(dǎo)出。如圖5所示,圖中為外部輸入的功率,對于熱梯度梁的任意微段d,在穩(wěn)態(tài)下能量輸入、能量耗散及能量傳導(dǎo)三者間具有平衡關(guān)系。對于采用的遲滯阻尼模型,Cremer等的研究表明時間周期內(nèi)對振動能量的平均耗散功率可表示為

    圖5 微元體的能量流平衡Fig.5 Energy flow balance in a differential element

    (42)

    由式(29)、式(34)和式(36)可知,時間和局部空間平均后勢能密度與動能密度相等。這一性質(zhì)在Navazi等進(jìn)行的試驗中也得到過驗證?;诖?式(42)可表示為

    (43)

    考慮圖5所示的微元體的能量平衡關(guān)系,可得

    (44)

    將式(41)代入式(44),可得

    (45)

    式(45)即為建立的能量密度控制方程。該方程與Wohlever和Bernhard建立的常溫環(huán)境下梁的能量密度控制方程具有一致的形式。但式(45) 的導(dǎo)出同時考慮了熱環(huán)境產(chǎn)生的熱應(yīng)力以及材料的溫度相關(guān)性,因此式(45)中的關(guān)鍵參數(shù)群速度將與未考慮熱環(huán)境影響的群速度不同。

    2.3 能量密度控制方程的求解

    在高頻分析時,結(jié)構(gòu)的輸入點導(dǎo)納一般可在忽略邊界條件的情況下使用無限結(jié)構(gòu)的輸入導(dǎo)納進(jìn)行近似。Chen等使用波動法精確求解了預(yù)應(yīng)力梁在兩端簡支條件下的輸入導(dǎo)納,結(jié)果表明考慮具體邊界條件得到的輸入導(dǎo)納在無限結(jié)構(gòu)導(dǎo)納上下隨頻率振蕩,且二者隨激勵頻率的增加逐漸接近。高頻分析一般在倍頻程下以頻率平均意義進(jìn)行,因此在高頻時使用無限結(jié)構(gòu)的導(dǎo)納可以較好地反映在頻率平均時外界激勵對結(jié)構(gòu)統(tǒng)計意義下的輸入功率。

    圖6所示為忽略邊界時處于熱梯度環(huán)境中的梁及梁中的彈性波,圖中為外部集中形式的激振力,和分別為左側(cè)半梁和右側(cè)半梁的彎矩,和分別為左側(cè)半梁和右側(cè)半梁的剪力。由于忽略了邊界的存在,以激勵位置為分割點,兩側(cè)的梁均只有由外界力激起的遠(yuǎn)離激勵點的彈性波。忽略時間簡諧項,左側(cè)半梁(梁1)和右側(cè)半梁(梁2)中的彈性波位移可分別表示為

    圖6 無限尺寸熱梯度梁中的彈性波Fig.6 Elastic waves in infinite beam in thermal gradient environment

    (46)

    式中:和分別為左側(cè)半梁和右側(cè)半梁的橫向位移;和分別為左側(cè)半梁遠(yuǎn)場波和近場波的幅值;和分別為右側(cè)半梁遠(yuǎn)場波和近場波的幅值。

    式(46)共有4個未知數(shù),由激勵位置處的連續(xù)性條件可得到4個協(xié)調(diào)方程。

    在計算得到式(46)中的4個參數(shù)后,使用輸入導(dǎo)納的定義可得處于熱梯度環(huán)境下無限梁的輸入導(dǎo)納為

    (47)

    由式(47)給出的輸入點導(dǎo)納可得輸入功率:

    (48)

    將式(48)確定的外界激勵輸入功率流代入式(45) 便可得梁上的能量密度的空間分布。式(45) 在數(shù)學(xué)形式上為定義在空間域上的二階常微分方程,它的解可表示為

    (49)

    (50)

    式中:為復(fù)平面內(nèi)的積分變量。

    由式(50)可見式(49)中的積分在復(fù)平面內(nèi)有兩個留數(shù)點。由留數(shù)定理可得式(45)的特解為

    (51)

    將式(51)代入式(49)便可得熱梯度環(huán)境下梁上的能量密度分布。對于單梁,在梁的左右兩端應(yīng)不存在能量流的輸入輸出,因此在左右兩端的邊界條件可表示為

    (52)

    式(52)中的能量流可通過式(41)確定。由式(52) 中的兩個邊界條件便可確定兩個未知系數(shù)和。基于式(19)~式(52)得到的是能量密度控制方程的解析解,當(dāng)使用有限元方法進(jìn)行數(shù)值求解時,該模型也被稱為能量有限元方法。相比于數(shù)值解,解析形式的能量流模型在參數(shù)化分析和優(yōu)化時將具有優(yōu)勢。

    3 數(shù)值算例與討論

    對建立的能量流模型進(jìn)行驗證,并通過數(shù)值算例對處于非均勻熱環(huán)境中梁的高頻振動特性進(jìn)行研究。算例中采用的梁材料為Ti-6Al-4V。Ti-6Al-4V具有良好的低溫和高溫性能,可穩(wěn)定工作于400~500 ℃的環(huán)境下。溫度對這種材料彈性模量和熱膨脹系數(shù)的影響如表1所示,材料密度為4 439 kg/m,不隨溫度變化。

    表1 Ti-6Al-4V隨溫度變化的參數(shù)Table 1 Temperature-dependent coefficients of Ti-6Al-4V

    3.1 模型有效性驗證

    首先對建立的能量流模型的有效性進(jìn)行驗證。基于能量泛函和多項式形函數(shù)可得熱梯度環(huán)境下梁的單元矩陣,通過組裝可得有限元模型。在驗證對比時使用有限元方法在稠密網(wǎng)格時得到的結(jié)果作為參考對比值。為確保有限元結(jié)果的收斂性,在有限元中設(shè)定網(wǎng)格長度小于1/10波長,因此有限元計算所需的自由度和計算量遠(yuǎn)遠(yuǎn)大于能量流模型。

    能量流模型對結(jié)構(gòu)高頻振動響應(yīng)的預(yù)測需要首先實現(xiàn)對外界激勵功率輸入的估計。在式(48)中使用無限結(jié)構(gòu)的導(dǎo)納求取外界激勵的功率輸入,圖7使用有限元法并劃分稠密網(wǎng)格精確求解了幾種常見邊界條件下單位幅值簡諧力作用于梁中點時引起的功率流輸入,梁的長度為3 m,截面寬度和高度分別為5×10m和1×10m,材料阻尼=0.05且不隨溫度變化。梁的上表面處于溫度為500 K的熱環(huán)境中,下表面溫度為300 K,軸向熱膨脹約束彈簧剛度為1 000 N/m。

    由圖7可見由無限結(jié)構(gòu)導(dǎo)納法得到的輸入功率流在頻率內(nèi)較為平滑,且整體上隨激勵頻率升高而下降。而考慮邊界條件時得到的精確輸入功率流則在前者上下振蕩。這是由于無限結(jié)構(gòu)導(dǎo)納法不考慮邊界,也因此不存在由駐波效應(yīng)引起的共振和反共振現(xiàn)象。從圖7中還可看出雖然考慮具體邊界條件時輸入功率流在頻域內(nèi)有一定振蕩,但它們總體是在無限結(jié)構(gòu)輸入功率的上下波動。同時隨頻率升高考慮邊界條件得到的輸入功率逐漸接近無限結(jié)構(gòu)的輸入功率。這一現(xiàn)象可通過邊界對波傳遞的影響解釋。對于保守邊界,無論其具體形式,波在傳遞至邊界而反射時只產(chǎn)生相位的變化而不產(chǎn)生幅值的變化。隨波長變小,相位變化的影響逐漸變小,且這一影響往往可經(jīng)頻率平均進(jìn)一步進(jìn)行消除。因此,式(48)可較為合理地預(yù)測各種邊界條件下高頻激勵時外界的功率輸入。

    圖7 外界功率流隨激勵頻率的變化Fig.7 Variation of power input with excitation frequency

    圖8為波長與波數(shù)隨激勵頻率的變化曲線,可見隨激勵頻率增加,波長逐漸變短,而波數(shù)逐漸增大;隨波長變短,近場波的影響范圍將逐漸減?。浑S著波數(shù)的增加,波動在空間中的振蕩頻率將更快,此時進(jìn)行局部空間平均以消除能量密度和能量流中與波數(shù)同頻的振蕩項也將是合理的。

    圖8 波長與波數(shù)隨激勵頻率變化Fig.8 Variations of wavenumber and wavelength with excitation frequency

    圖9為具有單位幅值的簡諧力作用于梁中點時,在不同激勵頻率下能量密度的空間分布。在圖9中,使用具有稠密網(wǎng)格的有限元結(jié)果作為對比參考值。有限元計算在以激勵頻率為中心頻率的1/3倍頻程內(nèi)進(jìn)行,并取1/3倍頻程內(nèi)100個頻率對應(yīng)結(jié)果的均值為基準(zhǔn)值。由圖9可見能量流模型預(yù)測的能量密度結(jié)果在空間內(nèi)較為平滑,而基于確定性有限元方法所得的結(jié)果在空間中呈現(xiàn)振蕩。能量流模型的結(jié)果可較為準(zhǔn)確地反映1/3倍頻程內(nèi)響應(yīng)的空間整體變化趨勢,且與基準(zhǔn)值具有較好的吻合性。能量流模型的預(yù)測結(jié)果在激勵位置及邊界處與有限元結(jié)果具有一定的差別,這一現(xiàn)象是由于激勵位置和邊界處存在近場波。由圖8可知隨激勵頻率增加波長逐漸減小,從而近場波存在的范圍也越來越小,這一現(xiàn)象也可在圖9中觀察到。在圖9中隨激勵頻率增加,邊界處振蕩的范圍與幅度都趨向減小。

    圖9 不同激勵頻率下梁振動能量密度分布情況Fig.9 Energy density distributions of beam at different excitation frequencies

    圖10為不同激勵頻率下梁中的功率流。由圖7可知隨激勵頻率增加,由簡諧激勵力引起的功率輸入逐漸減小。由于激勵處于對稱位置,因此梁上的功率流也呈對稱性分布,從激勵位置處向邊界傳播,并在傳播過程中衰減。

    圖10 不同激勵頻率f下梁中的功率流Fig.10 Power flow in beam at different excitation frequency f

    由式(43)中遲滯阻尼模型可知,在相同的功率輸入下,總體的振動能量密度將反比于激勵頻率。因此對比圖9中3種激勵頻率下梁的能量密度響應(yīng)可見,隨著激勵頻率的增加梁的整體振動能量響應(yīng)減小。同時由于波長減小,輸入功率流向邊界時將在空間內(nèi)經(jīng)歷更多的波長周期衰減,因此在圖9中隨激勵頻率的增加能量密度在空間中的衰減趨勢也更為明顯。

    圖11為激勵位置在1/3梁長度(激勵點坐標(biāo)=1,坐標(biāo)系原點位于梁左端點處)、激勵頻率為4 kHz時梁上的能量密度分布。其他參數(shù)與圖7中的算例相同。由圖11可見對于非對稱激勵,能量流模型給出的結(jié)果仍能與參考值較好吻合。在靠近激勵位置的左側(cè)邊界由于近場波效應(yīng)明顯,基準(zhǔn)值具有較大的振蕩性。但總體來說,能量流模型的結(jié)果較好地反映了振動能量的空間分布變化趨勢。能量流模型的準(zhǔn)確性依賴于對外界功率輸入的準(zhǔn)確估計。對于高頻激勵,式(48)基于無限結(jié)構(gòu)導(dǎo)納給出的功率輸入計算公式在激勵位置距離邊界數(shù)個波長時便較為準(zhǔn)確,而波長在高頻時較短,因此對于梁上的絕大多數(shù)位置式(48)均能給出合理的輸入功率流。對于圖11 中的算例,雖然激勵處在非對稱位置,但其距離邊界位置仍有若干波長,因此能量流模型也較準(zhǔn)確地得到了梁上的振動能量分布特性。

    圖11 非對稱激勵時梁中的能量密度Fig.11 Energy density in beam with asymmetrical excitation

    由上述對比可見,能量流模型的結(jié)果較為準(zhǔn)確地反映了處于熱梯度環(huán)境下的梁在承受高頻激勵時的振動能量分布特性。

    3.2 熱環(huán)境的影響

    高速飛行器等設(shè)備在運行中可能處于溫度不斷變化的熱環(huán)境中。溫度的變化將直接引起結(jié)構(gòu)剛度及熱應(yīng)力的變化。圖12為梁的上表面溫度變化時,截面剛度與軸向熱應(yīng)力的變化趨勢,其中下表面溫度被設(shè)定為=0.6,激勵頻率為4 kHz。由圖12可見隨上表面溫度增加,截面剛度逐漸減小而軸向熱應(yīng)力變大。本質(zhì)上這兩種變化趨勢都將降低系統(tǒng)的固有頻率,也即使得梁的振動更趨向于高頻振動。

    圖12 溫度對彎曲剛度和軸向力的影響Fig.12 Influences of temperature on bending stiffness and axial force

    圖13為上表面溫度變化時波長和輸入功率的變化曲線。由圖13可見隨溫度增加,輸入功率也增加,而波長將變短。由圖13還可看出,對于不同的軸向約束彈簧剛度系數(shù),波長與輸入功率基本相同。這表明此時軸向力的大小對高頻振動特性的影響較小。對于細(xì)長梁這類薄壁結(jié)構(gòu)而言,其臨界軸向壓應(yīng)力一般較小。Chen等的研究表明,較小的軸向力對高階固有頻率影響不明顯,因此對高頻響應(yīng)的影響也較小。

    圖13 溫度對波長和輸入功率的影響Fig.13 Influences of temperature on wavelength and power input

    圖14為簡支梁在僅考慮軸向熱應(yīng)力和僅材料變化時時固有頻率與常溫下的梁的固有頻率平方之比。固有頻率使用經(jīng)典的分離變量法由式(13) 得到。由圖14可知,當(dāng)僅考慮軸向力時,低頻模態(tài)相比于常溫下的梁具有明顯的改變,熱梯度環(huán)境引起的軸向力較大幅度地降低了低階固有頻率,但其對高頻模態(tài)的影響并不顯著。由于高頻響應(yīng)一般由激勵頻率倍頻程內(nèi)的若干高階模態(tài)決定,因此軸向熱應(yīng)力的影響將不顯著。而非均勻熱環(huán)境造成的材料變化引起了剛度的降低,從整體上降低了固有頻率,因此其對全頻段的響應(yīng)均將有明顯影響。圖14中的趨勢雖是基于簡支邊界條件,但高頻時結(jié)構(gòu)響應(yīng)對邊界并不敏感,因此上述討論對于任意邊界具有一般性。

    圖14 材料變化與熱應(yīng)力對固有頻率的影響Fig.14 Influences of material variation and thermal stress on natural frequencies

    圖15為不同上表面溫度時梁的振動能量密度分布。分析頻率為4 kHz,=1 000 N/m,上表面溫度分別設(shè)置為30 ℃和500 ℃。

    圖15 溫度對梁上能量密度的影響Fig.15 Influence of temperature on energy distribution of beam

    由圖15可見溫度的增加總體上提高了梁上的振動能量,同時在溫度較高時梁上的能量密度在空間中的衰減更為明顯。這一現(xiàn)象對應(yīng)于圖13 所得結(jié)論。隨溫度升高外界激勵力對梁的功率輸入增大,因此梁上的總體振動能量水平增高,而較小的波長使振動能量在空間中的衰減更為明顯。

    4 結(jié) 論

    通過求解穩(wěn)態(tài)熱傳導(dǎo)方程得到了熱梯度下梁的內(nèi)部溫度場,考慮溫度場對材料物性參數(shù)的改變及產(chǎn)生的熱應(yīng)力,引入物理中性層概念消除軸向拉伸與橫向彎曲變形之間的耦合,進(jìn)而基于歐拉梁理論建立了運動控制方程,并詳細(xì)推導(dǎo)了能量流模型及其解析解。通過數(shù)值算例將能量流模型與基準(zhǔn)解進(jìn)行對比,結(jié)果表明:

    1) 基于無限結(jié)構(gòu)導(dǎo)納得到的功率輸入能較好地評估外部高頻激勵對熱梯度梁的功率流輸入。

    2) 與基準(zhǔn)解的對比表明,建立的能量流模型能以較小的計算量準(zhǔn)確預(yù)測熱梯度環(huán)境中的梁在高頻激勵下的振動能量分布情況。

    3) 熱梯度環(huán)境對梁高頻振動特性的改變主要是由于結(jié)構(gòu)的材料屬性在非均勻溫度場中產(chǎn)生了變化,而熱應(yīng)力對高頻振動產(chǎn)生的影響較小。

    4) 隨外界溫度升高,非均勻熱環(huán)境整體上使熱梯度梁中彎曲波波長變短,外界激勵力輸入功率增加,從而振動能量具有更明顯的空間衰減。

    猜你喜歡
    熱應(yīng)力軸向密度
    大型立式單級引黃離心泵軸向力平衡的研究
    『密度』知識鞏固
    密度在身邊 應(yīng)用隨處見
    WNS型鍋爐煙管管端熱應(yīng)力裂紋原因分析
    “玩轉(zhuǎn)”密度
    密度應(yīng)用知多少
    荒銑加工軸向切深識別方法
    采用單元基光滑點插值法的高溫管道熱應(yīng)力分析
    微小型薄底零件的軸向車銑實驗研究
    基于流熱固耦合的核電蒸汽發(fā)生器傳熱管熱應(yīng)力數(shù)值模擬
    香蕉精品网在线| 日韩av在线免费看完整版不卡| kizo精华| 亚洲,欧美,日韩| 中文资源天堂在线| 亚洲精品第二区| 99久国产av精品国产电影| 久久精品夜色国产| 国产人妻一区二区三区在| 香蕉精品网在线| 欧美日韩亚洲高清精品| 午夜精品国产一区二区电影 | 日本黄色片子视频| 超碰97精品在线观看| 午夜激情福利司机影院| 乱系列少妇在线播放| 91精品伊人久久大香线蕉| 好男人在线观看高清免费视频| 校园人妻丝袜中文字幕| 色吧在线观看| 一本色道久久久久久精品综合| www.色视频.com| 精品午夜福利在线看| 亚洲av免费在线观看| 菩萨蛮人人尽说江南好唐韦庄| 国产成年人精品一区二区| 一级a做视频免费观看| 国产成人精品婷婷| 男女国产视频网站| 大片电影免费在线观看免费| 一级a做视频免费观看| 日本wwww免费看| 午夜福利在线观看免费完整高清在| 国产探花极品一区二区| 国产一区二区三区综合在线观看 | 性插视频无遮挡在线免费观看| 一区二区三区乱码不卡18| 国产精品一区www在线观看| 欧美日韩亚洲高清精品| 91精品国产九色| 性色av一级| 麻豆乱淫一区二区| 成年免费大片在线观看| 精品人妻一区二区三区麻豆| 国产欧美日韩精品一区二区| 青春草视频在线免费观看| 美女被艹到高潮喷水动态| 日韩在线高清观看一区二区三区| 搞女人的毛片| 久久久久久久国产电影| 韩国av在线不卡| 91久久精品国产一区二区成人| 97在线视频观看| 亚洲av福利一区| 免费电影在线观看免费观看| 精品久久久久久电影网| 亚洲av成人精品一二三区| 亚洲天堂av无毛| 国产精品一区二区性色av| 岛国毛片在线播放| 欧美激情国产日韩精品一区| 午夜视频国产福利| 男女那种视频在线观看| 97在线人人人人妻| 亚洲不卡免费看| 中文字幕亚洲精品专区| 欧美成人精品欧美一级黄| 亚洲欧美日韩另类电影网站 | 成人毛片60女人毛片免费| 免费看光身美女| av在线亚洲专区| 美女xxoo啪啪120秒动态图| 搡老乐熟女国产| 黄色怎么调成土黄色| 色播亚洲综合网| 亚洲一区二区三区欧美精品 | 青春草亚洲视频在线观看| 噜噜噜噜噜久久久久久91| 嫩草影院精品99| 97超视频在线观看视频| 亚洲精品,欧美精品| 日韩欧美 国产精品| 亚洲成人av在线免费| 国产成人一区二区在线| 日韩,欧美,国产一区二区三区| 80岁老熟妇乱子伦牲交| 亚洲人成网站高清观看| 欧美极品一区二区三区四区| 在线 av 中文字幕| 一区二区三区免费毛片| 精品少妇久久久久久888优播| 激情 狠狠 欧美| 久久久国产一区二区| 国产国拍精品亚洲av在线观看| 亚洲综合色惰| 哪个播放器可以免费观看大片| 日日啪夜夜撸| 日韩人妻高清精品专区| 午夜日本视频在线| 亚洲精品亚洲一区二区| 日韩av不卡免费在线播放| 精品国产露脸久久av麻豆| 日日啪夜夜撸| 哪个播放器可以免费观看大片| 久久久久久伊人网av| 亚洲欧洲国产日韩| 日本午夜av视频| 九九久久精品国产亚洲av麻豆| 交换朋友夫妻互换小说| av.在线天堂| 亚洲经典国产精华液单| 欧美区成人在线视频| 国产高清不卡午夜福利| 黄色配什么色好看| 国产探花极品一区二区| 晚上一个人看的免费电影| 建设人人有责人人尽责人人享有的 | 成人国产麻豆网| 日本色播在线视频| 国产女主播在线喷水免费视频网站| 国产老妇女一区| 在线a可以看的网站| 色网站视频免费| 永久网站在线| 一级毛片黄色毛片免费观看视频| 高清午夜精品一区二区三区| 久久女婷五月综合色啪小说 | 欧美精品一区二区大全| 免费大片黄手机在线观看| 97超碰精品成人国产| 精品少妇黑人巨大在线播放| 波野结衣二区三区在线| 国产黄片美女视频| 99视频精品全部免费 在线| 国内精品美女久久久久久| 99久久九九国产精品国产免费| 内射极品少妇av片p| 免费av观看视频| 国产在线男女| 国产在线男女| 国产午夜精品一二区理论片| 日本熟妇午夜| 在线精品无人区一区二区三 | 久久精品夜色国产| 麻豆乱淫一区二区| 高清毛片免费看| 亚洲av中文字字幕乱码综合| 中文资源天堂在线| 网址你懂的国产日韩在线| 特级一级黄色大片| 国产在线男女| 内射极品少妇av片p| 精品少妇黑人巨大在线播放| 最近中文字幕2019免费版| 高清午夜精品一区二区三区| 国产永久视频网站| 自拍偷自拍亚洲精品老妇| 婷婷色综合大香蕉| 国产一区二区在线观看日韩| 99热网站在线观看| 99热国产这里只有精品6| 国产免费一级a男人的天堂| 亚洲欧美精品自产自拍| 日本色播在线视频| 欧美日韩在线观看h| 成人黄色视频免费在线看| 99热国产这里只有精品6| 菩萨蛮人人尽说江南好唐韦庄| 免费播放大片免费观看视频在线观看| 国产精品久久久久久精品电影| 国产精品国产三级国产av玫瑰| xxx大片免费视频| 亚洲av免费在线观看| 日本wwww免费看| 中文精品一卡2卡3卡4更新| 国国产精品蜜臀av免费| 久久99热这里只频精品6学生| 日韩成人av中文字幕在线观看| 国模一区二区三区四区视频| 国产美女午夜福利| 国产精品99久久99久久久不卡 | 亚洲一级一片aⅴ在线观看| 精品久久久久久电影网| 看非洲黑人一级黄片| 久久久成人免费电影| 国产男人的电影天堂91| 男女边吃奶边做爰视频| 嫩草影院新地址| 日本一二三区视频观看| 国产欧美日韩一区二区三区在线 | 99热这里只有是精品在线观看| 国产精品久久久久久精品古装| 九九爱精品视频在线观看| 国产在线一区二区三区精| 午夜精品一区二区三区免费看| 精品久久久久久久人妻蜜臀av| 日韩成人伦理影院| 成人高潮视频无遮挡免费网站| 白带黄色成豆腐渣| 精品人妻偷拍中文字幕| 国产成年人精品一区二区| 亚洲精品国产成人久久av| 嫩草影院新地址| 欧美xxxx黑人xx丫x性爽| 在现免费观看毛片| 青春草国产在线视频| 如何舔出高潮| 涩涩av久久男人的天堂| 欧美97在线视频| 在线播放无遮挡| 国产精品一区二区三区四区免费观看| 99精国产麻豆久久婷婷| 日韩欧美 国产精品| 少妇人妻一区二区三区视频| 久久久久久伊人网av| 精品久久久久久久久av| 天美传媒精品一区二区| 一本久久精品| 欧美xxxx性猛交bbbb| 久久久久久久大尺度免费视频| 中文在线观看免费www的网站| 精品久久久久久久人妻蜜臀av| 免费看光身美女| 免费观看av网站的网址| 午夜福利在线在线| 亚洲av成人精品一二三区| 超碰av人人做人人爽久久| 春色校园在线视频观看| 午夜精品国产一区二区电影 | 岛国毛片在线播放| 久久99精品国语久久久| 国产欧美另类精品又又久久亚洲欧美| 国产一区二区亚洲精品在线观看| av国产久精品久网站免费入址| 秋霞在线观看毛片| 国产色爽女视频免费观看| 白带黄色成豆腐渣| 亚洲熟女精品中文字幕| 色吧在线观看| 日本爱情动作片www.在线观看| 丝袜脚勾引网站| 亚洲不卡免费看| 免费看不卡的av| 亚洲成人av在线免费| 国产亚洲5aaaaa淫片| 中文精品一卡2卡3卡4更新| 亚洲婷婷狠狠爱综合网| 国产黄色视频一区二区在线观看| 亚洲不卡免费看| 成人午夜精彩视频在线观看| 久热这里只有精品99| 精品国产一区二区三区久久久樱花 | 亚洲天堂av无毛| 五月玫瑰六月丁香| 天美传媒精品一区二区| 丝袜脚勾引网站| 久久久久久久国产电影| 亚洲精品一二三| 制服丝袜香蕉在线| 亚洲av成人精品一二三区| 亚洲精品国产色婷婷电影| 亚洲婷婷狠狠爱综合网| 国产中年淑女户外野战色| 亚洲精品一区蜜桃| 中文字幕人妻熟人妻熟丝袜美| tube8黄色片| 国产高清不卡午夜福利| 日产精品乱码卡一卡2卡三| 国产欧美另类精品又又久久亚洲欧美| 69av精品久久久久久| 国产欧美日韩一区二区三区在线 | 国产日韩欧美在线精品| 波野结衣二区三区在线| 日本午夜av视频| 久久久久九九精品影院| 99热这里只有是精品50| 男人舔奶头视频| 激情五月婷婷亚洲| 国产精品一区www在线观看| 欧美高清成人免费视频www| av在线观看视频网站免费| 国产综合精华液| 最近的中文字幕免费完整| 黄色欧美视频在线观看| 欧美精品一区二区大全| 久久影院123| 成人亚洲精品av一区二区| 精品酒店卫生间| 精品久久久久久电影网| 大香蕉97超碰在线| 国产伦精品一区二区三区四那| 亚洲成人精品中文字幕电影| av国产精品久久久久影院| a级毛色黄片| 80岁老熟妇乱子伦牲交| 永久网站在线| 国产精品久久久久久精品电影| 91久久精品国产一区二区成人| a级毛色黄片| 在线播放无遮挡| 中文字幕av成人在线电影| 激情五月婷婷亚洲| 简卡轻食公司| 精品亚洲乱码少妇综合久久| 91精品一卡2卡3卡4卡| 纵有疾风起免费观看全集完整版| 亚洲国产av新网站| 久久精品久久精品一区二区三区| 狂野欧美白嫩少妇大欣赏| 中文字幕人妻熟人妻熟丝袜美| 看免费成人av毛片| 日产精品乱码卡一卡2卡三| 欧美一区二区亚洲| 人妻系列 视频| 午夜免费观看性视频| 两个人的视频大全免费| 菩萨蛮人人尽说江南好唐韦庄| 亚洲色图综合在线观看| 国产美女午夜福利| 免费人成在线观看视频色| 一级毛片黄色毛片免费观看视频| 黄片wwwwww| 青青草视频在线视频观看| 最近的中文字幕免费完整| 成人综合一区亚洲| 你懂的网址亚洲精品在线观看| 黄片wwwwww| 亚洲高清免费不卡视频| 国产色爽女视频免费观看| 在线看a的网站| 国产av不卡久久| 免费在线观看成人毛片| 亚洲精华国产精华液的使用体验| 汤姆久久久久久久影院中文字幕| 欧美少妇被猛烈插入视频| 香蕉精品网在线| 午夜福利在线观看免费完整高清在| 中国国产av一级| 免费黄网站久久成人精品| 亚洲av欧美aⅴ国产| 一级片'在线观看视频| 五月天丁香电影| 欧美一级a爱片免费观看看| 国产免费一级a男人的天堂| 日韩av免费高清视频| 丝袜脚勾引网站| 国产成人午夜福利电影在线观看| 午夜福利在线观看免费完整高清在| 爱豆传媒免费全集在线观看| 春色校园在线视频观看| tube8黄色片| 日本-黄色视频高清免费观看| 久久久久久伊人网av| 性插视频无遮挡在线免费观看| 中文字幕免费在线视频6| 国产毛片a区久久久久| 少妇丰满av| 久久韩国三级中文字幕| 国产91av在线免费观看| 久久女婷五月综合色啪小说 | 18禁裸乳无遮挡动漫免费视频 | av福利片在线观看| 国产成人91sexporn| 一级毛片我不卡| 少妇的逼水好多| 欧美日韩视频精品一区| 欧美成人一区二区免费高清观看| 久久精品熟女亚洲av麻豆精品| 亚洲精品视频女| 一二三四中文在线观看免费高清| 亚洲国产精品国产精品| 亚洲欧美精品自产自拍| 精品午夜福利在线看| 最近手机中文字幕大全| 日日啪夜夜爽| 亚洲欧美一区二区三区黑人 | 各种免费的搞黄视频| 人妻少妇偷人精品九色| 亚洲精华国产精华液的使用体验| 亚洲成人中文字幕在线播放| 国产69精品久久久久777片| 丝瓜视频免费看黄片| 少妇高潮的动态图| 亚洲精品乱久久久久久| 亚洲综合精品二区| 国产乱人偷精品视频| 我的老师免费观看完整版| 有码 亚洲区| 啦啦啦啦在线视频资源| 亚洲成色77777| 欧美精品国产亚洲| 亚洲精品成人av观看孕妇| 久久亚洲国产成人精品v| 又大又黄又爽视频免费| 亚洲国产最新在线播放| 日韩电影二区| 亚洲真实伦在线观看| 国产视频首页在线观看| 极品教师在线视频| 国产精品久久久久久久电影| av免费在线看不卡| 国内精品宾馆在线| 久热这里只有精品99| 亚洲精品自拍成人| 听说在线观看完整版免费高清| 女的被弄到高潮叫床怎么办| 波野结衣二区三区在线| 麻豆国产97在线/欧美| 简卡轻食公司| 99re6热这里在线精品视频| 交换朋友夫妻互换小说| 国产国拍精品亚洲av在线观看| 亚洲精品影视一区二区三区av| 午夜激情福利司机影院| 免费大片黄手机在线观看| 一级毛片 在线播放| 国产精品久久久久久精品电影小说 | 一级a做视频免费观看| 欧美3d第一页| 欧美区成人在线视频| 一级毛片电影观看| 国产欧美日韩一区二区三区在线 | 亚洲经典国产精华液单| 草草在线视频免费看| 成年人午夜在线观看视频| 欧美老熟妇乱子伦牲交| 亚洲av.av天堂| 美女视频免费永久观看网站| 三级男女做爰猛烈吃奶摸视频| 亚洲国产欧美人成| 成人亚洲精品av一区二区| 免费av不卡在线播放| 搡女人真爽免费视频火全软件| 亚洲av免费高清在线观看| 亚洲精品一二三| 少妇丰满av| 99热6这里只有精品| 欧美人与善性xxx| 亚洲av免费在线观看| 最近的中文字幕免费完整| 亚洲色图av天堂| 爱豆传媒免费全集在线观看| 久热久热在线精品观看| 两个人的视频大全免费| 免费大片18禁| 亚洲av在线观看美女高潮| 99精国产麻豆久久婷婷| 一级毛片电影观看| 如何舔出高潮| 插逼视频在线观看| 日韩强制内射视频| 久久精品国产亚洲av天美| 边亲边吃奶的免费视频| 黄色怎么调成土黄色| 老司机影院成人| 波野结衣二区三区在线| 午夜激情久久久久久久| 国产精品一区www在线观看| 亚洲精品一二三| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 亚洲av免费高清在线观看| 少妇高潮的动态图| 黑人高潮一二区| 亚洲国产成人一精品久久久| 人妻夜夜爽99麻豆av| 久久人人爽av亚洲精品天堂 | 热re99久久精品国产66热6| 99精国产麻豆久久婷婷| 国产亚洲最大av| 观看免费一级毛片| 草草在线视频免费看| 成年人午夜在线观看视频| 国产精品一区二区性色av| 色婷婷久久久亚洲欧美| 欧美xxⅹ黑人| 大香蕉97超碰在线| 国内精品美女久久久久久| 高清欧美精品videossex| 免费大片黄手机在线观看| 久久97久久精品| 国产精品一区二区三区四区免费观看| 69人妻影院| 亚洲欧美日韩东京热| 色视频在线一区二区三区| 亚洲性久久影院| 久久综合国产亚洲精品| av专区在线播放| 日韩av免费高清视频| 久热这里只有精品99| 久久ye,这里只有精品| 国产精品国产av在线观看| 亚洲精品中文字幕在线视频 | 国产av国产精品国产| 欧美潮喷喷水| 国产一区亚洲一区在线观看| 久久久a久久爽久久v久久| 亚洲av男天堂| 欧美成人午夜免费资源| 中文欧美无线码| 亚洲精品成人av观看孕妇| 国产午夜福利久久久久久| 王馨瑶露胸无遮挡在线观看| 网址你懂的国产日韩在线| 少妇裸体淫交视频免费看高清| 黄色欧美视频在线观看| 免费播放大片免费观看视频在线观看| 秋霞在线观看毛片| 我的老师免费观看完整版| 午夜福利视频1000在线观看| av在线亚洲专区| 男女边吃奶边做爰视频| 国产色爽女视频免费观看| 日韩不卡一区二区三区视频在线| 欧美高清性xxxxhd video| 精品久久久久久久久亚洲| 久久99热6这里只有精品| 蜜臀久久99精品久久宅男| 国产精品久久久久久精品电影小说 | 久久精品国产亚洲网站| 交换朋友夫妻互换小说| 春色校园在线视频观看| 麻豆久久精品国产亚洲av| 禁无遮挡网站| 欧美日韩一区二区视频在线观看视频在线 | 成人特级av手机在线观看| 男女无遮挡免费网站观看| 成人国产av品久久久| 国产成人精品一,二区| 久久99热这里只有精品18| 中文字幕亚洲精品专区| 日韩中字成人| 亚洲精品第二区| 永久免费av网站大全| av国产久精品久网站免费入址| 一个人看视频在线观看www免费| 久久久精品94久久精品| 免费少妇av软件| 亚洲av.av天堂| 亚洲精品一区蜜桃| 日韩不卡一区二区三区视频在线| 成人二区视频| 日日啪夜夜撸| 少妇人妻精品综合一区二区| av在线蜜桃| 国产亚洲精品久久久com| 成人亚洲精品一区在线观看 | 三级国产精品欧美在线观看| 18禁在线无遮挡免费观看视频| 日本熟妇午夜| 黄色怎么调成土黄色| 国产伦理片在线播放av一区| 男女边吃奶边做爰视频| 日韩视频在线欧美| 最近最新中文字幕免费大全7| 亚洲人成网站在线播| 国产一区二区三区av在线| 三级国产精品片| 精品人妻偷拍中文字幕| 国产片特级美女逼逼视频| 全区人妻精品视频| 国模一区二区三区四区视频| 一个人看的www免费观看视频| 亚洲aⅴ乱码一区二区在线播放| 亚洲最大成人中文| 深爱激情五月婷婷| 久久精品国产亚洲av涩爱| 日日摸夜夜添夜夜爱| 国产成人一区二区在线| 国产黄片视频在线免费观看| 亚洲欧美日韩东京热| 免费观看性生交大片5| 日本熟妇午夜| 亚洲色图综合在线观看| 亚洲av二区三区四区| 国产免费一区二区三区四区乱码| 免费av不卡在线播放| 日韩欧美精品免费久久| 国产av码专区亚洲av| 午夜福利视频1000在线观看| av在线蜜桃| 日韩亚洲欧美综合| 国内精品美女久久久久久| 高清欧美精品videossex| 国产熟女欧美一区二区| 国产亚洲精品久久久com| 亚洲精品aⅴ在线观看| 久久久久国产精品人妻一区二区| 久久久亚洲精品成人影院| 在线免费观看不下载黄p国产| 蜜臀久久99精品久久宅男| 一级av片app| 99久久精品热视频| 国产大屁股一区二区在线视频| 精品久久久久久电影网| 男女那种视频在线观看| 一级a做视频免费观看| 欧美高清成人免费视频www| 国产精品精品国产色婷婷| 亚洲在久久综合| 亚洲欧美日韩卡通动漫| 女人久久www免费人成看片| 午夜福利视频精品| 夜夜看夜夜爽夜夜摸| 在线 av 中文字幕| 久久久午夜欧美精品| 成人二区视频| 久久99热这里只有精品18| 亚洲在线观看片| 亚洲精品国产色婷婷电影| 特级一级黄色大片| 欧美成人a在线观看| 国产熟女欧美一区二区| 国产淫片久久久久久久久| 中文字幕av成人在线电影| 韩国av在线不卡|