• <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ì)算
    在线观看人妻少妇| 一本—道久久a久久精品蜜桃钙片| 街头女战士在线观看网站| 各种免费的搞黄视频| 国产国拍精品亚洲av在线观看| 啦啦啦在线观看免费高清www| 少妇猛男粗大的猛烈进出视频| 欧美成人一区二区免费高清观看| 亚洲欧美中文字幕日韩二区| 久久久久国产精品人妻一区二区| 久久久久久久久久久丰满| 免费看光身美女| 国产极品天堂在线| 91久久精品国产一区二区三区| 五月天丁香电影| 成年免费大片在线观看| 乱码一卡2卡4卡精品| 五月开心婷婷网| 22中文网久久字幕| 国产精品福利在线免费观看| 欧美一区二区亚洲| 卡戴珊不雅视频在线播放| 国产亚洲一区二区精品| 中国美白少妇内射xxxbb| 日韩av在线免费看完整版不卡| 国产精品一及| 一级毛片我不卡| 高清毛片免费看| 欧美激情极品国产一区二区三区 | 国产精品久久久久成人av| 亚洲不卡免费看| www.av在线官网国产| 亚洲欧美精品专区久久| 人人妻人人添人人爽欧美一区卜 | 三级国产精品欧美在线观看| 一级毛片电影观看| 国产中年淑女户外野战色| 一级毛片aaaaaa免费看小| 日本黄大片高清| 精品午夜福利在线看| 日韩免费高清中文字幕av| 51国产日韩欧美| 中文字幕久久专区| 国产永久视频网站| 91久久精品国产一区二区三区| 狂野欧美白嫩少妇大欣赏| 国产淫语在线视频| 亚洲一区二区三区欧美精品| 啦啦啦中文免费视频观看日本| 成人二区视频| 黑丝袜美女国产一区| 一区二区av电影网| 最后的刺客免费高清国语| 亚洲图色成人| 人妻夜夜爽99麻豆av| 国产精品人妻久久久影院| 日本一二三区视频观看| 黄色欧美视频在线观看| 午夜福利高清视频| 国产精品一区二区在线不卡| 成年免费大片在线观看| 国产精品一区二区三区四区免费观看| 99久久精品热视频| 97热精品久久久久久| 久久精品国产鲁丝片午夜精品| 三级经典国产精品| 日日撸夜夜添| 国产精品熟女久久久久浪| 国产成人a区在线观看| 欧美xxⅹ黑人| 久久ye,这里只有精品| 亚洲精品成人av观看孕妇| 另类亚洲欧美激情| 国产一区二区三区av在线| 校园人妻丝袜中文字幕| 久久精品国产亚洲av天美| 人人妻人人澡人人爽人人夜夜| 日韩中字成人| 国产精品欧美亚洲77777| 黄色怎么调成土黄色| 美女福利国产在线 | 国产精品久久久久久久久免| 精品亚洲成a人片在线观看 | 夫妻午夜视频| 啦啦啦视频在线资源免费观看| 大香蕉97超碰在线| 18禁裸乳无遮挡免费网站照片| 丰满少妇做爰视频| 少妇 在线观看| 亚洲av不卡在线观看| h视频一区二区三区| 看免费成人av毛片| 色婷婷av一区二区三区视频| 人妻少妇偷人精品九色| 色5月婷婷丁香| 中文欧美无线码| 久久97久久精品| 99久久人妻综合| 国产有黄有色有爽视频| 亚洲国产欧美人成| 人妻系列 视频| 亚洲性久久影院| 日本午夜av视频| 高清不卡的av网站| 又黄又爽又刺激的免费视频.| 大话2 男鬼变身卡| 欧美一区二区亚洲| 国产爽快片一区二区三区| av天堂中文字幕网| 一级片'在线观看视频| 乱系列少妇在线播放| 午夜免费男女啪啪视频观看| 欧美bdsm另类| 最后的刺客免费高清国语| 在线观看一区二区三区激情| 午夜福利在线在线| 六月丁香七月| 人妻 亚洲 视频| 国产极品天堂在线| 免费观看无遮挡的男女| tube8黄色片| 国产免费一级a男人的天堂| 97在线视频观看| 天堂俺去俺来也www色官网| 国产成人freesex在线| 亚洲欧美成人综合另类久久久| 婷婷色综合www| 亚洲三级黄色毛片| 街头女战士在线观看网站| 啦啦啦视频在线资源免费观看| 国产亚洲av片在线观看秒播厂| 国产高潮美女av| 国产亚洲最大av| 超碰97精品在线观看| 女人久久www免费人成看片| 免费不卡的大黄色大毛片视频在线观看| 下体分泌物呈黄色| 亚洲欧美日韩东京热| 九九久久精品国产亚洲av麻豆| 中文字幕精品免费在线观看视频 | 欧美一区二区亚洲| 精品亚洲成国产av| 午夜福利在线在线| 性色avwww在线观看| 国产大屁股一区二区在线视频| 寂寞人妻少妇视频99o| 欧美少妇被猛烈插入视频| 亚洲av综合色区一区| 少妇高潮的动态图| 啦啦啦在线观看免费高清www| 三级国产精品片| 五月开心婷婷网| 精品午夜福利在线看| 韩国av在线不卡| 成人高潮视频无遮挡免费网站| 精品一区二区免费观看| 最近的中文字幕免费完整| 91狼人影院| 亚洲伊人久久精品综合| 97精品久久久久久久久久精品| 欧美激情国产日韩精品一区| 国产成人精品婷婷| 亚洲精品国产av蜜桃| 日本-黄色视频高清免费观看| 亚洲精品国产成人久久av| 特大巨黑吊av在线直播| 国产免费一级a男人的天堂| 国产综合精华液| av在线app专区| 寂寞人妻少妇视频99o| 深爱激情五月婷婷| 国产色爽女视频免费观看| 国产乱人视频| 午夜老司机福利剧场| 80岁老熟妇乱子伦牲交| 亚洲怡红院男人天堂| 国产老妇伦熟女老妇高清| 国产毛片在线视频| 日韩三级伦理在线观看| 极品教师在线视频| 亚洲电影在线观看av| 亚洲精品国产色婷婷电影| 蜜桃久久精品国产亚洲av| 国产熟女欧美一区二区| 精品国产乱码久久久久久小说| 久久久久精品久久久久真实原创| 久热这里只有精品99| 一级毛片黄色毛片免费观看视频| av免费观看日本| 爱豆传媒免费全集在线观看| 在线观看美女被高潮喷水网站| 一区二区三区免费毛片| 啦啦啦啦在线视频资源| 国产精品国产av在线观看| 亚洲av欧美aⅴ国产| xxx大片免费视频| 夜夜骑夜夜射夜夜干| 精品人妻一区二区三区麻豆| 成人亚洲精品一区在线观看 | 免费久久久久久久精品成人欧美视频 | 在线观看免费日韩欧美大片 | 九九在线视频观看精品| 99九九线精品视频在线观看视频| 免费少妇av软件| 亚洲天堂av无毛| 热99国产精品久久久久久7| 大又大粗又爽又黄少妇毛片口| 在线看a的网站| 高清午夜精品一区二区三区| 亚洲人成网站在线观看播放| 精品熟女少妇av免费看| 一区二区av电影网| 亚洲高清免费不卡视频| 亚洲美女搞黄在线观看| 少妇高潮的动态图| 最新中文字幕久久久久| 卡戴珊不雅视频在线播放| 国产爽快片一区二区三区| 国产成人免费无遮挡视频| 国产欧美日韩一区二区三区在线 | 亚洲,一卡二卡三卡| 久久精品国产鲁丝片午夜精品| 成人亚洲欧美一区二区av| 国产男女内射视频| 亚洲成人一二三区av| 国产精品.久久久| 草草在线视频免费看| 亚洲婷婷狠狠爱综合网| 97在线人人人人妻| 一个人看的www免费观看视频| 哪个播放器可以免费观看大片| 在线天堂最新版资源| 丰满人妻一区二区三区视频av| 亚洲高清免费不卡视频| 中文字幕亚洲精品专区| 大陆偷拍与自拍| 又粗又硬又长又爽又黄的视频| 亚洲欧美日韩另类电影网站 | 欧美xxⅹ黑人| 大片电影免费在线观看免费| 大香蕉97超碰在线| 国产伦精品一区二区三区四那| 久久久亚洲精品成人影院| 伦理电影免费视频| 嘟嘟电影网在线观看| 成人亚洲精品一区在线观看 | 狂野欧美白嫩少妇大欣赏| 欧美人与善性xxx| 日韩中字成人| 亚洲国产最新在线播放| 少妇精品久久久久久久| 成人毛片a级毛片在线播放| 欧美少妇被猛烈插入视频| 天堂俺去俺来也www色官网| 内射极品少妇av片p| 在线观看一区二区三区| av卡一久久| 亚洲不卡免费看| 亚洲av二区三区四区| 日韩三级伦理在线观看| 国产精品久久久久成人av| 免费人成在线观看视频色| 国产女主播在线喷水免费视频网站| 一边亲一边摸免费视频| 国产美女午夜福利| 亚洲精华国产精华液的使用体验| 亚洲一区二区三区欧美精品| 啦啦啦在线观看免费高清www| 久久久久网色| 国产伦精品一区二区三区四那| 男女边摸边吃奶| 伦理电影免费视频| 精品人妻熟女av久视频| 联通29元200g的流量卡| 我的老师免费观看完整版| 久久久国产一区二区| 国产美女午夜福利| 亚洲性久久影院| 免费观看a级毛片全部| 久久精品久久精品一区二区三区| 熟妇人妻不卡中文字幕| 波野结衣二区三区在线| 丰满人妻一区二区三区视频av| 精品人妻熟女av久视频| 91狼人影院| 少妇被粗大猛烈的视频| 美女内射精品一级片tv| 一级爰片在线观看| 天天躁夜夜躁狠狠久久av| 青青草视频在线视频观看| 老师上课跳d突然被开到最大视频| .国产精品久久| 十分钟在线观看高清视频www | 搡老乐熟女国产| 18禁动态无遮挡网站| 亚洲电影在线观看av| 精品久久国产蜜桃| a级毛色黄片| 日韩,欧美,国产一区二区三区| 26uuu在线亚洲综合色| 亚洲高清免费不卡视频| 国产成人a区在线观看| 久久精品人妻少妇| 成人漫画全彩无遮挡| 18禁动态无遮挡网站| 肉色欧美久久久久久久蜜桃| 高清av免费在线| 亚洲精品456在线播放app| 国产熟女欧美一区二区| 国产男女内射视频| 欧美丝袜亚洲另类| 久久国产精品大桥未久av | 久久精品久久精品一区二区三区| 又黄又爽又刺激的免费视频.| 中文欧美无线码| 热99国产精品久久久久久7| 内射极品少妇av片p| 亚洲av男天堂| 最近手机中文字幕大全| 人人妻人人看人人澡| 日韩制服骚丝袜av| 亚洲四区av| 日韩成人伦理影院| 亚洲av成人精品一二三区| 夜夜骑夜夜射夜夜干| 精品一品国产午夜福利视频| 亚洲欧美精品专区久久| 一级毛片我不卡| 亚洲av中文字字幕乱码综合| 日本午夜av视频| 老司机影院成人| 亚洲国产精品999| 国产精品一及| 日韩一区二区三区影片| 黄色日韩在线| 国产亚洲最大av| 中文字幕免费在线视频6| 夜夜看夜夜爽夜夜摸| 久久人人爽人人爽人人片va| 特大巨黑吊av在线直播| 两个人的视频大全免费| xxx大片免费视频| 亚洲精品日本国产第一区| 久久毛片免费看一区二区三区| 日韩伦理黄色片| 18禁在线播放成人免费| 亚洲四区av| 在线播放无遮挡| 久久久久久九九精品二区国产| 最近最新中文字幕免费大全7| 热99国产精品久久久久久7| 亚洲性久久影院| 亚洲欧美日韩另类电影网站 | av网站免费在线观看视频| 男女国产视频网站| 内地一区二区视频在线| 亚洲国产高清在线一区二区三| 日韩强制内射视频| 99久久人妻综合| 午夜福利在线观看免费完整高清在| 激情 狠狠 欧美| 久久久久久人妻| 婷婷色麻豆天堂久久| 美女内射精品一级片tv| 国产欧美日韩精品一区二区| 街头女战士在线观看网站| 日本一二三区视频观看| 91狼人影院| 亚洲精品色激情综合| 国产精品一二三区在线看| av天堂中文字幕网| 舔av片在线| 一个人免费看片子| 波野结衣二区三区在线| 精品久久国产蜜桃| 欧美极品一区二区三区四区| 日韩欧美 国产精品| 国产成人免费无遮挡视频| 人人妻人人看人人澡| 日韩大片免费观看网站| 99热6这里只有精品| 免费观看在线日韩| 嫩草影院新地址| 亚洲丝袜综合中文字幕| 久久精品国产自在天天线| 插阴视频在线观看视频| 国产精品不卡视频一区二区| 国产精品国产三级国产av玫瑰| 日韩免费高清中文字幕av| 国产成人午夜福利电影在线观看| 中国三级夫妇交换| 美女视频免费永久观看网站| 纵有疾风起免费观看全集完整版| 六月丁香七月| kizo精华| 在线观看三级黄色| 亚洲人成网站在线播| 日韩av在线免费看完整版不卡| 好男人视频免费观看在线| 97超碰精品成人国产| 高清日韩中文字幕在线| 超碰97精品在线观看| av.在线天堂| 久久青草综合色| 午夜视频国产福利| 高清日韩中文字幕在线| 国产精品久久久久久精品古装| 欧美精品一区二区大全| 亚洲国产最新在线播放| 国产成人一区二区在线| 欧美精品国产亚洲| 欧美极品一区二区三区四区| 精品人妻偷拍中文字幕| 亚洲成人手机| 午夜免费观看性视频| 日本wwww免费看| 99热这里只有是精品在线观看| 午夜福利在线观看免费完整高清在| 亚洲aⅴ乱码一区二区在线播放| 成人美女网站在线观看视频| 3wmmmm亚洲av在线观看| 狂野欧美白嫩少妇大欣赏| 插逼视频在线观看| 最近中文字幕高清免费大全6| 精品一区二区三卡| 视频区图区小说| 热99国产精品久久久久久7| 网址你懂的国产日韩在线| 色视频在线一区二区三区| av视频免费观看在线观看| 一级毛片 在线播放| 啦啦啦啦在线视频资源| 日本黄大片高清| 六月丁香七月| 免费大片黄手机在线观看| 欧美少妇被猛烈插入视频| 涩涩av久久男人的天堂| av国产免费在线观看| 麻豆乱淫一区二区| 亚洲av欧美aⅴ国产| 18禁裸乳无遮挡免费网站照片| 国产成人精品福利久久| 亚洲国产精品一区三区| 国产精品麻豆人妻色哟哟久久| 最黄视频免费看| 国产片特级美女逼逼视频| 久久人妻熟女aⅴ| 毛片一级片免费看久久久久| 97超碰精品成人国产| 国产精品三级大全| 亚洲第一av免费看| 欧美亚洲 丝袜 人妻 在线| 久久久色成人| 你懂的网址亚洲精品在线观看| 啦啦啦视频在线资源免费观看| 亚洲久久久国产精品| 在线观看国产h片| 80岁老熟妇乱子伦牲交| 黑人高潮一二区| 亚洲不卡免费看| 精华霜和精华液先用哪个| 成年免费大片在线观看| 在线观看av片永久免费下载| 久久国内精品自在自线图片| 成人亚洲精品一区在线观看 | 国产精品久久久久久精品古装| 日韩av不卡免费在线播放| 99久国产av精品国产电影| 五月开心婷婷网| 精品久久久久久久久av| 免费在线观看成人毛片| av.在线天堂| 国模一区二区三区四区视频| 岛国毛片在线播放| 91精品一卡2卡3卡4卡| 国产午夜精品久久久久久一区二区三区| 日韩免费高清中文字幕av| 免费大片黄手机在线观看| 久久6这里有精品| 岛国毛片在线播放| av在线播放精品| 亚洲欧美日韩东京热| 18禁在线播放成人免费| 91精品一卡2卡3卡4卡| 久久ye,这里只有精品| 国产伦精品一区二区三区四那| 在线观看三级黄色| 国产 一区 欧美 日韩| 国产精品av视频在线免费观看| 欧美 日韩 精品 国产| 在线观看免费日韩欧美大片 | 久久久久国产精品人妻一区二区| 在线观看av片永久免费下载| 色婷婷久久久亚洲欧美| 大陆偷拍与自拍| 18禁在线无遮挡免费观看视频| 亚洲在久久综合| 日本av手机在线免费观看| 亚洲性久久影院| 激情 狠狠 欧美| 日韩成人伦理影院| 激情 狠狠 欧美| 高清视频免费观看一区二区| 久久影院123| 尾随美女入室| 超碰av人人做人人爽久久| 国产伦在线观看视频一区| 青青草视频在线视频观看| 亚洲国产精品999| 观看美女的网站| 国产伦在线观看视频一区| 精品少妇久久久久久888优播| 在线精品无人区一区二区三 | 久久鲁丝午夜福利片| 一区二区三区四区激情视频| 精品少妇黑人巨大在线播放| 一级毛片久久久久久久久女| 国产毛片在线视频| 欧美一区二区亚洲| 亚洲久久久国产精品| 精品少妇久久久久久888优播| 免费看光身美女| 天美传媒精品一区二区| 成人高潮视频无遮挡免费网站| 高清不卡的av网站| 一个人看的www免费观看视频| 大片免费播放器 马上看| 亚洲色图av天堂| 久久午夜福利片| 插逼视频在线观看| 我要看黄色一级片免费的| 欧美xxⅹ黑人| 亚洲欧美一区二区三区国产| 婷婷色麻豆天堂久久| 建设人人有责人人尽责人人享有的 | av.在线天堂| 一本色道久久久久久精品综合| 成人美女网站在线观看视频| 2018国产大陆天天弄谢| 日韩电影二区| 亚洲精品国产成人久久av| 人人妻人人看人人澡| 中国国产av一级| 免费看光身美女| 99久久综合免费| 成人亚洲欧美一区二区av| 中文字幕精品免费在线观看视频 | 午夜激情久久久久久久| 国产欧美日韩精品一区二区| 久久人人爽人人爽人人片va| 国产伦精品一区二区三区四那| 久久影院123| 免费看光身美女| 一区二区三区四区激情视频| 国产男女超爽视频在线观看| 欧美成人一区二区免费高清观看| 久久人妻熟女aⅴ| 欧美精品亚洲一区二区| 亚洲精品视频女| 性色av一级| 看非洲黑人一级黄片| 欧美xxxx黑人xx丫x性爽| 一级毛片黄色毛片免费观看视频| 男人和女人高潮做爰伦理| 欧美精品一区二区大全| 韩国av在线不卡| 街头女战士在线观看网站| 免费黄网站久久成人精品| 色吧在线观看| 国产欧美日韩精品一区二区| 狂野欧美激情性xxxx在线观看| 一个人看视频在线观看www免费| 美女福利国产在线 | 日日撸夜夜添| 大又大粗又爽又黄少妇毛片口| 午夜激情久久久久久久| 在线播放无遮挡| 91久久精品国产一区二区成人| 波野结衣二区三区在线| 色婷婷久久久亚洲欧美| 亚洲精品亚洲一区二区| 国产精品蜜桃在线观看| 一本—道久久a久久精品蜜桃钙片| h视频一区二区三区| 激情 狠狠 欧美| 女性生殖器流出的白浆| 日韩免费高清中文字幕av| 欧美精品人与动牲交sv欧美| 欧美丝袜亚洲另类| 日日啪夜夜撸| 99热这里只有精品一区| 国产精品国产三级国产av玫瑰| 人妻 亚洲 视频| 久久精品久久精品一区二区三区| 99精国产麻豆久久婷婷| 国产精品一区www在线观看| 国产黄频视频在线观看| 蜜桃久久精品国产亚洲av| 亚州av有码| 18+在线观看网站| 国产精品精品国产色婷婷| 女人久久www免费人成看片| 亚洲精品国产av成人精品| 男女边摸边吃奶| 亚洲精品456在线播放app| 99热网站在线观看| 国产伦精品一区二区三区四那| 欧美日韩视频高清一区二区三区二| 国产成人精品福利久久| 国产亚洲精品久久久com| freevideosex欧美| 亚洲av综合色区一区| 菩萨蛮人人尽说江南好唐韦庄| 日本欧美视频一区| 国产高清国产精品国产三级 |