• <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)控制
    欧美乱妇无乱码| 我要看黄色一级片免费的| 在线观看66精品国产| 免费黄频网站在线观看国产| 在线天堂中文资源库| 国产伦人伦偷精品视频| 精品国产乱码久久久久久小说| 国产成人啪精品午夜网站| xxxhd国产人妻xxx| 中文亚洲av片在线观看爽 | 国产又爽黄色视频| 国产免费福利视频在线观看| 十分钟在线观看高清视频www| 丁香六月欧美| 久久久久久人人人人人| 欧美激情 高清一区二区三区| 满18在线观看网站| 可以免费在线观看a视频的电影网站| 窝窝影院91人妻| 麻豆成人av在线观看| 国产在线观看jvid| 欧美日韩亚洲综合一区二区三区_| 汤姆久久久久久久影院中文字幕| 91大片在线观看| 下体分泌物呈黄色| 美女午夜性视频免费| 欧美日韩亚洲综合一区二区三区_| 视频区图区小说| 亚洲色图综合在线观看| 18禁观看日本| 99在线人妻在线中文字幕 | 欧美大码av| 下体分泌物呈黄色| 久热这里只有精品99| 亚洲av第一区精品v没综合| e午夜精品久久久久久久| 日日夜夜操网爽| 建设人人有责人人尽责人人享有的| 日韩精品免费视频一区二区三区| 18在线观看网站| 老司机福利观看| 精品少妇内射三级| 悠悠久久av| 国产欧美日韩一区二区三| 久久精品国产亚洲av高清一级| 欧美日韩一级在线毛片| 亚洲成人手机| 国产精品影院久久| 国产有黄有色有爽视频| av网站在线播放免费| 老司机深夜福利视频在线观看| 人人妻,人人澡人人爽秒播| 在线天堂中文资源库| 久9热在线精品视频| 9色porny在线观看| 亚洲精品国产区一区二| 久久久水蜜桃国产精品网| 欧美精品一区二区免费开放| 青青草视频在线视频观看| 精品国产亚洲在线| 99在线人妻在线中文字幕 | 久久热在线av| 蜜桃国产av成人99| 久久av网站| 久久性视频一级片| 看免费av毛片| 欧美+亚洲+日韩+国产| av有码第一页| 日韩有码中文字幕| 一区二区三区乱码不卡18| 精品人妻1区二区| 国产精品自产拍在线观看55亚洲 | 在线av久久热| 老司机深夜福利视频在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 久久人妻熟女aⅴ| 一区在线观看完整版| 亚洲欧洲精品一区二区精品久久久| 97人妻天天添夜夜摸| av欧美777| 高清黄色对白视频在线免费看| 一级片'在线观看视频| 下体分泌物呈黄色| 国产精品二区激情视频| 一区二区av电影网| 久久国产亚洲av麻豆专区| 伊人久久大香线蕉亚洲五| 视频区欧美日本亚洲| 国产成人一区二区三区免费视频网站| 欧美人与性动交α欧美软件| av欧美777| 中亚洲国语对白在线视频| 夫妻午夜视频| 蜜桃国产av成人99| 90打野战视频偷拍视频| 91麻豆av在线| 亚洲成人免费电影在线观看| 亚洲熟妇熟女久久| 中文字幕人妻熟女乱码| 一本大道久久a久久精品| 国产精品久久久av美女十八| 男女下面插进去视频免费观看| 成年女人毛片免费观看观看9 | 黄色丝袜av网址大全| 九色亚洲精品在线播放| 成年动漫av网址| www.999成人在线观看| 嫁个100分男人电影在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 久久人人爽av亚洲精品天堂| 精品亚洲乱码少妇综合久久| 亚洲国产欧美一区二区综合| 国产一区二区 视频在线| 蜜桃国产av成人99| 中文字幕人妻丝袜制服| 美女扒开内裤让男人捅视频| 国产精品亚洲一级av第二区| 亚洲中文av在线| 亚洲伊人久久精品综合| 激情在线观看视频在线高清 | av超薄肉色丝袜交足视频| 欧美另类亚洲清纯唯美| 国产精品久久久久久精品古装| 人人妻人人爽人人添夜夜欢视频| aaaaa片日本免费| 美国免费a级毛片| 国产精品麻豆人妻色哟哟久久| 天堂中文最新版在线下载| 亚洲综合色网址| 免费日韩欧美在线观看| 亚洲精品久久成人aⅴ小说| 一级片免费观看大全| 国产片内射在线| 亚洲熟女毛片儿| 捣出白浆h1v1| 久久久久国产一级毛片高清牌| 欧美激情久久久久久爽电影 | 啦啦啦中文免费视频观看日本| 中文字幕制服av| 久久中文字幕一级| 久久久国产一区二区| 一级黄色大片毛片| 日韩视频一区二区在线观看| 欧美国产精品一级二级三级| 欧美成人午夜精品| 90打野战视频偷拍视频| 正在播放国产对白刺激| 天天躁夜夜躁狠狠躁躁| 最近最新免费中文字幕在线| 亚洲精品av麻豆狂野| 99久久人妻综合| 少妇裸体淫交视频免费看高清 | 天天躁夜夜躁狠狠躁躁| 精品福利永久在线观看| 黑人猛操日本美女一级片| 午夜两性在线视频| 大香蕉久久网| 欧美另类亚洲清纯唯美| 91字幕亚洲| 少妇裸体淫交视频免费看高清 | 麻豆av在线久日| 国产一区二区 视频在线| 久久久久久人人人人人| cao死你这个sao货| 久久久久久久久免费视频了| 国产精品麻豆人妻色哟哟久久| 最近最新免费中文字幕在线| 99九九在线精品视频| 国产真人三级小视频在线观看| 最近最新免费中文字幕在线| 国产1区2区3区精品| 久久久国产欧美日韩av| 久久性视频一级片| 热99国产精品久久久久久7| 欧美在线黄色| 69精品国产乱码久久久| 一边摸一边抽搐一进一出视频| 国产精品国产av在线观看| 久久午夜亚洲精品久久| 中文字幕精品免费在线观看视频| 亚洲久久久国产精品| 久久久国产成人免费| 午夜福利影视在线免费观看| 欧美午夜高清在线| 免费在线观看黄色视频的| 99国产精品一区二区三区| 国产精品国产av在线观看| 成人手机av| 女人爽到高潮嗷嗷叫在线视频| 久久狼人影院| 操美女的视频在线观看| 国产精品99久久99久久久不卡| 日本精品一区二区三区蜜桃| 欧美另类亚洲清纯唯美| 好男人电影高清在线观看| 最近最新中文字幕大全免费视频| 国产99久久九九免费精品| 18禁观看日本| 视频在线观看一区二区三区| 亚洲欧美一区二区三区久久| 男女床上黄色一级片免费看| 熟女少妇亚洲综合色aaa.| 免费在线观看影片大全网站| 午夜福利在线观看吧| 中文字幕av电影在线播放| 天天影视国产精品| 99久久国产精品久久久| 日韩视频一区二区在线观看| 国产高清国产精品国产三级| videos熟女内射| 母亲3免费完整高清在线观看| 天堂8中文在线网| 涩涩av久久男人的天堂| 丝袜人妻中文字幕| 亚洲专区字幕在线| 亚洲熟妇熟女久久| 高清在线国产一区| 亚洲国产成人一精品久久久| 天堂中文最新版在线下载| 50天的宝宝边吃奶边哭怎么回事| 成年动漫av网址| 精品国产乱码久久久久久小说| 岛国在线观看网站| 欧美激情极品国产一区二区三区| 最近最新中文字幕大全电影3 | 亚洲全国av大片| 国产亚洲精品久久久久5区| 国产在线精品亚洲第一网站| 午夜福利欧美成人| 菩萨蛮人人尽说江南好唐韦庄| 国产精品1区2区在线观看. | 国产精品 国内视频| 国产免费视频播放在线视频| 国产免费av片在线观看野外av| 99久久国产精品久久久| 日本黄色视频三级网站网址 | av片东京热男人的天堂| 久久久久久免费高清国产稀缺| 黑丝袜美女国产一区| 成人影院久久| 啦啦啦免费观看视频1| 一区在线观看完整版| 成年版毛片免费区| 中文字幕人妻丝袜制服| 99国产精品99久久久久| a级片在线免费高清观看视频| 无人区码免费观看不卡 | 亚洲一区二区三区欧美精品| 少妇猛男粗大的猛烈进出视频| 久热这里只有精品99| 久久久精品国产亚洲av高清涩受| 国产视频一区二区在线看| 99精国产麻豆久久婷婷| 久久99热这里只频精品6学生| 国产精品麻豆人妻色哟哟久久| 欧美人与性动交α欧美精品济南到| 老司机在亚洲福利影院| 肉色欧美久久久久久久蜜桃| 国产成人精品久久二区二区免费| 久热这里只有精品99| 极品人妻少妇av视频| 在线观看免费视频日本深夜| 狠狠婷婷综合久久久久久88av| 国产成人av激情在线播放| 日韩视频在线欧美| 成人国语在线视频| 久久婷婷成人综合色麻豆| 男女下面插进去视频免费观看| 中文字幕高清在线视频| 人妻一区二区av| 亚洲第一欧美日韩一区二区三区 | 麻豆av在线久日| 精品久久久精品久久久| 国产精品九九99| 国产高清videossex| 女人高潮潮喷娇喘18禁视频| 久久久精品国产亚洲av高清涩受| 超碰97精品在线观看| 亚洲精品中文字幕一二三四区 | 精品免费久久久久久久清纯 | 亚洲 国产 在线| 男女边摸边吃奶| 精品久久久精品久久久| 老司机影院毛片| 精品第一国产精品| 欧美午夜高清在线| 亚洲欧洲日产国产| 男女无遮挡免费网站观看| 国产成人av教育| 国产日韩欧美视频二区| 交换朋友夫妻互换小说| 亚洲va日本ⅴa欧美va伊人久久| 色94色欧美一区二区| 最新美女视频免费是黄的| 超碰97精品在线观看| 亚洲精品在线观看二区| 久久久国产成人免费| 一级黄色大片毛片| 亚洲精品自拍成人| 成年版毛片免费区| cao死你这个sao货| 色婷婷久久久亚洲欧美| 中文字幕人妻丝袜制服| 国产99久久九九免费精品| 国产av一区二区精品久久| aaaaa片日本免费| 丝袜在线中文字幕| 制服诱惑二区| 精品一区二区三区视频在线观看免费 | 99久久精品国产亚洲精品| 麻豆乱淫一区二区| 国产精品久久久久久人妻精品电影 | 久久毛片免费看一区二区三区| 99热网站在线观看| 男女午夜视频在线观看| 一区二区av电影网| 老熟女久久久| 女警被强在线播放| 黄色 视频免费看| 久久人妻熟女aⅴ| 亚洲成人手机| 老司机影院毛片| 啦啦啦中文免费视频观看日本| 久久av网站| 成人18禁在线播放| 亚洲精品久久成人aⅴ小说| 日韩中文字幕视频在线看片| 丁香六月天网| 国产成人欧美在线观看 | a级片在线免费高清观看视频| 亚洲avbb在线观看| 男女无遮挡免费网站观看| 国产97色在线日韩免费| 少妇精品久久久久久久| 99久久精品国产亚洲精品| 狠狠精品人妻久久久久久综合| 亚洲色图 男人天堂 中文字幕| 午夜激情av网站| 男女免费视频国产| av不卡在线播放| 国产成人精品在线电影| 午夜91福利影院| 亚洲中文av在线| 黄色怎么调成土黄色| 一本一本久久a久久精品综合妖精| 伦理电影免费视频| 在线永久观看黄色视频| 在线播放国产精品三级| 欧美变态另类bdsm刘玥| 国产人伦9x9x在线观看| 国产欧美日韩一区二区三| svipshipincom国产片| 国产一区有黄有色的免费视频| 最新在线观看一区二区三区| 精品亚洲成a人片在线观看| bbb黄色大片| 久久婷婷成人综合色麻豆| 香蕉国产在线看| 亚洲欧美激情在线| 制服诱惑二区| 成人av一区二区三区在线看| 在线观看免费午夜福利视频| 纯流量卡能插随身wifi吗| 亚洲天堂av无毛| 精品国产超薄肉色丝袜足j| 一边摸一边抽搐一进一小说 | 老熟妇仑乱视频hdxx| 欧美乱妇无乱码| 亚洲三区欧美一区| 亚洲黑人精品在线| 国产精品99久久99久久久不卡| 亚洲一码二码三码区别大吗| 国产精品国产高清国产av | 美女扒开内裤让男人捅视频| 一区福利在线观看| 成人影院久久| 高清视频免费观看一区二区| av国产精品久久久久影院| 亚洲久久久国产精品| 午夜免费鲁丝| 在线天堂中文资源库| 丰满迷人的少妇在线观看| 国产高清激情床上av| av网站在线播放免费| 日韩精品免费视频一区二区三区| 国产成人啪精品午夜网站| 欧美黑人精品巨大| 久久久精品94久久精品| 一进一出好大好爽视频| 两个人看的免费小视频| 一本一本久久a久久精品综合妖精| 人人澡人人妻人| 国产一卡二卡三卡精品| 精品少妇久久久久久888优播| 美女高潮喷水抽搐中文字幕| 亚洲色图综合在线观看| videosex国产| 国产在视频线精品| 一夜夜www| 757午夜福利合集在线观看| 满18在线观看网站| 亚洲精品成人av观看孕妇| 亚洲一区二区三区欧美精品| 国产一区二区三区视频了| 两性夫妻黄色片| 亚洲精品在线美女| 亚洲一卡2卡3卡4卡5卡精品中文| 国产在线一区二区三区精| 极品人妻少妇av视频| 丁香六月天网| 老司机福利观看| 亚洲国产看品久久| 正在播放国产对白刺激| 国产欧美日韩精品亚洲av| 日本av手机在线免费观看| 黄色视频不卡| 国产成人精品久久二区二区免费| 天堂中文最新版在线下载| 女性被躁到高潮视频| 欧美在线一区亚洲| 一边摸一边做爽爽视频免费| 99香蕉大伊视频| 狠狠精品人妻久久久久久综合| 亚洲伊人久久精品综合| 国产精品麻豆人妻色哟哟久久| 2018国产大陆天天弄谢| 蜜桃在线观看..| 高潮久久久久久久久久久不卡| 老汉色av国产亚洲站长工具| 精品卡一卡二卡四卡免费| 大陆偷拍与自拍| 夜夜骑夜夜射夜夜干| 18禁美女被吸乳视频| 下体分泌物呈黄色| 欧美日韩成人在线一区二区| 久久久久久人人人人人| 两人在一起打扑克的视频| 国产在线观看jvid| 免费av中文字幕在线| 亚洲第一青青草原| 啦啦啦中文免费视频观看日本| 国产精品一区二区在线观看99| 精品久久久精品久久久| 一级片免费观看大全| 日本黄色视频三级网站网址 | 窝窝影院91人妻| 一个人免费在线观看的高清视频| 午夜日韩欧美国产| 狠狠婷婷综合久久久久久88av| 在线十欧美十亚洲十日本专区| av网站免费在线观看视频| 亚洲欧美日韩另类电影网站| 国产欧美日韩一区二区精品| 久久人人97超碰香蕉20202| 亚洲免费av在线视频| 啦啦啦免费观看视频1| 久久久久精品国产欧美久久久| 中文字幕制服av| 国产国语露脸激情在线看| 91成人精品电影| 亚洲第一欧美日韩一区二区三区 | 国产成人欧美在线观看 | 色精品久久人妻99蜜桃| 999精品在线视频| 久久这里只有精品19| 新久久久久国产一级毛片| 午夜激情久久久久久久| 91麻豆精品激情在线观看国产 | 在线播放国产精品三级| 一二三四社区在线视频社区8| 免费在线观看日本一区| 99在线人妻在线中文字幕 | 91成人精品电影| 亚洲avbb在线观看| 日本精品一区二区三区蜜桃| 咕卡用的链子| 久久亚洲真实| 蜜桃国产av成人99| 国产成人免费观看mmmm| 久久亚洲精品不卡| 黄色怎么调成土黄色| 99精国产麻豆久久婷婷| 欧美精品亚洲一区二区| 亚洲人成电影免费在线| 国产精品麻豆人妻色哟哟久久| 正在播放国产对白刺激| 美女高潮到喷水免费观看| 免费在线观看视频国产中文字幕亚洲| 亚洲一区二区三区欧美精品| 男人舔女人的私密视频| 色综合欧美亚洲国产小说| 久久狼人影院| 黄色a级毛片大全视频| 亚洲av电影在线进入| 欧美性长视频在线观看| 大香蕉久久成人网| 久久精品国产a三级三级三级| 色综合欧美亚洲国产小说| 18禁国产床啪视频网站| 美女福利国产在线| 天堂中文最新版在线下载| 亚洲熟女毛片儿| 无人区码免费观看不卡 | xxxhd国产人妻xxx| 91av网站免费观看| 国产日韩欧美视频二区| 精品久久久久久电影网| 国产主播在线观看一区二区| 丁香六月欧美| 满18在线观看网站| 飞空精品影院首页| 精品欧美一区二区三区在线| 欧美黄色片欧美黄色片| 日韩免费高清中文字幕av| av超薄肉色丝袜交足视频| 久久久久久久大尺度免费视频| 亚洲精品乱久久久久久| 99精品欧美一区二区三区四区| 色综合婷婷激情| av超薄肉色丝袜交足视频| 每晚都被弄得嗷嗷叫到高潮| 日韩中文字幕视频在线看片| 午夜福利,免费看| 咕卡用的链子| 久久久久国内视频| 十八禁高潮呻吟视频| 午夜久久久在线观看| 岛国毛片在线播放| 曰老女人黄片| 69精品国产乱码久久久| av在线播放免费不卡| 怎么达到女性高潮| 最近最新中文字幕大全免费视频| 悠悠久久av| av又黄又爽大尺度在线免费看| 欧美黑人欧美精品刺激| 一本久久精品| 青草久久国产| 涩涩av久久男人的天堂| 久久影院123| 老司机影院毛片| 一级a爱视频在线免费观看| 国产极品粉嫩免费观看在线| 国产又色又爽无遮挡免费看| 啦啦啦视频在线资源免费观看| 亚洲av欧美aⅴ国产| 亚洲国产成人一精品久久久| 法律面前人人平等表现在哪些方面| 99在线人妻在线中文字幕 | h视频一区二区三区| 免费在线观看日本一区| 最近最新中文字幕大全免费视频| 成人av一区二区三区在线看| 自拍欧美九色日韩亚洲蝌蚪91| 99精品欧美一区二区三区四区| 国产成+人综合+亚洲专区| 在线观看免费日韩欧美大片| 亚洲美女黄片视频| 高潮久久久久久久久久久不卡| 免费观看av网站的网址| 欧美日韩黄片免| 成人国语在线视频| 十八禁网站网址无遮挡| 在线 av 中文字幕| 99久久国产精品久久久| 国产免费av片在线观看野外av| netflix在线观看网站| 久久久久国内视频| 涩涩av久久男人的天堂| 我的亚洲天堂| 亚洲一码二码三码区别大吗| 狠狠精品人妻久久久久久综合| 99国产极品粉嫩在线观看| 成人18禁高潮啪啪吃奶动态图| 99精国产麻豆久久婷婷| 热re99久久精品国产66热6| 精品午夜福利视频在线观看一区 | 久久国产精品人妻蜜桃| 亚洲九九香蕉| 飞空精品影院首页| 国产激情久久老熟女| 亚洲人成77777在线视频| 国产精品亚洲一级av第二区| 一本色道久久久久久精品综合| 人妻一区二区av| 国产老妇伦熟女老妇高清| 欧美人与性动交α欧美精品济南到| 日韩精品免费视频一区二区三区| 91精品国产国语对白视频| 欧美精品人与动牲交sv欧美| 人人妻,人人澡人人爽秒播| 国产aⅴ精品一区二区三区波| 一夜夜www| 久久久久久久大尺度免费视频| 久久久精品94久久精品| 大陆偷拍与自拍| 老汉色∧v一级毛片| 国产男女内射视频| 黄色视频,在线免费观看| 麻豆乱淫一区二区| 精品午夜福利视频在线观看一区 | 精品第一国产精品| 亚洲欧美一区二区三区黑人| 午夜视频精品福利| 搡老乐熟女国产| 久久中文字幕一级| 久久久久久人人人人人| 高清在线国产一区| 性色av乱码一区二区三区2| 久久婷婷成人综合色麻豆| 免费久久久久久久精品成人欧美视频| 日韩视频一区二区在线观看| 热re99久久精品国产66热6| 欧美中文综合在线视频| 亚洲中文av在线|