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

    基于伴隨方程和自由變形技術(shù)的跨聲速機(jī)翼氣動(dòng)設(shè)計(jì)方法研究

    2014-05-04 09:40:32白俊強(qiáng)陳頌華俊孫智偉黃江濤
    關(guān)鍵詞:機(jī)翼外形氣動(dòng)

    白俊強(qiáng),陳頌,華俊,孫智偉,黃江濤

    (1.西北工業(yè)大學(xué) 航空學(xué)院,陜西 西安 710072;2.中國(guó)航空研究院,北京 100012;3.中國(guó)空氣動(dòng)力研究與發(fā)展中心,四川 綿陽(yáng) 621000)

    基于伴隨方程和自由變形技術(shù)的跨聲速機(jī)翼氣動(dòng)設(shè)計(jì)方法研究

    白俊強(qiáng)1,陳頌1,華俊2,孫智偉1,黃江濤3

    (1.西北工業(yè)大學(xué) 航空學(xué)院,陜西 西安 710072;2.中國(guó)航空研究院,北京 100012;3.中國(guó)空氣動(dòng)力研究與發(fā)展中心,四川 綿陽(yáng) 621000)

    將連續(xù)伴隨方程法與自由變形技術(shù)(Eree Eorm Deform-FFD)相結(jié)合開(kāi)展了跨聲速機(jī)翼氣動(dòng)外形優(yōu)化設(shè)計(jì)方法研究。采用Bernstein基函數(shù)建立了空間FFD參數(shù)化方法,并應(yīng)用基于控制理論的連續(xù)伴隨方程方法建立了目標(biāo)函數(shù)對(duì)于待優(yōu)化幾何外形的梯度求解模式,將幾何外形參數(shù)化方法、連續(xù)伴隨方法以及CED數(shù)值模擬技術(shù)相結(jié)合,研究、構(gòu)建了適合跨聲速機(jī)翼的氣動(dòng)外形優(yōu)化設(shè)計(jì)系統(tǒng)。利用該系統(tǒng)對(duì)ONERA M6機(jī)翼及某型民用客機(jī)機(jī)翼進(jìn)行了氣動(dòng)減阻設(shè)計(jì),算例驗(yàn)證表明該方法應(yīng)用于跨聲速機(jī)翼氣動(dòng)減阻設(shè)計(jì)效果明顯,并且能較好的保持幾何表面連續(xù)性和光滑性。

    氣動(dòng)設(shè)計(jì);伴隨方程;自由變形;跨聲速機(jī)翼

    0 引 言

    飛行器的氣動(dòng)設(shè)計(jì)直接影響了飛行器的氣動(dòng)性能以及飛行品質(zhì),是飛行器設(shè)計(jì)中的重要內(nèi)容,一直以來(lái),飛行器設(shè)計(jì)師以及航空學(xué)者們都在為研究開(kāi)發(fā)高效的氣動(dòng)外形設(shè)計(jì)方法和設(shè)計(jì)工具而不斷努力。

    在氣動(dòng)設(shè)計(jì)方法中,Jameson提出的基于控制理論的連續(xù)伴隨方程法[1]以其極高的效率成為當(dāng)前國(guó)內(nèi)外航空界廣泛研究并應(yīng)用的一種設(shè)計(jì)方法。該方法以偏微分方程系統(tǒng)的最優(yōu)控制理論為基礎(chǔ),將待優(yōu)化的氣動(dòng)外形幾何邊界作為控制函數(shù),把流動(dòng)控制方程作為約束條件在目標(biāo)函數(shù)中引入,將約束問(wèn)題轉(zhuǎn)化為無(wú)約束問(wèn)題,設(shè)計(jì)問(wèn)題轉(zhuǎn)化為控制問(wèn)題。通過(guò)流場(chǎng)控制方程及其伴隨方程同時(shí)求解目標(biāo)函數(shù)相對(duì)于所有設(shè)計(jì)變量的梯度,每一次梯度求解的計(jì)算量只相當(dāng)于約兩倍的流場(chǎng)求解計(jì)算量[2],并且?guī)缀跖c設(shè)計(jì)變量個(gè)數(shù)的多少無(wú)關(guān),應(yīng)用這種方法可以大大提高氣動(dòng)外形設(shè)計(jì)的效率。1988年,Jameson等人首先將連續(xù)伴隨方程法應(yīng)用到了跨聲速翼型氣動(dòng)優(yōu)化設(shè)計(jì)中[3-6],從此該方法在跨聲速氣動(dòng)外形優(yōu)化設(shè)計(jì)中獲得了廣泛的應(yīng)用。該方法最初是基于結(jié)構(gòu)網(wǎng)格推導(dǎo)的,1997年Anderson和Venkatakrishnan等人將其推廣到了非結(jié)構(gòu)網(wǎng)格上[7]。國(guó)內(nèi)對(duì)于伴隨方程法的研究開(kāi)始較晚,喬志德、楊旭東等對(duì)結(jié)合控制理論與NS方程的氣動(dòng)反設(shè)計(jì)方法進(jìn)行了研究[8],唐智禮等應(yīng)用特征不變量分析法處理了伴隨方程的物面和遠(yuǎn)場(chǎng)邊界條件,取得了較好效果[9],陳作斌、周鑄、黃勇等將連續(xù)伴隨方程法結(jié)合B樣條參數(shù)化方法發(fā)展了翼型氣動(dòng)優(yōu)化方法[10-11]。

    應(yīng)用伴隨方程法進(jìn)行氣動(dòng)優(yōu)化需要求解目標(biāo)函數(shù)相對(duì)于設(shè)計(jì)變量的梯度,并且根據(jù)梯度信息進(jìn)行設(shè)計(jì)變量的更新,從而得到新的幾何外形,因此設(shè)計(jì)變量的選擇對(duì)整個(gè)氣動(dòng)優(yōu)化過(guò)程和結(jié)果有著決定性的影響。選擇合適的設(shè)計(jì)變量不僅能大大減少優(yōu)化設(shè)計(jì)時(shí)間,而且也是優(yōu)化成功與否的關(guān)鍵[12]。在早期采用伴隨方程法進(jìn)行的氣動(dòng)外形設(shè)計(jì)工作中,每一個(gè)表面網(wǎng)格點(diǎn)都作為設(shè)計(jì)變量,在優(yōu)化三維構(gòu)型的情況下,設(shè)計(jì)變量可能達(dá)到成千上萬(wàn)個(gè),并且容易造成氣動(dòng)外形不光順的情況[12]。R.Hicks[13]和其他人則用一族形函數(shù)對(duì)設(shè)計(jì)空間進(jìn)行了參數(shù)化,這族函數(shù)既可用以產(chǎn)生新的幾何外形,也可用以對(duì)原有的幾何外形產(chǎn)生擾動(dòng),但這些形函數(shù)存在的問(wèn)題是不能形成完備的空間,有可能得不到設(shè)計(jì)問(wèn)題的最優(yōu)解。另一類(lèi)方法是使用B樣條控制點(diǎn)作為設(shè)計(jì)變量,雖然也可以減少設(shè)計(jì)變量的個(gè)數(shù),還可以進(jìn)行局部控制,但在實(shí)際中發(fā)現(xiàn),當(dāng)控制點(diǎn)較多時(shí),B樣條容易產(chǎn)生波動(dòng)[12],不適合于三維氣動(dòng)外形的設(shè)計(jì)。

    自由變形方法(Free-Form Deform FFD)由Sederberg和Parry于1986年首次提出[14],已經(jīng)廣泛應(yīng)用在圖形圖像處理以及動(dòng)畫(huà)設(shè)計(jì)領(lǐng)域,近幾年來(lái)開(kāi)始在飛行器氣動(dòng)外形設(shè)計(jì)中得到初步應(yīng)用。該方法具有變形空間大,不需要對(duì)初始外形進(jìn)行擬合,操作簡(jiǎn)單等優(yōu)點(diǎn),并且可以很好的保持原有幾何的連續(xù)性、光滑性,因此它可以作為一種應(yīng)用于三維氣動(dòng)外形設(shè)計(jì)的參數(shù)化方法。F.Palacios等采用FFD方法對(duì)某超聲速客機(jī)構(gòu)型進(jìn)行了氣動(dòng)優(yōu)化設(shè)計(jì),顯著減小了激波阻力[15],高正紅、黃江濤[16]以及范召林、馬曉永[17,18]等也分別采用FFD技術(shù)對(duì)翼型以及機(jī)翼進(jìn)行了氣動(dòng)優(yōu)化設(shè)計(jì),取得了較好的效果。

    本文采用自由變形方法對(duì)待優(yōu)化幾何外形進(jìn)行參數(shù)化處理,選取FFD控制體頂點(diǎn)的位移作為設(shè)計(jì)變量,并應(yīng)用基于三維Euler方程的連續(xù)伴隨方程方法建立目標(biāo)函數(shù)對(duì)設(shè)計(jì)變量的梯度求解模式,建立了針對(duì)跨聲速機(jī)翼的氣動(dòng)設(shè)計(jì)系統(tǒng),并對(duì)M6機(jī)翼以及某客機(jī)構(gòu)型機(jī)翼進(jìn)行了氣動(dòng)減阻設(shè)計(jì),算例結(jié)果表明該方法優(yōu)化效果明顯,效率較高,并且能夠較好的保持原有幾何外形的連續(xù)性和光滑性。

    1 流動(dòng)控制方程及數(shù)值求解方法

    本文采用三維Euler方程作為流場(chǎng)控制方程,其守恒形式如式(1-1):

    其中x1,x2,x3和u1,u2,u3為笛卡爾坐標(biāo)系下的坐標(biāo)和速度分量,守恒變量w和fi表示為:

    采用有限體積法求解Euler方程,空間離散格式為Jameson-Schmidt-Turkel(JST)中心格式,時(shí)間推進(jìn)方法采用隱式歐拉法。

    2 連續(xù)伴隨方程法

    氣動(dòng)設(shè)計(jì)問(wèn)題是以求解外形變化對(duì)氣動(dòng)特能的影響為基礎(chǔ)而進(jìn)行的,為便于描述連續(xù)伴隨方程法,選定計(jì)算坐標(biāo)系(ξ1,ξ2,ξ3),一般情況下氣動(dòng)特性由目標(biāo)函數(shù)表述為:

    其中dBξ和dDξ分別為計(jì)算空間中的表面與空間積分單元,M與P取決于守恒變量w及計(jì)算空間矩陣S,分別表示目標(biāo)函數(shù)表達(dá)式中的邊界面積分部分和空間體積分部分。在滿足流動(dòng)控制方程約束條件下,氣動(dòng)外形的變化將導(dǎo)致流動(dòng)變量變分δw與矩陣變分δS,因此目標(biāo)函數(shù)的變化可表示為:

    定常狀態(tài)下,對(duì)應(yīng)于氣動(dòng)外形變化的Euler方程可表述為:

    其中Fi為流場(chǎng)守恒量,δFi也可寫(xiě)成與δw、δS相關(guān)的貢獻(xiàn),即:

    引入共軛變量Ψ=(Ψ1,Φ1,Φ2,Φ3,θ)T,在整個(gè)計(jì)算空間求積,則有:

    假定Ψ可微,方程可進(jìn)一步寫(xiě)成:

    由高斯定理可得:

    聯(lián)系式(4)與式(5),把流動(dòng)控制方程變分加入到目標(biāo)函數(shù)變分中,則目標(biāo)函數(shù)I的變分可寫(xiě)為:

    由于假定共軛變量Ψ是任意一個(gè)可微函數(shù),因而ψ就可選取為使δI中不再出現(xiàn)守恒變量變分δw的量,這樣目標(biāo)函數(shù)梯度就只與變換矩陣變分δS有關(guān),而不必重新去求由于設(shè)計(jì)變量的擾動(dòng)而引起的流場(chǎng)物理量變分δw。

    把式(10)空間積分項(xiàng)中流動(dòng)物理量變分δw的系數(shù)項(xiàng)組合在一起,則可得伴隨方程:

    對(duì)應(yīng)的伴隨方程邊界條件可由式(10)邊界積分中流動(dòng)物理量變分δw的系數(shù)項(xiàng)來(lái)得到,即:

    這樣目標(biāo)函數(shù)變分項(xiàng)不再與流動(dòng)變量變分δw相關(guān),剩余項(xiàng)為目標(biāo)函數(shù)變分δI的最終形式:

    在氣動(dòng)減阻設(shè)計(jì)中,目標(biāo)函數(shù)為[5]:

    其中:Sx和Sy是兩個(gè)坐標(biāo)方向上的表面投影面積,Sref是參考面積,dξ1和dξ2是計(jì)算坐標(biāo)系下的所求表面的微分。于是設(shè)計(jì)問(wèn)題轉(zhuǎn)換為了控制問(wèn)題,控制函數(shù)即為幾何外形。選擇最優(yōu)的幾何外形以使目標(biāo)函數(shù)I達(dá)到最小,在給定升力系數(shù)下額外的約束條件為:δCL=0。目標(biāo)函數(shù)的梯度最終可以寫(xiě)作[5]:

    其中:ψT為共軛變量,Qij為相應(yīng)的雅可比矩陣。

    3 自由變形參數(shù)化方法

    自由變形方法(Free Form Deform-FFD)由Sederberg和Parry于1986年首次提出[14],該方法仿照彈性物體受外力后發(fā)生相應(yīng)的變形這一物理現(xiàn)象,將研究對(duì)象置于控制體當(dāng)中,給控制體施加外力,則控制體內(nèi)的所有幾何發(fā)生相應(yīng)變形,處于其中的研究對(duì)象形狀也即發(fā)生變化;這一過(guò)程用數(shù)學(xué)關(guān)系表示即為:移動(dòng)控制體中的控制點(diǎn)(這使得控制體形狀改變,模仿了控制體受力后形狀改變的現(xiàn)象),研究對(duì)象的每一點(diǎn)隨控制頂點(diǎn)的移動(dòng),按照某一函數(shù)關(guān)系相應(yīng)的移動(dòng),其幾何形狀隨之改變。

    本文采用六面體形狀控制體,以控制體的三條邊的方向建立局部坐標(biāo)系,S、T、U是局部坐標(biāo)系軸矢量,沿三個(gè)邊的方向分別均勻布置l+1、m+1、n+1個(gè)控制節(jié)點(diǎn),控制點(diǎn)沿各邊均勻分布,令Pi,j,k為控制點(diǎn)的全局坐標(biāo)值;

    其中,i=0,1,…l;j=0,1,…m;k=0,1,…n。

    令X0(O)為局部坐標(biāo)系的原點(diǎn)在全局坐標(biāo)系下的坐標(biāo);X是控制框架內(nèi)任意一點(diǎn)的全局坐標(biāo),s、t、u分別為該點(diǎn)的局部坐標(biāo),0≤s,t,u≤1,則有:

    其中:

    本文采用Bernstein基函數(shù)建立控制體各頂點(diǎn)坐標(biāo)與控制體內(nèi)任意一點(diǎn)坐標(biāo)之間的映射關(guān)系,

    其中:Xffd為控制框架內(nèi)任意一點(diǎn)的全局坐標(biāo)值;Pi,j,k為控制點(diǎn)的全局坐標(biāo)值;Bil(s)、Bjm(t)、Bkm(u)分別為l、m、n次Bernstein多項(xiàng)式基函數(shù);s、t、u分別為該點(diǎn)在局部坐標(biāo)系O’STU中沿S、T、U方向的坐標(biāo)值。當(dāng)控制點(diǎn)的位置變化,Pi,j,k亦隨之改變,導(dǎo)致控制體內(nèi)任意一點(diǎn)坐標(biāo)受此擾動(dòng)后也發(fā)生改變。

    當(dāng)構(gòu)造較為復(fù)雜的形體時(shí),需要對(duì)變形對(duì)象施加兩個(gè)或者更多的FFD變形,對(duì)于局部變形來(lái)說(shuō),為保持切矢以及曲率連續(xù),對(duì)FFD空間拓?fù)湟约翱刂埔蟾鼑?yán)格。對(duì)于局部參數(shù)為(v,m)的曲面可以表示為:

    任意FFD空間組合X0(s,t,u),X1(s,t,u)共享邊界s0=s1=0,曲面變形后沿v,m方向一階導(dǎo)矢為:

    4 FFD技術(shù)與連續(xù)伴隨方程相結(jié)合的氣動(dòng)外形優(yōu)化設(shè)計(jì)系統(tǒng)

    為了克服以每一個(gè)網(wǎng)格點(diǎn)作為設(shè)計(jì)變量而造成的設(shè)計(jì)外形不光順或者采用B樣條參數(shù)化方法造成的波動(dòng)等問(wèn)題,本文將基于連續(xù)伴隨方程的梯度求解模式與FFD方法相結(jié)合,針對(duì)跨聲速機(jī)翼外形進(jìn)行了氣動(dòng)減阻設(shè)計(jì)。

    采用自由變形參數(shù)化方法將待優(yōu)化幾何外形進(jìn)行參數(shù)化,以每個(gè)FFD控制頂點(diǎn)的位移作為設(shè)計(jì)變量,求解目標(biāo)函數(shù)相對(duì)于設(shè)計(jì)變量的梯度。根據(jù)梯度信息更新設(shè)計(jì)變量,在設(shè)計(jì)變量的值更新后,采用FFD方法進(jìn)行幾何外形的變形,不需要進(jìn)行網(wǎng)格光順即可保持原有幾何屬性如連續(xù)性、光滑性等。目標(biāo)函數(shù)相對(duì)于每個(gè)設(shè)計(jì)變量的梯度可以通過(guò)式(22)求得,其中為FFD控制頂點(diǎn)的位置改變:

    采用FFD方法可以處理機(jī)翼氣動(dòng)外形設(shè)計(jì)在工程實(shí)際應(yīng)用中需要考慮的容積約束問(wèn)題。由FFD變形所引起的任意實(shí)體微元的變化可由FFD的雅可比行列式得到,若FFD變形由下式給出:

    則其雅可比行列式為:

    若變形前微元的體積為dxdydz,則變形后此微元的體積變?yōu)镴(E)dxdydz,變形后整個(gè)實(shí)體的體積是雅可比行列式在變形前實(shí)體上的三重積分,因此,若已知J(E)在變形區(qū)域上的上下界,當(dāng)J(E)由三變量的Bernstein多項(xiàng)式表示時(shí),則可以估計(jì)出其上下界。此時(shí),多項(xiàng)式的最大和最小系數(shù)即是體積變化的上下界[19]。

    結(jié)合FFD技術(shù)和連續(xù)伴隨方程法構(gòu)建機(jī)翼氣動(dòng)外形優(yōu)化設(shè)計(jì)系統(tǒng),優(yōu)化流程如圖1所示。以數(shù)值計(jì)算中的表面網(wǎng)格點(diǎn)定義待優(yōu)化的幾何外形,首先在待優(yōu)化外形周?chē)贾肍FD控制體,并求解其在控制體中的局部坐標(biāo)(即式(18)中的s、t、u),建立從幾何空間到局部空間的映射關(guān)系。以FFD控制頂點(diǎn)的位移作為設(shè)計(jì)變量,實(shí)現(xiàn)對(duì)待優(yōu)化外形的幾何參數(shù)化。而后由CED求解器采用三維歐拉方程(式(1))對(duì)初始外形流場(chǎng)進(jìn)行數(shù)值計(jì)算,待計(jì)算收斂后,采用空間流場(chǎng)信息作為輸入數(shù)據(jù),求解連續(xù)伴隨方程(式(11)),得到優(yōu)化目標(biāo)函數(shù)對(duì)于設(shè)計(jì)變量的梯度,根據(jù)梯度信息采用優(yōu)化算法更新設(shè)計(jì)變量的值,由新的設(shè)計(jì)變量應(yīng)用自由變形方法對(duì)待優(yōu)化幾何外形進(jìn)行變形,并且根據(jù)變形后的FFD控制體估算機(jī)翼的容積,若容積減小則按照外點(diǎn)罰函數(shù)方法進(jìn)行懲罰。根據(jù)變形后的幾何外形表面網(wǎng)格,采用動(dòng)網(wǎng)格方法進(jìn)行相應(yīng)的空間網(wǎng)格變形,生成與新的外形相對(duì)應(yīng)的空間網(wǎng)格。而后對(duì)新的外形進(jìn)行流場(chǎng)數(shù)值計(jì)算,并判斷優(yōu)化是否完成,如未完成,繼續(xù)求解梯度并且更新設(shè)計(jì)變量直至優(yōu)化收斂。本文中動(dòng)網(wǎng)格技術(shù)采用了Torsional Spring動(dòng)網(wǎng)格方法[20]。

    圖1 氣動(dòng)設(shè)計(jì)方法流程圖Fig.1 Flow chart of aerodynamic design based on continuous adjoint method

    5 驗(yàn)證算例

    本文對(duì)于ONERA M6機(jī)翼以及某型民用客機(jī)機(jī)翼外形進(jìn)行了氣動(dòng)減阻設(shè)計(jì)。

    5.1 ONERA M6機(jī)翼氣動(dòng)減阻設(shè)計(jì)

    在M6機(jī)翼周?chē)FD控制框,如圖2所示FFD控制點(diǎn)共計(jì)126個(gè),控制點(diǎn)沿各邊均勻分布。以消除機(jī)翼上表面明顯的λ形激波為氣動(dòng)減阻設(shè)計(jì)的主要目的,把控制體上表面各控制點(diǎn)垂直水平面方向的位移作為設(shè)計(jì)變量,則設(shè)計(jì)變量共63個(gè)。約束條件為升力系數(shù)以及機(jī)翼內(nèi)部容積相對(duì)于原始構(gòu)型不減小,優(yōu)化算法為最速下降法。整個(gè)優(yōu)化過(guò)程共調(diào)用CED求解器及伴隨方程求解器各23次,共花費(fèi)計(jì)算機(jī)機(jī)時(shí)2小時(shí)(8核PC機(jī)),從初始外形得到設(shè)計(jì)外形。由于FFD方法的特性,優(yōu)化后的幾何外形保持了原有的連續(xù)性和光滑性。優(yōu)化前后機(jī)翼幾何外形以及FFD控制點(diǎn)分布的對(duì)比如圖2所示。設(shè)計(jì)前后表面壓力分布對(duì)比如圖3所示。從壓力分布的對(duì)比可見(jiàn),經(jīng)過(guò)氣動(dòng)優(yōu)化設(shè)計(jì),機(jī)翼上表面的激波顯著減弱,優(yōu)化前后機(jī)翼阻力系數(shù)的對(duì)比如表1所示。優(yōu)化設(shè)計(jì)的收斂過(guò)程如圖4所示。經(jīng)過(guò)優(yōu)化設(shè)計(jì)后機(jī)翼阻力系數(shù)相比原始構(gòu)型從0.0123減小到0.0093,相比原始構(gòu)型阻力減小約24.39%,升阻比從21.38提高到28.17,相比原始構(gòu)型提高了31.76%,同時(shí)升力系數(shù)下降了0.001,相對(duì)容積只減少了1.3%。

    圖2 M6機(jī)翼優(yōu)化前后表面壓力分布以及控制點(diǎn)分布對(duì)比Fig.2 Surface pressuredistribution and thedistribution of control points of M6 wing before and after optimization

    圖3 M6機(jī)翼優(yōu)化前后表面壓力分布對(duì)比Fig.3 Surface pressure isobardistribution of M6 wing before and after optimization

    圖4 M6機(jī)翼優(yōu)化收斂歷史Fig.4 Convergence history of thedesign of M6

    表1 M6機(jī)翼優(yōu)化前后氣動(dòng)特性對(duì)比Table 1 Comparison of aerodynamic performance of the M6 wing before and after opimization

    5.2 某型民用客機(jī)機(jī)翼氣動(dòng)減阻設(shè)計(jì)

    針對(duì)某民用客機(jī)機(jī)翼外形進(jìn)行了第二個(gè)驗(yàn)證算例,機(jī)翼幾何外形以及FFD控制頂點(diǎn)分布如圖5所示,F(xiàn)FD控制點(diǎn)共198個(gè),取FFD控制體上表面的控制頂點(diǎn)垂直水平面的位移量作為設(shè)計(jì)變量,則設(shè)計(jì)變量共計(jì)99個(gè)。約束條件為升力系數(shù)和相對(duì)容積相比原始構(gòu)型不減小,優(yōu)化算法為最速下降法。整個(gè)優(yōu)化過(guò)程共調(diào)用CED求解器及伴隨方程求解器各40次,共花費(fèi)計(jì)算機(jī)機(jī)時(shí)20小時(shí)(8核PC機(jī))。優(yōu)化后機(jī)翼表面壓力分布以及FFD控制頂點(diǎn)分布如圖6所示。優(yōu)化前后機(jī)翼的氣動(dòng)特性參數(shù)對(duì)比如表2所示,優(yōu)化后機(jī)翼阻力系數(shù)從0.0261減小到0.0212,優(yōu)化后阻力相對(duì)原始構(gòu)型減小18.77%,升阻比從20.07提高到了25.94,提高29.3%,同時(shí)優(yōu)化前后相對(duì)容積只減小了1.9%,升力系數(shù)不變,并且優(yōu)化后機(jī)翼表面幾何外形很好的保持了原有外形的連續(xù)性和光滑性。優(yōu)化前后機(jī)翼表面壓力分布對(duì)比如圖7、圖8所示,優(yōu)化設(shè)計(jì)的收斂過(guò)程如圖9所示。可見(jiàn)經(jīng)過(guò)設(shè)計(jì)之后,機(jī)翼內(nèi)外翼的激波均明顯減弱,減阻效果顯著。

    圖5 某型客機(jī)初始機(jī)翼表面壓力分布以及控制點(diǎn)分布Fig.5 Initial surface pressure distribution and the distribution of control points

    圖6 某型客機(jī)優(yōu)化后機(jī)翼表面壓力分布以及控制點(diǎn)分布Fig.6 Optimum surface pressure distribution and the distribution of control points

    圖7 某型客機(jī)初始機(jī)翼表面壓力分布Fig.7 Initial surface pressure distribution of civil jet wing

    圖8 某型客機(jī)優(yōu)化后機(jī)翼表面壓力分布Fig.8 Optimum surface pressure distribution of wing

    圖9 某型客機(jī)優(yōu)化后機(jī)翼表面壓力分布Fig.9 Convergence history of the design of M6

    表2 某型客機(jī)機(jī)翼優(yōu)化前后氣動(dòng)特性對(duì)比Table 2 Comparison of aerodynamic performance of the aero transport’s wing before and after optimization

    6 結(jié) 論

    (1)將FFD技術(shù)與連續(xù)伴隨方程法相結(jié)合,以FFD控制點(diǎn)位移代替每個(gè)表面網(wǎng)格點(diǎn)的位移作為設(shè)計(jì)變量,可以使設(shè)計(jì)外形很好的保持原有外形的幾何連續(xù)和光滑性。

    (2)采用FFD技術(shù)進(jìn)行氣動(dòng)優(yōu)化設(shè)計(jì),可以有效估計(jì)優(yōu)化過(guò)程中所產(chǎn)生的每個(gè)結(jié)果的幾何容積。

    (3)兩個(gè)驗(yàn)證算例中構(gòu)型的氣動(dòng)減阻設(shè)計(jì)在8核PC機(jī)的計(jì)算平臺(tái)下均在一天之內(nèi)完成(分別為3小時(shí)和20小時(shí)),從設(shè)計(jì)結(jié)果可見(jiàn),F(xiàn)FD技術(shù)和連續(xù)伴隨方程相結(jié)合的跨聲速機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)系統(tǒng)對(duì)于跨聲速機(jī)翼的氣動(dòng)減阻設(shè)計(jì)效果明顯且效率較高。

    [1]JAMESON A.Aerodynamicdesign via control theory[J].Journal of Scientific Computing,1988,3:233-260.

    [2] WIDHALM M,RONZHEIMER A,HEPPERLE M.Comparison between gradient free and adjoint based aerodynamic optimization of a flying wing transport aircraft in the preliminarydesign[R].AIAA 2007-4060.

    [3]JAMESON A.Optimum aerodynamicdesign using CED and control theory[R].AIAA 95-1729-CP,1995.

    [4]JAMESON A.Control theory based airfoildesign using Euler equation[R].AIAA 94-472-CP,1994.

    [5]REUTHER J,JAMESON A.Control theory based airfoildesign for potential flow and a finite volumediscretization[R].AIAA 91-499,1994.

    [6]REUTHER J,JAMESON A.Aerodynamic shape optimization of complex aircraft configuration via an adjoint formulation[R].AIAA 96-0094.1996.

    [7]ANDERSON W K,VENKATAKRISHNAN V.Aerodynamicdesign optimization on unstructured grid with a continuous adjoint formulation[R].AIAA 97-0643,1997.

    [8]YANG X D,QIAO Z D,ZHU B.Aerodynamicdesign method based on control theory and Navier-Stokes equations[J].ACTA Aerodynamica Sinica,2005,23(1):46-52.(in Chinese)

    楊旭東,喬志德,朱兵.基于控制理論和NS方程的氣動(dòng)設(shè)計(jì)方法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2005,23(1):46-52.

    [9]TANG Z L,HUANG M K.Control theory based airfoildesign using Euler equations[J].ACTA Aerod ynamica Sinica,2001,19(3):262-270.(in Chinese)

    唐智禮,黃明?。诳刂评碚摰腅uler方程翼型減阻優(yōu)化設(shè)計(jì)[J].空氣動(dòng)力學(xué)學(xué)報(bào),2001,19(3):262-270.

    [10]HUANG Y,CHEN Z B,LIU G.An investigation of aerodynamic optimizationdesign for airfoil based on adjoint formulation[J].ACTA Aerod ynamica Sinica,1999,17(04):413-422.(in Chinese)

    黃勇,陳作斌,劉剛.基于伴隨方程的翼型數(shù)值優(yōu)化設(shè)計(jì)方法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),1999,17(04):413-422.

    [11]ZHOU Z,CHEN Z B.Airfoil aerodynamicdesign and optimization based on Navier-Stokes equations[J].ACTA Aerodynamica Sinica,2002,20(02):141-149

    周鑄,陳作斌.基于N-S方程的翼型氣動(dòng)優(yōu)化設(shè)計(jì)[J].空氣動(dòng)力學(xué)學(xué)報(bào),2002,20(02):141-149

    [12]XIONG J T.Aerodynamic shape optimizationdesign based on adjoint method and Navier-Stokes equation[D].Xi′an:Northwestern Polytechnical University.2007

    熊俊濤.基于伴隨方法和Navier-Stokes方程的氣動(dòng)優(yōu)化設(shè)計(jì)[D].西安:西北工業(yè)大學(xué),2007

    [13]HICKS R M,HENNE P A.Wingdesign by numerical optimization[R].AIAA 79-0080.1979.

    [14]SCOTT T S,PARRY R.Eree-formdeformation of solid geometric models[J].Computer Graphics,86,20(4):150-161.(SIGGRAPH 86.)

    [15]PALACIOS E,ALONSO J J,COLONO M,et al.Adjointbased method for supersonic aircraftdesign using equivalent areadistribution[R].AIAA 2012-0269,2012.

    [16]HUANG J T,GAO Z H,BAIJ Q,et al.Laminar airfoil aerodynamic shape optimizationdesign based on Delaunay graph mapping and FFD technique[J].Acta Aeronautica et Astronautica,2012,33(10):1817-1826

    黃江濤,高正紅,白俊強(qiáng),等.應(yīng)用Delaunay圖映射與FFD技術(shù)的層流翼型氣動(dòng)優(yōu)化設(shè)計(jì)[J].航空學(xué)報(bào),2012,33(10):1817-1826.

    [17]MA X Y,WU W H,EAN Z L.Ereedeformation method of wing based on NURBS[J].Journal of Sichuan University(Engineering Science Edition),2010,42(增2):195-197.

    馬曉永,吳文華,范召林.基于NURBS的機(jī)翼自由變形方法[J].四川大學(xué)學(xué)報(bào),2010,42(增2):195-197.

    [18]MA X Y,EAN Z L,WU W H,et al.Aerodynamic shape optimization based on FFD[C].13rdSystem Simulation Technology&Application,2011.

    馬曉永,范召林,吳文華,等.基于FFD技術(shù)的機(jī)翼氣動(dòng)外形優(yōu)化仿真[C].第13屆中國(guó)系統(tǒng)仿真技術(shù)及其應(yīng)用學(xué)術(shù)年會(huì)(13卷),2011.

    [19]朱心雄.自由曲線曲面造型技術(shù)[M].北京:科學(xué)出版社,2000.

    [20]DEGAND C,EARHAT C.A three-dimensional torsional spring analogy method for unstructureddynamic meshes[J].Computers&Structures,2002,80(3):305-316.

    Transonic wing aerodynamicdesign based on
    continuous adjoint method and free formdeform technique

    BAI Junqiang1,CHEN Song1,HUA Jun2,SUN Zhiwei1,HUANG Jiangtao1
    (1.Northwest Polytechnical University,Xi’an 710072,China;2.China Aeronautical Establishment,Beijing 100012,China)

    The approach of transonic wing aerodynamicdesign based on the combination of continuous adjoint method with free formdeform(FFD)technique has been presented.Bernstein basis function has been chosen to establish the FFD parameterization of arbitrary spatial property to overcome thedisadvantages ofdiscontinuity or unsmoothness by the bump function method or the B-spline approach.Besides,the scheme to compute the sensitivity of aerodynamic objective function respecting to geometrical shape has been established using the continuous adjoint method based on control theory,in which the sensitivity computation of all thedesign variables would cost one flow field computation and one adjoint field computation.Combinedwith torsional spring meshdeformation technique,F(xiàn)FD geometric parameterization,continuous adjoint method and computational fluiddynamics have been integrated to establish adesign system for transonic wing aerodynamic configuration aimmed atdrag reduction,which is good at both efficiency and robustness.This system has been testified with optimizationdesign cases of ONERA M6 wing and a typical large civil jet wing,the results of which have shown that this method is of good feasibility and efficiency in aerodynamic optimizationdesign fordrag reduction.Moreover,the results have shown that this method is of good capability of conserving the geometry continuity and smoothness.

    aerodynamicdesign;continuous adjoint method;free formdeform;transonic wing

    V211.3;TP273

    Adoi:10.7638/kqdlxxb-2012.0198

    0258-1825(2014)06-0820-08

    2012-12-03;

    2013-03-10

    973計(jì)劃大型客機(jī)減阻機(jī)理和方法研究(2014CB744804)

    白俊強(qiáng)(1970-),男,河南新鄉(xiāng)人,教授,主要研究方向:飛行器設(shè)計(jì).

    白俊強(qiáng),陳頌,華俊,等.基于伴隨方程和自由變形技術(shù)的跨聲速機(jī)翼氣動(dòng)設(shè)計(jì)方法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2014,32(6):820-826.

    10.7638/kqdlxxb-2012.0198 BAI J Q,CHEN S,HUA J,et al.Transonic wing aerodynamicdesign based on continuous adjoint method and free formdeform technique[J].ACTA Aerodynamica Sinica,2014,32(6):820-826.

    猜你喜歡
    機(jī)翼外形氣動(dòng)
    中寰氣動(dòng)執(zhí)行機(jī)構(gòu)
    適盒A4BOX 多功能料理鍋
    基于NACA0030的波紋狀翼型氣動(dòng)特性探索
    變時(shí)滯間隙非線性機(jī)翼顫振主動(dòng)控制方法
    基于反饋線性化的RLV氣動(dòng)控制一體化設(shè)計(jì)
    論袁牧之“外形的演技”
    足趾移植再造手指術(shù)后外形的整形
    機(jī)翼跨聲速抖振研究進(jìn)展
    KJH101-127型氣動(dòng)司控道岔的改造
    基于模糊自適應(yīng)的高超聲速機(jī)翼顫振的主動(dòng)控制
    少妇熟女欧美另类| 国产成人精品福利久久| freevideosex欧美| 草草在线视频免费看| 国产视频首页在线观看| 简卡轻食公司| 国产伦精品一区二区三区四那| 成年女人在线观看亚洲视频 | 精品久久久久久电影网| 视频中文字幕在线观看| 久久久久性生活片| 午夜视频国产福利| 亚洲精品日韩在线中文字幕| 国产亚洲最大av| 成年人午夜在线观看视频| 日本三级黄在线观看| 人妻一区二区av| 亚洲激情五月婷婷啪啪| 国产高清国产精品国产三级 | 观看免费一级毛片| 国产毛片a区久久久久| 街头女战士在线观看网站| 男女国产视频网站| 国产精品一区www在线观看| 亚洲人成网站高清观看| 肉色欧美久久久久久久蜜桃 | 国产一区有黄有色的免费视频| 国产毛片在线视频| 亚洲精品,欧美精品| 激情五月婷婷亚洲| 建设人人有责人人尽责人人享有的 | 久久午夜福利片| 香蕉精品网在线| 人人妻人人爽人人添夜夜欢视频 | 秋霞伦理黄片| 欧美精品一区二区大全| 高清av免费在线| 五月天丁香电影| 国产亚洲av片在线观看秒播厂| 国产探花极品一区二区| 日韩伦理黄色片| 久久99热这里只频精品6学生| 国产老妇伦熟女老妇高清| 中文在线观看免费www的网站| 亚洲在久久综合| 久久99精品国语久久久| 精品少妇久久久久久888优播| 色婷婷久久久亚洲欧美| 赤兔流量卡办理| 午夜福利高清视频| 国产精品精品国产色婷婷| 久久人人爽人人片av| 精品一区二区三卡| 男女边摸边吃奶| 有码 亚洲区| 在线 av 中文字幕| 看黄色毛片网站| 亚洲精品亚洲一区二区| 日韩成人av中文字幕在线观看| 精品人妻熟女av久视频| 午夜视频国产福利| 男人舔奶头视频| 久久精品综合一区二区三区| 国产美女午夜福利| 国产精品不卡视频一区二区| 国产极品天堂在线| 久久97久久精品| 亚洲精品国产av蜜桃| 色综合色国产| 免费不卡的大黄色大毛片视频在线观看| 午夜福利视频精品| 亚洲精品久久久久久婷婷小说| 五月玫瑰六月丁香| 精品人妻熟女av久视频| 亚洲欧美成人精品一区二区| av播播在线观看一区| 一级毛片久久久久久久久女| 亚洲久久久久久中文字幕| 成人免费观看视频高清| 日韩不卡一区二区三区视频在线| av国产久精品久网站免费入址| 99热这里只有是精品在线观看| 免费人成在线观看视频色| 身体一侧抽搐| 在线观看国产h片| 国产亚洲一区二区精品| 高清视频免费观看一区二区| 一区二区av电影网| 国产淫片久久久久久久久| 成年女人在线观看亚洲视频 | 久久精品久久久久久久性| 尾随美女入室| 成人黄色视频免费在线看| 日本av手机在线免费观看| 国产v大片淫在线免费观看| 欧美激情国产日韩精品一区| 国产 一区精品| 久久97久久精品| 久久精品国产a三级三级三级| 亚洲最大成人中文| 亚洲av男天堂| 亚洲欧美一区二区三区国产| 国产毛片a区久久久久| 精品一区二区免费观看| 毛片女人毛片| 亚洲国产精品999| 深夜a级毛片| 色播亚洲综合网| 又爽又黄a免费视频| 亚洲欧美一区二区三区黑人 | 亚洲美女视频黄频| 男插女下体视频免费在线播放| 国产日韩欧美亚洲二区| 国产乱人视频| 欧美性猛交╳xxx乱大交人| 婷婷色av中文字幕| 美女脱内裤让男人舔精品视频| 日韩伦理黄色片| 日韩欧美精品免费久久| 小蜜桃在线观看免费完整版高清| 自拍偷自拍亚洲精品老妇| 免费大片18禁| 内地一区二区视频在线| 边亲边吃奶的免费视频| 成人特级av手机在线观看| 美女高潮的动态| 亚洲最大成人中文| 丰满人妻一区二区三区视频av| 国语对白做爰xxxⅹ性视频网站| 国产久久久一区二区三区| 大陆偷拍与自拍| 777米奇影视久久| 最近2019中文字幕mv第一页| 色综合色国产| 有码 亚洲区| 啦啦啦啦在线视频资源| av福利片在线观看| 成人特级av手机在线观看| 美女视频免费永久观看网站| 插逼视频在线观看| 欧美性感艳星| 中文在线观看免费www的网站| 免费av不卡在线播放| 精品人妻偷拍中文字幕| 国产精品女同一区二区软件| av黄色大香蕉| 亚洲精品亚洲一区二区| 毛片一级片免费看久久久久| 国产一区有黄有色的免费视频| 夜夜看夜夜爽夜夜摸| 一级片'在线观看视频| 久久人人爽人人片av| 黄色一级大片看看| 麻豆国产97在线/欧美| 色视频www国产| 亚洲色图av天堂| 一级毛片 在线播放| 亚洲人与动物交配视频| 国产精品一区二区在线观看99| 精品久久久久久久久亚洲| av国产久精品久网站免费入址| videos熟女内射| 亚洲第一区二区三区不卡| 亚洲精品国产成人久久av| 极品教师在线视频| 高清午夜精品一区二区三区| 国产黄片视频在线免费观看| 乱码一卡2卡4卡精品| 日韩欧美 国产精品| 欧美性猛交╳xxx乱大交人| 国产久久久一区二区三区| 边亲边吃奶的免费视频| 精品久久久久久久末码| 久热这里只有精品99| av免费观看日本| 一本久久精品| 欧美激情国产日韩精品一区| 少妇的逼好多水| 成人二区视频| 在线观看av片永久免费下载| 亚洲,一卡二卡三卡| 国产乱人偷精品视频| 午夜福利高清视频| 插阴视频在线观看视频| 99re6热这里在线精品视频| www.色视频.com| 欧美xxxx性猛交bbbb| 免费看光身美女| 日韩大片免费观看网站| 亚洲精品456在线播放app| 免费观看av网站的网址| 欧美激情在线99| 天堂网av新在线| 街头女战士在线观看网站| 插阴视频在线观看视频| 日韩在线高清观看一区二区三区| 免费观看a级毛片全部| 婷婷色综合www| 亚洲精华国产精华液的使用体验| 国内揄拍国产精品人妻在线| 国产精品国产三级国产专区5o| 亚洲国产日韩一区二区| a级毛片免费高清观看在线播放| 久久精品久久久久久久性| 国产精品久久久久久久电影| 噜噜噜噜噜久久久久久91| 熟妇人妻不卡中文字幕| 久久99精品国语久久久| 丝袜脚勾引网站| 看免费成人av毛片| 综合色av麻豆| 伊人久久国产一区二区| 久久国产乱子免费精品| 亚洲人与动物交配视频| 插逼视频在线观看| 色视频在线一区二区三区| 久久精品国产亚洲av天美| 22中文网久久字幕| 日本一二三区视频观看| 亚洲av二区三区四区| 国产一区二区亚洲精品在线观看| 在线亚洲精品国产二区图片欧美 | 三级国产精品片| 三级男女做爰猛烈吃奶摸视频| 777米奇影视久久| 久久久久久久大尺度免费视频| 在线天堂最新版资源| 中国国产av一级| 国产 一区 欧美 日韩| 免费高清在线观看视频在线观看| 国产精品久久久久久av不卡| 不卡视频在线观看欧美| av一本久久久久| 久久精品国产a三级三级三级| 男的添女的下面高潮视频| 丝袜美腿在线中文| 少妇的逼水好多| 国产成人91sexporn| 韩国高清视频一区二区三区| 男女啪啪激烈高潮av片| 婷婷色综合www| 亚洲最大成人中文| 成人漫画全彩无遮挡| 免费黄频网站在线观看国产| 日本午夜av视频| 国产片特级美女逼逼视频| 亚洲色图av天堂| 在线观看国产h片| 国产成人a区在线观看| 哪个播放器可以免费观看大片| 亚洲综合精品二区| 久久久久久久久久成人| 国产高清国产精品国产三级 | 插逼视频在线观看| 男男h啪啪无遮挡| 在线播放无遮挡| 在线观看三级黄色| 新久久久久国产一级毛片| 2021天堂中文幕一二区在线观| 大片免费播放器 马上看| 欧美高清成人免费视频www| 菩萨蛮人人尽说江南好唐韦庄| 少妇 在线观看| 亚洲国产高清在线一区二区三| 亚洲最大成人av| 丝袜喷水一区| 中文欧美无线码| 国产淫片久久久久久久久| 人人妻人人澡人人爽人人夜夜| 亚洲内射少妇av| 尤物成人国产欧美一区二区三区| 国产永久视频网站| 日韩欧美精品免费久久| 日本熟妇午夜| 又粗又硬又长又爽又黄的视频| 男人狂女人下面高潮的视频| 天天躁夜夜躁狠狠久久av| 日韩欧美精品免费久久| 亚洲国产日韩一区二区| 日韩强制内射视频| 国产精品麻豆人妻色哟哟久久| 欧美日韩在线观看h| 麻豆成人午夜福利视频| 亚洲av二区三区四区| 午夜亚洲福利在线播放| a级一级毛片免费在线观看| 97超视频在线观看视频| 91精品国产九色| 91午夜精品亚洲一区二区三区| 亚洲精品影视一区二区三区av| a级一级毛片免费在线观看| 一本色道久久久久久精品综合| 韩国高清视频一区二区三区| 午夜福利高清视频| 午夜福利视频1000在线观看| 大香蕉久久网| 中国国产av一级| 美女被艹到高潮喷水动态| 久久久久久久久大av| 香蕉精品网在线| 免费看a级黄色片| 麻豆成人午夜福利视频| 国产有黄有色有爽视频| 国产综合精华液| 另类亚洲欧美激情| 啦啦啦在线观看免费高清www| 久久久久久久国产电影| 中文欧美无线码| 免费黄色在线免费观看| 身体一侧抽搐| 插阴视频在线观看视频| 久久久久久久久久久免费av| 久久综合国产亚洲精品| 午夜免费男女啪啪视频观看| 又大又黄又爽视频免费| 国产 精品1| 午夜免费观看性视频| 亚洲成人精品中文字幕电影| 亚洲av电影在线观看一区二区三区 | av在线老鸭窝| av网站免费在线观看视频| 综合色av麻豆| 国产精品国产三级国产专区5o| 成人免费观看视频高清| 亚洲最大成人av| 日韩欧美一区视频在线观看 | 日韩欧美精品v在线| 水蜜桃什么品种好| av福利片在线观看| 亚洲国产欧美人成| 嘟嘟电影网在线观看| 日韩av在线免费看完整版不卡| 一本久久精品| 夫妻性生交免费视频一级片| 国产av国产精品国产| 亚洲国产精品999| 亚洲精品国产成人久久av| 国产欧美日韩精品一区二区| 一本色道久久久久久精品综合| 国产爽快片一区二区三区| 亚洲精品成人久久久久久| 性色av一级| 久久久午夜欧美精品| 久久精品夜色国产| 国产女主播在线喷水免费视频网站| 欧美成人一区二区免费高清观看| 精品少妇黑人巨大在线播放| 青春草亚洲视频在线观看| 国产午夜精品一二区理论片| 一级a做视频免费观看| 国产淫语在线视频| av又黄又爽大尺度在线免费看| 国产91av在线免费观看| 日韩中字成人| 成人二区视频| 亚洲欧美日韩无卡精品| 性色av一级| 国产精品成人在线| 久久久久性生活片| 久久久午夜欧美精品| 亚洲精品成人av观看孕妇| 22中文网久久字幕| 国产美女午夜福利| 久久精品人妻少妇| 少妇 在线观看| 亚洲一区二区三区欧美精品 | 深夜a级毛片| 婷婷色av中文字幕| 欧美日韩视频精品一区| 男人狂女人下面高潮的视频| 久久国产乱子免费精品| 国产精品久久久久久久久免| 精品人妻视频免费看| 亚洲精品第二区| 亚洲av国产av综合av卡| 狂野欧美白嫩少妇大欣赏| 亚洲精品aⅴ在线观看| 亚洲av免费高清在线观看| 寂寞人妻少妇视频99o| 2021少妇久久久久久久久久久| 免费av不卡在线播放| 亚洲婷婷狠狠爱综合网| 亚洲成色77777| 欧美xxxx黑人xx丫x性爽| 成人亚洲欧美一区二区av| 永久免费av网站大全| 国产在视频线精品| 插逼视频在线观看| 久久国内精品自在自线图片| 国产免费福利视频在线观看| 精品一区二区三区视频在线| 久久久久久九九精品二区国产| 婷婷色综合www| 制服丝袜香蕉在线| 在线免费十八禁| 制服丝袜香蕉在线| 亚洲国产精品成人久久小说| 22中文网久久字幕| 国产中年淑女户外野战色| 久久久午夜欧美精品| 天堂俺去俺来也www色官网| 免费看a级黄色片| 超碰97精品在线观看| 成人黄色视频免费在线看| 欧美极品一区二区三区四区| 成人国产麻豆网| 中文字幕亚洲精品专区| 国产黄频视频在线观看| 欧美3d第一页| 哪个播放器可以免费观看大片| 热re99久久精品国产66热6| 国产精品成人在线| 欧美日韩亚洲高清精品| 亚洲欧洲国产日韩| 国产乱人偷精品视频| 丰满人妻一区二区三区视频av| 又爽又黄a免费视频| 亚洲av免费在线观看| 亚洲三级黄色毛片| 亚洲国产欧美人成| 中文资源天堂在线| 校园人妻丝袜中文字幕| 国产黄片美女视频| 亚洲欧美清纯卡通| 亚洲av日韩在线播放| 亚洲欧美一区二区三区国产| 亚洲精品,欧美精品| 欧美日韩亚洲高清精品| 亚洲国产av新网站| videos熟女内射| 亚洲av福利一区| 欧美变态另类bdsm刘玥| 久久综合国产亚洲精品| 色综合色国产| 婷婷色综合www| 男的添女的下面高潮视频| 国产乱来视频区| 哪个播放器可以免费观看大片| 蜜桃亚洲精品一区二区三区| 婷婷色综合大香蕉| 日本三级黄在线观看| 亚洲精品乱久久久久久| 人妻 亚洲 视频| 久久久精品欧美日韩精品| 日韩,欧美,国产一区二区三区| 校园人妻丝袜中文字幕| 夫妻性生交免费视频一级片| 热99国产精品久久久久久7| av网站免费在线观看视频| 韩国高清视频一区二区三区| 最新中文字幕久久久久| 国产日韩欧美亚洲二区| 精品久久久久久久久av| 熟女av电影| 啦啦啦中文免费视频观看日本| 久久精品国产亚洲av涩爱| 91久久精品国产一区二区成人| 久久女婷五月综合色啪小说 | www.av在线官网国产| 国产免费一区二区三区四区乱码| 国产av码专区亚洲av| 亚洲欧美中文字幕日韩二区| 一级毛片 在线播放| 免费黄频网站在线观看国产| 国产精品爽爽va在线观看网站| 日日摸夜夜添夜夜添av毛片| 在线观看免费高清a一片| 国产av码专区亚洲av| 黄色视频在线播放观看不卡| av免费观看日本| 中文在线观看免费www的网站| 黑人高潮一二区| 国产又色又爽无遮挡免| 日日摸夜夜添夜夜爱| 99视频精品全部免费 在线| 在线免费十八禁| 日韩欧美精品v在线| 亚洲精品一区蜜桃| 国产男女内射视频| 欧美激情国产日韩精品一区| 国产在视频线精品| 国产精品国产三级专区第一集| 美女视频免费永久观看网站| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 深夜a级毛片| 日日啪夜夜爽| 男女边吃奶边做爰视频| 久久鲁丝午夜福利片| 高清欧美精品videossex| 午夜福利在线在线| 99热这里只有精品一区| 九色成人免费人妻av| 国产成人一区二区在线| 国模一区二区三区四区视频| 亚洲精品亚洲一区二区| 国产久久久一区二区三区| 国产黄a三级三级三级人| 国产探花极品一区二区| 2021少妇久久久久久久久久久| 日韩人妻高清精品专区| 国产一级毛片在线| 18禁裸乳无遮挡动漫免费视频 | 精品少妇久久久久久888优播| 三级经典国产精品| 水蜜桃什么品种好| 最新中文字幕久久久久| 美女被艹到高潮喷水动态| 亚洲精品久久午夜乱码| 欧美三级亚洲精品| 色视频www国产| 高清日韩中文字幕在线| 国产精品不卡视频一区二区| 最后的刺客免费高清国语| 婷婷色麻豆天堂久久| 亚洲人成网站在线播| 国产国拍精品亚洲av在线观看| 免费av观看视频| av卡一久久| 久久久久久九九精品二区国产| 国产69精品久久久久777片| 亚洲综合色惰| 看黄色毛片网站| 亚洲精品影视一区二区三区av| 国产成人精品久久久久久| 天堂俺去俺来也www色官网| 七月丁香在线播放| 色5月婷婷丁香| 亚洲欧美清纯卡通| 亚洲成人精品中文字幕电影| 国产一区有黄有色的免费视频| 亚洲精品国产色婷婷电影| 免费大片黄手机在线观看| 免费观看的影片在线观看| 国产精品国产三级专区第一集| 精品99又大又爽又粗少妇毛片| 日韩三级伦理在线观看| 国产成人freesex在线| 免费观看的影片在线观看| 国产午夜精品久久久久久一区二区三区| 亚洲美女搞黄在线观看| 综合色av麻豆| 人体艺术视频欧美日本| 三级国产精品片| 国产精品国产三级国产av玫瑰| 成人午夜精彩视频在线观看| 亚洲婷婷狠狠爱综合网| 免费高清在线观看视频在线观看| 亚洲美女搞黄在线观看| 91久久精品国产一区二区三区| 免费大片黄手机在线观看| 丝袜脚勾引网站| 日韩不卡一区二区三区视频在线| 国产淫语在线视频| 国产成人a∨麻豆精品| 亚洲精品自拍成人| videos熟女内射| 欧美日韩视频高清一区二区三区二| 人人妻人人澡人人爽人人夜夜| 简卡轻食公司| 九九久久精品国产亚洲av麻豆| 欧美97在线视频| 国产av国产精品国产| 国产爱豆传媒在线观看| 国产爽快片一区二区三区| 久久久久久久午夜电影| 国产老妇女一区| 啦啦啦啦在线视频资源| 一级片'在线观看视频| 成年版毛片免费区| 秋霞在线观看毛片| 高清日韩中文字幕在线| 国产探花极品一区二区| 亚洲成人久久爱视频| 一区二区av电影网| 久久亚洲国产成人精品v| 大香蕉97超碰在线| 少妇丰满av| 国产精品人妻久久久久久| 婷婷色综合www| 亚洲精品成人久久久久久| a级一级毛片免费在线观看| 欧美成人一区二区免费高清观看| 男人舔奶头视频| 在线观看人妻少妇| 免费观看性生交大片5| 中文字幕av成人在线电影| 国产日韩欧美亚洲二区| 精品少妇久久久久久888优播| 国产免费福利视频在线观看| 精品99又大又爽又粗少妇毛片| 日本猛色少妇xxxxx猛交久久| 大香蕉久久网| 麻豆乱淫一区二区| 国产一区亚洲一区在线观看| 成年版毛片免费区| 国产精品一二三区在线看| 国产一区亚洲一区在线观看| 伦理电影大哥的女人| 韩国av在线不卡| 国产男人的电影天堂91| 网址你懂的国产日韩在线| 久久久久性生活片| 国产亚洲91精品色在线| h日本视频在线播放| 国产69精品久久久久777片| 少妇的逼好多水| 九九在线视频观看精品| 国产v大片淫在线免费观看| 激情五月婷婷亚洲| 午夜福利在线观看免费完整高清在| 蜜桃亚洲精品一区二区三区| 久久久精品94久久精品| 欧美成人一区二区免费高清观看| 亚洲精品久久久久久婷婷小说| 亚洲天堂国产精品一区在线|