亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        考慮靜氣動(dòng)彈性影響的客機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)研究

        2017-09-04 02:29:07楊體浩白俊強(qiáng)孫智偉史亞云西北工業(yè)大學(xué)航空學(xué)院陜西西安710072
        關(guān)鍵詞:氣動(dòng)彈性機(jī)翼構(gòu)型

        楊體浩, 白俊強(qiáng), 辛 亮, 孫智偉, 史亞云(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072)

        考慮靜氣動(dòng)彈性影響的客機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)研究

        楊體浩, 白俊強(qiáng)*, 辛 亮, 孫智偉, 史亞云
        (西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072)

        針對(duì)跨聲速客機(jī)氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)問(wèn)題,建立了考慮靜氣動(dòng)彈性影響的氣動(dòng)/結(jié)構(gòu)一體化優(yōu)化設(shè)計(jì)方法,并針對(duì)現(xiàn)代跨聲速民用客機(jī)開展了氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)研究。數(shù)值評(píng)估選擇全速勢(shì)方程加附面層修正,氣彈分析采用基于RBF插值技術(shù)的松耦合分析方法,優(yōu)化方法使用改進(jìn)的微分進(jìn)化算法。通過(guò)對(duì)CRM和DLR-F6標(biāo)模進(jìn)行計(jì)算并與實(shí)驗(yàn)數(shù)據(jù)對(duì)比,驗(yàn)證了采用的氣動(dòng)數(shù)值評(píng)估手段和靜氣動(dòng)彈性分析方法可靠性。利用建立的優(yōu)化設(shè)計(jì)方法對(duì)跨聲速客機(jī)機(jī)翼進(jìn)行了分別以扭轉(zhuǎn)角分布和剖面翼型為設(shè)計(jì)變量的考慮靜氣動(dòng)彈性影響的氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì),航程分別提高了5.63%和3.05%。航程的提高主要得益于機(jī)翼的載荷分布和結(jié)構(gòu)厚度分布的改變,以扭轉(zhuǎn)角分布為設(shè)計(jì)變量的優(yōu)化設(shè)計(jì)以2.56%的結(jié)構(gòu)重量損失獲得了6.53%的升阻比的提高,以剖面翼型外形為設(shè)計(jì)變量的優(yōu)化重量減小了3.56%同時(shí)升阻比提高了1.53%。

        多學(xué)科;空氣動(dòng)力學(xué);結(jié)構(gòu);全速勢(shì)方程;微分進(jìn)化算法;靜氣動(dòng)彈性

        0 引 言

        隨著航空燃油價(jià)格的不斷上升,改善燃油經(jīng)濟(jì)性是現(xiàn)代飛機(jī)設(shè)計(jì)追求的目標(biāo),采用展弦比更大的薄機(jī)翼有利于進(jìn)一步提高飛機(jī)的氣動(dòng)特性,然而機(jī)翼柔性的增大將會(huì)造成更為明顯的氣動(dòng)彈性變形,對(duì)飛機(jī)氣動(dòng)特性造成顯著的影響[1]。但是傳統(tǒng)的設(shè)計(jì)方法在進(jìn)行氣動(dòng)外形設(shè)計(jì)時(shí),并沒(méi)有很好的考慮靜氣動(dòng)彈性的影響,為了維持有利的巡航外形,往往采用增大機(jī)翼剛度的方法來(lái)保證飛機(jī)具有足夠的氣動(dòng)彈性穩(wěn)定性,結(jié)構(gòu)強(qiáng)度以及避免不同飛行狀態(tài)和機(jī)動(dòng)操縱時(shí)機(jī)翼發(fā)生不利的扭轉(zhuǎn)變形。而這種方式將會(huì)顯著增加機(jī)翼的結(jié)構(gòu)重量。因此隨著復(fù)合材料等新材料的大量使用,以及飛機(jī)柔性的增大,需要在初步設(shè)計(jì)階段就考慮柔性機(jī)翼的氣動(dòng)/結(jié)構(gòu)耦合問(wèn)題。

        對(duì)于柔性較大的飛機(jī),不同的飛行狀態(tài)對(duì)應(yīng)不同的載荷分布,進(jìn)而影響飛機(jī)的靜氣動(dòng)彈性變形,而靜氣動(dòng)彈性變形產(chǎn)生的扭轉(zhuǎn)角分布變化以及機(jī)翼明顯的上彎擾度變形所形成的弧形機(jī)翼[2-3]對(duì)飛機(jī)的氣動(dòng)特性會(huì)帶來(lái)明顯的影響。柔性機(jī)翼氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)充分利用了機(jī)翼的柔性特點(diǎn),結(jié)合氣動(dòng)彈性分析技術(shù),通過(guò)合理的氣動(dòng)、結(jié)構(gòu)耦合設(shè)計(jì)使得在多種飛行狀態(tài)下機(jī)翼都能夠獲得期望的氣動(dòng)彈性變形,以提高飛機(jī)的綜合性能。

        國(guó)內(nèi)外學(xué)者在氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)方面都做了大量的研究工作。董波等[4]將全速勢(shì)方程CFD數(shù)值評(píng)估方法與基于工程梁理論的機(jī)翼結(jié)構(gòu)設(shè)計(jì)方法相結(jié)合,進(jìn)行了氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì),但是沒(méi)有考慮靜氣動(dòng)彈性的影響。余雄慶[5]等利用代理模型預(yù)測(cè)氣動(dòng)力系數(shù)和翼尖位移、扭轉(zhuǎn)變形,以機(jī)翼剖面翼型扭轉(zhuǎn)角和結(jié)構(gòu)尺寸為設(shè)計(jì)變量進(jìn)行考慮靜氣動(dòng)彈性影響的氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)的研究,并對(duì)設(shè)計(jì)結(jié)果進(jìn)行了簡(jiǎn)單的描述。代理模型的采用在一定程度上提高了計(jì)算效率。陳海昕[6]等利用智能算法采用“三步”多目標(biāo)氣動(dòng)/結(jié)構(gòu)優(yōu)化設(shè)計(jì)策略對(duì)機(jī)翼平面形狀、機(jī)翼扭轉(zhuǎn)角、翼型剖面外形和結(jié)構(gòu)尺寸參數(shù)進(jìn)行了氣動(dòng)/結(jié)構(gòu)多學(xué)科優(yōu)化設(shè)計(jì)研究?!叭健眱?yōu)化策略在一定程度上將氣動(dòng),結(jié)構(gòu)進(jìn)行解耦,減少了具有大規(guī)模設(shè)計(jì)變量的復(fù)雜氣動(dòng)/結(jié)構(gòu)多學(xué)科優(yōu)化問(wèn)題的計(jì)算量。但是,這種優(yōu)化策略也在一定程度上削弱了氣動(dòng)、結(jié)構(gòu)之間的耦合作用。左英桃[7]利用離散伴隨方法將基于歐拉方程的氣動(dòng)力模型和有限元結(jié)構(gòu)求解器相結(jié)合,進(jìn)行了以剖面翼型和扭轉(zhuǎn)角為設(shè)計(jì)變量的,考慮了靜氣動(dòng)彈性影響的氣動(dòng)優(yōu)化設(shè)計(jì)研究。伴隨方法的使用顯著地提高了優(yōu)化設(shè)計(jì)的效率。經(jīng)過(guò)優(yōu)化設(shè)計(jì),全機(jī)升阻比提高了4.06%,但是整個(gè)優(yōu)化設(shè)計(jì)過(guò)程并沒(méi)考慮結(jié)構(gòu)參數(shù)的變化。

        在國(guó)外,以美國(guó)Martins團(tuán)隊(duì)為代表的課題組在以伴隨方法為基礎(chǔ)的氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)方法方面的研究較為成熟。Martins[8-9]利用伴隨方法建立了基于RANS方程和TACS (the Analysis of Commposite Structures)結(jié)構(gòu)求解器的針對(duì)復(fù)雜構(gòu)型大設(shè)計(jì)變量的氣動(dòng)/結(jié)構(gòu)優(yōu)化設(shè)計(jì)框架,并成功地應(yīng)用到了大展弦比跨聲速客機(jī)機(jī)翼的設(shè)計(jì)研究中。

        雖然目前國(guó)內(nèi)外在氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)方面有大量的研究可尋。但是,大多數(shù)研究工作都是圍繞著以提高氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)方法的設(shè)計(jì)效率、計(jì)算精度和計(jì)算復(fù)雜度開展。少有論文從應(yīng)用角度出發(fā),將氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)與單學(xué)科的優(yōu)化設(shè)計(jì)結(jié)果進(jìn)行對(duì)比研究,較詳細(xì)地分析考慮靜氣動(dòng)彈性影響之后,對(duì)機(jī)翼氣動(dòng)、結(jié)構(gòu)設(shè)計(jì)帶來(lái)的影響。

        因此,本文從應(yīng)用角度出發(fā),針對(duì)現(xiàn)代大型跨聲速民用客機(jī)開展氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)問(wèn)題的研究。在設(shè)計(jì)方法上,為了兼顧計(jì)算效率和計(jì)算復(fù)雜度之間的矛盾,采用基于全速勢(shì)方程加附面層修正的CFD數(shù)值評(píng)估方法,板殼有限元模型和松耦合靜氣動(dòng)彈性分析流程的氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)方法。針對(duì)0.85馬赫數(shù)跨聲速客機(jī)的機(jī)翼設(shè)計(jì),利用建立的優(yōu)化設(shè)計(jì)方法進(jìn)行了分別以機(jī)翼扭轉(zhuǎn)角分布和剖面翼型為設(shè)計(jì)變量的單點(diǎn)優(yōu)化設(shè)計(jì)研究。并將氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)結(jié)果與單學(xué)科的氣動(dòng)設(shè)計(jì)結(jié)果進(jìn)行對(duì)比,初步研究了考慮靜氣動(dòng)彈性影響后,對(duì)機(jī)翼氣動(dòng)載荷、翼型厚度和彎度分布,以及結(jié)構(gòu)尺寸厚度分布設(shè)計(jì)等方面帶來(lái)的影響。

        1 全速勢(shì)方程加附面層修正的氣動(dòng)數(shù)值評(píng)估方法

        隨著CFD技術(shù)在飛機(jī)氣動(dòng)外形設(shè)計(jì)中的廣泛應(yīng)用,數(shù)值評(píng)估方法的可靠性、準(zhǔn)確性及計(jì)算速度決定了設(shè)計(jì)階段的效率。研究[10-12]表明飛機(jī)在巡航狀態(tài)沒(méi)有強(qiáng)激波及大分離區(qū)存在時(shí),是可以采用全速勢(shì)方程進(jìn)行流場(chǎng)的快速評(píng)估的,而對(duì)于更為復(fù)雜的流動(dòng)現(xiàn)象,需要更高精度的數(shù)值模擬方法。本文采用的全速勢(shì)流場(chǎng)評(píng)估方法通過(guò)在邊界層理論范圍內(nèi)建立有黏-無(wú)黏耦合迭代進(jìn)行流場(chǎng)求解。首先由全速勢(shì)方程求解全流場(chǎng),得到流場(chǎng)的無(wú)黏解,再以無(wú)黏解為邊界條件,由可壓縮附面層方程計(jì)算附面層流動(dòng),得到附面層流動(dòng),以及附面層內(nèi)各流動(dòng)參數(shù)。然后再次計(jì)算全速勢(shì)方程,如此進(jìn)行有黏無(wú)黏迭代求解直到收斂。

        通過(guò)可壓縮附面層方程和有限差分格式進(jìn)行黏性附面層求解。

        通過(guò)將全速勢(shì)方程加附面層修正的數(shù)值方法應(yīng)用于標(biāo)模(Common Research Model,CRM)[13]在設(shè)計(jì)點(diǎn)(Ma=0.85,Cl=0.5,Re=1×106)的計(jì)算,并與實(shí)驗(yàn)數(shù)據(jù)和RANS方程(采用SST湍流模型,Roe格式,網(wǎng)格量800萬(wàn))計(jì)算結(jié)果進(jìn)行對(duì)比,驗(yàn)證所采用的數(shù)值方法的可靠性。CRM見圖1。

        圖1 CRM標(biāo)模剖面Cp位置Fig.1 The position of sectional Cp on CRM

        由圖2可以看出,在沿展向的6個(gè)剖面位置處,全速勢(shì)方程加附面層修正的計(jì)算結(jié)果與RANS的計(jì)算結(jié)果總體比較吻合,并與實(shí)驗(yàn)數(shù)據(jù)趨勢(shì)一致,尤其頭部峰值和下翼面的壓力分布重合度比較高,兩種數(shù)值計(jì)算方法得到的壓力恢復(fù)梯度、激波位置和激波強(qiáng)度略有差異,但是整個(gè)機(jī)翼的壓力云圖(圖3所示)形態(tài)一致,在工程誤差范圍之內(nèi)。因此,在初始機(jī)翼氣動(dòng)力設(shè)計(jì)階段,在沒(méi)有強(qiáng)激波和大分離流動(dòng)時(shí),完全可以利用全速勢(shì)方程加附面層修正的方法代替精度較高的RANS方法,作為馬赫數(shù)0.85的高馬赫數(shù)超臨界機(jī)翼的氣動(dòng)評(píng)估手段,以顯著提高設(shè)計(jì)效。

        (a)η=0.1306 (b)η=0.02828 (c)η=0.3971

        (d)η=0.5024 (e)η=0.7268 (f)η=0.95

        圖2 全速勢(shì)方程、RANS方程計(jì)算結(jié)果與實(shí)驗(yàn)剖面翼型壓力分布對(duì)比
        Fig.2 Comparison of the sectional pressure among the result of full-potential, RANS and experiment data

        圖3 全速勢(shì)方程,RANS方程計(jì)算結(jié)果表面壓力云圖對(duì)比Fig.3 Comparison of the surface pressure contour between the result of full-potential and RANS

        2 靜氣動(dòng)彈性計(jì)算方法及算例驗(yàn)證

        2.1 靜氣動(dòng)彈性分析流程

        求解氣動(dòng)彈性問(wèn)題的耦合方法通??梢苑譃榫o耦合和松耦合兩種。松耦合方法以模塊化的形式進(jìn)行氣動(dòng)彈性計(jì)算,克服了緊耦合方法面臨的計(jì)算求解困難的問(wèn)題[14]。松耦合方法氣動(dòng)力模型和結(jié)構(gòu)模型不斷耦合,交換數(shù)據(jù),當(dāng)兩者達(dá)到收斂時(shí),得到一個(gè)最終的靜氣動(dòng)彈性分析結(jié)果。本文采用C2型基函數(shù)的RBF插值方法[15-16]作為氣動(dòng)與結(jié)構(gòu)求解器之間的數(shù)據(jù)交換技術(shù),全速勢(shì)方程加附面層修正的方法作為氣動(dòng)力數(shù)值評(píng)估方法。整個(gè)靜氣動(dòng)彈性分析流程如圖4所示。

        圖4 靜氣動(dòng)彈性分析流程示意圖Fig.4 The process of static aeroelasticity analysis

        2.2 DLR F6翼身組合體驗(yàn)證算例

        選取DLR-F6翼身組合體構(gòu)型進(jìn)行靜氣動(dòng)彈性分析程序的驗(yàn)證。計(jì)算狀態(tài)為Ma=0.75,Re=3×106,定升力系數(shù)CL=0.5。機(jī)翼有限元模型采用實(shí)體模擬機(jī)翼,并挖空其中的翼盒部分,通過(guò)MSC PATRAN進(jìn)行有限元建模。材料彈性模量E=2.05×1011,泊松比υ=0.3。

        圖5和圖6分別為本文所建立的靜氣動(dòng)彈性分析方法與用Workbench商用軟件計(jì)算得到的沿展向Z方向后緣撓度變形、扭轉(zhuǎn)變形與實(shí)驗(yàn)數(shù)據(jù)[17]的對(duì)比圖??梢钥闯鰞煞N不同計(jì)算方法得到的Z向后緣撓度變形和扭轉(zhuǎn)角變形基本貼合,僅在翼梢附近具有一定明顯差別。與實(shí)驗(yàn)數(shù)據(jù)相比,兩種不同的計(jì)算方法得到的變形結(jié)果整體趨勢(shì)與實(shí)驗(yàn)數(shù)據(jù)較為吻合,但仍存在一定誤差。

        對(duì)于靜氣動(dòng)彈性,結(jié)構(gòu)變形與時(shí)間無(wú)關(guān),因此靜氣動(dòng)彈性分析取決于結(jié)構(gòu)剛度分布的準(zhǔn)確性、定常氣動(dòng)載荷計(jì)算的精度和數(shù)值交換的準(zhǔn)確性。圖7為分別采用全速勢(shì)方程加附面層修正的數(shù)值方法與RANS方程計(jì)算(采用SST湍流模型,Roe格式,網(wǎng)格量800萬(wàn))F6標(biāo)模得到的沿展向4個(gè)站位處的剖面壓力分布結(jié)果以及與實(shí)驗(yàn)數(shù)據(jù)對(duì)比圖,從圖中可以看出兩種CFD數(shù)值評(píng)估方法計(jì)算得到的壓力分布除在翼梢附近及頭部峰值外基本重合,并與實(shí)驗(yàn)數(shù)據(jù)貼合相對(duì)較好,說(shuō)明所采用的全速勢(shì)方程加附面層修正的數(shù)值評(píng)估方法能夠比較準(zhǔn)確的給出機(jī)翼的載荷分布和剖面壓力分布。由于無(wú)法得到完全符合實(shí)驗(yàn)?zāi)P偷慕Y(jié)構(gòu)有限元模型,因此兩種不同計(jì)算方法得到的計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的誤差主要來(lái)源于結(jié)構(gòu)有限元模型與實(shí)驗(yàn)?zāi)P椭g的差異。計(jì)算結(jié)果表明所建立的靜氣動(dòng)彈性分析方法是可靠的。

        圖5 計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)擾度變形對(duì)比Fig.5 Comparison of displacement between the calculation result and the experiment data

        圖6 計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)扭轉(zhuǎn)角變形對(duì)比Fig.6 Comparison of twist deformation between the calculation result and the experiment data

        (a) y/b=0.23

        (b) y/b=0.37

        (c) y/b=0.51

        (d) y/b=0.85

        3 多學(xué)科優(yōu)化設(shè)計(jì)框架

        整個(gè)系統(tǒng)優(yōu)化采用改進(jìn)的微分進(jìn)化算法。1996 年微分進(jìn)化算法參加了首屆IEEE進(jìn)化算法大賽,在所有參賽的進(jìn)化算法中,DE被證明為最優(yōu)的進(jìn)化算法,隨后該算法在各個(gè)領(lǐng)域都得到了廣泛的應(yīng)用[19]。DE算法的基本操作包括變異、交叉和選擇,標(biāo)準(zhǔn)微分進(jìn)化算法的進(jìn)本操作計(jì)算如下:

        變異算子:

        圖8 優(yōu)化框架流程圖Fig.8 Optimization framework flowchart

        交叉算子:

        選擇算子:

        通過(guò)將交叉算子中交叉對(duì)象按適應(yīng)值排序,并將次好的個(gè)體和差的個(gè)體分別與最好的個(gè)體進(jìn)行交叉,同時(shí)求取兩個(gè)交叉算子所得的平均值作為最終交叉結(jié)果對(duì)基本微分進(jìn)化算法進(jìn)行改進(jìn),在保證算法全局性的情況下,提高算法的收斂速度。改進(jìn)后的變異算子為:

        整個(gè)系統(tǒng)級(jí)的優(yōu)化以航程最大為目標(biāo),為了保證氣動(dòng)學(xué)科計(jì)算時(shí)升力系數(shù)不變,將機(jī)翼重量的減少量增加為燃油重量,如公式(8)所示。最后,整個(gè)系統(tǒng)以計(jì)算得到的航程為依據(jù),將氣動(dòng)和結(jié)構(gòu)設(shè)計(jì)結(jié)果進(jìn)行綜合評(píng)估,引導(dǎo)整個(gè)系統(tǒng)的優(yōu)化方向。

        其中,K為巡航升阻比,qe為發(fā)動(dòng)機(jī)耗油率,W0為巡航起始階段飛機(jī)重量,W1為零油重量減機(jī)翼結(jié)構(gòu)重量,W為機(jī)翼結(jié)構(gòu)重量,表達(dá)式為:

        其中,W1為直接通過(guò)有限元模型計(jì)算得到的結(jié)構(gòu)重量;c為放大系數(shù)[20],用于考慮由于有限元模型的簡(jiǎn)化而忽略掉的其它部件(不包括增升裝置)以及部件之間的鏈接構(gòu)件等的重量;W2為忽略掉的前后緣增升裝置及其滑軌等相關(guān)部件的重量,因?yàn)槲闹衅矫嫘螤钜约霸錾b置的布置形式不變,因此為了簡(jiǎn)化處理,可近似認(rèn)為這部分重量不發(fā)生變化。

        4 優(yōu)化算例

        4.1 考慮靜氣動(dòng)彈性影響的機(jī)翼扭轉(zhuǎn)角優(yōu)化設(shè)計(jì)

        選取現(xiàn)代大型寬體客機(jī)(起飛總重為218000 kg)翼身組合體進(jìn)行考慮靜氣動(dòng)彈性影響的機(jī)翼氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)。其中氣動(dòng)設(shè)計(jì)僅僅考慮扭轉(zhuǎn)角變化,從翼根(機(jī)翼安裝角不變)到翼梢選取9個(gè)控制剖面,結(jié)構(gòu)設(shè)計(jì)僅僅進(jìn)行承力結(jié)構(gòu)的尺寸優(yōu)化。結(jié)構(gòu)有限元模型只考慮了組成承力翼盒的蒙皮、前后梁腹板和翼肋,長(zhǎng)桁和梁緣條厚度采用文獻(xiàn)[21]和文獻(xiàn)[22]中的方法將其“打扁”計(jì)入蒙皮厚度,以減小建模復(fù)雜度。機(jī)翼結(jié)構(gòu)有限元模型如圖9所示,材料為鋁材。每?jī)蓚€(gè)翼肋作為一個(gè)設(shè)計(jì)變量,前后梁以及上下蒙皮被翼肋劃分為多塊,其中梁每?jī)蓧K為一個(gè)設(shè)計(jì)變量。上、下蒙皮每一塊為一個(gè)設(shè)計(jì)變量??偣灿?25個(gè)設(shè)計(jì)變量,其中氣動(dòng)設(shè)計(jì)變量9個(gè),結(jié)構(gòu)設(shè)計(jì)變量316個(gè),如表1所示。氣動(dòng)計(jì)算狀態(tài)為Ma=0.85,CL=0.515,結(jié)構(gòu)設(shè)計(jì)變量的取值范圍為,后梁厚度范圍為10~50 mm,前梁厚度范圍5~40 mm,上、下蒙皮厚度變化范圍5~50 mm,翼肋厚度變化范圍5~30 mm。優(yōu)化目標(biāo)為航程最大,結(jié)構(gòu)約束包括應(yīng)力屈服極限(考慮2.5倍的過(guò)載系數(shù))約束以及翼梢彎曲和扭轉(zhuǎn)變形約束。為了對(duì)比研究考慮靜氣動(dòng)彈性影響的多學(xué)科優(yōu)化設(shè)計(jì)對(duì)設(shè)計(jì)結(jié)果產(chǎn)生的影響,進(jìn)行了以扭轉(zhuǎn)角為設(shè)計(jì)變量的單學(xué)科的氣動(dòng)優(yōu)化設(shè)計(jì),并將初始構(gòu)型(Ori_aeroelasicity)、考慮靜氣動(dòng)彈性影響的多學(xué)科設(shè)計(jì)結(jié)果(Opt_aeroelasicity)與單學(xué)科氣動(dòng)設(shè)計(jì)優(yōu)化結(jié)果(Opt_aerodynamic)利用RANS方程對(duì)設(shè)計(jì)結(jié)果進(jìn)行校核。

        圖9 結(jié)構(gòu)有限元模型及設(shè)計(jì)變量Fig.9 The structure finite element model and design variables

        AerodynamicVariablesDescriptionNumberStructuralVariablesDescriptionNumberTwist9UpperSkin126LowSkin126Spars42Ribs22

        經(jīng)過(guò)多學(xué)科優(yōu)化,Opt_aeroelascity構(gòu)型相比于初始構(gòu)型而言航程提高了5.63%,如表2所示。其中升阻比提高了6.53%,機(jī)翼重量增加了2.56%。初始構(gòu)型靜氣動(dòng)彈性變形后翼梢撓度為1.34 m,優(yōu)化設(shè)計(jì)結(jié)果的翼梢撓度為1.61 m(圖10)。從初始構(gòu)型與優(yōu)化結(jié)果的表面壓力云圖(圖11)對(duì)比可以看出,優(yōu)化設(shè)計(jì)結(jié)果相比于初始構(gòu)型載荷適當(dāng)外移,顯著的改善了機(jī)翼的氣動(dòng)特性,同時(shí)載荷的外移使得機(jī)翼結(jié)構(gòu)重量有所增加。由于建立的有限元模型將梁緣條以及長(zhǎng)桁的厚度計(jì)入蒙皮中,因此彎矩主要由上下蒙皮承受,扭矩主要由梁腹板以及翼盒承受,因此優(yōu)化結(jié)果顯示蒙皮厚度普遍較厚。從圖12至圖13可以看出,優(yōu)化后外翼段載荷較大使得中段翼承受更大的彎矩,因此中段翼蒙皮厚度增加明顯,同時(shí)載荷的外移使得翼根彎矩增大,因此翼根處蒙皮厚度也有所增大。由于采用的有限元模型沒(méi)有考慮開口以及由起落架、發(fā)動(dòng)機(jī)引起的集中載荷的影響,同時(shí)氣動(dòng)僅僅考慮了巡航點(diǎn)的氣動(dòng)性能,因此優(yōu)化結(jié)果與實(shí)際機(jī)翼結(jié)構(gòu)存在差異。

        表2 優(yōu)化設(shè)計(jì)結(jié)果Table 2 Optimization result

        圖10 優(yōu)化結(jié)果與初始構(gòu)型機(jī)翼變形對(duì)比Fig.10 Comparison of the deformation between the optimization result and the original result

        圖11 表面壓力云圖對(duì)比Fig.11 Comparison of the surface pressure contour

        圖14為單學(xué)科氣動(dòng)優(yōu)化設(shè)計(jì)(藍(lán)色虛線)與多學(xué)科設(shè)計(jì)(綠色實(shí)線)以及初始構(gòu)型(紅色實(shí)線)的機(jī)翼升力系數(shù)和載荷分布對(duì)比圖。單學(xué)科氣動(dòng)優(yōu)化設(shè)計(jì)結(jié)果載荷分布更加貼近橢圓形載荷分布,因此具有更好的升阻特性。而考慮靜氣動(dòng)彈性影響的設(shè)計(jì)結(jié)果載荷分布介于初始構(gòu)型與單學(xué)科氣動(dòng)設(shè)計(jì)之間,為典型的三角形分布。

        圖12 優(yōu)化設(shè)計(jì)結(jié)果結(jié)構(gòu)厚度分布Fig.12 Structural thickness distribution of the optimization result

        圖13 初始構(gòu)型結(jié)構(gòu)厚度分布Fig.13 Structural thickness distribution of the original result

        (a)

        (b)

        圖15為5個(gè)典型剖面的壓力分布對(duì)比圖,圖16為扭轉(zhuǎn)角分布對(duì)比??梢钥闯?,單學(xué)科氣動(dòng)優(yōu)化設(shè)計(jì)為了獲得較高的升阻特性,減小了外翼段的負(fù)扭轉(zhuǎn),且整個(gè)扭轉(zhuǎn)角分布成波浪形。而氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)結(jié)果一方面適當(dāng)?shù)臏p小了外翼段負(fù)扭轉(zhuǎn),同時(shí)有效的控制了沿展向的扭轉(zhuǎn)角變化,避免出現(xiàn)明顯的連續(xù)波浪形,除了展向11m位置處略有起伏外,總體上分布比較光滑??梢娋C合性能最好的構(gòu)型并不是單學(xué)科氣動(dòng)特性最優(yōu)的構(gòu)型,而是氣動(dòng)特性和結(jié)構(gòu)特性之間的一種折衷。

        圖15 剖面壓力分布對(duì)比圖Fig.15 Comparison of sectional pressure distribution

        圖16 機(jī)翼展向扭轉(zhuǎn)角分布對(duì)比圖Fig.16 Comparison of wing twist distribution

        4.2 考慮靜氣動(dòng)彈性影響的機(jī)翼剖面翼型優(yōu)化設(shè)計(jì)

        選取4.1節(jié)中的初始構(gòu)型作為本節(jié)優(yōu)化算例的初始構(gòu)型,進(jìn)行考慮靜氣動(dòng)彈性影響的機(jī)翼剖面翼型優(yōu)化設(shè)計(jì)研究,其中扭轉(zhuǎn)角分布不變。如圖17所示,沿著機(jī)翼選取9個(gè)控制剖面,利用FFD方法,每個(gè)剖面有12個(gè)設(shè)計(jì)變量,結(jié)構(gòu)設(shè)計(jì)僅僅進(jìn)行承力結(jié)構(gòu)的尺寸優(yōu)化,設(shè)計(jì)變量描述如表3所示。優(yōu)化目標(biāo)為航程最大,設(shè)計(jì)約束包括油箱容積不減,應(yīng)力屈服極限(考慮2.5倍的過(guò)載系數(shù))約束以及翼梢彎曲和扭轉(zhuǎn)變形約束。為了對(duì)比研究考慮靜氣動(dòng)彈性影響之后對(duì)機(jī)翼剖面翼型以及結(jié)構(gòu)設(shè)計(jì)產(chǎn)生的影響,進(jìn)行了以剖面翼型為設(shè)計(jì)變量的單學(xué)科的氣動(dòng)優(yōu)化設(shè)計(jì)。最終,利用RANS方程對(duì)初始構(gòu)型(Ori_aeroelasicity)、考慮靜氣動(dòng)彈性影響的設(shè)計(jì)結(jié)果(Opt_aeroelasicity)以及單學(xué)科氣動(dòng)設(shè)計(jì)優(yōu)化結(jié)果(Opt_aerodynamic)進(jìn)行校核、對(duì)比分析。

        圖17 剖面翼型位置Fig.17 Position of airfoil sections

        AerodynamicVariablesDescriptionNumberStructuralVariablesDescriptionNumberFoil108UpperSkin126LowSkin126Spars42Ribs22

        如表4所示,經(jīng)過(guò)優(yōu)化航程提高了3.05%,其中升阻比提高了1.53%,機(jī)翼重量減少了3.56%。優(yōu)化后的構(gòu)型翼梢撓度為1.28 m(圖18)。從表面壓力云圖(圖19)和載荷分布(圖20)對(duì)比可以看出,考慮靜氣動(dòng)彈性影響的設(shè)計(jì)結(jié)果相比于初始構(gòu)型載荷適當(dāng)向中段翼移動(dòng),內(nèi)翼段和翼稍附近載荷相對(duì)較小。載荷分布的這種變化有利于在控制機(jī)翼翼根彎矩的情況下減小誘導(dǎo)阻力。相比之下,單學(xué)科的氣動(dòng)優(yōu)化設(shè)計(jì)結(jié)果則表現(xiàn)為載荷的明顯外移。

        表4 優(yōu)化設(shè)計(jì)結(jié)果Table 4 Optimization result

        圖18 優(yōu)化結(jié)果與初始構(gòu)型機(jī)翼變形對(duì)比Fig.18 Comparison of the deformation between the optimization result and the original result

        圖19 表面壓力云圖對(duì)比Fig.19 Comparison of the surface pressure contour

        (a)

        機(jī)翼剖面翼型的厚度、彎度、前加載和后加載對(duì)剖面當(dāng)?shù)剌d荷有著重要的影響。圖21顯示考慮靜氣動(dòng)彈性影響的設(shè)計(jì)結(jié)果相比于初始構(gòu)型,中段翼最大厚度分布明顯較小,同時(shí)圖23顯示該處上翼面彎度基本沒(méi)有減小,下翼面具有一定的前加載,而翼稍以及內(nèi)翼段呈現(xiàn)出相反的趨勢(shì),因此載荷向中段翼移動(dòng)。剖面壓力分布對(duì)比圖(圖24)顯示,通過(guò)減小下翼面的厚度以及采用前加載或者后加載,可以使下翼面產(chǎn)生更多的升力,進(jìn)而避免上翼面提供過(guò)大的升力,以此控制激波的強(qiáng)度。單學(xué)科的氣動(dòng)優(yōu)化設(shè)計(jì)結(jié)果在外翼段翼型厚度減小,出現(xiàn)了明顯的前加載和后加載,上翼面激波強(qiáng)度被顯著削弱。機(jī)翼翼盒面積分布(圖22)顯示,為了滿足油箱容積約束,無(wú)論是多學(xué)科的優(yōu)化設(shè)計(jì)還是單學(xué)科的氣動(dòng)優(yōu)化設(shè)計(jì)結(jié)果,翼根附近的翼盒面積都有所增大,以此彌補(bǔ)其它剖面翼盒面積的減小造成的油箱容積的損失,并且可以適當(dāng)減小內(nèi)翼段的載荷,使載荷外移。

        圖21 最大厚度分布對(duì)比圖Fig.21 Comparison of the distribution of maximum thickness

        從圖13和圖25可以看出,相比于初始構(gòu)型,多學(xué)科優(yōu)化結(jié)果除了kink附近蒙皮厚度有所增大外,整個(gè)外翼段蒙皮以及后梁厚度都明顯減小。這主要是由于多學(xué)科優(yōu)化設(shè)計(jì)結(jié)果的整個(gè)外翼段的翼型相對(duì)厚度明顯增大(如圖23),使得機(jī)翼翼盒承載能力增強(qiáng),因此可以適當(dāng)?shù)販p小外翼段蒙皮以及后梁厚度。在中段翼,雖然剖面翼型最大厚度減小,在一定程度上削弱了翼盒承載能力,但是后梁附近的機(jī)翼厚度并沒(méi)有明顯減小甚至有所增大,因此中段翼可以以較小的結(jié)構(gòu)重量代價(jià)來(lái)維持機(jī)翼的承載能力。

        圖22 翼盒面積分布對(duì)比圖Fig.22 Comparison of the distribution of wing-box area

        圖23 剖面翼型對(duì)比圖Fig.23 Comparison of airfoil sections

        圖24 剖面壓力分布對(duì)比圖Fig.24 Comparison of sectional pressure distribution

        圖25 優(yōu)化設(shè)計(jì)結(jié)果結(jié)構(gòu)厚度分布Fig.25 Structural thickness distribution of the optimization result

        5 結(jié) 論

        本文兼顧計(jì)算效率和計(jì)算復(fù)雜度之間的矛盾,建立了考慮靜氣動(dòng)彈性影響的氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)方法,并針對(duì)高馬赫數(shù)跨聲速客機(jī),進(jìn)行了考慮轉(zhuǎn)角分布、剖面翼型和結(jié)構(gòu)尺寸厚度影響的機(jī)翼多學(xué)科優(yōu)化設(shè)計(jì)研究,得到的結(jié)論主要如下:

        1) 通過(guò)進(jìn)行CRM標(biāo)模計(jì)算并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比,說(shuō)明基于全速式方程的數(shù)值方法可以作為初始設(shè)計(jì)階段的數(shù)值評(píng)估手段,以顯著地提高整個(gè)設(shè)計(jì)階段的設(shè)計(jì)效率。同時(shí),對(duì)F6模型進(jìn)行的靜氣動(dòng)彈性驗(yàn)證算例顯示,建立的松耦合靜氣動(dòng)彈性分析方法是可信的,可以用于考慮靜氣動(dòng)彈性影響的優(yōu)化設(shè)計(jì)研究中。

        2) 在以扭轉(zhuǎn)角為設(shè)計(jì)變量的算例中,相比于初始構(gòu)型,多學(xué)科優(yōu)化設(shè)計(jì)結(jié)果航程提高了5.63%。這主要得益于外翼段負(fù)扭轉(zhuǎn)的減小,使得載荷向外翼段移動(dòng),升阻特性得到明顯改善。相比于單學(xué)科的氣動(dòng)優(yōu)化設(shè)計(jì),考慮靜氣動(dòng)彈性影響的多學(xué)科優(yōu)化設(shè)計(jì)有效地控制了沿展向的扭轉(zhuǎn)角分布,避免載荷分布過(guò)分外移。最終,以2.56%的結(jié)構(gòu)重量損失代價(jià)獲得了6.53%的升阻比的提高。

        3) 在以剖面翼型為設(shè)計(jì)變量的算例中,相比于初始構(gòu)型,多學(xué)科優(yōu)化設(shè)計(jì)結(jié)果航程提高了3.05%。這主要得益于在保證油箱容積約束的條件下,通過(guò)調(diào)整機(jī)翼厚度分布以及前、后加載,在避免機(jī)翼激波強(qiáng)度明顯增大的情況下,將載荷向中段翼移動(dòng)。在此基礎(chǔ)上,保證后梁所處位置附近的翼型厚度不減,并適當(dāng)增大外翼段的翼盒面積,使得結(jié)構(gòu)重量減少了3.56%,同時(shí)使升阻比略有提高(提高了1.53%)。

        [1] Ma T L, Ma D L, Zhang H, et al. Aerodynamic characteristic analysis of high-aspect ratio elastic wing[J]. Journal of Beijing University of Aeronautics and Astronautics, 2007, 33(7): 781-784 (in Chinese)馬鐵林, 馬東立, 張華, 等. 大展弦比柔性機(jī)翼氣動(dòng)特性分析[J]. 北京航空航天大學(xué)學(xué)報(bào), 2007, 33(7): 781-784

        [2] Nguyen N. NASA innovation fund 2010 project: elastically shaped future air vehicle concept[R]. 2010

        [3] Cone C D. The theory of induced lift and minimum induced drag of nonplanar lifting systems[R]. National Aeronautics and Space Administration, 1962

        [4] Dong B, Zhang X D, Li Z N, et al. Integrated aerodynamic/structural design optimization for wing oftrunk liner[J]. Journal of Beijing University of Aeronautics and Astronautics, 2002, 28(4): 435-437. (in Chinese)董波, 張曉東, 酈正能, 等. 干線客機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)綜合設(shè)計(jì)研究[J]. 北京航空航天大學(xué)學(xué)報(bào), 2002, 28(4): 435-437

        [5] Hu J, Wang R H, Wang W J, et al. Multidisciplinary optimization of transport wing aerodynamic/structural integrated design[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2012, 44(4): 458-463. (in Chinese)胡婕, 王如華, 王穩(wěn)江, 等. 客機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)多學(xué)科優(yōu)化方法[J]. 南京航空航天大學(xué)學(xué)報(bào), 2012, 28(4): 458-463[6] Zhao T, Zhang Y F, Chen H X. Multi-objective aerodynamic-structural optimization of supercritical wing of wide body aircraft[C]//54thAIAA Aerospace Sciences Meeting, 2016

        [7] Zuo Y T, Lu J J, Chen G, et al. Efficient multidisciplinary aerodynamic optimization design based on discrete adjoint method[C]//54thAIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2013

        [8] Martins J, Kennedy G, Kenway G K. High aspect ratio wing design: optimal aerostructural tradeoffs for the next generation of materials[C]//Aerospace Sciences Meeting. 2014

        [9] Kenway G K W, Martins J R R A. Multi-point high-fidelity aerostructural optimization of a transport aircraft configuration[J]. Journal of Aircraft, 2014, 51(1): 144-160

        [10] Yang T H. Application of optimized design in wide-body aircraft aerodynamic design[D]. Xi’an: School of Aeronautics, Northwestern Polytechnical University, 2014. (in Chinese)楊體浩. 優(yōu)化設(shè)計(jì)在寬體客機(jī)氣動(dòng)設(shè)計(jì)中的應(yīng)用研究[D]. 西安: 西北工業(yè)大學(xué)航空學(xué)院, 2014

        [11] Chen A W, Curtin M, Carlson R B, et al. TRANAIR applications to engine/airframe integration[J]. Journal of Aircraft,

        1990, 27(8): 716-721

        [12] Zhou T, Zhang M, Li Y L, et al. Supercritical airfoil design based on full potential equations[J]. Aeronautical Computing Technique, 2009, 39(4): 58-64. (in Chinese)周濤, 張淼, 李亞林, 等. 基于全速勢(shì)方程的超臨界翼型設(shè)計(jì)[J]. 航空計(jì)算技術(shù), 2009, 39(4): 58-64

        [13] Brodersen O, Crippa S. RANS-based aerodynamic drag and pitching moment predictions for the common research model[M]. New Results in Numerical and Experimental Fluid Mechanics IX. Springer International Publishing, 2014: 485-493

        [14] Cui P, Han J L. Investigation of nonlinear aeroelastic analysis using CFD/CSD[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(3): 482-486. (in Chinese)崔鵬, 韓景龍. 基于CFD/CSD 的非線性氣動(dòng)彈性分析方法[J]. 航空學(xué)報(bào), 2010, 31(3): 482-486

        [15] Allen C B, Rendall T C S. Unified approach to CFD-CSD interpolation and mesh motion using radial basis functions[C]//25th AIAA Applied Aerodynamics Conference. Miami: AIAA 2007-3804

        [16] Yang G W. Recent progress on computational aeroelasticity[J]. Advances In Mechanics, 2009, 39(4): 406-420. (in Chinese)楊國(guó)偉. 計(jì)算氣動(dòng)彈性若干研究進(jìn)展[J]. 力學(xué)進(jìn)展, 2009, 39(4): 406-420

        [17] Burner A W, Goad W K, Massey E A. Wing deformation measurements of the DLR-F6 transport configuration in the national transonic facility[C]//AIAA Applied Aerodynamics Conference, 2009

        [18] Martins J R R A, Lambe A B. Multidisciplinary design optimization: a survey of architectures[J]. AIAA Journal, 2013, 51(9): 2049-2075

        [19] Zhao G Q. Differential evolution algorithm with greedy strategy and its applications[D]. Harbin: Department of Automatic Test and Control, Harbin Institute of Technology, 2007. (in Chinese)趙光權(quán). 基于貪婪策略的微分進(jìn)化算法及其應(yīng)用研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué)自動(dòng)化測(cè)試與控制系, 2007

        [20] Kennedy G J, Kenway G K, Martins J. High aspect ratio wing design optimal aerostructural tradeoffs for the next generation of materials[C]//Aerospace Sciences Meeting, 2014

        [21] Zhang K S, Han Z H, Li W J, et al. Multidisciplinary aerodynamic/structural design optimization for high subsonic transport wing using approximation technique[J]. Acta Aeronautica et Astronautica Sinica. 2006, 27(5): 810-815. (in Chinese)張科施, 韓忠華, 李為吉, 等. 基于近似技術(shù)的高亞音速運(yùn)輸機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)優(yōu)化設(shè)計(jì)[J]. 航空學(xué)報(bào), 2006, 27(5): 810-815

        [22] Kim Y, Kim J, Jeon Y, et al. Multidisciplinary design optimization of supersonic fighter wing using response surface methodology[C]//40th AIAA Aerospace Sciences Meeting & Exhibit, 2002.

        Aerodynamic/structural integrated design for aircraft wing with static aeroelasticity effect

        YANG Tihao, BAI Junqiang*, XIN Liang, SUN Zhiwei, SHI Yayun
        (College of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China)

        For the problem of aerodynamic/structural integrated design about a transonic aircraft, an aerodynamic/structural integrated design method is developed in consideration of the influence of static aeroelasticity. The numerical evaluation method in the integrated design is established on the basis of full potential equations and boundary layer correction. A loose coupled method based on three dimension radial basis functions (RBF) interpolation is adopted to analyze the static aeroelasticity. An improved differential evolution algorithm is chosen as the optimization algorithm for the design. The comparisons between the numerical results and the experiment data for benchmark models CRM and DLR-F6 show that the numerical technology and the loose coupled static aeroelasticity analysis method are reliable. By respectively choosing twist angle distributions and airfoil sections as the design parameter, together with the static aeroelasticity effect taken into account, the optimization method established is used in aerodynamic/structural integrated design for a transonic aircraft. Due to these two optimizations, the flying ranges are increased by 5.63% and 3.05%, respectively. The improvement in range owes to the distribution changes in the wing load and the structural thickness. In the twist distribution design case, 6.53% improvement in lift-drag ratio is obtained at the expense of 3.2% increase in the structural weight. However, in the airfoil sections design case, the lift-drag ratio is increased by 1.53%, and the structural weight is decreased by 3.56%. The design results show that the integrated optimization design method constructed in this paper is reasonable and practical.

        multidisciplinary; aerodynamics; structure; full potential equation; differential evolution algorithm; static aeroelasticity

        0258-1825(2017)04-0598-12

        2017-04-02;

        2017-06-25

        楊體浩(1989-),男,四川人,博士,研究方向:飛行器氣動(dòng)設(shè)計(jì),多學(xué)科優(yōu)化設(shè)計(jì). E-mail:xiaoyaoyangtihao@163.com

        白俊強(qiáng)(1971-),男,河南人,教授,博士,研究方向:飛行器總體及多學(xué)科優(yōu)化設(shè)計(jì)、計(jì)算流體力學(xué)、飛行力學(xué)、氣動(dòng)彈性與氣動(dòng)噪聲. E-mail:junqiang@nwpu.ecu.cn

        楊體浩, 白俊強(qiáng), 辛亮, 等. 考慮靜氣動(dòng)彈性影響的客機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)一體化設(shè)計(jì)研究[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2017, 35(4): 598-609.

        10.7638/kqdlxxb-2017.0082 YANG T H, BAI J Q, XIN L, et al. Aerodynamic/structural integrated design for aircraft wing with static aeroelasticity effect[J]. Acta Aerodynamica Sinica, 2017, 35(4): 598-609.

        V211.3

        A doi: 10.7638/kqdlxxb-2017.0082

        猜你喜歡
        氣動(dòng)彈性機(jī)翼構(gòu)型
        分子和離子立體構(gòu)型的判定
        變時(shí)滯間隙非線性機(jī)翼顫振主動(dòng)控制方法
        航天器受迫繞飛構(gòu)型設(shè)計(jì)與控制
        飛翼無(wú)人機(jī)嗡鳴氣動(dòng)彈性響應(yīng)分析
        模態(tài)選取對(duì)靜氣動(dòng)彈性分析的影響
        機(jī)翼跨聲速抖振研究進(jìn)展
        直升機(jī)的氣動(dòng)彈性問(wèn)題
        大型風(fēng)力機(jī)整機(jī)氣動(dòng)彈性響應(yīng)計(jì)算
        遙感衛(wèi)星平臺(tái)與載荷一體化構(gòu)型
        基于模糊自適應(yīng)的高超聲速機(jī)翼顫振的主動(dòng)控制
        色一情一乱一乱一区99av| 少妇精品揄拍高潮少妇桃花岛| 亚洲乱码中文字幕在线| 亚洲av中文无码乱人伦在线播放| 欧美一片二片午夜福利在线快| 国产亚洲无码1024| 日韩精品视频免费在线观看网站| 成人免费无遮挡在线播放| 亚洲欧美综合在线天堂| 欧美成人网视频| 国产一区二区三区最新地址| 男人的天堂免费a级毛片无码| 日韩精品无码久久久久久| 亚洲色偷偷偷综合网另类小说| 中文字幕人妻久久久中出| av无码av天天av天天爽| 成年人黄视频大全| 国产精品亚洲一区二区三区正片 | 不卡视频在线观看网站| 国产电影一区二区三区| 狠狠人妻久久久久久综合| 日本女优在线观看一区二区三区| 国产日产一区二区三区四区五区| 精品三级av无码一区| 人妻熟妇乱系列| 麻豆三级视频网站在线观看 | 亚洲国产av精品一区二区蜜芽| 四虎在线播放免费永久视频| 精品久久一区二区av| 日本少妇高潮喷水视频| www国产精品内射熟女| 免费a级毛片无码a∨免费| 青青草手机在线观看视频在线观看 | 久久国产精品男人的天堂av | 色婷婷激情在线一区二区三区| 水蜜桃精品视频在线观看| 国产成人aaaaa级毛片| 综合色天天久久| 国产诱惑人的视频在线观看| 亚洲av无码专区在线观看下载| 国产一级特黄无码免费视频|