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

    低速流場中柔性懸臂板的后顫振響應(yīng)

    2017-11-20 01:53:45陳濤徐敏謝丹安效民
    航空學(xué)報(bào) 2017年3期
    關(guān)鍵詞:氣動(dòng)彈性氣動(dòng)力懸臂

    陳濤, 徐敏, 謝丹, 安效民

    西北工業(yè)大學(xué) 航天學(xué)院, 西安 710072

    低速流場中柔性懸臂板的后顫振響應(yīng)

    陳濤, 徐敏*, 謝丹, 安效民

    西北工業(yè)大學(xué) 航天學(xué)院, 西安 710072

    建立了一個(gè)新的非線性氣動(dòng)彈性模型,對低速流場中柔性懸臂板的后顫振響應(yīng)特性進(jìn)行了分析。建模中考慮了結(jié)構(gòu)幾何非線性、氣動(dòng)力非線性以及兩者之間的強(qiáng)耦合效應(yīng)。通過實(shí)驗(yàn)數(shù)據(jù)對所建立的氣動(dòng)彈性模型進(jìn)行了驗(yàn)證。發(fā)現(xiàn)在低速流場中柔性懸臂板可能會(huì)以周期加倍的方式進(jìn)入混沌運(yùn)動(dòng)。結(jié)構(gòu)幾何非線性效應(yīng)和翼尖渦引起的非定常氣動(dòng)力效應(yīng)對柔性懸臂板的結(jié)構(gòu)響應(yīng)有顯著影響,而尾渦變形引起的非定常氣動(dòng)力對結(jié)構(gòu)運(yùn)動(dòng)的影響較小。還研究了不同耦合算法的差異,給出了小展弦比大柔性結(jié)構(gòu)非線性氣動(dòng)彈性數(shù)值仿真時(shí)耦合策略的選擇依據(jù)。

    低速流場; 柔性懸臂板; 后顫振; 混沌運(yùn)動(dòng); 強(qiáng)耦合

    低速流場中懸臂板的氣動(dòng)彈性行為可以用來解釋現(xiàn)實(shí)中的許多現(xiàn)象,如機(jī)翼顫振、旗幟飄動(dòng)和打鼾等。最近, Dunnmon[1]和Doaré[2]等甚至開始研究利用柔性懸壁板的顫振來收集能量用于發(fā)電。因此,低速流場中懸臂板的非線性氣動(dòng)彈性特性研究是一項(xiàng)很有意義的工作。

    關(guān)于低速流場中的非線性氣動(dòng)彈性問題,早期的研究可以追溯到二維情形,Lee等[3]對此做過詳細(xì)的綜述,二維問題中極限環(huán)振蕩(LCO)、混沌運(yùn)動(dòng)等現(xiàn)象都已通過數(shù)值計(jì)算或?qū)嶒?yàn)等手段發(fā)現(xiàn)。對于三維問題,非線性效應(yīng)比二維情況要復(fù)雜的多。Tang[4]和Attar[5]等通過耦合非定常渦格法和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)的方法建立了一個(gè)非線性氣動(dòng)彈性模型,并研究了切尖三角翼和小展弦比長方形機(jī)翼在低速流場中的極限環(huán)振蕩特性。研究發(fā)現(xiàn)這些非線性氣動(dòng)彈性現(xiàn)象主要受結(jié)構(gòu)幾何非線性效應(yīng)的主導(dǎo)。Cavallaro等[6]采用有限元和勢流理論方法研究了柔性連接翼的非線性氣動(dòng)彈性特性。Arena等[7]采用勢流理論方法研究了大柔性機(jī)翼的后顫振響應(yīng)特性,他們發(fā)現(xiàn)流場的非定常效應(yīng)會(huì)對機(jī)翼運(yùn)動(dòng)產(chǎn)生重要影響。此外,Relvas和Suleman[8]基于結(jié)構(gòu)有限元方法和非定常渦格法也建立了一個(gè)非線性氣動(dòng)彈性模型并研究了柔性懸臂板的極限環(huán)振蕩現(xiàn)象。國內(nèi)尹維龍和田東奎[9]基于勢流理論方法研究了二維柔性翼型的氣動(dòng)彈性特性。謝長川等[10]基于勢流理論和結(jié)構(gòu)有限元方法研究了柔性機(jī)翼的低速顫振特性,結(jié)果發(fā)現(xiàn)結(jié)構(gòu)幾何非線性效應(yīng)對柔性翼的顫振邊界有重要影響。張健和向錦武[11]基于ONERA失速氣動(dòng)力模型和非線性梁模型研究了柔性飛機(jī)的非線性氣動(dòng)彈性特性。

    以上研究表明,低速流場中三維柔性結(jié)構(gòu)的氣動(dòng)彈性建模需要考慮3種非線性效應(yīng):① 由于結(jié)構(gòu)大變形引起的結(jié)構(gòu)幾何非線性效應(yīng);② 由于流場渦等引起的非定常氣動(dòng)力效應(yīng);③ 這兩者的耦合效應(yīng)。首先,對于結(jié)構(gòu)幾何非線性問題,Gordnier和Visbal[12]研究表明,氣動(dòng)彈性建模中遇到的主要問題是像von Karmann板理論等方法在求解幾何大變形問題時(shí)會(huì)產(chǎn)生“鎖住”效應(yīng)。其次,對于低速氣動(dòng)彈性問題,目前國際上廣泛采用勢流理論方法進(jìn)行非定常氣動(dòng)力建模。雖然有更精確的方法如CFD等,但在低速問題中柔性結(jié)構(gòu)響應(yīng)周期往往較長,采用CFD方法會(huì)消耗大量計(jì)算時(shí)間。Murua等[13]的研究表明非定常渦格法在計(jì)算準(zhǔn)確性和計(jì)算效率2個(gè)方面都能達(dá)到較滿意的效果,國內(nèi)許多學(xué)者也采用非定常渦格法研究低速飛行器的非定常氣動(dòng)力建模問題[14-15],該方法還常用于氣動(dòng)力非定常效應(yīng)很強(qiáng)的撲翼氣動(dòng)彈性現(xiàn)象研究[16-17]。最后,關(guān)于結(jié)構(gòu)幾何非線性與非定常氣動(dòng)力的耦合效應(yīng),Piperno和Farhat[18]的研究表明,氣動(dòng)彈性建模中采用松耦合策略會(huì)產(chǎn)生一定的數(shù)值誤差。對于柔性翼這樣非線性特性很強(qiáng)的氣動(dòng)彈性問題,松耦合策略可能會(huì)引起較大的誤差。

    如前所述,Attar和Relvas等只研究了柔性壁板的周期性極限環(huán)振蕩現(xiàn)象,雖然他們都考慮了結(jié)構(gòu)和氣動(dòng)力的非線性效應(yīng),但對于耦合策略的影響仍然有待研究。本文的工作是對Attar[5]和Relvas[8]等研究工作的進(jìn)一步擴(kuò)展。在本文的氣動(dòng)彈性建模中不僅詳細(xì)考慮了結(jié)構(gòu)幾何非線性效應(yīng)帶來的“鎖住”問題和非定常氣動(dòng)力問題,還考慮了結(jié)構(gòu)與氣動(dòng)力之間的強(qiáng)耦合問題。其中,結(jié)構(gòu)建模采用一種多變量有限元方法,非定常氣動(dòng)力建模采用非定常渦格法,采用強(qiáng)耦合策略實(shí)現(xiàn)結(jié)構(gòu)和氣動(dòng)力求解器的耦合。本文的結(jié)構(gòu)安排如下,首先給出新的非線性氣動(dòng)彈性建模過程,然后通過柔性三角翼極限環(huán)振蕩的仿真對所提出的氣動(dòng)彈性模型進(jìn)行驗(yàn)證,最后對柔性懸臂板后顫振響應(yīng)以及相關(guān)的非線性影響因素給出分析結(jié)果。

    1 結(jié)構(gòu)幾何非線性建模

    柔性懸臂板的結(jié)構(gòu)建模采用多變量實(shí)體殼單元。這種單元的有限元列式過程簡要介紹如下。

    實(shí)體殼單元是一種多變量的有限元單元,這種單元的位移場、應(yīng)變場和應(yīng)力場是相互獨(dú)立的,可以分別假設(shè)獨(dú)立的插值模式。其有限元列式基于三場FHW變分原理[19]:

    (1)

    對于上述處理方法,筆者早期研究工作已經(jīng)證明這是一種解決“鎖住”效應(yīng)的有效手段[20]。當(dāng)給定了位移、應(yīng)變和應(yīng)力的插值形函數(shù)以后,就可以對式(1)進(jìn)行離散。本文的算法中采用一種無條件穩(wěn)定的能量-動(dòng)量守恒算法作為時(shí)間推進(jìn)算法,采用Newton-Raphson算法對式(1)進(jìn)行線性化:

    (2)

    (3)

    2 非定常氣動(dòng)力建模

    與Attar[5]和Relvas[8]等相同,本文采用非定常渦格法進(jìn)行非定常氣動(dòng)力建模。該方法假設(shè)結(jié)構(gòu)周圍的流動(dòng)是無旋附著流,并且不考慮黏性效應(yīng)。對于本文的懸臂板問題,非定常渦格法的非線性特性主要體現(xiàn)在以下3個(gè)方面:① 結(jié)構(gòu)表面的渦格跟隨結(jié)構(gòu)變形而非線性地改變形狀和位置。② 機(jī)翼尾部拖出的尾渦可以自由變形。③ 隨著結(jié)構(gòu)的變形,可能會(huì)有翼尖渦不斷拖出并卷起。在本文的研究中,翼尖渦的拖出方式與尾渦的拖出方式相同。這種翼尖渦模型在三角翼的相關(guān)研究中已被證明是一種有效的翼尖渦模擬方法[21]。在后面的數(shù)值仿真中會(huì)發(fā)現(xiàn),翼尖渦效應(yīng)其實(shí)對結(jié)構(gòu)的氣動(dòng)彈性響應(yīng)有很大影響。但是在Relvas等[8]的研究中并沒有考慮翼尖渦效應(yīng),因此,本文是對其研究工作的拓展。下面簡要介紹非定常渦格法的求解流程。

    定常渦格法的求解過程主要是為了得到結(jié)構(gòu)表面的渦量強(qiáng)度分布,這些渦量強(qiáng)度用來計(jì)算結(jié)構(gòu)表面的氣動(dòng)力載荷。首先,將結(jié)構(gòu)表面用一系列渦格離散,在每一個(gè)渦格中,渦的強(qiáng)度為常量。每一個(gè)渦格片段的渦強(qiáng)度由共享該片段的相鄰2個(gè)渦格所決定。通過Biot-Savart定律可以求解每一個(gè)渦格片段的誘導(dǎo)速度場。隨著時(shí)間的推移,結(jié)構(gòu)表面的渦會(huì)從結(jié)構(gòu)后緣拖出形成尾渦。剛拖出的尾渦強(qiáng)度與相鄰的結(jié)構(gòu)表面渦強(qiáng)度相同,當(dāng)結(jié)構(gòu)表面的渦拖出形成尾渦后,渦的強(qiáng)度保持不變。非定常渦格法通過滿足物面無穿透邊界條件計(jì)算得到物面的渦強(qiáng)分布:

    [V∞+Vbound(r,t)+Vwake(r,t)]·n=0

    AΓbound=Rhs

    (4)

    式中:Vbound和Vwake分別為結(jié)構(gòu)表面和尾渦中的誘導(dǎo)速度;r為位置向量;t為時(shí)間;n為結(jié)構(gòu)表面的單位法向量;A為氣動(dòng)力影響系數(shù)矩陣;Γbound為待求解的未知渦量向量;Rhs為殘差向量。當(dāng)物面渦量求得以后,結(jié)構(gòu)表面的空氣壓強(qiáng)可由Bernoulli方程求解得到:

    (5)

    式中:Φ為速度勢函數(shù);p為壓強(qiáng);V∞和p∞分別為遠(yuǎn)場處的速度和壓強(qiáng)。

    結(jié)構(gòu)上下表面的壓強(qiáng)差表達(dá)式為

    Δp=pu-pl=ρ∞[Vu·Vu-Vl·Vl+

    (6)

    式中:Δp為壓強(qiáng)差;pu和pl分別為結(jié)構(gòu)上下表面的壓強(qiáng);Vu和Vl分別為結(jié)構(gòu)上下表面的速度;ru和rl分別為結(jié)構(gòu)上下表面的位置向量。

    由壓強(qiáng)差便可求出對應(yīng)的氣動(dòng)力為

    Fi=ΔpiAini

    (7)

    3 非線性氣動(dòng)彈性建模

    基于上述非線性有限元求解方法和非定常氣動(dòng)力求解方法可以建立柔性壁板的氣動(dòng)彈性模型。下面簡要介紹建模中所采用的耦合策略和結(jié)構(gòu)/流場求解器之間的數(shù)據(jù)交換方法。

    3.1 耦合策略

    眾所周知,傳統(tǒng)的松耦合策略會(huì)產(chǎn)生時(shí)間滯后效應(yīng)。這會(huì)導(dǎo)致在每一個(gè)時(shí)間步計(jì)算得到的結(jié)構(gòu)/流體交界面的位置不準(zhǔn)確。然而,在引言中所提到的所有關(guān)于柔性懸臂板的研究都是采用松耦合算法。本文的氣動(dòng)彈性建模采用耦合精度更高的強(qiáng)耦合策略(見圖1),在每一個(gè)時(shí)間步內(nèi)引入預(yù)估-校正方法來提高結(jié)構(gòu)/流體交界面位置的計(jì)算精度。圖1中UVLM表示非定常渦格法氣動(dòng)力求解器,F(xiàn)EM表示有限元求解器。

    內(nèi)迭代的收斂標(biāo)準(zhǔn)給定為相鄰兩次求得的結(jié)構(gòu)位移向量的相對誤差小于1%。此外,在實(shí)際的數(shù)值算法實(shí)現(xiàn)中還應(yīng)給定最大內(nèi)迭代步數(shù),以防止出現(xiàn)永遠(yuǎn)不收斂的情況,不過在本文的研究中沒有出現(xiàn)這種不收斂的問題。關(guān)于上述的強(qiáng)耦合算法需要指出的是當(dāng)最大內(nèi)迭代步數(shù)給定為1時(shí),強(qiáng)耦合算法將自動(dòng)退化為松耦合算法。本文強(qiáng)耦合算法的最大內(nèi)迭代步數(shù)給定為10,實(shí)際中發(fā)現(xiàn)所有內(nèi)迭代計(jì)算均會(huì)在10步以內(nèi)達(dá)到收斂。

    3.2 數(shù)據(jù)交換方法

    在非線性氣動(dòng)彈性建模中,結(jié)構(gòu)和流體求解器之間的數(shù)據(jù)交換方法不僅影響計(jì)算速度,更影響計(jì)算精度,特別是當(dāng)結(jié)構(gòu)變形很大時(shí)(例如極限環(huán)振蕩等),這一問題更加突出。本文采用徑向基函數(shù)方法實(shí)現(xiàn)氣動(dòng)網(wǎng)格的變形[20],這一方法非常適合于大變形過程的動(dòng)網(wǎng)格插值。氣動(dòng)力的插值與結(jié)構(gòu)位移場插值不同,它需要保證插值前后力場的等效,即總的力和力矩需要保持守恒。本文采用有限元形狀函數(shù)插值方法實(shí)現(xiàn)氣動(dòng)力的插值。

    4 數(shù)值計(jì)算

    本節(jié)首先對上述所建立的非線性氣動(dòng)彈性模型進(jìn)行驗(yàn)證,在此基礎(chǔ)上對柔性懸臂板的后顫振特性進(jìn)行數(shù)值仿真并給出分析結(jié)果。

    4.1 非線性氣動(dòng)彈性模型的驗(yàn)證

    通過一個(gè)三角翼極限環(huán)振蕩(LCO)的數(shù)值仿真與實(shí)驗(yàn)結(jié)果的對比對所提出的非線性氣動(dòng)彈性模型進(jìn)行驗(yàn)證。該模型是一個(gè)45° 的直角三角形平板,厚度為1.6 mm,翼根中部60%的部位受到固支約束。彈性模量為3.3 GPa,密度為1 299 kg/m3[6]。圖2為該三角翼模型的示意圖,V∞為來流速度。

    經(jīng)過網(wǎng)格收斂性分析后,氣動(dòng)彈性仿真的結(jié)構(gòu)有限元網(wǎng)格共有819個(gè)實(shí)體殼單元。氣動(dòng)網(wǎng)格的大小為60×20(弦向×展向)。尾渦的長度取為弦長的4倍。

    圖3給出了在1° 迎角時(shí)機(jī)翼翼尖面外方向極限環(huán)振蕩的最大速度對比結(jié)果。計(jì)算所得該三角翼的顫振速度為23.0 m/s,實(shí)驗(yàn)值為24.0 m/s。

    圖3表明在來流速度小于26.0 m/s時(shí),本文的計(jì)算結(jié)果比Attar等[5]的計(jì)算結(jié)果要大,兩者之間的差別隨著流速的增大而逐漸減小。在較高的來流速度時(shí),計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果非常接近,并且與Attar等[5]采用高精度有限元方法的計(jì)算結(jié)果也非常接近。在較低的速度范圍時(shí),計(jì)算的顫振點(diǎn)和極限環(huán)振蕩結(jié)果與Cavallaro等[6]的結(jié)果比較接近。需要指出的是在較低的來流速度范圍時(shí),本文計(jì)算結(jié)果與其他兩個(gè)計(jì)算結(jié)果的差別主要是由于對結(jié)構(gòu)阻尼的不同處理導(dǎo)致的。Attar[5]和Cavallaro[6]等分別采用了不同的結(jié)構(gòu)阻尼形式,由于非線性氣動(dòng)彈性計(jì)算中結(jié)構(gòu)阻尼的影響是一個(gè)比較復(fù)雜的問題,目前國際上仍沒有定論,因此本文沒有考慮結(jié)構(gòu)阻尼。從計(jì)算結(jié)果的對比來看,本文所提出的非線性氣動(dòng)彈性模型能較好地模擬大幅值的極限環(huán)運(yùn)動(dòng)。

    需要指出的圖3中本文的計(jì)算結(jié)果采用的是強(qiáng)耦合策略,這主要是為了驗(yàn)證前文所建立的強(qiáng)耦合策略的正確性。強(qiáng)耦合策略由于帶有內(nèi)迭代過程而比較耗時(shí),對于這里的三角翼模型,傳統(tǒng)松耦合策略與強(qiáng)耦合策略的計(jì)算結(jié)果也比較接近,見圖4。

    由于本文的氣動(dòng)彈性建模借鑒并改進(jìn)了Relvas等[8]的方法,因此有必要將本文方法與其進(jìn)行對比。文獻(xiàn)[8]中研究的方形懸臂板模型展長(弦長)為0.3 m,板厚度為1 mm,材料為鋁合金,計(jì)算初始迎角為0°。圖5給出了采用Nastran線性分析方法得到的顫振V-g圖,圖6給出了本文的計(jì)算結(jié)果。由圖5可知,Nastran計(jì)算的顫振速度為53.7 m/s,由圖6可知,本文的計(jì)算結(jié)果為57.5 m/s。文獻(xiàn)[8]的計(jì)算結(jié)果為42.0 m/s,對比可知,本文結(jié)果與Nastran的相對誤差為7.08%,文獻(xiàn)[8]的結(jié)果與Nastran的相對誤差為21.79%。計(jì)算結(jié)果表明,本文的方法用于方形懸臂板后顫振響應(yīng)分析是可行的。

    4.2 柔性懸臂板的非線性氣動(dòng)彈性特性

    經(jīng)過實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證以后,將所建立的氣動(dòng)彈性模型用于柔性懸臂板(圖7)的后顫振響應(yīng)分析。板的厚度為0.5 mm,材料屬性如下:彈性模量為71 GPa,密度為2 800 kg/m3,泊松比為 0.33。需要說明的是文獻(xiàn)[8]中的方形懸臂板模型偏硬,不適合于較大幅值的極限環(huán)振蕩研究,而本文第2個(gè)算例主要是為了研究較大幅度的極限環(huán)振蕩,因此本文研究的懸臂板模型尺寸在文獻(xiàn)[8]的基礎(chǔ)上進(jìn)行了調(diào)整。在數(shù)值仿真之前,通過網(wǎng)格收斂性分析確定采用的計(jì)算網(wǎng)格。最終使用的結(jié)構(gòu)網(wǎng)格數(shù)量為30×30(弦向×展向),對應(yīng)的氣動(dòng)力計(jì)算網(wǎng)格數(shù)量為40×40。計(jì)算過程中尾渦的長度為弦長的4倍。

    對柔性懸臂板進(jìn)行氣動(dòng)彈性數(shù)值仿真,初始迎角為0°。表1給出了本文計(jì)算的結(jié)構(gòu)固有頻率和解析解[22](雙向梁函數(shù)組合級數(shù)逼近法)以及Nastran計(jì)算結(jié)果的對比。圖8給出了Nastran的顫振分析結(jié)果。圖9給出了來流速度分別為38.3、46、54、58、62、66 m/s時(shí)前緣翼尖的結(jié)構(gòu)位移響應(yīng)。

    由表1可知,本文采用實(shí)體殼單元計(jì)算的結(jié)構(gòu)固有頻率與解析解和Nastran的計(jì)算結(jié)果非常接近,與解析解的最大誤差僅為1.49%(二階)。由圖8可知,Nastran分析的結(jié)構(gòu)顫振邊界為35.7 m/s,顫振發(fā)生的機(jī)理為二階扭轉(zhuǎn)模態(tài)發(fā)散。圖9(a)表明本文基于有限元和渦格法計(jì)算的顫振邊界為38.3 m/s,兩者之間的相對誤差為6.79%。

    Method1st2nd3rd4th5thAnalysisresult10.53326.04464.56782.32194.088Nastran10.59925.61364.61082.70293.531Currentwork10.60525.65564.84583.11593.946

    結(jié)構(gòu)響應(yīng)曲線圖9表明,當(dāng)來流速度高于顫振速度時(shí),結(jié)構(gòu)會(huì)發(fā)生極限環(huán)振蕩。當(dāng)來流速度高于54 m/s時(shí),結(jié)構(gòu)響應(yīng)變得非常復(fù)雜。為了更好地分析結(jié)構(gòu)響應(yīng)隨來流的變化趨勢,圖10給出了相應(yīng)的相平面圖。

    相平面圖表明,當(dāng)來流速度為46 m/s時(shí),結(jié)構(gòu)響應(yīng)形式為周期性的極限環(huán)振蕩。當(dāng)來流速度為54 m/s時(shí),結(jié)構(gòu)響應(yīng)形式為準(zhǔn)周期的極限環(huán)振蕩。當(dāng)來流速度為58 m/s時(shí),結(jié)構(gòu)響應(yīng)出現(xiàn)周期加倍現(xiàn)象。當(dāng)來流速度增大到62 m/s和66 m/s時(shí),結(jié)構(gòu)響應(yīng)的規(guī)律變得非常復(fù)雜。筆者早期關(guān)于超聲速流中壁板顫振的研究表明這種運(yùn)動(dòng)形式的變化是一種進(jìn)入混沌運(yùn)動(dòng)的路徑[23]。因此,本文首次發(fā)現(xiàn)了低速流場中三維柔性懸臂板很有可能也會(huì)以周期加倍的方式進(jìn)入混沌運(yùn)動(dòng)。不過由于本文計(jì)算得到的時(shí)域響應(yīng)結(jié)果步數(shù)有限,尚不能斷定一定發(fā)生了混沌運(yùn)動(dòng),定量的證明將在后續(xù)工作中采用降階的方法進(jìn)行研究。

    從位移響應(yīng)曲線(圖9)還可以看出,在較大來流速度下,結(jié)構(gòu)首先經(jīng)歷一個(gè)較大幅值的瞬態(tài)響應(yīng)過程,然后收斂到穩(wěn)態(tài)響應(yīng)。需要指出的是,本文在氣動(dòng)彈性建模中沒有考慮結(jié)構(gòu)阻尼。另外,根據(jù)國際權(quán)威學(xué)者Hallissy和Cesnik[24]的研究,勢流理論與高精度的氣動(dòng)力計(jì)算方法(如CFD方法等)相比氣動(dòng)阻尼偏小,這可能導(dǎo)致結(jié)構(gòu)的瞬態(tài)響應(yīng)幅值偏大。在Cavallaro等[6]關(guān)于柔性連接翼非線性氣動(dòng)彈性特性的研究中也發(fā)現(xiàn)了類似的瞬態(tài)響應(yīng)幅值很大的情況,可見這是采用渦格法計(jì)算時(shí)的一種普遍現(xiàn)象。由于本文主要研究的是結(jié)構(gòu)穩(wěn)態(tài)響應(yīng)特性,因此有必要考察結(jié)構(gòu)瞬態(tài)響應(yīng)對穩(wěn)態(tài)響應(yīng)的影響。

    圖11給出了分別采用結(jié)構(gòu)一階模態(tài)和二階模態(tài)作為有限元計(jì)算的初始速度擾動(dòng)時(shí)結(jié)構(gòu)位移響應(yīng)的對比。由圖可知,結(jié)構(gòu)在不同初始擾動(dòng)激勵(lì)下經(jīng)歷了不同的瞬態(tài)響應(yīng)過程,但最終的穩(wěn)態(tài)響應(yīng)非常接近。因此說明結(jié)構(gòu)瞬態(tài)響應(yīng)不會(huì)對穩(wěn)態(tài)響應(yīng)造成很大影響。

    以上研究表明低速流場中柔性懸臂板的后顫振運(yùn)動(dòng)存在很強(qiáng)的非線性特性,因此有必要對其中的非線性因素進(jìn)行深入研究。下面將從結(jié)構(gòu)幾何非線性、氣動(dòng)力非線性以及耦合方式的影響3個(gè)方面進(jìn)行分析。

    圖12給出了懸臂板在來流速度為46 m/s時(shí)考慮和不考慮結(jié)構(gòu)幾何非線性效應(yīng)的計(jì)算結(jié)果對比,由圖可知在懸臂板的后顫振響應(yīng)中結(jié)構(gòu)幾何非線性效應(yīng)的影響非常大。

    為了進(jìn)一步分析結(jié)構(gòu)幾何非線性效應(yīng)對后顫振響應(yīng)的影響,表2給出了在穩(wěn)態(tài)響應(yīng)階段結(jié)構(gòu)最大位移處計(jì)算的模態(tài)頻率??梢钥闯?,隨著來流速度的增加,結(jié)構(gòu)在最大位移處的一階模態(tài)頻率變化較小,二階及高階模態(tài)頻率變化較大。在來流速度為54 m/s時(shí),結(jié)構(gòu)極限環(huán)振蕩最大位移達(dá)到了板厚度的4倍左右,此時(shí)結(jié)構(gòu)二階模態(tài)頻率增大了約10%。圖12中采用線性結(jié)構(gòu)模型的時(shí)域響應(yīng)在t=0.3 s時(shí)結(jié)構(gòu)最大位移已達(dá)到板厚度的5倍左右,此時(shí)結(jié)構(gòu)幾何非線性效應(yīng)已經(jīng)非常明顯。

    幾何非線性效應(yīng)主要是由結(jié)構(gòu)內(nèi)部應(yīng)力對結(jié)構(gòu)剛度的影響造成的,圖13給出了表2中各來流速度對應(yīng)的結(jié)構(gòu)Von Mises應(yīng)力分布。由圖可知,當(dāng)來流速度從38.3 m/s增大到46 m/s時(shí),結(jié)構(gòu)內(nèi)部應(yīng)力最大值由0.06 MPa增大到了2.4 MPa,此時(shí)結(jié)構(gòu)剛度產(chǎn)生了較大變化,幾何非線性效應(yīng)已非常明顯。

    除了結(jié)構(gòu)幾何非線性,氣動(dòng)力非線性可能也會(huì)對柔性懸臂板的后顫振響應(yīng)產(chǎn)生重要影響。下面主要考察尾渦和翼尖渦對非線性氣動(dòng)力的貢獻(xiàn)。圖14給出了考慮翼尖渦和不考慮翼尖渦的計(jì)算結(jié)果對比。圖15給出了在來流速度為66 m/s不考慮翼尖渦時(shí)結(jié)構(gòu)穩(wěn)態(tài)響應(yīng)的變形情況。圖16 給出了分別采用自由變形尾渦模型和平面尾渦模型的計(jì)算結(jié)果。需要指出的是,此時(shí)為了單獨(dú)考察尾渦變形的影響,數(shù)值計(jì)算時(shí)沒有考慮翼尖渦模型。

    當(dāng)不考慮翼尖渦時(shí)(圖14虛線),結(jié)構(gòu)的穩(wěn)態(tài)響應(yīng)形式為動(dòng)態(tài)屈曲。雖然初始迎角為0°,但在瞬態(tài)響應(yīng)過程中,結(jié)構(gòu)產(chǎn)生了扭轉(zhuǎn)變形(圖15),翼面上當(dāng)?shù)赜遣粸?°,所產(chǎn)生的升力使得結(jié)構(gòu)最終在某一平衡位置附近做幅值非常小的動(dòng)態(tài)屈曲運(yùn)動(dòng),這種現(xiàn)象在筆者早期關(guān)于超聲速氣流中懸臂板的混沌運(yùn)動(dòng)研究中發(fā)現(xiàn)過(見文獻(xiàn)[23])。但考慮翼尖渦后,結(jié)構(gòu)響應(yīng)形式完全不同,圖14表明是否考慮翼尖渦會(huì)對計(jì)算結(jié)果產(chǎn)生重要影響。圖16表明是否考慮尾渦的變形對結(jié)構(gòu)響應(yīng)幾乎沒有影響。由以上分析可知,翼尖渦對非線性氣動(dòng)力的貢獻(xiàn)較大,而尾渦的自由變形引起的氣動(dòng)力非線性效應(yīng)較小。

    表2 變形最大位移處結(jié)構(gòu)模態(tài)頻率

    Table 2 Structural modal frequencies at position of maximum deformation Hz

    Inflowvelocity/(m·s-1)1st2nd3rd4th5th010.60525.65564.84583.11593.94638.310.60625.66364.92983.19294.0564610.62626.98265.33083.54095.0025410.65428.25965.77484.80295.8795811.13441.06678.044108.315119.752

    圖17給出了在來流速度為54 m/s和66 m/s時(shí)分別采用松耦合、改進(jìn)的松耦合以及強(qiáng)耦合策略的計(jì)算結(jié)果對比。這里的改進(jìn)松耦合策略為Relvas和Suleman[8]在研究相似柔性懸臂板非線性氣動(dòng)彈性問題時(shí)所采用的耦合方式。由圖可知,在較低的來流速度54 m/s時(shí),改進(jìn)的松耦合方式與強(qiáng)耦合方式的計(jì)算結(jié)果非常接近,而松耦合算法與強(qiáng)耦合算法的差異隨著時(shí)間的推移會(huì)逐漸增大。在較高的來流速度66 m/s時(shí),不管是松耦合還是改進(jìn)的松耦合方式與強(qiáng)耦合方式的差別都是非常大的。

    仔細(xì)分析還發(fā)現(xiàn),在數(shù)值仿真過程中前述柔性三角翼與柔性懸臂板在最大來流速度時(shí)結(jié)構(gòu)的位移響應(yīng)均達(dá)到了板厚度的20倍左右(圖4和圖9(f)),因此,2個(gè)模型的幾何非線性效應(yīng)都是很強(qiáng)的。但在柔性三角翼的研究中,當(dāng)來流速度較高時(shí),不同的耦合策略計(jì)算得到的位移響應(yīng)差別不大。而對于柔性懸臂板模型,不同的耦合策略計(jì)算得到的位移響應(yīng)相差卻很大。2個(gè)模型的主要區(qū)別在于氣動(dòng)力非線性效應(yīng)的強(qiáng)弱,對于柔性三角翼,由于不存在翼尖渦引起的非線性氣動(dòng)力,因此,其氣動(dòng)力非線性效應(yīng)要弱于柔性懸臂板。據(jù)此可知,在小展弦比大柔性結(jié)構(gòu)的非線性氣動(dòng)彈性數(shù)值仿真時(shí),當(dāng)只有結(jié)構(gòu)幾何非線性效應(yīng)占主導(dǎo)地位時(shí),不同耦合策略的計(jì)算結(jié)果差別可能不是很大。當(dāng)結(jié)構(gòu)幾何非線性與氣動(dòng)力非線性效應(yīng)均較強(qiáng)時(shí)松耦合策略和改進(jìn)的松耦合策略都會(huì)引起很大的誤差。

    5 結(jié) 論

    建立一個(gè)新的非線性氣動(dòng)彈性模型對低速流場中柔性懸臂板的后顫振特性進(jìn)行了數(shù)值仿真。結(jié)構(gòu)幾何非線性問題采用多變量有限元方法求解,非定常氣動(dòng)力采用非定常渦格法求解,實(shí)現(xiàn)了結(jié)構(gòu)與氣動(dòng)力求解器之間的強(qiáng)耦合算法。采用柔性三角翼模型對所提出的氣動(dòng)彈性模型進(jìn)行了驗(yàn)證,結(jié)果表明,對于大柔性結(jié)構(gòu)的非線性氣動(dòng)彈性問題,所提出的算法可以給出較好的仿真結(jié)果。

    首次發(fā)現(xiàn)了在低速流場中三維柔性懸臂板可能會(huì)以位移響應(yīng)周期加倍的形式進(jìn)入混沌運(yùn)動(dòng)。數(shù)值仿真結(jié)果表明,在柔性懸臂板的后顫振響應(yīng)中,結(jié)構(gòu)幾何非線性效應(yīng)和氣動(dòng)力非線性效應(yīng)對結(jié)構(gòu)運(yùn)動(dòng)都有很大影響。其中,氣動(dòng)力非線性效應(yīng)主要由翼尖渦引起,而尾渦的變形對非線性氣動(dòng)力的貢獻(xiàn)相對較小。在小展弦比大柔性結(jié)構(gòu)的非線性氣動(dòng)彈性數(shù)值仿真時(shí),當(dāng)只有結(jié)構(gòu)幾何非線性效應(yīng)占主導(dǎo)地位時(shí),松耦合策略與強(qiáng)耦合策略的計(jì)算結(jié)果可能差別不是很大,此時(shí)可以采用松耦合策略或改進(jìn)的松耦合策略。當(dāng)結(jié)構(gòu)幾何非線性與氣動(dòng)力非線性效應(yīng)均較強(qiáng)時(shí)必須采用強(qiáng)耦合策略。

    需要說明的是本文目前只是初步研究了結(jié)構(gòu)幾何非線性、氣動(dòng)力非線性和耦合策略對柔性結(jié)構(gòu)非線性氣動(dòng)彈性仿真的影響。下一步的工作將對其他非線性因素如結(jié)構(gòu)阻尼效應(yīng)等進(jìn)行研究。另外,還將開展相關(guān)的風(fēng)洞試驗(yàn)研究。

    [1] DUNNMON J A, STANTON S C, MANN B P, et al. Power extraction from aeroelastic limit cycle oscillations[J]. Journal of Fluids and Structures, 2011, 27(8): 1182-1198.

    [2] DOARé O, MICHELIN S. Piezoelectric coupling in energy-harvesting fluttering flexible plates: Linear stability analysis and conversion efficiency[J]. Journal of Fluids and Structures, 2011, 27(8): 1357-1375.

    [3] LEE B H K, PRICE S J, WONG Y S. Nonlinear aeroelastic analysis of airfoils: Bifurcation and chaos[J]. Progress in Aerospace Sciences, 1999, 35(3): 205-334.

    [4] TANG D, DOWELL E H, HALL K C. Limit cycle oscillations of a cantilevered wing in low subsonic flow[J]. AIAA Journal, 1999, 37(3): 364-371.

    [5] ATTAR P J, DOWELL E H, WHITE J R. Modeling the LCO of a delta wing using a high fidelity structural model:AIAA-2004-1692[R]. Reston: AIAA, 2004.

    [6] CAVALLARO R, IANNELLI A, DEMASI L, et al. Phenomenology of nonlinear aeroelastic responses of highly deformable joined-wings configurations: AIAA-2014-1199[R]. Reston: AIAA, 2014.

    [7] ARENA A, LACARBONARA W, MARZOCCA P. Nonlinear aeroelastic formulation and postflutter analysis of flexible high aspect ratio wings[J]. Journal of Aircraft, 2013, 50(6): 1748-1764.

    [8] RELVAS A, SULEMAN A. Fluid-structure interaction modelling of nonlinear aeroelastic structures using the finite element corotational theory[J]. Journal of Fluids and Structures, 2006, 22(1): 59-75.

    [9] 尹維龍, 田東奎. 柔性翼型的氣動(dòng)彈性建模與顫振特性分析[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 2012, 44(9): 69-72. YIN W L, TIAN D K. Aeroelastic modeling and flutter characteristics of flexible aerofoil[J]. Journal of Harbin Institute of Technology, 2012, 44(9): 69-72 (in Chinese).

    [10] 謝長川, 吳志剛, 楊超. 大展弦比柔性機(jī)翼的氣動(dòng)彈性分析[J]. 北京航空航天大學(xué)學(xué)報(bào), 2003, 29(12): 1087-1090. XIE C C, WU Z G, YANG C. Aeroelastic analysis of flexible large aspect ratio wing[J]. Journal of Beijing University of Aeronautics and Astronautics, 2003, 29(12): 1087-1090 (in Chinese).

    [11] 張健, 向錦武. 柔性飛機(jī)非線性氣動(dòng)彈性與飛行動(dòng)力學(xué)耦合靜、動(dòng)態(tài)特性[J]. 航空學(xué)報(bào), 2011, 32(9): 1569-1582. ZHANG J, XIANG J W. Static and dynamic characteristics of coupled nonlinear aeroelasticity and flight dynamics of flexible aircraft[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(9): 1569-1582 (in Chinese).

    [12] GORDNIER R E, VISBAL M R. Development of a three-dimensional viscous aeroelastic solver for nonlinear panel flutter[J]. Journal of Fluids and Structures, 2002, 16(4): 497-527.

    [13] MURUA J, PALACIOS R, GRAHAM J M R. Applications of the unsteady vortex-lattice method in aircraft aeroelasticity and flight dynamics[J]. Progress in Aerospace Sciences, 2012, 55(5): 46-72.

    [14] 劉燚, 謝長川, 王立波, 等. 柔性飛機(jī)大變形曲面氣動(dòng)力計(jì)算及配平分析[J]. 工程力學(xué), 2015, 32(10): 239-249. LIU Y, XIE C C, WANG L B, et al. Nonplanar aerodynamic computation and trim analysis under large deflection of flexible aircraft[J]. Engineering Mechanics, 2015, 32(10):239-249 (in Chinese).

    [15] 宋磊, 楊華, 解靜峰, 等. 基于改進(jìn)渦格法的飛翼布局飛機(jī)穩(wěn)定性導(dǎo)數(shù)計(jì)算[J]. 南京航空航天大學(xué)學(xué)報(bào), 2014, 46(3): 457-462. SONG L, YANG H, XIE J F, et al. Predicting stability derivatives of flying wing aircraft based on improved vortex lattice method[J]. Journal of Nanjing University of Aeronautics and Astronautics, 2014, 46(3): 457-462 (in Chinese).

    [16] ROCCIA B A, PREIDIKMAN S, MASSA J C, et al. Modified unsteady vortex-lattice method to study flapping wings in hover flight[J]. AIAA Journal, 2013, 51(11): 2628-2642.

    [17] 賀紅林, 周翔. 柔性撲翼非定常渦格法氣動(dòng)力計(jì)算的改進(jìn)與實(shí)現(xiàn)[J]. 航空學(xué)報(bào), 2010, 31(6): 1121-1126. HE H L, ZHOU X. Implementation of an improved unsteady vortex lattice method for flexible flapping-wing aerodynamic computation[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(6): 1121-1126 (in Chinese).

    [18] PIPERNO S, FARHAT C. Partitioned procedures for the transient solution of coupled aeroelastic problems-Part II: Energy transfer analysis and three-dimensional applications[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(24-25): 3147-3170.

    [19] VU-QUOC L, TAN X G. Optimal solid shells for non-linear analysis of multilayer composites. 1. Static[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(9-10): 975-1016.

    [20] CHEN T, XU M, XIE L. Aeroelastic modeling using geometrically nonlinear solid-shell element[J]. AIAA Journal, 2014, 52(9): 1980-1993.

    [21] LEVIN D, KATZ J. A vortex-lattice method for the calculation of the nonsteady separated flow over delta wings: AIAA-1980-1803[R]. Reston: AIAA, 1980.

    [22] 曹志遠(yuǎn). 板殼振動(dòng)理論[M]. 北京:中國鐵道出版社, 1989: 455-457. CAO Z Y. Vibration theory of plates and shells[M]. Beijing: China Railway Press, 1989: 455-457 (in Chinese).

    [23] XIE D, XU M, DAI H, et al. Observation and evolution of chaos for a cantilever plate in supersonic flow[J]. Journal of Fluids and Structures, 2014, 50(6): 271-291.

    [24] HALLISSY P B, CESNIK C E S. High-fidelity aeroelastic analysis of very flexible aircraft: AIAA-2011-1914[R]. Reston: AIAA, 2011.

    (責(zé)任編輯:李明敏)

    *Corresponding author. E-mail: xumin@nwpu.edu.cn

    Post flutter response of flexible cantilever plate inlow speed flow field

    CHEN Tao, XU Min*, XIE Dan, AN Xiaomin

    SchoolofAstronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China

    A new nonlinear aeroelastic model is established for the analysis of post flutter behavior of the flexible cantilever plate in low speed flows. Geometric nonlinearity, aerodynamic nonlinearity as well as the strong coupling effect between geometric and aerodynamic nonlinearities are considered in the modeling. The aeroelastic model is verified with experimental data. The study shows that the flexible cantilever plate in low speed flows can enter into chaotic movement by means of periodic doubling in structural response. It is found that geometric nonlinearity and the unsteady aerodynamic effect caused by tip vortexes have significant impact on the structural response, while the unsteady aerodynamic force caused by deformation of wake vortexes has slight impact on the structural response. Different coupling strategies are also investigated to provide a reference for the selection of coupling strategy in nonlinear aeroelastic simulations of low-aspect-ratio flexible structures in low speed flows.

    low speed flow field; flexible cantilever plate; post flutter; chaotic movement; strong coupling

    2016-03-11; Revised:2016-10-10; Accepted:2016-11-14; Published online:2016-11-22 10:56

    URL:www.cnki.net/kcms/detail/11.1929.V.20161122.1056.002.html

    s:National Natural Science Foundation of China (11202165, 11402212)

    http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

    10.7527/S1000-6893.2016.0296

    2016-03-11; 退修日期:2016-10-10; 錄用日期:2016-11-14; 網(wǎng)絡(luò)出版時(shí)間:2016-11-22 10:56

    www.cnki.net/kcms/detail/11.1929.V.20161122.1056.002.html

    國家自然科學(xué)基金 (11202165, 11402212)

    *通訊作者.E-mail: xumin@nwpu.edu.cn

    陳濤, 徐敏, 謝丹, 等. 低速流場中柔性懸臂板的后顫振響應(yīng)[J]. 航空學(xué)報(bào), 2017, 38(3): 120215. CHEN T, XU M, XIE D, et al. Post flutter response of flexible cantilever plate in low speed flow field[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(3): 120215.

    V211.47

    A

    1000-6893(2017)03-120215-14

    猜你喜歡
    氣動(dòng)彈性氣動(dòng)力懸臂
    飛行載荷外部氣動(dòng)力的二次規(guī)劃等效映射方法
    懸臂式硫化罐的開發(fā)設(shè)計(jì)
    當(dāng)液壓遇上懸臂云臺(tái) 捷信GHFG1液壓懸臂云臺(tái)試用
    側(cè)風(fēng)對拍動(dòng)翅氣動(dòng)力的影響
    探討掛籃懸臂灌注連梁的施工
    LF爐懸臂爐蓋防傾覆結(jié)構(gòu)
    飛翼無人機(jī)嗡鳴氣動(dòng)彈性響應(yīng)分析
    模態(tài)選取對靜氣動(dòng)彈性分析的影響
    直升機(jī)的氣動(dòng)彈性問題
    大型風(fēng)力機(jī)整機(jī)氣動(dòng)彈性響應(yīng)計(jì)算
    亚洲精品一二三| 操美女的视频在线观看| 51午夜福利影视在线观看| 国产黄色免费在线视频| 又大又黄又爽视频免费| bbb黄色大片| 亚洲欧洲国产日韩| 国产男人的电影天堂91| 青春草亚洲视频在线观看| 国产成人av激情在线播放| 亚洲av男天堂| 成年女人毛片免费观看观看9 | 在线av久久热| 又粗又硬又长又爽又黄的视频| 一区二区av电影网| 考比视频在线观看| kizo精华| 欧美97在线视频| 免费黄频网站在线观看国产| 午夜福利一区二区在线看| 在线天堂中文资源库| 天天操日日干夜夜撸| 成人黄色视频免费在线看| 日韩伦理黄色片| 在线观看免费午夜福利视频| 午夜激情久久久久久久| 看免费av毛片| 欧美日韩av久久| av电影中文网址| 欧美国产精品va在线观看不卡| 超碰成人久久| 在线观看www视频免费| 纯流量卡能插随身wifi吗| 王馨瑶露胸无遮挡在线观看| 国产精品99久久99久久久不卡| 中文字幕制服av| 精品亚洲乱码少妇综合久久| 精品人妻在线不人妻| 老熟女久久久| 精品一品国产午夜福利视频| 捣出白浆h1v1| 亚洲七黄色美女视频| 视频区欧美日本亚洲| 美女午夜性视频免费| 免费黄频网站在线观看国产| 蜜桃在线观看..| 欧美激情 高清一区二区三区| 久久久久精品人妻al黑| 你懂的网址亚洲精品在线观看| 一级片免费观看大全| 精品国产乱码久久久久久男人| 久久久国产精品麻豆| 男的添女的下面高潮视频| 一区在线观看完整版| 亚洲综合色网址| 最近最新中文字幕大全免费视频 | 一级毛片黄色毛片免费观看视频| 亚洲精品国产区一区二| 99热全是精品| 一级a爱视频在线免费观看| 天堂俺去俺来也www色官网| 国产亚洲精品第一综合不卡| 亚洲久久久国产精品| 久热这里只有精品99| 蜜桃国产av成人99| 国产成人系列免费观看| 国产日韩欧美亚洲二区| 日韩电影二区| 亚洲欧美色中文字幕在线| 精品少妇黑人巨大在线播放| 久久久欧美国产精品| 午夜福利在线免费观看网站| 美女主播在线视频| 国产视频首页在线观看| 脱女人内裤的视频| 欧美日韩av久久| 黑人巨大精品欧美一区二区蜜桃| 18禁国产床啪视频网站| 国产男女超爽视频在线观看| 亚洲成人免费av在线播放| 男女午夜视频在线观看| 国产欧美日韩一区二区三区在线| 日韩欧美一区视频在线观看| 狠狠婷婷综合久久久久久88av| 日韩电影二区| 深夜精品福利| 在线看a的网站| 天天操日日干夜夜撸| 久久天躁狠狠躁夜夜2o2o | 国产一区二区 视频在线| 午夜91福利影院| 久久久精品区二区三区| 国产真人三级小视频在线观看| 岛国毛片在线播放| 成人手机av| www.精华液| 亚洲一区中文字幕在线| 新久久久久国产一级毛片| 人体艺术视频欧美日本| 如日韩欧美国产精品一区二区三区| www.精华液| 欧美性长视频在线观看| 欧美 日韩 精品 国产| 黄色视频不卡| 日本午夜av视频| 搡老乐熟女国产| 国产精品熟女久久久久浪| 久久国产精品男人的天堂亚洲| 国产高清不卡午夜福利| 免费在线观看黄色视频的| 青青草视频在线视频观看| 又紧又爽又黄一区二区| 悠悠久久av| 性高湖久久久久久久久免费观看| 香蕉国产在线看| 黑人猛操日本美女一级片| 色婷婷久久久亚洲欧美| 亚洲国产毛片av蜜桃av| 亚洲av日韩精品久久久久久密 | 久久人妻福利社区极品人妻图片 | 国产亚洲精品第一综合不卡| 人妻一区二区av| 亚洲人成电影观看| 一级毛片女人18水好多 | 99国产精品一区二区蜜桃av | 一区二区三区四区激情视频| 亚洲成人免费电影在线观看 | 色视频在线一区二区三区| 97人妻天天添夜夜摸| 亚洲成人免费av在线播放| 人人妻人人澡人人爽人人夜夜| 久久99精品国语久久久| 日日爽夜夜爽网站| 免费看av在线观看网站| 人妻一区二区av| 亚洲人成电影观看| 丰满迷人的少妇在线观看| 婷婷成人精品国产| 亚洲久久久国产精品| 在线观看免费日韩欧美大片| 国产片内射在线| 免费在线观看视频国产中文字幕亚洲 | 日韩av不卡免费在线播放| 色综合欧美亚洲国产小说| 午夜91福利影院| 韩国精品一区二区三区| 99国产精品免费福利视频| 亚洲欧美一区二区三区黑人| 国产高清videossex| 一级a爱视频在线免费观看| 考比视频在线观看| 久久精品成人免费网站| 大片电影免费在线观看免费| 97精品久久久久久久久久精品| 婷婷成人精品国产| 亚洲精品国产一区二区精华液| 99久久99久久久精品蜜桃| 国产亚洲av片在线观看秒播厂| 黄色一级大片看看| bbb黄色大片| 91精品国产国语对白视频| 国产精品一二三区在线看| 亚洲人成网站在线观看播放| 电影成人av| 国产精品免费大片| 亚洲欧美中文字幕日韩二区| 在线精品无人区一区二区三| 后天国语完整版免费观看| 人人澡人人妻人| 天堂俺去俺来也www色官网| 美女高潮到喷水免费观看| 99热国产这里只有精品6| 大片电影免费在线观看免费| 少妇 在线观看| 精品人妻在线不人妻| 色精品久久人妻99蜜桃| 自拍欧美九色日韩亚洲蝌蚪91| 一边摸一边做爽爽视频免费| 一区二区日韩欧美中文字幕| 亚洲精品自拍成人| www.精华液| 老司机靠b影院| 两个人免费观看高清视频| 亚洲国产欧美一区二区综合| 国产亚洲欧美精品永久| 丰满饥渴人妻一区二区三| 人人妻人人爽人人添夜夜欢视频| 青草久久国产| 精品少妇久久久久久888优播| 久久久久国产精品人妻一区二区| 日韩,欧美,国产一区二区三区| 69精品国产乱码久久久| 啦啦啦中文免费视频观看日本| 在线亚洲精品国产二区图片欧美| 成人亚洲精品一区在线观看| xxx大片免费视频| 在线 av 中文字幕| 国产av一区二区精品久久| 亚洲国产欧美一区二区综合| 国产日韩一区二区三区精品不卡| 亚洲精品国产区一区二| 国产欧美亚洲国产| 亚洲美女黄色视频免费看| 精品少妇久久久久久888优播| 婷婷成人精品国产| 19禁男女啪啪无遮挡网站| 久9热在线精品视频| 看免费av毛片| 一级黄片播放器| 久久久欧美国产精品| 国产成人精品久久二区二区91| 亚洲成国产人片在线观看| 我要看黄色一级片免费的| 午夜福利视频在线观看免费| 国产免费福利视频在线观看| 美女扒开内裤让男人捅视频| 久久久久国产精品人妻一区二区| 99re6热这里在线精品视频| 麻豆国产av国片精品| 精品少妇久久久久久888优播| 欧美 日韩 精品 国产| 国产亚洲欧美在线一区二区| 在线av久久热| 九色亚洲精品在线播放| 蜜桃在线观看..| 18在线观看网站| 国产精品秋霞免费鲁丝片| 中文精品一卡2卡3卡4更新| 少妇的丰满在线观看| 人人妻人人添人人爽欧美一区卜| 夜夜骑夜夜射夜夜干| 亚洲一区中文字幕在线| 国产无遮挡羞羞视频在线观看| 乱人伦中国视频| 免费看av在线观看网站| 亚洲精品一区蜜桃| 在线观看人妻少妇| 一区二区日韩欧美中文字幕| 久久精品人人爽人人爽视色| 99国产精品一区二区蜜桃av | 波多野结衣av一区二区av| tube8黄色片| 精品人妻熟女毛片av久久网站| 欧美激情高清一区二区三区| 亚洲黑人精品在线| 久久精品久久久久久噜噜老黄| 国产在线视频一区二区| 国产又爽黄色视频| 久久国产精品大桥未久av| 操美女的视频在线观看| av不卡在线播放| 老汉色∧v一级毛片| 国产日韩欧美在线精品| 国产一卡二卡三卡精品| 欧美精品一区二区大全| 国产精品亚洲av一区麻豆| 国产日韩一区二区三区精品不卡| 夜夜骑夜夜射夜夜干| 韩国精品一区二区三区| 丝瓜视频免费看黄片| 亚洲欧美色中文字幕在线| 五月开心婷婷网| 欧美人与善性xxx| 一区二区日韩欧美中文字幕| 90打野战视频偷拍视频| 成在线人永久免费视频| 亚洲中文av在线| 日本黄色日本黄色录像| 国产在线免费精品| 在线看a的网站| 国产精品成人在线| 黄色片一级片一级黄色片| 中文字幕最新亚洲高清| 精品第一国产精品| 大香蕉久久成人网| 一区二区三区精品91| 精品久久蜜臀av无| 国产色视频综合| 中文字幕色久视频| 欧美久久黑人一区二区| 国产成人啪精品午夜网站| 一区二区三区四区激情视频| av天堂在线播放| 交换朋友夫妻互换小说| 丰满迷人的少妇在线观看| 十八禁高潮呻吟视频| 狠狠婷婷综合久久久久久88av| av在线播放精品| 国产真人三级小视频在线观看| 精品免费久久久久久久清纯 | 欧美人与性动交α欧美精品济南到| 欧美精品av麻豆av| 久久国产精品人妻蜜桃| 女人高潮潮喷娇喘18禁视频| 夫妻性生交免费视频一级片| 国产一区亚洲一区在线观看| 在现免费观看毛片| 色网站视频免费| 一本色道久久久久久精品综合| 亚洲专区国产一区二区| 十八禁网站网址无遮挡| 亚洲国产欧美在线一区| 97精品久久久久久久久久精品| 另类精品久久| 亚洲专区国产一区二区| 亚洲精品美女久久久久99蜜臀 | 69精品国产乱码久久久| 高清av免费在线| xxxhd国产人妻xxx| 精品福利观看| 色播在线永久视频| 在线观看www视频免费| 久久精品久久久久久久性| 国产日韩欧美在线精品| 精品国产一区二区三区四区第35| 午夜两性在线视频| www.精华液| 黄色a级毛片大全视频| 国语对白做爰xxxⅹ性视频网站| 黄色片一级片一级黄色片| 亚洲美女黄色视频免费看| 18在线观看网站| 中文字幕色久视频| 色综合欧美亚洲国产小说| 99久久综合免费| 999精品在线视频| 国产欧美亚洲国产| 久久99一区二区三区| 亚洲欧洲国产日韩| 国语对白做爰xxxⅹ性视频网站| 美国免费a级毛片| 欧美日韩亚洲国产一区二区在线观看 | 欧美精品av麻豆av| 在现免费观看毛片| 天天躁狠狠躁夜夜躁狠狠躁| 欧美精品人与动牲交sv欧美| 国产成人欧美在线观看 | 尾随美女入室| 日韩一本色道免费dvd| 亚洲成国产人片在线观看| 久久精品国产综合久久久| 日韩制服丝袜自拍偷拍| 成人亚洲欧美一区二区av| 黄色片一级片一级黄色片| 国产高清videossex| 午夜激情av网站| 国产日韩欧美在线精品| 人体艺术视频欧美日本| 中文字幕人妻熟女乱码| 搡老乐熟女国产| 色视频在线一区二区三区| 少妇粗大呻吟视频| 少妇被粗大的猛进出69影院| 久久免费观看电影| av片东京热男人的天堂| 中文字幕色久视频| 大码成人一级视频| 婷婷色麻豆天堂久久| 一级毛片电影观看| 在线 av 中文字幕| 亚洲伊人久久精品综合| 黄色a级毛片大全视频| 日韩中文字幕欧美一区二区 | 国产免费又黄又爽又色| 秋霞在线观看毛片| 无遮挡黄片免费观看| 一二三四社区在线视频社区8| 一级,二级,三级黄色视频| 亚洲国产中文字幕在线视频| 午夜免费鲁丝| 欧美日韩精品网址| 亚洲精品一卡2卡三卡4卡5卡 | 国产视频一区二区在线看| 亚洲精品中文字幕在线视频| 少妇被粗大的猛进出69影院| 熟女少妇亚洲综合色aaa.| 免费高清在线观看日韩| 19禁男女啪啪无遮挡网站| 黑丝袜美女国产一区| 日本色播在线视频| 高清黄色对白视频在线免费看| 亚洲精品日韩在线中文字幕| 国产成人av教育| 国产97色在线日韩免费| 久久国产精品人妻蜜桃| 中文字幕精品免费在线观看视频| a级毛片黄视频| 精品国产超薄肉色丝袜足j| 人妻一区二区av| 成在线人永久免费视频| 亚洲av国产av综合av卡| 国产精品秋霞免费鲁丝片| 日日爽夜夜爽网站| 我要看黄色一级片免费的| 亚洲第一av免费看| 亚洲av日韩在线播放| 亚洲国产av影院在线观看| 国产成人精品久久二区二区免费| 黑人巨大精品欧美一区二区蜜桃| 精品视频人人做人人爽| 亚洲av美国av| 欧美xxⅹ黑人| 别揉我奶头~嗯~啊~动态视频 | 97在线人人人人妻| 蜜桃在线观看..| 女警被强在线播放| 高清欧美精品videossex| 91麻豆精品激情在线观看国产 | 观看av在线不卡| 国产99久久九九免费精品| 中文字幕亚洲精品专区| 女人久久www免费人成看片| 爱豆传媒免费全集在线观看| 我要看黄色一级片免费的| 91成人精品电影| 美女大奶头黄色视频| 久久ye,这里只有精品| 国产精品一区二区精品视频观看| 国产男女超爽视频在线观看| 搡老乐熟女国产| 巨乳人妻的诱惑在线观看| a级毛片在线看网站| 丝袜美足系列| 高清视频免费观看一区二区| 国产精品免费大片| 99久久99久久久精品蜜桃| 欧美激情极品国产一区二区三区| 色婷婷av一区二区三区视频| 日韩伦理黄色片| 国产一区二区在线观看av| 国产片内射在线| 天天躁日日躁夜夜躁夜夜| 亚洲精品国产av成人精品| 丝袜脚勾引网站| 亚洲国产看品久久| 久久精品成人免费网站| 国产精品三级大全| 亚洲一区中文字幕在线| 97在线人人人人妻| 久久久久精品人妻al黑| 交换朋友夫妻互换小说| 天天躁夜夜躁狠狠久久av| 午夜老司机福利片| 久久国产精品人妻蜜桃| a级毛片黄视频| 嫩草影视91久久| 国产av一区二区精品久久| 伊人亚洲综合成人网| 18禁黄网站禁片午夜丰满| 午夜91福利影院| 91成人精品电影| 国产亚洲欧美在线一区二区| 最新在线观看一区二区三区 | 日韩av免费高清视频| 久久毛片免费看一区二区三区| 久久免费观看电影| 久久99热这里只频精品6学生| 18禁观看日本| 男女免费视频国产| www.精华液| 欧美 日韩 精品 国产| 99精国产麻豆久久婷婷| 日韩中文字幕视频在线看片| 9热在线视频观看99| 久久毛片免费看一区二区三区| 日本色播在线视频| 国产免费视频播放在线视频| 午夜精品国产一区二区电影| 中文字幕人妻丝袜制服| 精品一区在线观看国产| 国产精品久久久久久精品古装| 国产一区二区激情短视频 | 美女高潮到喷水免费观看| 日本wwww免费看| 国产老妇伦熟女老妇高清| 日韩电影二区| 亚洲天堂av无毛| 午夜91福利影院| 宅男免费午夜| 欧美日韩一级在线毛片| 男人爽女人下面视频在线观看| 久久精品国产亚洲av涩爱| 可以免费在线观看a视频的电影网站| 久久久久久免费高清国产稀缺| 中文字幕另类日韩欧美亚洲嫩草| 亚洲国产欧美网| 1024香蕉在线观看| xxx大片免费视频| 欧美日韩亚洲高清精品| 久久99一区二区三区| 久久 成人 亚洲| 国产午夜精品一二区理论片| 国产国语露脸激情在线看| 免费在线观看视频国产中文字幕亚洲 | av网站免费在线观看视频| 青春草亚洲视频在线观看| 十八禁高潮呻吟视频| 男女边摸边吃奶| 啦啦啦视频在线资源免费观看| 日韩精品免费视频一区二区三区| 国产精品偷伦视频观看了| 咕卡用的链子| 男女床上黄色一级片免费看| 少妇的丰满在线观看| 三上悠亚av全集在线观看| 母亲3免费完整高清在线观看| 免费一级毛片在线播放高清视频 | 日韩电影二区| 欧美成人精品欧美一级黄| 亚洲av欧美aⅴ国产| 欧美日韩国产mv在线观看视频| 国产精品久久久人人做人人爽| 久久久精品国产亚洲av高清涩受| 激情视频va一区二区三区| 欧美日本中文国产一区发布| 国产亚洲欧美精品永久| 乱人伦中国视频| 欧美在线一区亚洲| 男女免费视频国产| 亚洲伊人久久精品综合| 99久久综合免费| 99国产精品免费福利视频| 激情五月婷婷亚洲| 精品亚洲乱码少妇综合久久| 日韩av在线免费看完整版不卡| 久久 成人 亚洲| 欧美黑人精品巨大| 麻豆国产av国片精品| 欧美精品啪啪一区二区三区 | 国产91精品成人一区二区三区 | av网站免费在线观看视频| 亚洲av片天天在线观看| 国产在线视频一区二区| 久久久久精品国产欧美久久久 | 亚洲欧美色中文字幕在线| 十八禁人妻一区二区| 免费日韩欧美在线观看| 国产欧美日韩精品亚洲av| 一边亲一边摸免费视频| 国产欧美日韩一区二区三区在线| 亚洲一区二区三区欧美精品| 真人做人爱边吃奶动态| avwww免费| e午夜精品久久久久久久| 少妇粗大呻吟视频| 1024视频免费在线观看| 一区二区三区乱码不卡18| 热99久久久久精品小说推荐| 免费久久久久久久精品成人欧美视频| 国产高清视频在线播放一区 | 久久国产精品人妻蜜桃| 日本五十路高清| 香蕉国产在线看| 一二三四在线观看免费中文在| 五月开心婷婷网| 久久影院123| 欧美精品高潮呻吟av久久| 母亲3免费完整高清在线观看| 一级毛片我不卡| 亚洲成av片中文字幕在线观看| 国产男女内射视频| 男女午夜视频在线观看| 国产精品香港三级国产av潘金莲 | 观看av在线不卡| 国产99久久九九免费精品| 日韩欧美一区视频在线观看| 精品一区二区三区四区五区乱码 | 男人爽女人下面视频在线观看| 国产在线视频一区二区| 操美女的视频在线观看| 亚洲国产欧美日韩在线播放| 国产免费现黄频在线看| 久久精品亚洲av国产电影网| 国产成人免费观看mmmm| 亚洲国产看品久久| 久久人人爽av亚洲精品天堂| 老司机在亚洲福利影院| 免费在线观看完整版高清| 亚洲,一卡二卡三卡| 久久女婷五月综合色啪小说| 高清欧美精品videossex| 国产女主播在线喷水免费视频网站| 97精品久久久久久久久久精品| 亚洲色图 男人天堂 中文字幕| 伊人久久大香线蕉亚洲五| 国产一区亚洲一区在线观看| 午夜福利一区二区在线看| 亚洲精品久久午夜乱码| 免费高清在线观看日韩| 美女国产高潮福利片在线看| 丝袜人妻中文字幕| 午夜精品国产一区二区电影| 十分钟在线观看高清视频www| 国产主播在线观看一区二区 | 天天躁狠狠躁夜夜躁狠狠躁| 久久久久精品国产欧美久久久 | 国产一卡二卡三卡精品| 免费在线观看日本一区| 天天添夜夜摸| 妹子高潮喷水视频| 亚洲国产中文字幕在线视频| 精品国产一区二区三区久久久樱花| 久久久久国产精品人妻一区二区| 国产精品亚洲av一区麻豆| 19禁男女啪啪无遮挡网站| 亚洲专区国产一区二区| 各种免费的搞黄视频| 五月开心婷婷网| 国产亚洲一区二区精品| 久久久国产一区二区| 国产高清不卡午夜福利| 国产成人欧美| 欧美中文综合在线视频| 亚洲九九香蕉| 国产视频一区二区在线看|