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

    旋翼翼型高維多目標(biāo)氣動(dòng)優(yōu)化設(shè)計(jì)

    2022-02-17 12:00:54宋超周鑄李偉斌羅驍
    關(guān)鍵詞:高維旋翼氣動(dòng)

    宋超,周鑄,李偉斌,羅驍

    (中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所,綿陽 621000)

    旋翼是直升機(jī)的核心氣動(dòng)部件,旋翼氣動(dòng)設(shè)計(jì)技術(shù)為直升機(jī)平臺(tái)性能提供堅(jiān)實(shí)的支撐[1]。旋翼翼型是旋翼氣動(dòng)設(shè)計(jì)的基礎(chǔ),直接影響直升機(jī)的飛行包線、載荷、噪聲等關(guān)鍵指標(biāo),其復(fù)雜的工作狀態(tài)決定了翼型設(shè)計(jì)是一項(xiàng)十分復(fù)雜的工作。旋翼槳葉之間存在強(qiáng)烈的氣動(dòng)干擾,且展向剖面翼型相對氣流速度受高速旋轉(zhuǎn)而差別巨大。在直升機(jī)飛行的不同工況下,前飛速度與旋轉(zhuǎn)速度疊加,使得旋翼在不同方位角的工況也劇烈變化。例如,在前飛時(shí),前行槳葉槳尖區(qū)域上存在跨聲速流動(dòng),伴隨激波的產(chǎn)生;在機(jī)動(dòng)飛行時(shí),旋翼后行槳葉上存在分離或動(dòng)態(tài)失速現(xiàn)象。因此,先進(jìn)的旋翼翼型要求兼顧高速低阻、低速高升力、低零升力矩等特性[2],其設(shè)計(jì)是典型的多設(shè)計(jì)點(diǎn)、多目標(biāo)和強(qiáng)約束優(yōu)化問題。

    基于數(shù)值優(yōu)化方法的翼型設(shè)計(jì)技術(shù)已經(jīng)取得了很大進(jìn)步,但如何滿足多設(shè)計(jì)點(diǎn)、多目標(biāo)、多約束的設(shè)計(jì)要求,仍是亟需突破的難點(diǎn)。對于多目標(biāo)優(yōu)化問題,常見的做法是采用目標(biāo)加權(quán)方法將多目標(biāo)統(tǒng)一為單目標(biāo)。例如,楊慧等[3]在開展旋翼翼型多目標(biāo)、多約束設(shè)計(jì)時(shí),將懸停狀態(tài)阻力、機(jī)動(dòng)狀態(tài)升力、阻力發(fā)散馬赫數(shù)特性等多設(shè)計(jì)目標(biāo)加權(quán)統(tǒng)一為單目標(biāo)問題。王清和招啟軍[4]提出了適用于中型運(yùn)輸直升機(jī)旋翼翼型設(shè)計(jì)的目標(biāo)及約束條件,利用遺傳算法對SC1095翼型進(jìn)行了優(yōu)化,也將多個(gè)設(shè)計(jì)目標(biāo)組合為單一目標(biāo)。

    目標(biāo)加權(quán)方法的優(yōu)化結(jié)果依賴于權(quán)因子的選取,當(dāng)優(yōu)化目標(biāo)較多時(shí),選取合適的權(quán)因子是十分困難的。采用基于Pareto非支配關(guān)系的多目標(biāo)進(jìn)化算法,如NSGA-Ⅱ,可以得到整個(gè)Pareto最優(yōu)前沿,以供設(shè)計(jì)者決策。Wang和Zhao[5]采用NSGA-Ⅱ算法進(jìn)行了考慮2個(gè)設(shè)計(jì)點(diǎn)的旋翼翼型設(shè)計(jì),旨在提高翼型升阻特性。Massaro和Benini[6]以SC1095翼型為基準(zhǔn),利用進(jìn)化算法與代理模型進(jìn)行了旋翼翼型兩目標(biāo)的設(shè)計(jì),得到了最優(yōu)Pareto前沿。然而,進(jìn)化多目標(biāo)算法解決高維多目標(biāo)設(shè)計(jì)問題(目標(biāo)數(shù)大于3)時(shí)面臨難以收斂的挑戰(zhàn)。具體而言,隨著目標(biāo)數(shù)的增多,Pareto占優(yōu)關(guān)系難以區(qū)分進(jìn)化個(gè)體,進(jìn)化動(dòng)力不足,描述Pareto前沿的解集數(shù)目也呈指數(shù)級增長,算法收斂緩慢,甚至難以收斂。為了避免高維多目標(biāo)優(yōu)化問題面臨的收斂困難問題,Zhao等[7]利用主成分分析(principal component analysis,PCA)降維方法,將6個(gè)目標(biāo)的旋翼設(shè)計(jì)問題轉(zhuǎn)化為2個(gè)目標(biāo)設(shè)計(jì)問題。然而,多目標(biāo)自身相關(guān)的嚴(yán)苛要求和降維對結(jié)果的不確定影響,使得降維方法的應(yīng)用受到了一定限制。

    在進(jìn)化算法研究領(lǐng)域,解決高維多目標(biāo)優(yōu)化算法難以收斂的問題也是熱點(diǎn)之一。Deb和Jain[8]利用聚類算子替換了NSGA-Ⅱ中的擁擠距離算子,提出了NSGA-Ⅲ算法,能夠有效應(yīng)用于高維多目標(biāo)優(yōu)化問題。Zhang和Li[9]提出了一種基于分解的多目標(biāo)優(yōu)化算法(multi-objective evolutionary algorithm based on decomposition,MOEA/D),其優(yōu)越的性能受到了學(xué)者們的廣泛關(guān)注。利用高維多目標(biāo)優(yōu)化算法,能夠充分考慮復(fù)雜氣動(dòng)優(yōu)化所涉及的多個(gè)目標(biāo),有助于提高綜合氣動(dòng)性能。然而,該類算法對復(fù)雜氣動(dòng)優(yōu)化設(shè)計(jì)問題的適應(yīng)性還未得到驗(yàn)證。

    高維多目標(biāo)優(yōu)化設(shè)計(jì)面臨的另一個(gè)挑戰(zhàn)是目標(biāo)空間的可視化問題。由于目標(biāo)空間維數(shù)的增加,Pareto前沿難以在二維或三維坐標(biāo)系中直觀展現(xiàn)。高維目標(biāo)空間的可視化也是亟待解決的問題與研究的熱點(diǎn)。高維目標(biāo)空間可視化技術(shù)可分為2類:①在平行坐標(biāo)系中表示目標(biāo)的解,如熱圖[10]和平行坐標(biāo)系[11];②構(gòu)建高維到二維的映射關(guān)系,并保留解之間的距離信息,如自組織圖映射(self-organizing mapping,SOM)[12]、徑向坐標(biāo)可視化[13]、雷達(dá)圖[14]等。SOM方法憑借其直觀的解集性質(zhì)表達(dá)形式,成為了目前流行的可視化方法。Kumano等[15]采用多目標(biāo)遺傳算法進(jìn)行了涉及4個(gè)目標(biāo)的公務(wù)機(jī)機(jī)翼氣動(dòng)/結(jié)構(gòu)多學(xué)科優(yōu)化,并利用SOM 對最優(yōu)Pareto前沿進(jìn)行了可視化分析。王超[16]利用SOM對旋翼翼型與戰(zhàn)斗機(jī)翼型最優(yōu)Pareto前沿進(jìn)行了聚類分析,為設(shè)計(jì)決策提供了全局視角。

    為實(shí)現(xiàn)旋翼高維多目標(biāo)氣動(dòng)優(yōu)化設(shè)計(jì),提高旋翼翼型綜合氣動(dòng)特性,基于MOEA/D算法建立了考慮高低速升阻特性、力矩特性、阻力發(fā)散特性等的旋翼翼型高維多目標(biāo)優(yōu)化設(shè)計(jì)方法,并采用高精度kriging模型以提高優(yōu)化設(shè)計(jì)效率。采用SOM方法對最優(yōu)解集進(jìn)行了聚類分析,并對典型最優(yōu)翼型進(jìn)行了驗(yàn)證與分析。

    1 高維多目標(biāo)優(yōu)化設(shè)計(jì)框架

    針對旋翼翼型的高維多目標(biāo)優(yōu)化設(shè)計(jì)問題,本文采用MOEA/D算法進(jìn)行問題求解。算法在優(yōu)化過程中需要計(jì)算大量的種群,對利用高精度CFD分析的優(yōu)化問題,需要耗費(fèi)巨大的計(jì)算資源。本文利用kriging模型,建立設(shè)計(jì)變量與氣動(dòng)特性之間較為精確的模型,以代替耗時(shí)的高精度CFD分析。本節(jié)簡略介紹MOEA/D算法與kriging代理模型。

    1.1 MOEA/D算法

    MOEA/D算法并不直接近似Pareto前沿,而是將多目標(biāo)優(yōu)化問題分解成一定數(shù)量的單目標(biāo)優(yōu)化子問題,然后利用進(jìn)化算法同時(shí)求解這些單目標(biāo)子問題。這些單目標(biāo)最優(yōu)解集的集合是真實(shí)Pareto前沿面的一個(gè)良好近似。每個(gè)單目標(biāo)子問題通過初始權(quán)向量之間的距離組成鄰域,子問題在對應(yīng)的鄰域內(nèi)與其他子問題進(jìn)行協(xié)同優(yōu)化。一般可以采用加權(quán)函數(shù)、Tchebycheff函數(shù)、基于懲罰的邊界交叉(penalty-based boundary intersection,PBI)等方法將多目標(biāo)轉(zhuǎn)換為單目標(biāo)求解。本文采用Tchebycheff函數(shù)聚合方法,其數(shù)學(xué)表達(dá)式為

    式中:x為決策變量;z*=為參考點(diǎn);為權(quán)重向量,共有N組權(quán)重向量。

    采用Das-Dennis方法[17]生成指定數(shù)量的相對均勻的權(quán)重向量,它們同時(shí)滿足如下條件:

    式中:H為用戶定義的正整數(shù);權(quán)重向量個(gè)數(shù)(種群個(gè)數(shù))滿足。

    下面簡單描述MOEA/D算法步驟:

    步驟1算法初始化。

    構(gòu)造權(quán)重向量λ1,λ2,…,λN,并計(jì)算任意2個(gè)權(quán)重向量之間的歐氏距離,為每個(gè)權(quán)重向量選出最近的T個(gè)向量作為它的鄰域B(i)={i1,i2,…,iT}。

    初始化種群x1,x2,…,xN,設(shè)FVi=F(xi),i=1,2,…,N。

    初始化z=(z1,z2,…,zk),zj=),其中j=1,2,…,k。

    步驟2更新。

    繁殖:從B(i)中選擇2個(gè)個(gè)體r1、r2,利用xr1、xr2交叉生成新的個(gè)體ˉy,再以一定的變異概率生成y。

    修正:在y的基礎(chǔ)上應(yīng)用特定的修正或啟發(fā)式改進(jìn)策略,得到新的個(gè)體y′;若fj(y′)<zj,則zj=fj(y′),j=1,2,…,k。

    步驟3停止判斷。若滿足停止條件,算法停止,輸出Pareto前沿,否則執(zhí)行步驟2。

    1.2 kriging模型

    kriging模型是一種基于隨機(jī)過程的估計(jì)方差最小無偏估計(jì)模型,適用于擬合具有高度非線性、多峰值問題,在氣動(dòng)優(yōu)化設(shè)計(jì)領(lǐng)域中有廣泛的應(yīng)用。隨機(jī)函數(shù)Y(x)定義為

    式中:fj(x)為基函數(shù),一般可選為簡單的多項(xiàng)式;βj為基函數(shù)對應(yīng)的系數(shù);代表Y(x)的數(shù)學(xué)期望值;Z(x)為均值為0、方差為的靜態(tài)隨機(jī)過程。對于設(shè)計(jì)空間中的任意2個(gè)位置x與x′,隨機(jī)量的協(xié)方差可表述為

    式中:相關(guān)函數(shù)R(x,x′)只與空間距離有關(guān),代表了不同位置隨機(jī)變量之間的相關(guān)性,相關(guān)函數(shù)在兩位置距離無窮大時(shí)R=0,距離為零時(shí)R=1;R隨距離的增大而減小。

    本文選用常用的高斯函數(shù):

    式中:θ為kriging超參數(shù),本文采用最大似然估計(jì)方法訓(xùn)練得到θ。

    采用交叉驗(yàn)證(leave-one-out)方法來驗(yàn)證kriging模型精度,定義為

    式中:y(xi)為第i個(gè)樣本點(diǎn)的響應(yīng)值;y-(i)(xi)為xi的預(yù)測值,特殊之處在于該預(yù)測值是由利用排除第i個(gè)樣本點(diǎn)之后的樣本集所建立模型的預(yù)測值。

    2 求解器精度驗(yàn)證

    進(jìn)行旋翼翼型氣動(dòng)優(yōu)化設(shè)計(jì)研究前,先驗(yàn)證求解器的精度。氣動(dòng)特性求解采用PMB3D[18]。求解RANS方程所采用的時(shí)間推進(jìn)方法為LU-SGS,空間離散采用JST格式。湍流模型為Spalart-Allmaras模型。計(jì)算網(wǎng)格為C型結(jié)構(gòu)網(wǎng)格,網(wǎng)格示意圖如圖1所示。

    圖1 OA309翼型計(jì)算網(wǎng)格Fig.1 Computational grids for OA309 airfoil

    為排除計(jì)算網(wǎng)格對結(jié)果的影響,開展網(wǎng)格無關(guān)性檢驗(yàn),選取3套網(wǎng)格作為檢驗(yàn)對象。表1為不同網(wǎng)格量對應(yīng)的升力系數(shù)Cl和阻力系數(shù)Cd。計(jì)算狀態(tài)為Ma=0.6,Re=2.5×106,Cl=0.50。從表中可以看出,3組不同網(wǎng)格量計(jì)算得到的升力系數(shù)幾乎一致,阻力系數(shù)隨網(wǎng)格量的增大而減小,其中網(wǎng)格2與網(wǎng)格3的阻力系數(shù)僅相差0.67%。因此,兼顧網(wǎng)格規(guī)模和計(jì)算精度考慮,選用網(wǎng)格2進(jìn)行流場分析。

    表1 OA309翼型氣動(dòng)特性隨網(wǎng)格量的變化Table 1 Aerodynamic performance of OA309 airfoil with different amounts of computational grids

    為驗(yàn)證求解器的可靠性,進(jìn)行了OA309翼型不同狀態(tài)下計(jì)算值與實(shí)驗(yàn)值的對比。實(shí)驗(yàn)值來自西北工業(yè)大學(xué)翼型葉柵空氣動(dòng)力學(xué)國防科技重點(diǎn)實(shí)驗(yàn)室NF-6風(fēng)洞,實(shí)驗(yàn)雷諾數(shù)表示為Re=Ma×4.2×106。圖2為OA309翼型表面壓力分布計(jì)算值與實(shí)驗(yàn)值對比(Ma=0.6,Re=2.5×106,Cl=0.6),計(jì)算值在壓力峰值處略小于實(shí)驗(yàn)值,兩者總體吻合良好(Cp為表面壓力系數(shù))。圖3為相同狀態(tài)下OA309翼型極曲線的計(jì)算值和實(shí)驗(yàn)值對比,可以看出兩者的阻力值吻合良好。在Cl=0.5時(shí),阻力計(jì)算值相比實(shí)驗(yàn)值的計(jì)算誤差約為6.4%。由此可以看出,PMB3D計(jì)算結(jié)果能夠很好反映翼型的氣動(dòng)特性。

    圖2 OA309翼型表面壓力系數(shù)分布計(jì)算值與實(shí)驗(yàn)值對比(Ma=0.6,Re=2.5×106,Cl=0.6)Fig.2 Comparison of pressure coefficient at surface between numerical results and experimental data for OA309 airfoil(Ma=0.6,Re=2.5×106,Cl=0.6)

    圖3 OA309翼型極曲線計(jì)算值與實(shí)驗(yàn)值對比(Ma=0.6,Re=2.5×106)Fig.3 Comparison of polars between numerical results and experimental data for OA309 airfoil(Ma=0.6,Re=2.5×106)

    3 旋翼翼型高維多目標(biāo)算例

    旋翼不同站位的翼型流場環(huán)境隨旋翼工作狀態(tài)和旋翼方位角的變化而劇烈變化,因此,旋翼翼型的設(shè)計(jì)需要充分考慮其站位的不同。一般將旋翼翼型劃分為外段翼型(90%R~100%R,R為旋翼半徑)、中段翼型(80%R~90%R)和內(nèi)段翼型(75%R~80%R)。外段翼型強(qiáng)調(diào)高速性能;內(nèi)段翼型強(qiáng)調(diào)高升力性能;中段翼型則兼顧兩方面性能。然而,旋翼翼型設(shè)計(jì)指標(biāo)之間往往相互矛盾、相互制約,需要根據(jù)直升機(jī)性能需求確定設(shè)計(jì)目標(biāo),通過仔細(xì)地權(quán)衡和折中,在上述要求中取得一個(gè)較好的平衡。

    本文充分考慮旋翼翼型的各項(xiàng)性能指標(biāo),提煉出高維多目標(biāo)設(shè)計(jì)優(yōu)化問題,避免了設(shè)計(jì)目標(biāo)難以權(quán)衡的問題。下面分別以布置在中段的9%厚度翼型和內(nèi)段的12%厚度翼型為例,進(jìn)行旋翼翼型高維多目標(biāo)優(yōu)化。

    3.1 9%厚度旋翼翼型高維多目標(biāo)算例

    以O(shè)A309翼型為初始翼型,翼型參數(shù)化采用CST(class function/shape function transformation)方法[19],翼型上下表面型函數(shù)階數(shù)均為6,設(shè)計(jì)變量數(shù)目為14。以O(shè)A309翼型的CST參數(shù)為基準(zhǔn),設(shè)計(jì)變量變化范圍為±20%。

    旋翼設(shè)計(jì)同時(shí)考慮高升力能力、懸停效率和阻力發(fā)散特性,具體的設(shè)計(jì)目標(biāo)和約束條件如表2所示。Ma=0.4時(shí)的最大升力系數(shù)是衡量旋翼翼型性能的重要參數(shù),同時(shí)需要考慮Ma=0.4時(shí)的力矩特性,對應(yīng)于表中的前2個(gè)設(shè)計(jì)目標(biāo)。Ma=0.6時(shí),機(jī)動(dòng)狀態(tài)需要較大的最大升力系數(shù),對應(yīng)表2中設(shè)計(jì)目標(biāo)3;在Ma=0.6時(shí),為懸停狀態(tài),為了提高懸停效率,需要減小翼型阻力,對應(yīng)設(shè)計(jì)目標(biāo)4。最后,考慮阻力發(fā)散特性,同時(shí)約束力矩系數(shù)的絕對值。本設(shè)計(jì)問題的設(shè)計(jì)目標(biāo)個(gè)數(shù)為5,幾何約束與力矩約束共5個(gè),涉及的計(jì)算狀態(tài)共7個(gè)。定義多目標(biāo)設(shè)計(jì)問題為

    表2 9%厚度旋翼翼型設(shè)計(jì)目標(biāo)與約束Table 2 Design objectives and constraints for a r otor airfoil with 9% thickness

    式中:對5個(gè)設(shè)計(jì)目標(biāo)進(jìn)行了歸一化處理,參考值為OA309翼型的性能;下標(biāo)數(shù)字對應(yīng)設(shè)計(jì)狀態(tài);設(shè)計(jì)狀態(tài)5、6、7對應(yīng)的馬赫數(shù)分別為0.80、0.81、0.82,Re=Ma×8×106;上標(biāo)0表示基準(zhǔn)翼型;A為翼型面積;t為翼型最大相對厚度;Cm為力矩系數(shù)。

    根據(jù)設(shè)計(jì)問題,采用進(jìn)化算法開展優(yōu)化設(shè)計(jì),需要大量評估翼型的氣動(dòng)性能,而高精度CFD流場分析耗時(shí)耗力。為提高優(yōu)化效率,采用kriging模型建立氣動(dòng)特性的響應(yīng)模型,替代優(yōu)化過程中的流場計(jì)算。首先,采用拉丁超立方試驗(yàn)設(shè)計(jì)方法(latin hypercube sampling,LHS)[20]生成700個(gè)樣本點(diǎn),在7個(gè)給定的設(shè)計(jì)狀態(tài)下進(jìn)行翼型氣動(dòng)特性分析。進(jìn)而采用交叉驗(yàn)證方法檢驗(yàn)?zāi)P途?。以設(shè)計(jì)狀態(tài)1為例,在700個(gè)樣本點(diǎn)中隨機(jī)抽取1個(gè)樣本點(diǎn)進(jìn)行交叉驗(yàn)證,重復(fù)進(jìn)行10次,得到的最大預(yù)測誤差與平均預(yù)測誤差如表3所示。從表中可以看出力矩系數(shù)的預(yù)測誤差相對較大,但最大預(yù)測誤差僅為2.596%,阻力系數(shù)預(yù)測誤差保持在0.2%以下,升力系數(shù)預(yù)測更為準(zhǔn)確。交叉驗(yàn)證結(jié)果表明,kriging模型具有較高的預(yù)測精度,可以替代高精度的CFD流場分析。

    表3 設(shè)計(jì)狀態(tài)1的kriging模型交叉驗(yàn)證結(jié)果Table 3 Results of cross validation for kriging model under design condition 1

    在建立了高精度kriging模型之后,利用kriging模型進(jìn)行5目標(biāo)的全局優(yōu)化設(shè)計(jì)。首先,采用MOEA/D算法,參數(shù)設(shè)置采用Das-Dennis方法:H=9,得到權(quán)重向量個(gè)數(shù)為715,鄰居子問題規(guī)模為20,交叉概率為0.5,最大迭代步數(shù)為200步。

    圖4給出了優(yōu)化迭代200步后的Pareto前沿。從圖中可以看出,f1與f4、f1與f5存在較為明顯的非支配關(guān)系,然而f1與f2、f3的支配關(guān)系不明顯。類似圖4的二維圖像不能很好表達(dá)設(shè)計(jì)目標(biāo)之間的關(guān)系,很難從中選擇出綜合性能良好的設(shè)計(jì)結(jié)果。

    圖4 旋翼翼型優(yōu)化最優(yōu)Pareto前沿Fig.4 Optimal Pareto front for optimization of rotor airfoil

    本文采用SOM處理最優(yōu)Pareto前沿,選取的SOM單元個(gè)數(shù)為6×6。SOM 聚類分析將每組目標(biāo)視為一個(gè)個(gè)體,根據(jù)個(gè)體之間的距離將715組5維目標(biāo)映射在36個(gè)單元上,處理后的結(jié)果如圖5所示。比較圖5(a)與圖5(d),圖5(a)中深色區(qū)域基本對應(yīng)圖5(b)中淺色區(qū)域,反映出目標(biāo)1與目標(biāo)4是相互矛盾的,這也與圖4中反映的趨勢一致,說明了SOM 聚類分析的正確性。圖5(f)為每個(gè)SOM 單元包含的個(gè)體數(shù),綜合圖5(a)~(e),折中選取圖5(f)亮黃色單元中的個(gè)體進(jìn)行分析。表4中給出了亮黃色單元中的8組目標(biāo),可以看出8組最優(yōu)目標(biāo)具有高度相似性,其中目標(biāo)1和目標(biāo)4與基本翼型基本保持一致,目標(biāo)3略有下降,目標(biāo)2和目標(biāo)5提升明顯。下面選取其中第1個(gè)個(gè)體,并利用CFD進(jìn)行氣動(dòng)特性分析。

    圖5 最優(yōu)Pareto前沿SOM可視化結(jié)果(9%厚度翼型)Fig.5 Visualization results of optimal Pareto front using SOM(airfoil with 9% thickness)

    表4 Pareto前沿中選取的8組最優(yōu)目標(biāo)Table 4 Eight groups of optimal objectives selected from Pareto front

    圖6為選取的優(yōu)化翼型形狀,表5為優(yōu)化翼型與OA309翼型的幾何信息,包括最大相對厚度與翼型面積,從表中可以看出,優(yōu)化翼型的相對厚度與無量綱面積均滿足約束條件。

    表5 OA309翼型與優(yōu)化翼型幾何信息Table 5 Geometry information of OA309 and optimized airfoils

    圖6 OA309翼型與優(yōu)化翼型形狀Fig.6 Geometry shape of OA309 and optimized airfoils

    優(yōu)化翼型與OA309翼型的升阻特性對比如圖7所示,可以看出,優(yōu)化翼型與OA309翼型在Ma=0.4時(shí)的升阻特性基本一致。圖8給出了Ma=0.4時(shí)的力矩特性??梢钥闯?,優(yōu)化翼型的零升力矩絕對值相比OA309翼型基本不變,在更大的升力系數(shù)范圍內(nèi),優(yōu)化翼型相比OA309翼型的力矩系數(shù)更接近0,力矩系數(shù)幅值(Cl≈0.8)減小約50.7%,力矩特性得到改善。在高馬赫數(shù)Ma=0.6時(shí),氣動(dòng)特性可以由圖9、圖10中看出,高馬赫數(shù)下的最大升力系數(shù)和最大升阻比均有明顯提高,其中最大升力系數(shù)提高約6.5%,最大升阻比提高約7.7%。對于力矩特性,優(yōu)化翼型在較大范圍具有更小的力矩幅值,力矩特性得到提升。

    圖7 OA309翼型與優(yōu)化翼型升阻特性對比(Ma=0.4,Re=3.2×106)Fig.7 Comparison of lift drag ratio curves between OA309 and optimized airfoils(Ma=0.4,Re=3.2×106)

    圖8 OA309翼型與優(yōu)化翼型力矩特性對比(Ma=0.4,Re=3.2×106)Fig.8 Comparison of moment coefficient curves between OA309 and optimized airfoils(Ma=0.4,Re=3.2×106)

    圖9 OA309翼型與優(yōu)化翼型升阻特性對比(Ma=0.6,Re=4.8×106)Fig.9 Comparison of lift drag ratio curves between OA309 and optimized airfoils(Ma=0.6,Re=4.8×106)

    圖10 OA309翼型與優(yōu)化翼型力矩特性對比(Ma=0.6,Re=4.8×106)Fig.10 Comparison of moment coefficient curves between OA309 and optimized airfoil(Ma=0.6,Re=4.8×106)

    圖11給出了優(yōu)化翼型與OA309翼型的阻力發(fā)散特性。從圖中可以看出,在一定速度范圍內(nèi),優(yōu)化翼型具有更低的阻力水平,在Ma=0.83之前,阻力系數(shù)增長更加緩慢。盡管之后阻力系數(shù)增長較快,但翼型很少工作在更高的馬赫數(shù)下。因此,優(yōu)化翼型具有更好的阻力發(fā)散特性。綜上所述,優(yōu)化翼型保持了低速的高升力特性,提高了高速時(shí)的阻力性能,同時(shí)提升了力矩特性。

    圖11 OA309翼型與優(yōu)化翼型阻力系數(shù)隨馬赫數(shù)變化特性對比(Cl=0,Re=Ma×8×106)Fig.11 Comparison of drag coefficient varying with Mach number between OA309 and optimized airfoils(Cl=0,Re=Ma×8×106)

    3.2 12%厚度旋翼翼型高維多目標(biāo)算例

    以O(shè)A312翼型為初始翼型,進(jìn)行旋翼內(nèi)段翼型優(yōu)化,設(shè)計(jì)狀態(tài)如表6所示。設(shè)計(jì)考慮低速時(shí)的最大升力特性、高速時(shí)的翼型阻力及阻力發(fā)散馬赫數(shù)。定義設(shè)計(jì)問題如下:

    表6 12%厚度旋翼翼型設(shè)計(jì)目標(biāo)與約束Table 6 Design objectives and constraints for a rotor airfoil with 12% thickness

    式中:下標(biāo)數(shù)字對應(yīng)設(shè)計(jì)狀態(tài);設(shè)計(jì)狀態(tài)5、6、7對應(yīng)的馬赫數(shù)分別為0.74、0.75、0.76;Re=Ma×8×106;其余符號含義同3.1節(jié)。

    翼型參數(shù)化方法、設(shè)計(jì)變量選取和優(yōu)化算法參數(shù)設(shè)置與3.1節(jié)相同。采用SOM方法進(jìn)行優(yōu)化結(jié)果后的處理。SOM 可視化結(jié)果如圖12所示。從圖中明顯看出,目標(biāo)1與目標(biāo)3、4是相互矛盾的,而目標(biāo)1與目標(biāo)5則有較好的正相關(guān)關(guān)系。OA312翼型具有較高的低速最大升力系數(shù),因此,在確保升力系數(shù)不損失的情況下,盡可能降低高速工況阻力。最終,選取的一組解在圖12(f)中標(biāo)記出來,相應(yīng)的目標(biāo)值如表7所示。選取表中第2個(gè)最優(yōu)解,利用CFD進(jìn)行詳細(xì)性能分析。

    圖12 最優(yōu)Pareto前沿SOM可視化結(jié)果(12%厚度翼型)Fig.12 Visualization results of optimal Pareto front using SOM(airfoil with 12% thickness)

    表7 Pareto前沿中選取的5組最優(yōu)目標(biāo)Table 7 Five gr oups of optimal objectives selected from Pareto front

    圖13為優(yōu)化翼型與OA312翼型的幾何形狀對比,表8中為翼型的幾何信息,包括最大相對厚度與翼型面積,從表中看出,優(yōu)化翼型的相對厚度與無量綱面積均滿足約束條件。

    表8 OA312翼型與優(yōu)化翼型幾何信息Table 8 Geometr y information of OA312 and optimized airfoils

    圖13 OA312翼型與優(yōu)化翼型形狀Fig.13 Geometry shape of OA312 and optimized airfoils

    圖14給出了翼型在Ma=0.4時(shí)的升阻特性。從圖中看出,在升力線性段,兩者無明顯差別,優(yōu)化翼型的最大升力系數(shù)略有下降。圖15給出了Ma=0.4時(shí)的翼型力矩特性,優(yōu)化翼型相比OA312翼型的力矩特性有明顯改善。在Ma=0.6時(shí),優(yōu)化翼型的升阻特性略優(yōu)于OA312翼型(見圖16)。同時(shí)力矩特性改善明顯,如圖17所示,的最大值減小約11%。

    圖14 OA312翼型與優(yōu)化翼型升阻比對比(Ma=0.4,Re=3.2×106)Fig.14 Comparison of lift-to-drag ratio curves between OA312 and optimized airfoils(Ma=0.4,Re=3.2×106)

    圖15 OA312翼型與優(yōu)化翼型力矩特性對比(Ma=0.4,Re=3.2×106)Fig.15 Comparison of moment coefficient curves between OA312 and optimized airfoils(Ma=0.4,Re=3.2×106)

    圖16 OA312翼型與優(yōu)化翼型升阻特性對比(Ma=0.6,Re=4.8×106)Fig.16 Comparison of lift-to-drag ratio curves between OA312 and optimized airfoils(Ma=0.6,Re=4.8×106)

    圖17 OA312翼型與優(yōu)化翼型力矩特性對比(Ma=0.6,Re=4.8×106)Fig.17 Comparison of moment coefficient curves between OA312 and optimized airfoils(Ma=0.6,Re=4.8×106)

    圖18給出了優(yōu)化翼型與OA312翼型的阻力發(fā)散特性。從圖中可以明顯看出,優(yōu)化翼型在較寬的馬赫數(shù)范圍內(nèi)保持著更低的阻力系數(shù),其中在Ma=0.76處減阻效果為6.1%。

    圖18 OA312翼型與優(yōu)化翼型阻力系數(shù)隨馬赫數(shù)變化特性對比(Cl=0,Re=Ma×8×106)Fig.18 Comparison of drag coefficient varying with Mach number between OA312 and optimized airfoils(Cl=0,Re=Ma×8×106)

    綜上,旋翼內(nèi)段翼型的優(yōu)化取得了良好效果,在保持低速高升力特性的同時(shí),明顯提高了高速時(shí)的阻力性能,同時(shí)提升了力矩特性。

    4 結(jié)論

    基于MOEA/D算法,本文提出了旋翼翼型高維多目標(biāo)優(yōu)化設(shè)計(jì)方法,采用kriging代理模型提高了求解效率,并應(yīng)用SOM方法對最優(yōu)Pareto解集進(jìn)行了可視化分析。旋翼翼型5目標(biāo)優(yōu)化設(shè)計(jì)的結(jié)果表明:

    1)MOEA/D算法能夠有效求解高維多目標(biāo)氣動(dòng)優(yōu)化設(shè)計(jì)問題;結(jié)合高精度代理模型,能夠高效解決旋翼翼型復(fù)雜氣動(dòng)優(yōu)化設(shè)計(jì)問題。

    2)高維多目標(biāo)優(yōu)化設(shè)計(jì)結(jié)果難以通過二維圖表直觀地展示最優(yōu)Pareto前沿。SOM 聚類分析能夠?yàn)樵O(shè)計(jì)者的決策提供較好的全局可視化結(jié)果,有助于從Pareto前沿中選取綜合性能良好的設(shè)計(jì)結(jié)果。

    3)旋翼翼型高維多目標(biāo)設(shè)計(jì)避免了設(shè)計(jì)目標(biāo)的加權(quán)或折中處理,利用高維多目標(biāo)進(jìn)化算法實(shí)現(xiàn)了翼型阻力特性、力矩特性和阻力發(fā)散特性的綜合提升,能夠提高旋翼懸停、機(jī)動(dòng)性能和操縱特性。

    猜你喜歡
    高維旋翼氣動(dòng)
    中寰氣動(dòng)執(zhí)行機(jī)構(gòu)
    改進(jìn)型自抗擾四旋翼無人機(jī)控制系統(tǒng)設(shè)計(jì)與實(shí)現(xiàn)
    基于NACA0030的波紋狀翼型氣動(dòng)特性探索
    大載重長航時(shí)油動(dòng)多旋翼無人機(jī)
    基于STM32的四旋翼飛行器的設(shè)計(jì)
    電子制作(2019年9期)2019-05-30 09:41:48
    一種改進(jìn)的GP-CLIQUE自適應(yīng)高維子空間聚類算法
    基于反饋線性化的RLV氣動(dòng)控制一體化設(shè)計(jì)
    基于加權(quán)自學(xué)習(xí)散列的高維數(shù)據(jù)最近鄰查詢算法
    四旋翼無人機(jī)動(dòng)態(tài)面控制
    一般非齊次非線性擴(kuò)散方程的等價(jià)變換和高維不變子空間
    国产老妇伦熟女老妇高清| 国产免费视频播放在线视频 | 我要看日韩黄色一级片| 综合色丁香网| 日本爱情动作片www.在线观看| 少妇的逼好多水| 免费观看性生交大片5| 一本一本综合久久| 又爽又黄a免费视频| 国产午夜精品一二区理论片| 午夜福利网站1000一区二区三区| 日韩欧美精品v在线| 国产老妇女一区| 国产一区二区在线观看日韩| 久久6这里有精品| 人妻系列 视频| 啦啦啦观看免费观看视频高清| 国产精品1区2区在线观看.| 你懂的网址亚洲精品在线观看 | 我要看日韩黄色一级片| 国语对白做爰xxxⅹ性视频网站| 女人被狂操c到高潮| 精品久久久噜噜| 成人一区二区视频在线观看| 不卡视频在线观看欧美| 国产伦一二天堂av在线观看| 国产 一区 欧美 日韩| 黄色欧美视频在线观看| av国产久精品久网站免费入址| 亚洲av成人精品一二三区| 亚洲自偷自拍三级| 国产黄片视频在线免费观看| 日本免费a在线| 免费观看的影片在线观看| a级毛色黄片| 国产精品蜜桃在线观看| 直男gayav资源| 成人欧美大片| 国产成人freesex在线| 久久久久久久久久成人| 国产精品一区二区三区四区免费观看| 日本色播在线视频| 超碰av人人做人人爽久久| 国产精品一区www在线观看| 又黄又爽又刺激的免费视频.| 黄色日韩在线| 日韩精品有码人妻一区| 18禁在线无遮挡免费观看视频| 91精品国产九色| 国产精品不卡视频一区二区| 中文字幕av成人在线电影| 五月玫瑰六月丁香| 国产麻豆成人av免费视频| 亚洲国产成人一精品久久久| 欧美3d第一页| 三级男女做爰猛烈吃奶摸视频| 高清视频免费观看一区二区 | 日本一本二区三区精品| 国产精品久久久久久久电影| 高清av免费在线| 在现免费观看毛片| 美女被艹到高潮喷水动态| 亚洲欧美精品自产自拍| 国产私拍福利视频在线观看| 国产成人a∨麻豆精品| 老司机福利观看| 两性午夜刺激爽爽歪歪视频在线观看| 国产不卡一卡二| 一个人免费在线观看电影| 女人被狂操c到高潮| 亚洲,欧美,日韩| 国产一区二区三区av在线| 成人漫画全彩无遮挡| 99热全是精品| 亚洲成av人片在线播放无| 国内少妇人妻偷人精品xxx网站| 美女国产视频在线观看| 欧美性感艳星| 日韩人妻高清精品专区| 99热这里只有精品一区| 嘟嘟电影网在线观看| 老师上课跳d突然被开到最大视频| 国产黄片美女视频| 久久久久久伊人网av| 一级毛片我不卡| 亚洲精品一区蜜桃| 国产免费福利视频在线观看| 九九热线精品视视频播放| 久久精品夜色国产| 国产伦在线观看视频一区| 亚洲精品,欧美精品| 欧美一区二区亚洲| 99热这里只有精品一区| 一个人看的www免费观看视频| 免费观看在线日韩| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 免费av观看视频| 欧美三级亚洲精品| 婷婷六月久久综合丁香| 免费看美女性在线毛片视频| 中文字幕免费在线视频6| 少妇高潮的动态图| 高清日韩中文字幕在线| 毛片女人毛片| 国产 一区 欧美 日韩| 久久精品91蜜桃| 日韩三级伦理在线观看| 午夜激情福利司机影院| 伊人久久精品亚洲午夜| 亚洲av男天堂| 乱系列少妇在线播放| a级一级毛片免费在线观看| 精品人妻偷拍中文字幕| 日韩强制内射视频| 国产v大片淫在线免费观看| 欧美激情国产日韩精品一区| 久久婷婷人人爽人人干人人爱| 丰满少妇做爰视频| 成人午夜精彩视频在线观看| 黄片无遮挡物在线观看| 一区二区三区乱码不卡18| 乱码一卡2卡4卡精品| 全区人妻精品视频| 青春草亚洲视频在线观看| 建设人人有责人人尽责人人享有的 | 午夜精品国产一区二区电影 | 久久久久久久国产电影| 久久精品91蜜桃| 高清av免费在线| 中文字幕精品亚洲无线码一区| 欧美高清成人免费视频www| 插逼视频在线观看| 亚洲在久久综合| 婷婷色综合大香蕉| 亚洲综合精品二区| 日韩欧美在线乱码| 亚洲欧美成人综合另类久久久 | 身体一侧抽搐| 狂野欧美激情性xxxx在线观看| 少妇熟女欧美另类| 一个人看的www免费观看视频| 七月丁香在线播放| 少妇人妻精品综合一区二区| 人妻制服诱惑在线中文字幕| 亚洲av男天堂| 人妻系列 视频| 亚洲欧美日韩东京热| 最新中文字幕久久久久| 国产片特级美女逼逼视频| 国产午夜精品一二区理论片| 草草在线视频免费看| 一边摸一边抽搐一进一小说| 中文字幕av在线有码专区| 日本与韩国留学比较| 日本av手机在线免费观看| 内地一区二区视频在线| 91久久精品电影网| 两个人视频免费观看高清| 又爽又黄a免费视频| 国产免费一级a男人的天堂| 淫秽高清视频在线观看| 久热久热在线精品观看| 午夜福利在线观看免费完整高清在| 国产黄a三级三级三级人| 三级国产精品片| 日韩国内少妇激情av| 国内揄拍国产精品人妻在线| 成人av在线播放网站| 国产亚洲精品av在线| 欧美精品一区二区大全| 久久婷婷人人爽人人干人人爱| 网址你懂的国产日韩在线| 小蜜桃在线观看免费完整版高清| 欧美激情久久久久久爽电影| 一区二区三区高清视频在线| 男女下面进入的视频免费午夜| 男女边吃奶边做爰视频| 性插视频无遮挡在线免费观看| 成人鲁丝片一二三区免费| 久久亚洲国产成人精品v| 男女啪啪激烈高潮av片| 亚洲精品影视一区二区三区av| 男女那种视频在线观看| 国产精品无大码| 一区二区三区免费毛片| 精品久久久久久久人妻蜜臀av| 深夜a级毛片| 一区二区三区乱码不卡18| 秋霞伦理黄片| 草草在线视频免费看| 日韩亚洲欧美综合| 国产亚洲精品av在线| 成人国产麻豆网| 久久久久网色| 亚洲成色77777| 国产成人精品婷婷| 免费看日本二区| 精品不卡国产一区二区三区| 少妇裸体淫交视频免费看高清| 亚洲国产精品国产精品| 日产精品乱码卡一卡2卡三| 亚洲国产精品成人综合色| 又黄又爽又刺激的免费视频.| 久久欧美精品欧美久久欧美| 春色校园在线视频观看| 少妇熟女aⅴ在线视频| 久久久久久伊人网av| 国产精品久久久久久精品电影小说 | 高清在线视频一区二区三区 | 亚洲精品国产av成人精品| 久久久久久久久中文| 国内揄拍国产精品人妻在线| 三级国产精品片| 色网站视频免费| 色综合站精品国产| 乱码一卡2卡4卡精品| 午夜福利在线在线| 午夜福利视频1000在线观看| 国产成人freesex在线| 国产精品人妻久久久影院| 三级国产精品片| 人体艺术视频欧美日本| 中国国产av一级| 99久久人妻综合| 91午夜精品亚洲一区二区三区| 日韩av在线大香蕉| kizo精华| 亚洲av成人av| 精品人妻一区二区三区麻豆| www日本黄色视频网| 亚洲av免费高清在线观看| 麻豆成人av视频| 日本wwww免费看| 久99久视频精品免费| 亚洲国产精品合色在线| 男女国产视频网站| 国产精品麻豆人妻色哟哟久久 | 一级毛片久久久久久久久女| 青春草视频在线免费观看| 日本午夜av视频| 国产片特级美女逼逼视频| 岛国毛片在线播放| 深爱激情五月婷婷| 夜夜爽夜夜爽视频| 哪个播放器可以免费观看大片| 亚洲最大成人手机在线| 伊人久久精品亚洲午夜| 激情 狠狠 欧美| 尤物成人国产欧美一区二区三区| 国产精品国产三级专区第一集| 一边摸一边抽搐一进一小说| www.av在线官网国产| 男女边吃奶边做爰视频| 国产日韩欧美在线精品| 日本wwww免费看| 麻豆精品久久久久久蜜桃| 美女脱内裤让男人舔精品视频| 一级毛片久久久久久久久女| 看黄色毛片网站| 免费电影在线观看免费观看| 女人被狂操c到高潮| 国产女主播在线喷水免费视频网站 | 嫩草影院入口| 久久久久久久久久成人| 噜噜噜噜噜久久久久久91| 精品人妻一区二区三区麻豆| 久久久久精品久久久久真实原创| 欧美丝袜亚洲另类| 日本黄色视频三级网站网址| av卡一久久| 18+在线观看网站| 97超视频在线观看视频| 精品一区二区三区视频在线| 免费无遮挡裸体视频| 亚洲欧洲日产国产| 亚洲人成网站高清观看| 国产精品久久久久久精品电影| 色视频www国产| 国产女主播在线喷水免费视频网站 | 一个人看视频在线观看www免费| 国产一区二区在线观看日韩| 久久韩国三级中文字幕| 在线播放国产精品三级| 国产熟女欧美一区二区| 久久精品综合一区二区三区| 免费看日本二区| 国产成人精品久久久久久| 亚洲自偷自拍三级| 韩国高清视频一区二区三区| 中文字幕av在线有码专区| 日本与韩国留学比较| 免费看美女性在线毛片视频| 亚洲久久久久久中文字幕| 精品不卡国产一区二区三区| 日本午夜av视频| 国产亚洲最大av| 大香蕉97超碰在线| 高清在线视频一区二区三区 | 久久国产乱子免费精品| 最近视频中文字幕2019在线8| 简卡轻食公司| 国产探花在线观看一区二区| 国产成人精品久久久久久| 免费不卡的大黄色大毛片视频在线观看 | 日本wwww免费看| 午夜免费激情av| videos熟女内射| 精品国产露脸久久av麻豆 | 亚洲第一区二区三区不卡| 国产成年人精品一区二区| 国产老妇伦熟女老妇高清| 欧美高清成人免费视频www| 91aial.com中文字幕在线观看| 欧美zozozo另类| 国产精品永久免费网站| 永久网站在线| 国产黄a三级三级三级人| 韩国av在线不卡| 国产精品.久久久| 亚洲精品乱久久久久久| 国产成人aa在线观看| 久久久久久久久久成人| 久久99热这里只有精品18| 最近手机中文字幕大全| 国产激情偷乱视频一区二区| 亚洲精品亚洲一区二区| 亚洲欧美精品专区久久| 亚洲精品色激情综合| 美女xxoo啪啪120秒动态图| 国产黄色视频一区二区在线观看 | 成人三级黄色视频| 69人妻影院| 国产色爽女视频免费观看| 国产午夜福利久久久久久| 欧美一区二区国产精品久久精品| 中文字幕制服av| 亚洲精品亚洲一区二区| av视频在线观看入口| 成人美女网站在线观看视频| 欧美又色又爽又黄视频| 中文字幕熟女人妻在线| 欧美另类亚洲清纯唯美| 国产精品乱码一区二三区的特点| 日韩精品有码人妻一区| 丝袜美腿在线中文| 国产v大片淫在线免费观看| 在线观看66精品国产| 亚洲av免费在线观看| 亚洲人与动物交配视频| 老女人水多毛片| 亚洲中文字幕日韩| 有码 亚洲区| 日韩 亚洲 欧美在线| 成人午夜精彩视频在线观看| 蜜臀久久99精品久久宅男| 校园人妻丝袜中文字幕| 国产高清国产精品国产三级 | 最新中文字幕久久久久| 国产在视频线精品| 国产伦精品一区二区三区四那| 亚洲乱码一区二区免费版| 一边摸一边抽搐一进一小说| 直男gayav资源| 国产高潮美女av| 日韩强制内射视频| 国产美女午夜福利| 97超碰精品成人国产| 汤姆久久久久久久影院中文字幕 | 日韩av不卡免费在线播放| 一区二区三区四区激情视频| 99热这里只有是精品50| 国产午夜精品久久久久久一区二区三区| 日韩欧美精品免费久久| 国产又色又爽无遮挡免| 26uuu在线亚洲综合色| 韩国高清视频一区二区三区| 美女大奶头视频| 婷婷色av中文字幕| 亚洲国产精品sss在线观看| 国产69精品久久久久777片| 只有这里有精品99| 亚洲一区高清亚洲精品| av免费观看日本| 嘟嘟电影网在线观看| 成人三级黄色视频| 久久久久久久国产电影| 精品一区二区免费观看| av在线老鸭窝| 少妇的逼好多水| 美女被艹到高潮喷水动态| 精品久久久久久久久av| 国产激情偷乱视频一区二区| 舔av片在线| 日韩一区二区视频免费看| 日韩精品有码人妻一区| 非洲黑人性xxxx精品又粗又长| 国产精品一区二区三区四区免费观看| 亚洲av.av天堂| 国产黄色视频一区二区在线观看 | 日本免费a在线| 精品国内亚洲2022精品成人| 热99在线观看视频| 99久久九九国产精品国产免费| 亚洲国产精品sss在线观看| 乱系列少妇在线播放| 麻豆成人av视频| 国产精品av视频在线免费观看| 国产亚洲av嫩草精品影院| 国产激情偷乱视频一区二区| 国产白丝娇喘喷水9色精品| 国产私拍福利视频在线观看| 简卡轻食公司| 欧美人与善性xxx| 婷婷六月久久综合丁香| 1024手机看黄色片| 亚洲在久久综合| 99热网站在线观看| 久久久国产成人精品二区| 午夜福利网站1000一区二区三区| 男人舔女人下体高潮全视频| 久久精品久久久久久噜噜老黄 | 亚洲天堂国产精品一区在线| 欧美日本亚洲视频在线播放| 久久韩国三级中文字幕| 日韩一区二区视频免费看| 色综合亚洲欧美另类图片| 久久午夜福利片| 伦精品一区二区三区| 亚洲精品456在线播放app| 日本熟妇午夜| 国产一区亚洲一区在线观看| 色视频www国产| 中文亚洲av片在线观看爽| 国产精品爽爽va在线观看网站| 成年av动漫网址| 两个人的视频大全免费| 日本爱情动作片www.在线观看| 真实男女啪啪啪动态图| 免费看a级黄色片| 一级毛片我不卡| 免费看光身美女| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 久久人人爽人人片av| 少妇熟女aⅴ在线视频| 97在线视频观看| 精品人妻一区二区三区麻豆| 欧美日韩在线观看h| 亚洲欧美一区二区三区国产| 成人毛片60女人毛片免费| 国产亚洲精品久久久com| 国产精品综合久久久久久久免费| 麻豆精品久久久久久蜜桃| 97人妻精品一区二区三区麻豆| 欧美极品一区二区三区四区| 七月丁香在线播放| 国产欧美日韩精品一区二区| 亚洲av一区综合| 看十八女毛片水多多多| 久久久色成人| 亚洲真实伦在线观看| 晚上一个人看的免费电影| 亚洲欧洲国产日韩| 十八禁国产超污无遮挡网站| 亚洲国产欧美人成| 伊人久久精品亚洲午夜| 久久久久九九精品影院| 又粗又爽又猛毛片免费看| 久久精品久久久久久久性| 精品国产三级普通话版| 日韩视频在线欧美| 免费大片18禁| 亚洲成人av在线免费| 日本欧美国产在线视频| 亚洲真实伦在线观看| 久久久久久久国产电影| av又黄又爽大尺度在线免费看 | 高清日韩中文字幕在线| 看片在线看免费视频| 日本免费a在线| 老司机影院毛片| 在线免费十八禁| 亚洲自拍偷在线| 久久久久久久久久久丰满| 国产亚洲91精品色在线| 欧美日本视频| 99久久无色码亚洲精品果冻| 麻豆久久精品国产亚洲av| 建设人人有责人人尽责人人享有的 | 九九热线精品视视频播放| 日韩国内少妇激情av| 一边摸一边抽搐一进一小说| 日韩av不卡免费在线播放| 一级av片app| 久久精品国产自在天天线| 国产高清三级在线| 欧美性猛交╳xxx乱大交人| 国产黄片美女视频| 国产精品一区www在线观看| 精品人妻视频免费看| 精品人妻一区二区三区麻豆| 成人特级av手机在线观看| 黄色配什么色好看| 国产亚洲精品久久久com| 成人毛片a级毛片在线播放| 男人的好看免费观看在线视频| 插逼视频在线观看| 久久精品国产亚洲网站| 国内揄拍国产精品人妻在线| 中文字幕久久专区| 日本免费a在线| av天堂中文字幕网| 午夜老司机福利剧场| 人人妻人人澡人人爽人人夜夜 | 99久久无色码亚洲精品果冻| 色综合色国产| 午夜精品在线福利| av播播在线观看一区| 别揉我奶头 嗯啊视频| 精品久久久久久久久av| 国内精品一区二区在线观看| 好男人在线观看高清免费视频| 亚洲自拍偷在线| 国产综合懂色| 舔av片在线| 夫妻性生交免费视频一级片| 少妇人妻精品综合一区二区| 日本wwww免费看| 亚洲av中文字字幕乱码综合| 久久6这里有精品| 身体一侧抽搐| 国产精品久久久久久久久免| 五月玫瑰六月丁香| 一个人看视频在线观看www免费| 日韩欧美精品免费久久| 日产精品乱码卡一卡2卡三| 亚洲精品一区蜜桃| 日韩强制内射视频| 成人亚洲欧美一区二区av| 精品人妻熟女av久视频| 久久精品夜夜夜夜夜久久蜜豆| 免费人成在线观看视频色| 两个人视频免费观看高清| 亚洲第一区二区三区不卡| 国产69精品久久久久777片| 免费电影在线观看免费观看| 色5月婷婷丁香| 插阴视频在线观看视频| 免费看a级黄色片| 国产大屁股一区二区在线视频| 乱人视频在线观看| 亚洲中文字幕日韩| 超碰97精品在线观看| 最近最新中文字幕免费大全7| 中国美白少妇内射xxxbb| 插阴视频在线观看视频| 村上凉子中文字幕在线| 亚洲人成网站高清观看| 日本爱情动作片www.在线观看| 寂寞人妻少妇视频99o| 欧美成人精品欧美一级黄| 国产精品人妻久久久影院| 两性午夜刺激爽爽歪歪视频在线观看| 男人和女人高潮做爰伦理| 国产一区二区三区av在线| 亚洲av成人精品一区久久| 日本av手机在线免费观看| 国产精品一区二区性色av| 中文字幕亚洲精品专区| 国产精品日韩av在线免费观看| 看十八女毛片水多多多| 非洲黑人性xxxx精品又粗又长| www.色视频.com| 日韩成人伦理影院| 2021少妇久久久久久久久久久| 美女大奶头视频| 精品国内亚洲2022精品成人| 卡戴珊不雅视频在线播放| 亚洲激情五月婷婷啪啪| 99久久成人亚洲精品观看| 国产高清国产精品国产三级 | 中文字幕制服av| 国产精品久久久久久精品电影| 日本一本二区三区精品| 一边亲一边摸免费视频| 国产又色又爽无遮挡免| 嫩草影院新地址| 熟女电影av网| 蜜臀久久99精品久久宅男| 超碰97精品在线观看| 亚洲欧美日韩无卡精品| 国产亚洲av片在线观看秒播厂 | 国内少妇人妻偷人精品xxx网站| 亚洲欧美中文字幕日韩二区| 亚洲乱码一区二区免费版| 别揉我奶头 嗯啊视频| 亚洲内射少妇av| 国产伦一二天堂av在线观看| 一夜夜www| 日本-黄色视频高清免费观看| 51国产日韩欧美| 最新中文字幕久久久久| 成人综合一区亚洲| 舔av片在线| 欧美不卡视频在线免费观看| 久久午夜福利片| 国内少妇人妻偷人精品xxx网站| 如何舔出高潮| 国产高清视频在线观看网站| 中国美白少妇内射xxxbb| 亚洲婷婷狠狠爱综合网| 亚洲国产高清在线一区二区三| 国产精品嫩草影院av在线观看| 国产乱人视频| 免费看日本二区| 亚洲av免费在线观看|