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

    考慮放寬靜穩(wěn)定度的民用客機(jī)氣動(dòng)優(yōu)化設(shè)計(jì)

    2017-11-20 01:20:25李立白俊強(qiáng)郭同彪陳頌
    航空學(xué)報(bào) 2017年9期
    關(guān)鍵詞:配平穩(wěn)定度升力

    李立, 白俊強(qiáng),*, 郭同彪, 陳頌

    1.西北工業(yè)大學(xué) 航空學(xué)院, 西安 710072 2.湖北航天飛行器研究所, 武漢 430040

    考慮放寬靜穩(wěn)定度的民用客機(jī)氣動(dòng)優(yōu)化設(shè)計(jì)

    李立1, 白俊強(qiáng)1,*, 郭同彪1, 陳頌2

    1.西北工業(yè)大學(xué) 航空學(xué)院, 西安 710072 2.湖北航天飛行器研究所, 武漢 430040

    為了更加有效地減小民用客機(jī)考慮配平約束后的阻力,針對(duì)典型跨聲速民用客機(jī)機(jī)翼-機(jī)身-平尾構(gòu)型研究了不同靜穩(wěn)定度下的氣動(dòng)優(yōu)化設(shè)計(jì),并總結(jié)出在民用客機(jī)的減阻設(shè)計(jì)中考慮放寬靜穩(wěn)定度具有較大的減阻潛力。通過(guò)自由型面變形(FFD)技術(shù)對(duì)全機(jī)外形進(jìn)行參數(shù)化,實(shí)現(xiàn)機(jī)翼型面的變形,進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì)并改變平尾的偏轉(zhuǎn)保證全機(jī)能夠力矩配平。采用基于雷諾平均Navier-Stokes(RANS)方程的離散伴隨方法求解目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的梯度,然后基于序列二次規(guī)劃算法進(jìn)行基于梯度的氣動(dòng)優(yōu)化設(shè)計(jì)?;贑RM(Common Research Model)構(gòu)型,針對(duì)不同參考重心位置進(jìn)行了考慮配平約束的減阻優(yōu)化設(shè)計(jì)研究,驗(yàn)證了優(yōu)化設(shè)計(jì)系統(tǒng)的有效性,算例結(jié)果表明,隨著重心位置后移即放寬靜穩(wěn)定度,優(yōu)化構(gòu)型配平阻力減小,外翼段前緣吸力峰值明顯降低且雙激波的強(qiáng)度得到有效減弱,此外機(jī)翼的升力系數(shù)分布更加貼合最佳升力系數(shù)分布。

    放寬靜穩(wěn)定度; 離散伴隨; 氣動(dòng)優(yōu)化設(shè)計(jì); 配平阻力; 升力系數(shù)分布

    提升飛機(jī)飛行性能并有效降低油耗是目前民用飛機(jī)設(shè)計(jì)研究的重要課題,而減小全機(jī)阻力是完成這些任務(wù)的最有效手段之一。具體來(lái)說(shuō),降低全機(jī)阻力可以從以下兩個(gè)方面考慮[1]。一方面是減小配平阻力。飛機(jī)在任何飛行狀態(tài)下都處于配平狀態(tài),即繞重心的力矩為0。通常來(lái)說(shuō),民用客機(jī)都是有靜穩(wěn)定裕度的,為了保證飛機(jī)配平,后置平尾會(huì)通過(guò)負(fù)偏轉(zhuǎn)產(chǎn)生一個(gè)負(fù)升力,以提供一定的抬頭力矩使得飛機(jī)配平,此時(shí)就會(huì)導(dǎo)致機(jī)翼上有一個(gè)正升力的增量,來(lái)保證全機(jī)總升力系數(shù)不變。這一正升力的增量會(huì)使得機(jī)翼產(chǎn)生一個(gè)升致阻力增量, 連同平尾產(chǎn)生的升致阻力增量, 一起被稱為飛機(jī)的配平阻力[2]。另一方面是降低翼身組合體的激波阻力和壓差阻力等。對(duì)于前者,通常是減小翼身組合體或者平尾的型阻來(lái)實(shí)現(xiàn)配平阻力的降低,但對(duì)后平尾飛機(jī)而言,通過(guò)向后移動(dòng)全機(jī)重心也可以實(shí)現(xiàn)配平阻力的減小,這樣的設(shè)計(jì)被稱為放寬靜穩(wěn)定度設(shè)計(jì);對(duì)于后者,可以通過(guò)合理優(yōu)化機(jī)翼型面和調(diào)整機(jī)翼的載荷分布來(lái)實(shí)現(xiàn)阻力的減小。

    國(guó)內(nèi)外對(duì)考慮靜穩(wěn)定度影響的飛行器設(shè)計(jì)已經(jīng)進(jìn)行了大量的研究:Lutze提出了不同重心位置對(duì)于全機(jī)氣動(dòng)力和幾何參數(shù)產(chǎn)生影響的解析表達(dá)式,研究了如何有效減小全機(jī)配平阻力[1];Sachs推導(dǎo)出機(jī)翼下洗對(duì)機(jī)翼/平尾干擾阻力的分析表達(dá)式,并基于此闡述了最小配平阻力和全機(jī)最優(yōu)重心的關(guān)系[3];王子方用工程估算的方法分析了放寬靜穩(wěn)定對(duì)改善飛機(jī)性能的影響及其氣動(dòng)原理[4];周慧鐘和李忠應(yīng)分析了放寬靜穩(wěn)定度可能對(duì)飛行器帶來(lái)的影響以及放寬靜穩(wěn)定度后飛行器本身的動(dòng)態(tài)特性[5];許維進(jìn)和劉志敏為了研究如何降低飛機(jī)配平狀態(tài)下阻力的問(wèn)題,推導(dǎo)了飛機(jī)升致阻力與重心位置間的關(guān)系式,并得到通過(guò)向后移動(dòng)重心位置實(shí)現(xiàn)放寬靜穩(wěn)定度應(yīng)在飛機(jī)設(shè)計(jì)初期就加以考慮的結(jié)論[2];王華友等討論了飛機(jī)放寬縱向靜穩(wěn)定度后, 重心后移量對(duì)飛機(jī)縱向氣動(dòng)參數(shù)的影響, 以及由此所引起的飛機(jī)配平迎角、平尾偏角及基本飛行性能的變化[6];馮小剛等介紹了放寬靜穩(wěn)定性技術(shù)、實(shí)現(xiàn)途徑和放寬程度分析,歸納了靜穩(wěn)定度放寬后對(duì)飛機(jī)帶來(lái)的影響[7]。但是,綜合國(guó)內(nèi)外來(lái)看,針對(duì)民用寬體客機(jī),均沒(méi)有實(shí)現(xiàn)考慮放寬靜穩(wěn)定度的機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì),無(wú)法充分挖掘設(shè)計(jì)結(jié)果在考慮放寬靜穩(wěn)定度后的全機(jī)減阻收益,也無(wú)法保證考慮放寬靜穩(wěn)定度后的構(gòu)型有較好的氣動(dòng)性能。

    優(yōu)化算法是氣動(dòng)外形優(yōu)化設(shè)計(jì)的頂層架構(gòu),由于其能控制整個(gè)優(yōu)化設(shè)計(jì)的運(yùn)行和方向,因而需要合理選擇。非梯度類優(yōu)化算法雖然開發(fā)難度低、魯棒性較好、理論上可以收斂到全局最優(yōu)解,但優(yōu)化周期較長(zhǎng),計(jì)算量代價(jià)大,且對(duì)于大規(guī)模設(shè)計(jì)變量處理效率并不高效[8]?;谔荻确ǖ臍鈩?dòng)優(yōu)化設(shè)計(jì)由于可高效迭代求解快速收斂到局部最優(yōu)解而受到很大關(guān)注,而且Lyu等的研究指出梯度法求解的不同局部最優(yōu)解之間非常接近,氣動(dòng)阻力相差在0.1個(gè)阻力單位以內(nèi)[9],此外,優(yōu)化過(guò)程中對(duì)目標(biāo)函數(shù)的調(diào)用次數(shù)比較少,處理大規(guī)模設(shè)計(jì)變量時(shí)其計(jì)算量與設(shè)計(jì)變量個(gè)數(shù)之間可實(shí)現(xiàn)基本解耦,并且能夠處理大規(guī)模約束條件。因此,應(yīng)用伴隨方程法是基于梯度算法進(jìn)行大規(guī)模設(shè)計(jì)變量氣動(dòng)外形優(yōu)化設(shè)計(jì)的較好選擇[10-11]。目前國(guó)外許多航空航天研究機(jī)構(gòu)以及高校都針對(duì)基于伴隨方程的梯度法氣動(dòng)優(yōu)化設(shè)計(jì)開展了大量研究工作,而且開發(fā)形成了一批實(shí)用的軟件工具,如美國(guó)斯坦福大學(xué)飛行器設(shè)計(jì)實(shí)驗(yàn)室開發(fā)的SU2開源程序軟件[12]、法國(guó)宇航院(ONERA)開發(fā)的elsA軟件[13]、德國(guó)宇航研究院(DLR)開發(fā)的Tau軟件[14]等。國(guó)內(nèi)在該方面的研究工作起步較晚,西北工業(yè)大學(xué)[15-16]、南京航空航天大學(xué)[17-18]、中國(guó)空氣動(dòng)力研究與發(fā)展中心[19-20]等開展了相關(guān)的工作,獲得了一定的研究成果。但尚未應(yīng)用于考慮放寬靜穩(wěn)定度的機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)中。

    本文針對(duì)大型民用寬體客機(jī)全機(jī)構(gòu)型考慮放寬靜穩(wěn)定度的優(yōu)化設(shè)計(jì)問(wèn)題,以CRM(Common Research Model)標(biāo)模作為研究對(duì)象,采用基于伴隨方程法的梯度優(yōu)化算法進(jìn)行了氣動(dòng)減阻優(yōu)化設(shè)計(jì)。通過(guò)改變?nèi)珯C(jī)的參考重心位置來(lái)改變?nèi)珯C(jī)靜穩(wěn)定度,針對(duì)不同靜穩(wěn)定度的情況,通過(guò)機(jī)翼型面的變形進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì)并利用平尾的偏轉(zhuǎn)最終保證本文實(shí)現(xiàn)的全機(jī)設(shè)計(jì)結(jié)果滿足配平約束條件。算例表明了優(yōu)化設(shè)計(jì)系統(tǒng)的可行性,并對(duì)比研究了不同靜穩(wěn)定度對(duì)全機(jī)氣動(dòng)特性的影響,證明了民機(jī)氣動(dòng)減阻設(shè)計(jì)考慮放寬靜穩(wěn)定度具有較大的減阻潛力。

    1 基于FFD技術(shù)的氣動(dòng)外形參數(shù)化

    氣動(dòng)外形參數(shù)化方法是將飛行器的氣動(dòng)外形通過(guò)一定的函數(shù)關(guān)系,映射為一組設(shè)計(jì)參數(shù)(即設(shè)計(jì)變量),設(shè)計(jì)者可以通過(guò)改變這些參數(shù)或變量的值,對(duì)氣動(dòng)外形進(jìn)行操作。實(shí)用的氣動(dòng)外形參數(shù)化方法要求能夠采用較少的參數(shù)來(lái)定義幾何外形,并且保證一定精度,同時(shí)提供充分的設(shè)計(jì)空間,能夠描述盡可能豐富的外形變化。

    本文采用自由型面變形(FFD)方法對(duì)研究對(duì)象進(jìn)行參數(shù)化建模,該方法以彈性體受力后變形的思想實(shí)現(xiàn)幾何參數(shù)化過(guò)程[21]。由于FFD方法參數(shù)化對(duì)象為幾何的空間變化量,因此不需要擬合研究對(duì)象的初始構(gòu)型。且該方法能保證幾何的光滑連續(xù)性,從而能準(zhǔn)確方便地描述研究對(duì)象的整體和局部幾何信息。

    對(duì)研究對(duì)象進(jìn)行FFD方法參數(shù)化建模的過(guò)程為:首先,在待變形的研究對(duì)象周圍選取空間控制點(diǎn),連接組成FFD控制體,如圖1所示。然后,建立研究對(duì)象與控制體頂點(diǎn)之間的映射關(guān)系,此過(guò)程研究對(duì)象和控制體處于同一個(gè)空間坐標(biāo)系,稱為全局坐標(biāo)系,而控制體上建立局部坐標(biāo)系。本文選取了基于Bernstein多項(xiàng)式的FFD方法,針對(duì)研究對(duì)象上任一點(diǎn),都是通過(guò)式(1)實(shí)現(xiàn)全局坐標(biāo)與局部坐標(biāo)之間的映射關(guān)系,表達(dá)式為

    XFFD(s,v,u)=

    (1)

    式中:XFFD(s,v,u)為研究對(duì)象上任一點(diǎn)在全局坐標(biāo)系下的坐標(biāo)值,s、v、u為該點(diǎn)在控制體內(nèi)的局部坐標(biāo)值;Pi,j,k為FFD控制點(diǎn)的全局坐標(biāo);Bil(s)、Bjm(v)、Bkn(u)分別為l、m、n次Bernstein多項(xiàng)式基函數(shù)。在變形過(guò)程中,即FFD控制點(diǎn)的全局坐標(biāo)產(chǎn)生一個(gè)位移量時(shí),研究對(duì)象任一點(diǎn)的局部坐標(biāo)保持不變,通過(guò)式(1)的映射關(guān)系得到該點(diǎn)變形后的全局坐標(biāo),即可得到變形后的幾何外形。

    本文采用FFD方法是為了實(shí)現(xiàn)機(jī)翼的氣動(dòng)外形參數(shù)化和平尾的整體偏轉(zhuǎn)。對(duì)于機(jī)翼,控制機(jī)翼型面的FFD控制點(diǎn)共有380個(gè)(圖1中藍(lán)色控制點(diǎn)),其中沿展向分布10個(gè),沿弦向分布19個(gè),垂直于弦向2個(gè)。對(duì)于平尾,其整體偏轉(zhuǎn)是通過(guò)讓包圍它的FFD控制體繞給定的轉(zhuǎn)軸(圖2中紅色線)旋轉(zhuǎn)實(shí)現(xiàn)的,其中給定轉(zhuǎn)軸的定義參考實(shí)際飛機(jī)的平尾偏轉(zhuǎn),為平尾翼根和翼梢40%弦向處連接而成的直線。圖2為平尾繞給定轉(zhuǎn)軸旋轉(zhuǎn)±5° 的對(duì)比示意圖,可見(jiàn)平尾FFD控制體的整體偏轉(zhuǎn)實(shí)現(xiàn)了其包圍的平尾表面網(wǎng)格相應(yīng)的轉(zhuǎn)動(dòng),轉(zhuǎn)動(dòng)后平尾翼根與機(jī)身的網(wǎng)格仍然正常連接,也說(shuō)明本文采取的FFD方法可有效實(shí)現(xiàn)平尾的偏轉(zhuǎn)從而進(jìn)行考慮配平約束的優(yōu)化設(shè)計(jì)。

    圖1 全機(jī)構(gòu)型的FFD控制體Fig.1 FFD control framework for aircraft

    圖2 FFD控制體對(duì)平尾實(shí)現(xiàn)偏轉(zhuǎn)Fig.2 Rotated tail by FFD control frame

    2 流場(chǎng)數(shù)值求解

    本文采用雷諾平均Navier-Stokes(RANS)方程作為流場(chǎng)控制方程進(jìn)行流場(chǎng)數(shù)值求解,其守恒形式的表達(dá)式為

    (2)

    式中:U為流場(chǎng)守恒變量;F1、F2、F3為無(wú)黏通量項(xiàng);G1、G2、G3為黏性通量項(xiàng)。流場(chǎng)數(shù)值計(jì)算采用全湍計(jì)算求解,湍流模型為S-A(Spalart-Allmaras)??臻g離散格式為二階中心差分格式,時(shí)間推進(jìn)采用LU-SGS (Lower-Upper Symmetric Gauss-Seidel) 方法隱式求解。

    本文選取的研究對(duì)象是NASA第四屆和第五屆阻力預(yù)測(cè)會(huì)議的CRM構(gòu)型,是典型的跨聲速民用寬體客機(jī)構(gòu)型,其力矩參考點(diǎn)在全機(jī)平均氣動(dòng)弦長(zhǎng)的25%處,坐標(biāo)為x= 33.677 9 m,y=0 m,z=4.519 93 m,機(jī)翼半展長(zhǎng)為29.4 m,參考面積為191.84 m2[22]。本文為驗(yàn)證所采用的流場(chǎng)數(shù)值求解方法的準(zhǔn)確性,以CRM標(biāo)模作為計(jì)算對(duì)象,以CRM構(gòu)型的典型試驗(yàn)狀態(tài)[23](馬赫數(shù)Ma=0.85,升力系數(shù)CL=0.485)進(jìn)行氣動(dòng)特性對(duì)比校驗(yàn)。在驗(yàn)證求解和之后的優(yōu)化設(shè)計(jì)過(guò)程中,采用的CFD計(jì)算網(wǎng)格規(guī)模均為600萬(wàn)。計(jì)算結(jié)果與試驗(yàn)結(jié)果的典型剖面壓力系數(shù)Cp分布對(duì)比如圖3所示,c為弦長(zhǎng),可看出計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合良好,說(shuō)明本文所采用的CFD流場(chǎng)數(shù)值求解方法準(zhǔn)確可靠,可用來(lái)進(jìn)行數(shù)值優(yōu)化設(shè)計(jì)。

    圖3 試驗(yàn)與本文方法計(jì)算的壓力系數(shù)分布對(duì)比Fig.3 Comparison of pressure coefficients distribution obtained by experiment and the present method calculation

    3 離散伴隨方程法

    離散伴隨方程法直接從已經(jīng)離散的目標(biāo)函數(shù)及流場(chǎng)控制方程出發(fā),可以獲得基于離散的目標(biāo)函數(shù)和流場(chǎng)解的數(shù)值精確導(dǎo)數(shù),且邊界條件無(wú)需特殊處理,因此本文采用離散伴隨方程求解氣動(dòng)目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的梯度。

    氣動(dòng)目標(biāo)函數(shù)W通常包括力系數(shù)、力矩系數(shù)等,可通過(guò)對(duì)表面網(wǎng)格的每一個(gè)單元積分得到。構(gòu)建目標(biāo)函數(shù)表達(dá)式為

    W=W(U(τ),G(τ))

    (3)

    式中:τ為一組關(guān)于氣動(dòng)外形的設(shè)計(jì)變量;G(τ)是由設(shè)計(jì)變量決定的CFD計(jì)算網(wǎng)格。當(dāng)給定設(shè)計(jì)變量τ時(shí),通過(guò)參數(shù)化方法可確定描述氣動(dòng)外形的表面網(wǎng)格,繼而根據(jù)動(dòng)網(wǎng)格算法可以求得對(duì)應(yīng)的空間網(wǎng)格,最終通過(guò)求解可獲得流場(chǎng)解向量U(τ),因此式(3)中氣動(dòng)目標(biāo)函數(shù)W寫成關(guān)于解向量U(τ)和計(jì)算網(wǎng)格G(τ)的形式。

    當(dāng)流場(chǎng)數(shù)值求解得到穩(wěn)定收斂的流場(chǎng)解向量時(shí),式(2)可寫為在空間網(wǎng)格離散后的流動(dòng)控制方程組,即

    R=R(U(τ),G(τ))=0

    (4)

    將式(3)和式(4)分別對(duì)τ求全導(dǎo)數(shù)可得

    (5)

    (6)

    將式(6)代入式(5)有

    (7)

    如果對(duì)式(7)每一項(xiàng)單獨(dú)求解,則會(huì)涉及到矩陣求逆項(xiàng)(?R/?U)-1,由于對(duì)大規(guī)模的矩陣求逆將會(huì)帶來(lái)巨大的計(jì)算消耗,所以設(shè)伴隨變量ψ,使得

    (8)

    此時(shí)可以將矩陣求逆項(xiàng)轉(zhuǎn)化為求解線性方程組,即轉(zhuǎn)化為求解伴隨變量:

    (9)

    此時(shí),式(7)可以寫為

    (10)

    以上推導(dǎo)即利用離散伴隨方程法求解氣動(dòng)目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的梯度的過(guò)程。式(10)中,?W/?G為目標(biāo)函數(shù)對(duì)計(jì)算網(wǎng)格的偏導(dǎo)數(shù),在收斂流場(chǎng)解的基礎(chǔ)上求偏導(dǎo)即可;?R/?G為流場(chǎng)殘差對(duì)計(jì)算網(wǎng)格的偏導(dǎo)數(shù),在流場(chǎng)計(jì)算中可求得該項(xiàng);dG/dτ為計(jì)算網(wǎng)格對(duì)設(shè)計(jì)變量的偏導(dǎo)數(shù),由參數(shù)化方法和動(dòng)網(wǎng)格變形算法確定。而由于伴隨方程式(9)本質(zhì)是一個(gè)大規(guī)模、高度稀疏的線性方程組,因此式(10)準(zhǔn)確高效求解的計(jì)算量主要體現(xiàn)在對(duì)該方程組的求解,本文采用廣義最小殘量(General Minimum RESidual, GMRES)算法[24]來(lái)進(jìn)行伴隨方程的求解。

    4 優(yōu)化設(shè)計(jì)系統(tǒng)

    本文建立的優(yōu)化設(shè)計(jì)系統(tǒng)主要由以下模塊組成:基于FFD技術(shù)的氣動(dòng)外形參數(shù)化模塊、網(wǎng)格變形模塊、流場(chǎng)數(shù)值求解模塊、伴隨方程模塊和梯度優(yōu)化算法模塊。其中梯度優(yōu)化算法模塊中的梯度優(yōu)化算法選用序列二次規(guī)劃(Sequential Quadratic Programming, SQP)算法[25],該算法能迅速收斂到局部最優(yōu)解,優(yōu)化過(guò)程中對(duì)目標(biāo)函數(shù)的調(diào)用次數(shù)比較少,且能夠處理大規(guī)模函數(shù)約束條件。具體優(yōu)化流程如圖4所示。

    圖4 優(yōu)化設(shè)計(jì)系統(tǒng)流程Fig.4 Flowchart of optimization design system

    5 考慮不同靜穩(wěn)定度的全機(jī)構(gòu)型優(yōu)化設(shè)計(jì)

    本文針對(duì)CRM機(jī)翼-機(jī)身-平尾構(gòu)型開展考慮靜穩(wěn)定度涉及多約束的氣動(dòng)減阻優(yōu)化設(shè)計(jì)。通過(guò)移動(dòng)飛機(jī)的參考重心(c.g.)位置來(lái)改變飛機(jī)的靜穩(wěn)定度,本文分別選取全機(jī)平均氣動(dòng)弦長(zhǎng)的15%、20%、25%、30%和35%處作為飛機(jī)參考重心,x坐標(biāo)具體位置分別為32.977 37 m,33.327 635 m,33.677 9 m,34.028 165 m和34.378 43 m,而y和z坐標(biāo)都保持為0 m和4.519 93 m。

    針對(duì)每一個(gè)參考重心位置,進(jìn)行此靜穩(wěn)定度下的氣動(dòng)優(yōu)化設(shè)計(jì),優(yōu)化的設(shè)計(jì)點(diǎn)狀態(tài)為馬赫數(shù)Ma=0.85,升力系數(shù)CL=0.5,雷諾數(shù)Re=40.0×106,優(yōu)化目標(biāo)為阻力系數(shù)CD最小,優(yōu)化設(shè)計(jì)變量包括來(lái)流迎角α、機(jī)翼FFD控制點(diǎn)的z向位移ω(共380個(gè)設(shè)計(jì)變量)以及偏轉(zhuǎn)平尾FFD控制體的扭轉(zhuǎn)角η。

    minCD_Ma0.85

    Vary:ω,α,η

    (11)

    5.1 優(yōu)化設(shè)計(jì)結(jié)果對(duì)比

    不同參考重心位置下,優(yōu)化前后設(shè)計(jì)點(diǎn)的機(jī)翼上表面壓力系數(shù)Cp分布對(duì)比如圖6所示,從圖中可看出,初始構(gòu)型機(jī)翼上表面均有明顯的激波,通過(guò)優(yōu)化,激波得到有效抑制并基本消除。

    圖5 初始構(gòu)型和優(yōu)化構(gòu)型氣動(dòng)特性對(duì)比(Cm=0)Fig.5 Comparison of aerodynamic performance of initial and optimized configurations (Cm=0)

    圖6 初始和優(yōu)化構(gòu)型壓力系數(shù)分布對(duì)比Fig.6 Comparison of pressure coefficients distribution of initial and optimized configurations

    5.2 考慮不同靜穩(wěn)定度優(yōu)化設(shè)計(jì)減阻分析

    本小節(jié)將討論考慮不同靜穩(wěn)定度優(yōu)化設(shè)計(jì)的具體減阻收益:針對(duì)初始和優(yōu)化后考慮配平的構(gòu)型,在不同參考重心位置下,對(duì)比分析機(jī)翼典型剖面的壓力分布以及對(duì)比研究機(jī)翼和平尾升力系數(shù)分布的影響。

    圖7為初始構(gòu)型在不同參考重心位置下機(jī)翼的剖面壓力系數(shù)分布對(duì)比圖,圖8為初始構(gòu)型在不同參考重心位置下機(jī)翼和平尾升力系數(shù)分布對(duì)比圖,其中橫坐標(biāo)y表示沿機(jī)翼展向的距離。

    圖7 初始構(gòu)型機(jī)翼剖面壓力系數(shù)分布對(duì)比Fig.7 Comparison of wing section pressure coefficients distribution among initial configurations

    由于此時(shí)初始構(gòu)型只能通過(guò)平尾旋轉(zhuǎn)實(shí)現(xiàn)配平,因此,圖7中機(jī)翼表面壓力分布趨勢(shì)基本一致,而從圖5和圖8中可以得出,當(dāng)重心位置后移時(shí),平尾負(fù)升力減小,機(jī)翼正升力也減小(維持全機(jī)總升力系數(shù)),導(dǎo)致機(jī)翼和平尾的升致阻力均減小,這就是此時(shí)配平阻力有效減小的原因,故可得出對(duì)于民用客機(jī)構(gòu)型,僅通過(guò)后移其參考重心位置即放寬其固定靜穩(wěn)定度就可有效實(shí)現(xiàn)減阻。

    圖9為在不同參考重心位置下優(yōu)化構(gòu)型機(jī)翼的剖面壓力系數(shù)分布對(duì)比。圖10為在不同參考重心位置下優(yōu)化構(gòu)型機(jī)翼和平尾升力系數(shù)分布對(duì)比。圖11~圖15為在每一個(gè)參考重心位置下優(yōu)化前后構(gòu)型機(jī)翼的剖面壓力系數(shù)分布對(duì)比。通過(guò)圖11~圖15首先可以看出在每一個(gè)參考重心位置下,優(yōu)化構(gòu)型均實(shí)現(xiàn)了不同程度的激波強(qiáng)度減弱,都表現(xiàn)出內(nèi)翼段激波強(qiáng)度明顯減弱,外翼段雙激波壓力分布形態(tài)得到改善,第一個(gè)激波的強(qiáng)度明顯減弱,即相對(duì)于初始構(gòu)型,優(yōu)化構(gòu)型有效減小了激波阻力;而對(duì)于不同的參考重心位置,通過(guò)圖9 對(duì)比可以看出,隨著重心位置后移,激波減弱程度也越來(lái)越大,特別是外翼段前緣吸力峰值明顯降低和雙激波的強(qiáng)度越來(lái)越弱,即激波阻力減小越來(lái)越多;此外,從圖10中可得出(圖中“Elliptical lift coefficient”為根據(jù)機(jī)翼橢圓形環(huán)量分布得到的最佳升力系數(shù)分布)隨著重心位置后移的優(yōu)化構(gòu)型,平尾負(fù)升力減小,機(jī)翼正升力也減小(圖5中的升力系數(shù)對(duì)比可看出具體數(shù)值變化),說(shuō)明配平阻力不僅有效減小而且減小程度越來(lái)越大,此外,機(jī)翼升力系數(shù)分布更加貼合最佳升力系數(shù)分布,這是因?yàn)殡S著重心后移,優(yōu)化為了滿足配平約束首先會(huì)使得機(jī)翼載荷外移減小重心后移帶來(lái)的抬頭力矩,進(jìn)而優(yōu)化會(huì)權(quán)衡配平阻力和機(jī)翼誘導(dǎo)阻力的變化,最終的優(yōu)化結(jié)果就是機(jī)翼載荷外移,有效減小機(jī)翼的誘導(dǎo)阻力,故可得出對(duì)于民用客機(jī)構(gòu)型,通過(guò)后移其參考重心位置即放寬其固定靜穩(wěn)定度然后再對(duì)機(jī)翼進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì)可有效改善機(jī)翼表面壓力分布形態(tài)和升力系數(shù)分布,并實(shí)現(xiàn)全機(jī)較大程度減阻的結(jié)論。

    圖8 初始構(gòu)型機(jī)翼和平尾沿展向升力系數(shù)分布對(duì)比Fig.8 Comparison of lift coefficients distributions of wing and tail along span among initial configurations

    圖9 優(yōu)化構(gòu)型機(jī)翼剖面壓力系數(shù)分布對(duì)比Fig.9 Comparison of wing section pressure coefficient distributions among optimized configurations

    圖10 優(yōu)化構(gòu)型機(jī)翼和平尾沿展向升力系數(shù)分布對(duì)比Fig.10 Comparison of lift coefficients distributions of wing and tail along span among optimized configurations

    圖參考重心時(shí)初始和優(yōu)化構(gòu)型機(jī)翼剖面壓力系數(shù)對(duì)比Fig.11 Comparison of wing section pressures coefficients of initial and optimized configurations (15% c.g. position)

    圖參考重心時(shí)初始和優(yōu)化構(gòu)型機(jī)翼剖面壓力系數(shù)對(duì)比Fig.12 Comparison of wing section pressure coefficients of initial and optimized configurations (20% c.g. position)

    圖參考重心時(shí)初始和優(yōu)化構(gòu)型機(jī)翼剖面壓力系數(shù)對(duì)比Fig.13 Comparison of wing section pressure coefficients of initial and optimized configurations (25% c.g. position)

    圖參考重心時(shí)初始和優(yōu)化構(gòu)型機(jī)翼剖面壓力系數(shù)對(duì)比Fig.14 Comparison of wing section pressure coefficients of initial and optimized configurations (30% c.g. position)

    圖參考重心時(shí)初始和優(yōu)化構(gòu)型機(jī)翼剖面壓力系數(shù)對(duì)比Fig.15 Comparison of wing section pressure coefficients of initial and optimized configurations (35% c.g.position)

    6 結(jié) 論

    1) 建立了一套基于離散伴隨方法的氣動(dòng)梯度優(yōu)化設(shè)計(jì)系統(tǒng),針對(duì)CRM機(jī)翼-機(jī)身-平尾構(gòu)型進(jìn)行了不同參考重心位置下的減阻優(yōu)化設(shè)計(jì),設(shè)計(jì)結(jié)果表明配平阻力有效減小,機(jī)翼激波強(qiáng)度減弱,表面壓力分布形態(tài)改善,全機(jī)阻力降低,驗(yàn)證了優(yōu)化系統(tǒng)的可行性。

    2) 對(duì)于初始構(gòu)型,僅后移其參考重心位置即放寬固定靜穩(wěn)定度,可通過(guò)改變機(jī)翼尾翼升力分布從而減小配平阻力來(lái)有效實(shí)現(xiàn)減阻。

    3) 針對(duì)不同參考重心位置下考慮全機(jī)配平約束的優(yōu)化設(shè)計(jì),所有優(yōu)化構(gòu)型均表現(xiàn)出內(nèi)翼段激波強(qiáng)度明顯減弱,外翼段雙激波壓力分布形態(tài)得到改善,明顯減弱了第一個(gè)激波的強(qiáng)度;而隨著重心位置后移,外翼段前緣吸力峰值明顯降低以及雙激波的強(qiáng)度有效減弱,且機(jī)翼的升力系數(shù)分布更加貼合最佳升力系數(shù)分布,最終綜合實(shí)現(xiàn)全機(jī)阻力的有效減小。

    4) 在民用客機(jī)的減阻設(shè)計(jì)中,考慮放寬靜穩(wěn)定度的氣動(dòng)優(yōu)化設(shè)計(jì)具有較大的減阻潛力。

    [1] LUTZE F H. Trimmed drag considerations[J]. Journal of Aircraft, 1977, 14(6): 544-546.

    [2] 許維進(jìn), 劉志敏. 重心位置對(duì)飛機(jī)阻力及其飛行性能的影響[J]. 飛行力學(xué), 1999, 17(1): 54-58.

    XU W J, LIU Z M. Effect of center-of-gravity position on aircraft drag and flight performance[J]. Flight Dynamics, 1999, 17(1): 54-58 (in Chinese).

    [3] SACHS G. Minimum trimmed drag and optimum c.g. position[J]. Journal of Aircraft, 1978, 15(8): 456-459.

    [4] 王子方. 放寬靜穩(wěn)定度對(duì)平衡阻力影響計(jì)算方法探討[J]. 飛行力學(xué), 1987, 2(3): 18-22.

    WANG Z F. Influence on the trimmed drag calculation method of relaxing static stability[J]. Flight Dynamics, 1987, 2(3): 18-22 (in Chinese).

    [5] 周慧鐘, 李忠應(yīng). 放寬飛行器縱向靜穩(wěn)定度問(wèn)題的分析研究[J]. 飛航導(dǎo)彈, 1998(2): 8-13.

    ZHOU H Z, LI Z Y. Study on relaxing longitudinal static stability of the flying vehicle[J]. Aerodynamic Missile Journal, 1998(2): 8-13 (in Chinese).

    [6] 王華友, 李振水, 高亞奎, 等. 某機(jī)放寬縱向靜穩(wěn)定度效益分析[J]. 飛行力學(xué), 2007, 25(1): 19-21.

    WANG H Y, LI Z S, GAO Y K, et al. Effect analysis on basic flight performance for a RLSS aircraft[J]. Flight Dynamics, 2007, 25(1): 19-21 (in Chinese).

    [7] 馮小剛, 李儼, 崔永青. 民機(jī)放寬靜穩(wěn)定性的研究[J]. 電子設(shè)計(jì)工程, 2013, 21(19): 118-122.

    FENF X G, LI Y, CUI Y Q. Research on relaxed static stability of civil aircraft[J]. Electronic Design Engineering, 2013, 21(19): 118-122 (in Chinese).

    [8] 熊俊濤, 喬志德, 楊旭東, 等. 基于黏性伴隨方法的跨聲速機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)[J]. 航空學(xué)報(bào), 2007, 28(2): 281-285.

    XIONG J T, QIAO Z D, YANG X D, et al. Optimum aerodynamic design of transonic wing based on viscous adjoint method[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(2): 281-285 (in Chinese).

    [9] LYU Z J, KENWAY G K W, MARTINS J R R A. Aerodynamic shape optimization investigations of the common research model wing benchmark[J]. AIAA Journal, 2015, 53(4): 968-985.

    [10] JAMESON A. Optimum aerodynamic design using CFD and control theory: AIAA-1995-1729-CP[R]. Reston, VA: AIAA, 1995.

    [11] REUTHER J, JAMESON A, FARMER J, et al. Aerodynamic shape optimization of complex aircraft configurations via an adjoint formulation: AIAA-1996-0094[R]. Reston, VA: AIAA, 1996.

    [12] ECONOMON T D, PALACIOS F, COPELAND S R, et al. SU2: An open-source suite for multiphysics simulation and design[J]. Journal of Aircraft, 2016, 54(3): 828-846.

    [13] CARRIER G, DESTARAC D, DUMONT A, et al. Gradient-based aerodynamic optimization with the elsA software[C]//Proceedings of 52nd Aerospace Sciences Meeting—AIAA SciTech 2014. Reston, VA: AIAA, 2014.

    [14] SCHWAMBORN D, GERHOLD T, HEINRICH R. The DLR TAU-code: Recent applications in research and industry[C]//Proceedings of the European Conference on Computational Fluid Dynamics, 2006.

    [15] 陳頌, 白俊強(qiáng), 史亞云, 等. 民用客機(jī)機(jī)翼/機(jī)身/平尾構(gòu)型氣動(dòng)外形優(yōu)化設(shè)計(jì)[J]. 航空學(xué)報(bào), 2015, 36(10): 3195-3207.

    CHEN S, BAI J Q, SHI Y Y, et al. Aerodynamics shape optimization design of civil jet wing-body-tail configuration[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(10): 3195-3207 (in Chinese).

    [16] 左英桃, 傅林, 高正紅, 等. 機(jī)翼-機(jī)身-短艙-掛架外形氣動(dòng)優(yōu)化設(shè)計(jì)方法[J]. 航空動(dòng)力學(xué)報(bào), 2013, 28(9): 2009-2015.

    ZUO Y T, FU L, GAO Z H, et al. Aerodynamic optimization design of wing-body-nacelle-pylon configuration[J]. Journal of Aerospace Power, 2013, 28(9): 2009-2015 (in Chinese).

    [17] 唐智禮. 約束最優(yōu)控制理論及其在氣動(dòng)優(yōu)化中的應(yīng)用[J]. 力學(xué)學(xué)報(bào), 2007, 39(2): 273-277.

    TANG Z L. Constrained optimum control theory: Application to aerodynamic design[J]. Chinese Journal of Theoretical Applied Mechanics, 2007, 39(2): 273-277 (in Chinese).

    [18] 徐兆可, 夏健, 高宜勝. 基于三維非結(jié)構(gòu)網(wǎng)格的連續(xù)伴隨優(yōu)化方法[J]. 南京航空航天大學(xué)學(xué)報(bào), 2015, 47(1): 145-152.

    XU Z K, XIA J, GAO Y S. Continuous adjoint approach to aerodynamic optimization on 3D unstructured grids[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2015, 47(1): 145-152 (in Chinese).

    [19] 吳文華, 范召林, 陳德華, 等. 基于伴隨算子的大飛機(jī)氣動(dòng)布局精細(xì)優(yōu)化設(shè)計(jì)[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2012, 30(6): 719-724.

    WU W H, FAN Z L, CHEN D H, et al. Adjoint based on high precise aerodynamic shape optimization for transonic civil aircraft[J]. Acta Aerodynamica Sinica, 2012, 30(6): 719-724 (in Chinese).

    [20] 李彬, 鄧有奇, 唐靜, 等. 基于三維非結(jié)構(gòu)混合網(wǎng)格的離散伴隨優(yōu)化方法[J]. 航空學(xué)報(bào), 2014, 35(3): 674-686.

    LI B, DENG Y Q, TANG J, et al. Discrete adjoint optimization method for 3D unstructured grid[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(3): 674-686 (in Chinese).

    [21] SEDERBERG T W, PARRY S R. Free-form deformation of solid geometric models[J]. Computer Graphics, 1986, 20(4): 151-160.

    [22] VASSBERG J C, DEHAAN M A, RIVERS S M, et al. Development of a common research model for applied CFD validation studies: AIAA-2008-6919[R]. Reston, VA: AIAA, 2008.

    [23] LEVY D W, LAFLIN K R, TINOCO E N, et al. Summary of data from the fifth computational fluid dynamics drag prediction workshop[J]. Journal of Aircraft, 2014, 51(4): 1194-1213.

    [24] SAAD Y, SCHULTZ M H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[J]. SIAM Journal on Scientific and Statistical Computing, 1986, 7(3): 856-869.

    [25] NOCEDAL J, WRIGHT S J. Numerical optimization[M]. New York: Springer, 1999: 526-572.

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

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

    Aerodynamic optimization design for civil aircraft considering relaxed static stability

    LI Li1, BAI Junqiang1,*, GUO Tongbiao1, CHEN Song2

    1.SchoolofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China2.HubeiSpaceVehicleResearchInstitute,Wuhan430040,China

    The study on the aerodynamic optimization design of the civil jet wing-body-tail configuration with different static stability is presented, which intends to efficiently reduce the aircraft’s total drag considering trimming. It is obtained that considering static stability in drag-reduction design can have great potential of reducing the drag. The Free-Form Deform (FFD) technique is adopted to parameterize the wing shape for aerodynamic optimization design and rotate horizontal tail for trimming the pitching moment of the whole aircraft. The discrete adjoint technique based on Reynolds Averaged Navier-Stokes (RANS) equations is used to solve the gradients of targets with regard to design variables, and sequential quadratic programming is used to conduct the gradient-based optimization design. Based on the Common Research Model (CRM), optimizations considering trimming constraint with different c.g. positions are carried out to reduce aerodynamic drag, and feasibility of the optimization system is confirmed. The results of optimization cases show that when the c.g. position moves backward, the optimization configuration has smaller trim drag, obviously lower negative pressure peak of the leading edge in the outer wing, and weaker shock wave. The distribution of wing spanwise lift coefficient is improved to achieve the optimized design result.

    relaxed static stability; discrete adjoint; aerodynamic optimization design; trim drag; lift coefficient distribution

    2017-01-07; Revised: 2017-01-19; Accepted: 2017-02-16; Published online: 2017-02-22 19:21

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

    National Basic Research Program of China (2014CB744804)

    V211.41+1

    A

    1000-6893(2017)09-121112-14

    2017-01-07; 退修日期: 2017-01-19; 錄用日期: 2017-02-16; 網(wǎng)絡(luò)出版時(shí)間: 2017-02-22 19:21

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

    國(guó)家“973”計(jì)劃(2014CB744804)

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

    李立, 白俊強(qiáng), 郭同彪, 等. 考慮放寬靜穩(wěn)定度的民用客機(jī)氣動(dòng)優(yōu)化設(shè)計(jì)[J]. 航空學(xué)報(bào), 2017, 38(9): 121112. LI L, BAI J Q, GUO T B, et al. Aerodynamic optimization design for civil aircraft considering relaxed static stability[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(9): 121112.

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

    10.7527/S1000-6893.2017.121112

    猜你喜歡
    配平穩(wěn)定度升力
    高速列車車頂–升力翼組合體氣動(dòng)特性
    高穩(wěn)晶振短期頻率穩(wěn)定度的仿真分析
    無(wú)人機(jī)升力測(cè)試裝置設(shè)計(jì)及誤差因素分析
    配平化學(xué)方程式小竅門——“單質(zhì)最后配平法”
    基于自適應(yīng)偽譜法的升力式飛行器火星進(jìn)入段快速軌跡優(yōu)化
    化學(xué)方程式的配平方法
    化合價(jià)歸零法配平復(fù)雜氧化還原反應(yīng)方程式
    B737NG飛機(jī)安定面配平非典型故障分析
    多MOSFET并聯(lián)均流的高穩(wěn)定度恒流源研究
    升力式再入飛行器體襟翼姿態(tài)控制方法
    老熟妇乱子伦视频在线观看| 久久久久久久久久黄片| 欧美最新免费一区二区三区 | 热99在线观看视频| 亚洲最大成人av| 欧美成狂野欧美在线观看| 精品99又大又爽又粗少妇毛片 | 成人av一区二区三区在线看| 丰满的人妻完整版| 欧美一区二区亚洲| 久久草成人影院| 午夜日韩欧美国产| 午夜福利高清视频| 女生性感内裤真人,穿戴方法视频| 最后的刺客免费高清国语| 国产亚洲欧美98| 欧美潮喷喷水| 日韩人妻高清精品专区| 天堂√8在线中文| 男人的好看免费观看在线视频| 内地一区二区视频在线| 在线观看66精品国产| www日本黄色视频网| 成年女人永久免费观看视频| 免费在线观看亚洲国产| 欧美日本视频| 深夜a级毛片| 舔av片在线| 熟女人妻精品中文字幕| 国产精品人妻久久久久久| 日韩欧美 国产精品| 久久久成人免费电影| 99国产精品一区二区蜜桃av| 精品午夜福利视频在线观看一区| 搞女人的毛片| 一区福利在线观看| 91在线精品国自产拍蜜月| 久久久久性生活片| 欧美激情国产日韩精品一区| 久久精品国产自在天天线| 少妇裸体淫交视频免费看高清| x7x7x7水蜜桃| 国产伦一二天堂av在线观看| 人人妻,人人澡人人爽秒播| 亚洲熟妇熟女久久| 国内精品久久久久久久电影| 国产成+人综合+亚洲专区| 亚洲 国产 在线| 亚洲内射少妇av| 久久午夜亚洲精品久久| 波多野结衣巨乳人妻| 国产高潮美女av| 亚洲内射少妇av| 简卡轻食公司| 亚洲av中文字字幕乱码综合| 午夜影院日韩av| 乱人视频在线观看| 欧美成人免费av一区二区三区| 搞女人的毛片| 一进一出好大好爽视频| 久久久色成人| 免费观看精品视频网站| 嫩草影视91久久| 久久精品综合一区二区三区| 欧美丝袜亚洲另类 | 日韩人妻高清精品专区| 午夜亚洲福利在线播放| 色尼玛亚洲综合影院| 亚洲经典国产精华液单 | 精品一区二区三区av网在线观看| 精品人妻偷拍中文字幕| 国产精品女同一区二区软件 | 欧美绝顶高潮抽搐喷水| 日本 av在线| 精品免费久久久久久久清纯| 成年女人毛片免费观看观看9| 搞女人的毛片| 免费电影在线观看免费观看| 熟女人妻精品中文字幕| 免费看光身美女| 搡老岳熟女国产| 99国产精品一区二区蜜桃av| 精品一区二区三区av网在线观看| av在线观看视频网站免费| av在线天堂中文字幕| 久久99热6这里只有精品| 亚洲无线观看免费| 无遮挡黄片免费观看| 观看美女的网站| 老司机午夜福利在线观看视频| 国产精品乱码一区二三区的特点| 亚洲综合色惰| 禁无遮挡网站| 少妇的逼水好多| 99国产极品粉嫩在线观看| 噜噜噜噜噜久久久久久91| 老司机午夜十八禁免费视频| 搡老熟女国产l中国老女人| 村上凉子中文字幕在线| 久久亚洲真实| 淫妇啪啪啪对白视频| 国产免费一级a男人的天堂| 九色成人免费人妻av| 能在线免费观看的黄片| 一本综合久久免费| 国产蜜桃级精品一区二区三区| 男人的好看免费观看在线视频| 高潮久久久久久久久久久不卡| 一进一出抽搐gif免费好疼| 欧美一区二区国产精品久久精品| 欧美极品一区二区三区四区| 变态另类成人亚洲欧美熟女| 亚洲内射少妇av| 免费在线观看亚洲国产| 亚洲精品一卡2卡三卡4卡5卡| 久久久久久久久大av| 午夜免费激情av| 亚洲在线观看片| 欧美一区二区精品小视频在线| 欧美日韩中文字幕国产精品一区二区三区| 老司机午夜十八禁免费视频| 欧美最新免费一区二区三区 | 国产精品永久免费网站| 91久久精品电影网| 超碰av人人做人人爽久久| 国产精品久久久久久久久免 | 日韩欧美免费精品| 黄色丝袜av网址大全| 亚洲电影在线观看av| 国产av不卡久久| 久久久久性生活片| 又紧又爽又黄一区二区| 久久精品综合一区二区三区| 欧美高清性xxxxhd video| 国产精品亚洲一级av第二区| 90打野战视频偷拍视频| 亚洲男人的天堂狠狠| 在线观看免费视频日本深夜| 午夜福利免费观看在线| 亚洲第一电影网av| 嫩草影视91久久| а√天堂www在线а√下载| 亚洲精品乱码久久久v下载方式| 国产精品久久久久久久电影| 亚洲性夜色夜夜综合| 欧美丝袜亚洲另类 | 99久久精品热视频| 日韩亚洲欧美综合| 国产探花极品一区二区| 美女高潮喷水抽搐中文字幕| 日本黄色片子视频| 一区二区三区激情视频| 国产伦人伦偷精品视频| 久久久久久久亚洲中文字幕 | 校园春色视频在线观看| 1000部很黄的大片| 日日干狠狠操夜夜爽| 精品人妻一区二区三区麻豆 | av专区在线播放| av在线观看视频网站免费| 一个人观看的视频www高清免费观看| 国产aⅴ精品一区二区三区波| 91久久精品电影网| 亚洲精品影视一区二区三区av| 国产高清三级在线| 亚洲欧美日韩无卡精品| 99riav亚洲国产免费| 国产精品野战在线观看| 日本 av在线| 久9热在线精品视频| 在线观看舔阴道视频| 亚洲精品久久国产高清桃花| 国产精品人妻久久久久久| av欧美777| 亚洲成av人片免费观看| 亚洲最大成人中文| 亚洲熟妇中文字幕五十中出| 老司机福利观看| 国产一级毛片七仙女欲春2| a级毛片a级免费在线| 婷婷丁香在线五月| 日韩中文字幕欧美一区二区| 精品久久国产蜜桃| 在线播放国产精品三级| 自拍偷自拍亚洲精品老妇| 一个人观看的视频www高清免费观看| 夜夜躁狠狠躁天天躁| 最好的美女福利视频网| 一本一本综合久久| 一个人看的www免费观看视频| 免费av毛片视频| 简卡轻食公司| 国产精品亚洲一级av第二区| 最近中文字幕高清免费大全6 | 午夜视频国产福利| 国产免费av片在线观看野外av| 精品人妻熟女av久视频| 又紧又爽又黄一区二区| 日本a在线网址| 国产精品一区二区免费欧美| 搡老岳熟女国产| 男人和女人高潮做爰伦理| 色5月婷婷丁香| 一进一出好大好爽视频| 99久国产av精品| 亚洲第一区二区三区不卡| 国内毛片毛片毛片毛片毛片| 啦啦啦韩国在线观看视频| 深夜a级毛片| 久久久色成人| 精品不卡国产一区二区三区| 嫩草影院新地址| 国产精品一区二区免费欧美| 一区福利在线观看| 首页视频小说图片口味搜索| 人妻夜夜爽99麻豆av| 欧美成人性av电影在线观看| 亚洲在线自拍视频| 久久精品久久久久久噜噜老黄 | 看片在线看免费视频| 少妇丰满av| 精品人妻熟女av久视频| 波野结衣二区三区在线| 在线观看美女被高潮喷水网站 | 黄色女人牲交| 亚洲成a人片在线一区二区| 小说图片视频综合网站| 精品免费久久久久久久清纯| 免费av观看视频| 亚洲片人在线观看| 色视频www国产| 成人特级黄色片久久久久久久| 亚洲性夜色夜夜综合| 国产色爽女视频免费观看| 亚洲美女黄片视频| 精品熟女少妇八av免费久了| 亚洲成人免费电影在线观看| 岛国在线免费视频观看| 亚洲国产日韩欧美精品在线观看| 长腿黑丝高跟| 午夜久久久久精精品| 亚洲真实伦在线观看| 男插女下体视频免费在线播放| 好看av亚洲va欧美ⅴa在| 日本精品一区二区三区蜜桃| 婷婷丁香在线五月| 最好的美女福利视频网| 午夜激情欧美在线| 亚洲电影在线观看av| а√天堂www在线а√下载| 国产欧美日韩精品一区二区| 日韩欧美三级三区| 国产久久久一区二区三区| 国产精品久久久久久亚洲av鲁大| 成人无遮挡网站| 成人永久免费在线观看视频| 欧美一区二区国产精品久久精品| 午夜久久久久精精品| 搡老岳熟女国产| 日韩 亚洲 欧美在线| 变态另类成人亚洲欧美熟女| 亚洲国产精品成人综合色| 波多野结衣巨乳人妻| 亚洲,欧美精品.| 国产精品久久久久久久久免 | 国产精品嫩草影院av在线观看 | 悠悠久久av| 午夜福利在线观看吧| 国产精品不卡视频一区二区 | 国产伦在线观看视频一区| 男女床上黄色一级片免费看| 国产又黄又爽又无遮挡在线| 色综合亚洲欧美另类图片| xxxwww97欧美| 一进一出抽搐gif免费好疼| 日韩欧美精品免费久久 | 午夜免费激情av| 成人高潮视频无遮挡免费网站| 国产一区二区亚洲精品在线观看| 啦啦啦观看免费观看视频高清| 黄色女人牲交| 亚洲美女搞黄在线观看 | 免费在线观看影片大全网站| 亚洲成av人片免费观看| 亚洲电影在线观看av| 毛片一级片免费看久久久久 | 99视频精品全部免费 在线| 一个人免费在线观看电影| 国产伦精品一区二区三区四那| 香蕉av资源在线| 亚洲av不卡在线观看| 亚洲av成人不卡在线观看播放网| 国产精品美女特级片免费视频播放器| 国产高潮美女av| 国产成人影院久久av| 天天一区二区日本电影三级| 国内毛片毛片毛片毛片毛片| 好看av亚洲va欧美ⅴa在| 久久久久久国产a免费观看| 中文字幕人妻熟人妻熟丝袜美| 国产一区二区在线观看日韩| 久久草成人影院| 久久久久亚洲av毛片大全| 午夜精品在线福利| 一区二区三区激情视频| 禁无遮挡网站| 亚洲国产高清在线一区二区三| 成人午夜高清在线视频| 两性午夜刺激爽爽歪歪视频在线观看| 极品教师在线免费播放| 好看av亚洲va欧美ⅴa在| 久久久久久国产a免费观看| 午夜精品一区二区三区免费看| 欧美午夜高清在线| www.熟女人妻精品国产| 2021天堂中文幕一二区在线观| 亚洲欧美日韩卡通动漫| 色综合欧美亚洲国产小说| 又爽又黄无遮挡网站| 久久久久免费精品人妻一区二区| 丁香欧美五月| 中亚洲国语对白在线视频| 国内精品久久久久精免费| 久久国产乱子伦精品免费另类| 亚洲欧美清纯卡通| 亚洲,欧美,日韩| 亚洲av成人av| 老司机午夜十八禁免费视频| 欧美精品啪啪一区二区三区| 黄色女人牲交| 老鸭窝网址在线观看| 国产在线精品亚洲第一网站| 国内少妇人妻偷人精品xxx网站| 日本免费a在线| 亚洲国产精品久久男人天堂| 日本免费一区二区三区高清不卡| 亚洲精品粉嫩美女一区| 中出人妻视频一区二区| 欧美日韩综合久久久久久 | 小说图片视频综合网站| 日韩欧美三级三区| 欧美黄色淫秽网站| av天堂中文字幕网| 有码 亚洲区| 欧美色欧美亚洲另类二区| 成人高潮视频无遮挡免费网站| 亚洲国产欧美人成| 午夜日韩欧美国产| 欧美一区二区亚洲| 亚洲va日本ⅴa欧美va伊人久久| 最近视频中文字幕2019在线8| 丁香六月欧美| 欧美一区二区国产精品久久精品| 午夜日韩欧美国产| 极品教师在线免费播放| 内射极品少妇av片p| 黄色视频,在线免费观看| 免费av观看视频| 90打野战视频偷拍视频| 日日摸夜夜添夜夜添av毛片 | 热99re8久久精品国产| 国产高潮美女av| av女优亚洲男人天堂| 国内精品美女久久久久久| 国产成人影院久久av| 亚洲人成网站在线播放欧美日韩| 亚洲成av人片免费观看| 老熟妇仑乱视频hdxx| 简卡轻食公司| 国产精品久久久久久亚洲av鲁大| 日本黄色片子视频| 18禁黄网站禁片免费观看直播| 国内精品久久久久久久电影| 国产综合懂色| 国产免费av片在线观看野外av| 别揉我奶头 嗯啊视频| 午夜福利欧美成人| 久久精品国产亚洲av天美| 男女做爰动态图高潮gif福利片| 精品99又大又爽又粗少妇毛片 | 最近最新免费中文字幕在线| 夜夜看夜夜爽夜夜摸| 午夜福利视频1000在线观看| 色哟哟哟哟哟哟| 久久午夜亚洲精品久久| 午夜福利高清视频| xxxwww97欧美| av国产免费在线观看| 精品午夜福利视频在线观看一区| 草草在线视频免费看| 欧洲精品卡2卡3卡4卡5卡区| 亚洲,欧美精品.| 天堂网av新在线| 哪里可以看免费的av片| 少妇的逼水好多| 黄色视频,在线免费观看| 亚洲一区二区三区不卡视频| а√天堂www在线а√下载| 黄色一级大片看看| 一边摸一边抽搐一进一小说| 变态另类成人亚洲欧美熟女| 久久久色成人| 久久欧美精品欧美久久欧美| 亚洲精品456在线播放app | 久久久久久久精品吃奶| 一个人看的www免费观看视频| 午夜福利在线观看吧| 日韩中文字幕欧美一区二区| 欧美日韩福利视频一区二区| av天堂在线播放| 桃色一区二区三区在线观看| 美女免费视频网站| 99久久成人亚洲精品观看| 欧美午夜高清在线| 啦啦啦韩国在线观看视频| 成人特级av手机在线观看| 中文资源天堂在线| 最近中文字幕高清免费大全6 | 精品国产亚洲在线| 午夜精品久久久久久毛片777| 成人国产一区最新在线观看| 首页视频小说图片口味搜索| 搡女人真爽免费视频火全软件 | 午夜福利视频1000在线观看| 国内少妇人妻偷人精品xxx网站| 中文字幕精品亚洲无线码一区| 日本五十路高清| 又爽又黄无遮挡网站| 久久精品久久久久久噜噜老黄 | 少妇熟女aⅴ在线视频| 69av精品久久久久久| 国产精品98久久久久久宅男小说| 18美女黄网站色大片免费观看| 亚洲国产高清在线一区二区三| 亚洲国产色片| 免费看a级黄色片| 国产v大片淫在线免费观看| 男人舔女人下体高潮全视频| 尤物成人国产欧美一区二区三区| 亚洲,欧美精品.| 美女被艹到高潮喷水动态| 国产 一区 欧美 日韩| 国产精品久久久久久精品电影| 国产色爽女视频免费观看| 丝袜美腿在线中文| 夜夜看夜夜爽夜夜摸| 麻豆成人午夜福利视频| 日本与韩国留学比较| 91麻豆av在线| 久久精品国产清高在天天线| a级一级毛片免费在线观看| 欧美一区二区亚洲| 国产精品1区2区在线观看.| 俺也久久电影网| 国产精品综合久久久久久久免费| 亚洲人成伊人成综合网2020| 最新中文字幕久久久久| 天堂网av新在线| or卡值多少钱| 国产乱人伦免费视频| 午夜免费激情av| 国内少妇人妻偷人精品xxx网站| 国产69精品久久久久777片| 精品国产三级普通话版| 日韩大尺度精品在线看网址| 超碰av人人做人人爽久久| h日本视频在线播放| 久久久色成人| 久久人人精品亚洲av| 中文字幕免费在线视频6| 国产91精品成人一区二区三区| 成人欧美大片| 天美传媒精品一区二区| 亚洲中文字幕日韩| 很黄的视频免费| 亚洲精品成人久久久久久| 国产aⅴ精品一区二区三区波| 麻豆国产97在线/欧美| 好看av亚洲va欧美ⅴa在| 日本 欧美在线| 国产免费一级a男人的天堂| 99热这里只有是精品在线观看 | 啦啦啦观看免费观看视频高清| 乱人视频在线观看| av天堂中文字幕网| 一进一出抽搐动态| av国产免费在线观看| 真实男女啪啪啪动态图| 亚洲无线在线观看| 精品福利观看| 久9热在线精品视频| 午夜福利18| 日本在线视频免费播放| 国产三级黄色录像| 757午夜福利合集在线观看| 日本一二三区视频观看| 国产精品影院久久| 国内久久婷婷六月综合欲色啪| 欧美精品国产亚洲| 成人永久免费在线观看视频| 久久久久九九精品影院| 亚洲真实伦在线观看| 色尼玛亚洲综合影院| 精品久久久久久久久av| 黄色配什么色好看| 国内精品久久久久久久电影| 男人舔奶头视频| 在线免费观看不下载黄p国产 | 亚洲最大成人av| 国产亚洲精品综合一区在线观看| 免费看美女性在线毛片视频| 亚洲中文字幕日韩| 最新在线观看一区二区三区| 国产私拍福利视频在线观看| 在线观看舔阴道视频| 国模一区二区三区四区视频| 此物有八面人人有两片| 日本三级黄在线观看| 亚洲国产欧美人成| 国产白丝娇喘喷水9色精品| 简卡轻食公司| 免费电影在线观看免费观看| 亚洲中文日韩欧美视频| 一个人观看的视频www高清免费观看| 两个人的视频大全免费| 国产精品伦人一区二区| 欧美黄色片欧美黄色片| 特级一级黄色大片| 美女cb高潮喷水在线观看| 国产男靠女视频免费网站| av黄色大香蕉| 国产精品爽爽va在线观看网站| 国产精品1区2区在线观看.| 欧美xxxx性猛交bbbb| 搡女人真爽免费视频火全软件 | bbb黄色大片| 亚洲片人在线观看| 波多野结衣巨乳人妻| 欧美成人a在线观看| 又粗又爽又猛毛片免费看| 日韩欧美在线二视频| 成年女人永久免费观看视频| 黄色日韩在线| 国内精品久久久久久久电影| 91麻豆精品激情在线观看国产| 亚洲中文日韩欧美视频| 亚洲精品久久国产高清桃花| 国产三级在线视频| 国产综合懂色| 一个人免费在线观看电影| 不卡一级毛片| 人妻久久中文字幕网| 亚洲精品色激情综合| 日韩国内少妇激情av| 成年人黄色毛片网站| av国产免费在线观看| 99精品在免费线老司机午夜| 欧美激情在线99| 日韩人妻高清精品专区| 丰满人妻熟妇乱又伦精品不卡| 色综合站精品国产| 精品不卡国产一区二区三区| 成年免费大片在线观看| 日本精品一区二区三区蜜桃| 黄色配什么色好看| 国产麻豆成人av免费视频| 国产精品av视频在线免费观看| 看片在线看免费视频| 欧美成人一区二区免费高清观看| 精品人妻1区二区| 十八禁国产超污无遮挡网站| 欧美日韩乱码在线| 一区二区三区高清视频在线| 久久精品国产亚洲av天美| 亚洲狠狠婷婷综合久久图片| 亚洲人成网站在线播放欧美日韩| 欧美中文日本在线观看视频| 天堂动漫精品| 久久久久国产精品人妻aⅴ院| bbb黄色大片| 久久久久久久久大av| 欧美激情国产日韩精品一区| 午夜精品久久久久久毛片777| 国产精品av视频在线免费观看| 亚洲午夜理论影院| 嫁个100分男人电影在线观看| 一级a爱片免费观看的视频| 国产探花在线观看一区二区| 免费人成视频x8x8入口观看| 人人妻,人人澡人人爽秒播| 国产精品一区二区免费欧美| 久久6这里有精品| 97碰自拍视频| 又黄又爽又刺激的免费视频.| 亚洲欧美日韩卡通动漫| 嫩草影视91久久| 久久人妻av系列| 搡老妇女老女人老熟妇| 99久久久亚洲精品蜜臀av| 亚洲男人的天堂狠狠| 乱人视频在线观看| 日本黄大片高清| 一卡2卡三卡四卡精品乱码亚洲| 欧美日韩中文字幕国产精品一区二区三区| 99视频精品全部免费 在线| 亚洲第一欧美日韩一区二区三区| 国产亚洲欧美在线一区二区| 99视频精品全部免费 在线| 国产v大片淫在线免费观看| 欧美一区二区国产精品久久精品| 在线天堂最新版资源| 老司机福利观看| 国产精品自产拍在线观看55亚洲| 三级国产精品欧美在线观看| 免费黄网站久久成人精品 |