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

    基于軸線馬赫數(shù)分布的噴管擴(kuò)張段無粘型面設(shè)計(jì)

    2016-07-05 12:53:35胡振震李震乾石義雷陳愛國(guó)
    實(shí)驗(yàn)流體力學(xué) 2016年4期
    關(guān)鍵詞:線網(wǎng)馬赫數(shù)軸線

    胡振震,李震乾,石義雷,陳愛國(guó)

    基于軸線馬赫數(shù)分布的噴管擴(kuò)張段無粘型面設(shè)計(jì)

    胡振震*,李震乾,石義雷,陳愛國(guó)

    (中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速空氣動(dòng)力研究所,四川綿陽 621000)

    針對(duì)高超聲速風(fēng)洞軸對(duì)稱噴管設(shè)計(jì)問題,開展了噴管擴(kuò)張段無粘型面設(shè)計(jì)研究。介紹了基于預(yù)設(shè)軸線馬赫數(shù)分布的直接設(shè)計(jì)方法,改進(jìn)了基于面積比的軸線馬赫數(shù)分布預(yù)設(shè)方法,提出了一種方便多點(diǎn)控制的軸線特征點(diǎn)分布方法。對(duì)設(shè)計(jì)噴管流場(chǎng)進(jìn)行特征線網(wǎng)三角化,與數(shù)值模擬結(jié)果進(jìn)行比較,并分析了影響噴管無粘型面的設(shè)計(jì)因素。表明:改進(jìn)的面積比方法可以保證軸線馬赫數(shù)分布預(yù)設(shè)的合理性;軸線馬赫數(shù)分布、軸線上特征點(diǎn)分布和邊界特征點(diǎn)數(shù)明顯影響噴管無粘型面。

    高超聲速風(fēng)洞;擴(kuò)張段無粘型面;直接設(shè)計(jì);軸線馬赫數(shù)分布

    0 引 言

    高超聲速風(fēng)洞噴管設(shè)計(jì)目前普遍采用無粘型面加邊界層修正的直接設(shè)計(jì)方法,相比于優(yōu)化設(shè)計(jì)[1-3],該方法原理簡(jiǎn)單計(jì)算便捷,在廣泛應(yīng)用中取得了成功。然而在實(shí)踐中有時(shí)也發(fā)現(xiàn)一個(gè)問題,即風(fēng)洞運(yùn)行的實(shí)際流場(chǎng)與噴管設(shè)計(jì)流場(chǎng)存在差異。該問題的產(chǎn)生因素有很多,包括來流非均勻非穩(wěn)態(tài)、非平衡膨脹、邊界層轉(zhuǎn)捩、熱化學(xué)非平衡、氣體組分冷凝、加工和組裝精度、噴管內(nèi)表面侵蝕和型面的力熱變型等,但噴管的無粘型面設(shè)計(jì)則是其中最重要的一個(gè)方面,因而分析其差異應(yīng)首先確認(rèn)噴管無粘型面設(shè)計(jì)方面的問題。

    無粘型面直接設(shè)計(jì)方法根據(jù)原理可分為3類:基于預(yù)設(shè)初始膨脹段型面的方法[4];基于泉流(或源流)假設(shè)的方法(如Foelsch[5]方法、Cresci[6]方法和Sivells[7-8]方法)和基于超聲速區(qū)預(yù)設(shè)軸線馬赫數(shù)分布的方法(如ASN[9]方法)。研究表明,第2類基于泉流假設(shè)的方法(特別是Sivells方法)在選擇合適的設(shè)計(jì)參數(shù)情況下可以得到良好的無粘噴管設(shè)計(jì),但該類方法設(shè)計(jì)的噴管長(zhǎng)度較長(zhǎng),調(diào)整不易。第1類和第3類方法均需計(jì)算完整的超聲速區(qū)特征線網(wǎng),區(qū)別在于第1類方法依據(jù)預(yù)設(shè)的初始膨脹段型面確定軸線上核心區(qū)的馬赫數(shù)分布,而第3類方法則根據(jù)預(yù)設(shè)的軸線馬赫數(shù)分布反求噴管型面。但在喉道附近利用特征線法從預(yù)設(shè)的噴管型面計(jì)算流場(chǎng),相比利用軸線上給定的速度分布反算能產(chǎn)生該預(yù)設(shè)速度分布的型面要困難。

    超聲速區(qū)軸線馬赫數(shù)分布的預(yù)設(shè)方法,最簡(jiǎn)單的是采用多項(xiàng)式分布描述,比如采用三次多項(xiàng)式[10]或五次多項(xiàng)式。多項(xiàng)式系數(shù)根據(jù)超聲速區(qū)起點(diǎn)和終點(diǎn)的馬赫數(shù)及其一階和二階導(dǎo)數(shù)條件確定,并可以通過附加項(xiàng)調(diào)整三次多項(xiàng)式,使得內(nèi)部分布可調(diào)。但隨著馬赫數(shù)升高或者噴管長(zhǎng)度增加,采用三次或五次多項(xiàng)式分布會(huì)使得軸線存在大于出口馬赫數(shù)的區(qū)域,導(dǎo)致軸線馬赫數(shù)分布出現(xiàn)問題。易仕和等人考慮采用Bezier曲線[9,11]和B樣條曲線[9,12]的軸線上速度分布方法直接避免該問題,而Tajfar等人[10]考慮利用泉流面積比和馬赫數(shù)的關(guān)系,首先由三次或五次多項(xiàng)式描述面積比分布,然后根據(jù)面積比反算馬赫數(shù)分布。采用了面積比的多項(xiàng)式分布后,軸線馬赫數(shù)不會(huì)超過設(shè)計(jì)出口馬赫數(shù),但仍可能產(chǎn)生分布不合理的現(xiàn)象,特別是在起點(diǎn)附近可能馬赫數(shù)增長(zhǎng)過快,進(jìn)而導(dǎo)致計(jì)算特征線網(wǎng)時(shí)特征線相交。為此本文在Tajfar方法的基礎(chǔ)上做出改進(jìn),解決局部速度分布不合理的問題。

    噴管直接設(shè)計(jì)方法應(yīng)用的基本技術(shù)是特征線網(wǎng)計(jì)算[4,13-14]和型面點(diǎn)確定技術(shù)。型面點(diǎn)的確定可采用流函數(shù)法和流量積分法,本文采用Moger和Ramsay[14]提出的拋物擬合的流量積分法。特征線網(wǎng)計(jì)算過程中分從前往后和從后往前2種推進(jìn)方法,在已知邊界情況下可利用左行特征線一條一條的往上游推進(jìn)或利用右行特征線往下游推進(jìn)(見圖1),本文采用向下游推進(jìn)的方法。為研究無粘型面設(shè)計(jì)和特征線網(wǎng)質(zhì)量的影響因素,本文將對(duì)設(shè)計(jì)噴管的特征線網(wǎng)進(jìn)行Delaunay三角化[15],并與數(shù)值模擬流場(chǎng)結(jié)果進(jìn)行比較分析。

    圖1 特征線推進(jìn)方法(a)往上游推進(jìn)(b)往下游推進(jìn)Fig.1 Characteristic marching method(a)march upstream(b)march downstream

    1 基于預(yù)設(shè)軸線馬赫數(shù)分布的直接設(shè)計(jì)方法

    1.1軸線馬赫數(shù)分布及面積比方法的改進(jìn)

    軸線上的馬赫數(shù)分布預(yù)設(shè)可分為分段預(yù)設(shè)和統(tǒng)一預(yù)設(shè)2種方法。Sivells方法在本文作者看來,可以理解為是一種分段預(yù)設(shè)方法。該方法將軸線分成3個(gè)部分,從聲速點(diǎn)I到點(diǎn)E是泉流開始的區(qū)域,速度分布用四次多項(xiàng)式描述,從E點(diǎn)到D點(diǎn)為泉流區(qū),從D點(diǎn)到F點(diǎn)是均勻流開始的區(qū)域,馬赫數(shù)分布用一個(gè)五次多項(xiàng)式描述(見圖2)。軸線分成3段,但結(jié)合起來就是一種軸線馬赫數(shù)分布,若給出喉部的一條初始線AJ(通過喉部跨聲速近似解[16-17]確定一條初值線并由MOC求得),再結(jié)合出口馬赫線FC,即可通過MOC計(jì)算噴管擴(kuò)張段完整特征線網(wǎng)(見圖3)進(jìn)而確定噴管型面,而無需分段計(jì)算特征線網(wǎng)(見圖2)。

    圖2 Sivells方法確定噴管型面Fig.2 Sivells method to determine the nozzle contour

    圖3 基于Sivells方法預(yù)設(shè)軸線馬赫數(shù)分布的特征線網(wǎng)Fig.3 Characteristic network based on the Sivells’axial Mach number distribution

    軸線上馬赫數(shù)分布的統(tǒng)一預(yù)設(shè)可采用多種方法,比如三次多項(xiàng)式、五次多項(xiàng)式、Bezier曲線、B樣條曲線,以及面積比的三次多項(xiàng)式等。預(yù)設(shè)馬赫數(shù)需滿足條件:J點(diǎn):Ma(xJ)=MaJ,Ma′(xJ)=Ma′J,Ma″(xJ)=Ma″J,F(xiàn)點(diǎn):Ma(xF)=MaF,Ma′(xF)=0,Ma″(xF)=0。J點(diǎn)信息在已知喉道曲率半徑和轉(zhuǎn)折角情況下通過跨聲速近似解或CFD解獲得,一階和二階導(dǎo)數(shù)通過差分方法計(jì)算。而F點(diǎn)的位置為:

    其中cF≥1為調(diào)整噴管長(zhǎng)度的系數(shù)。

    三次多項(xiàng)式分布為:

    改進(jìn)的三次多項(xiàng)式分布可通過a4系數(shù)進(jìn)行內(nèi)部微調(diào):

    五次多項(xiàng)式分布為:

    Bezier曲線和B樣條曲線分布公式相對(duì)復(fù)雜,2種曲線的定義、構(gòu)造和插值可參考文獻(xiàn)[18-19]。這里給出p次B樣條曲線的應(yīng)用。對(duì)于軸線上馬赫數(shù)分布問題,B樣條曲線是u(u∈[0,1])的函數(shù),由n+1個(gè)p階基函數(shù)構(gòu)成,基函數(shù)需要m=n+p+1個(gè)節(jié)點(diǎn)構(gòu)造,節(jié)點(diǎn)構(gòu)成的節(jié)點(diǎn)矢量ui(i=0,…,m)采用平均分布,控制點(diǎn)Pi(i=0,…,n)為(x,Ma)坐標(biāo),B樣條曲線描述的馬赫數(shù)分布由參數(shù)u的方程表示:

    控制點(diǎn)P0為:(x0,Ma0)=(xJ,MaJ)

    控制點(diǎn)P1由起點(diǎn)的一階導(dǎo)數(shù)確定:

    其中dL是可調(diào)參數(shù),常取控制點(diǎn)構(gòu)成的總弦長(zhǎng),當(dāng)?shù)玫降目刂泣c(diǎn)不合理時(shí),可適當(dāng)調(diào)整其大小。

    控制點(diǎn)P2由起點(diǎn)的二階導(dǎo)數(shù)確定:

    控制點(diǎn)Pn為:(xn,Man)=(xF,Maexit)

    控制點(diǎn)Pn-1由終點(diǎn)的一階導(dǎo)數(shù)確定:

    控制點(diǎn)Pn-2由終點(diǎn)的二階導(dǎo)數(shù)確定:

    上述6個(gè)控制點(diǎn)已經(jīng)足夠構(gòu)成B樣條曲線的馬赫數(shù)分布,若仍然需要進(jìn)行內(nèi)部微調(diào),可以在控制點(diǎn)2~(n-2)之間增加控制點(diǎn)。

    為解決軸線上局部區(qū)域馬赫數(shù)可能大于出口馬赫數(shù)的問題,Tajfar等人考慮利用泉流區(qū)面積比和馬赫數(shù)的關(guān)系,首先由三次多項(xiàng)式描述面積比分布,然后根據(jù)面積比反算馬赫數(shù)分布。泉流面積比Ar與馬赫數(shù)的關(guān)系為:

    而面積比由三次多項(xiàng)式表示

    其中J點(diǎn)面積及其導(dǎo)數(shù)信息可以由(10)式及其導(dǎo)數(shù)與J點(diǎn)馬赫數(shù)及其導(dǎo)數(shù)計(jì)算得到:

    圖4~5給出了噴管出口馬赫數(shù)為12,cF取3.0和5.57時(shí),不同預(yù)設(shè)方法的馬赫數(shù)分布比較。噴管喉道曲率半徑比等設(shè)計(jì)參數(shù)與Sivells方法選取的參數(shù)一致。cF=5.57是噴管長(zhǎng)度與Sivells噴管長(zhǎng)度一致時(shí)的取值。

    圖4 不同預(yù)設(shè)方法的軸線馬赫數(shù)分布(cF=3.0時(shí))Fig.4 Axial Mach number of different presetting methods with cF=3.0

    隨著cF取值增大,噴管長(zhǎng)度增加,三次或五次多項(xiàng)式描述的馬赫數(shù)會(huì)明顯超過出口馬赫數(shù)。圖6表明采用三次多項(xiàng)式分布,cF=4.0,軸線上特征點(diǎn)采用雙曲加密(兩端系數(shù)分別為0.99和0.50)時(shí)特征線出現(xiàn)相交。實(shí)踐表明:當(dāng)軸線上的馬赫數(shù)略超過出口馬赫數(shù)時(shí),特征線并不一定相交,但隨著超出程度增大,特征線出現(xiàn)相交。

    圖5 不同預(yù)設(shè)方法的軸線馬赫數(shù)分布(cF=5.57時(shí))Fig.5 Axial Mach number of different presetting methods with cF=5.57

    圖6 軸線上馬赫數(shù)大于出口馬赫數(shù)時(shí)的特征線相交Fig.6 Characteristic intersection when Mach number on axis larger then Maexit

    事實(shí)上,采用面積比的三次分布方法后,軸線馬赫數(shù)不會(huì)超過出口馬赫數(shù),但仍可能存在不合理的地方,特別是當(dāng)cF較小時(shí),采用面積比的三次多項(xiàng)式分布,其起點(diǎn)J附近馬赫數(shù)增長(zhǎng)明顯快過B樣條曲線和Sivells方法的結(jié)果(見圖4),而因?yàn)轳R赫數(shù)增長(zhǎng)過快,很可能導(dǎo)致計(jì)算特征線網(wǎng)時(shí)特征線相交。圖7給出了面積比三次多項(xiàng)式分布,cF=2.9,軸線點(diǎn)加密同圖6時(shí),起點(diǎn)J附近特征線網(wǎng)相交的情況。

    為了避免在起點(diǎn)附近出現(xiàn)速度增長(zhǎng)過快問題,本文在(10)式的基礎(chǔ)上做出改進(jìn),增加一個(gè)指數(shù)系數(shù)idx使起點(diǎn)附近馬赫數(shù)的增長(zhǎng)速度可調(diào),如式(13)所示:

    盡管式中增加了一個(gè)系數(shù)但利用相關(guān)函數(shù)控制馬赫數(shù)分布的本質(zhì)沒有變。當(dāng)idx取1.0時(shí),就是原面積比的方法,而當(dāng)需要減小馬赫數(shù)增長(zhǎng)速度時(shí),僅需調(diào)整系數(shù)idx取值小于1.0,反之亦然。圖4表明idx系數(shù)為0.8的面積比預(yù)設(shè)方法在起點(diǎn)附近馬赫數(shù)增長(zhǎng)速度明顯低于idx取1.0的情況。總的來說,面積比方法的改進(jìn),既滿足了超聲速區(qū)起點(diǎn)和終點(diǎn)馬赫數(shù)及其導(dǎo)數(shù)的連續(xù),又保證了軸線上的馬赫數(shù)不超過出口馬赫數(shù),避免過膨脹導(dǎo)致的冷凝問題,也保證了超聲速區(qū)起點(diǎn)和終點(diǎn)附近的馬赫數(shù)增長(zhǎng)的合理性。表1總結(jié)了幾種軸線馬赫數(shù)分布方法的特點(diǎn),表明改進(jìn)的面積比方法可以適用于任意長(zhǎng)度噴管設(shè)計(jì),且比采用B樣條曲線計(jì)算量小,也更易于理解。

    圖7 馬赫數(shù)增長(zhǎng)過快導(dǎo)致的特征線相交(J點(diǎn)附近)Fig.7 Characteristic intersection caused by too fast Mach number increase(near J point)

    表1 不同馬赫數(shù)分布方法比較Table 1 Comparison of different Mach number distribution methods

    1.2軸線特征點(diǎn)分布

    軸線上特征點(diǎn)的分布可采用指數(shù)函數(shù)加密[19]:

    其中p=1時(shí)為平均分布,p>1時(shí)為指數(shù)加密,Npt為特征點(diǎn)總數(shù),x1,x2為起點(diǎn)終點(diǎn)坐標(biāo)。也可以采用雙曲函數(shù)加密:

    其中b1,b2為加密參數(shù),b1控制首端,b2控制未端,取值范圍[0~1),取值越接近1,加密程度越高。

    上述2種加密方式,點(diǎn)的分布調(diào)整相對(duì)單調(diào)。文本提出一種多點(diǎn)控制的方法以實(shí)現(xiàn)更大的自由度。首先任意給定需要加密的位置(即控制點(diǎn)位置),并根據(jù)控制參數(shù)(統(tǒng)一給單側(cè)雙曲加密參數(shù),如(15)式中的b1并取b2=0)統(tǒng)一確定各位置的網(wǎng)格間距(該間距為第一個(gè)網(wǎng)格的間距),每個(gè)控制段(2個(gè)相鄰控制點(diǎn)之間)的各網(wǎng)格間距采用線性分布,則各控制段的所有網(wǎng)格確定,進(jìn)而確定各段的網(wǎng)格數(shù)。各控制段的網(wǎng)格數(shù)在總網(wǎng)格數(shù)的比例乘以實(shí)際的網(wǎng)格數(shù)就是最終各段的實(shí)際網(wǎng)格數(shù)。然后根據(jù)網(wǎng)格間距線性分布規(guī)律,自動(dòng)調(diào)整各段兩端的網(wǎng)格間距從而確定所有布點(diǎn)。

    考慮在[0,1)線段上布點(diǎn),特征點(diǎn)數(shù)為151,實(shí)際網(wǎng)格數(shù)為150,考慮控制點(diǎn)(0,0.05,0.15,0.6,1.0),控制參數(shù)(0.9999,0.99,0.95,0.90,0.1),計(jì)算過程如下:

    (1)由式(15)確定(式中x1=0,x2=1,b2=0)各控制點(diǎn)網(wǎng)格間距為:(0.6779×10-5,0.3585×10-3,0.1259×10-2,0.2076×10-2,0.6578×10-2)

    (2)各控制段的末端網(wǎng)格間距與首端比值rdi為:(52.8787,3.51314,1.64866,3.16826)

    (3)由(16)式計(jì)算各段應(yīng)布置網(wǎng)格數(shù)為:(273,123,269,92)

    其中i為第i個(gè)控制點(diǎn),Ni為控制點(diǎn)i~i+1段網(wǎng)格數(shù),di為i點(diǎn)的網(wǎng)格間距,rdi=di+1/di。

    (4)確定各段實(shí)際布置點(diǎn)數(shù)為:(54,24,53,19)

    (5)由(16)式確定各段實(shí)際的首端網(wǎng)格間距為:(0.3437×10-4,0.1846×10-2,0.6411×10-2,0.1010 ×10-1),末端網(wǎng)格間距和首端網(wǎng)格間距的比值不變,進(jìn)而根據(jù)網(wǎng)格間距線性分布確定所有的點(diǎn)的位置。圖8給出了上述3種不同布點(diǎn)方法的結(jié)果。

    2 馬赫12噴管的數(shù)值模擬分析

    噴管無粘型面的設(shè)計(jì)實(shí)踐表明,隨著噴管設(shè)計(jì)出口馬赫數(shù)的提高,噴管長(zhǎng)度增加,誤差會(huì)更為顯著,因此本文考慮設(shè)計(jì)馬赫12噴管,對(duì)無粘型面設(shè)計(jì)及其影響因素進(jìn)行分析。

    2.1軸線馬赫數(shù)分布的影響

    噴管設(shè)計(jì)考慮特征點(diǎn)數(shù):取AJ段81,JF段181,F(xiàn)C段251,轉(zhuǎn)折角12°,喉道曲率半徑比12。首先考慮與Sivells噴管長(zhǎng)度一致(cF取5.57),比較4種馬赫數(shù)分布:(Case1)Sivells分布,d Ma為D點(diǎn)馬赫數(shù)與D點(diǎn)最小可取值之差和F點(diǎn)馬赫數(shù)與該最小值之差的比值,取0.8;(Case2)面積比的三次多項(xiàng)式分布;(Case3)面積比的五次多項(xiàng)式分布;(Case4)B樣條分布,dL取0.2·(xF-xJ)。

    圖8 不同方法的布點(diǎn)比較Fig.8 Point distribution of different methods

    然后考慮稍短的噴管(cF取3.0),比較4種馬赫數(shù)分布的結(jié)果:(Case5)馬赫數(shù)的三次多項(xiàng)式;(Case6)面積比的三次多項(xiàng)式,idx=0.8;(Case7)面積比的五次多項(xiàng)式;(Case8)B樣條分布,dL取0.3 ·(xF-xJ)。

    圖9~12給出了Case1~4的噴管設(shè)計(jì)流場(chǎng)馬赫數(shù)等值線(特征線網(wǎng)三角化后的結(jié)果,紅虛線)與數(shù)值模擬結(jié)果(黑實(shí)線)的比較,4種不同軸線馬赫數(shù)分布的設(shè)計(jì)流場(chǎng)與數(shù)值模擬結(jié)果基本一致。實(shí)際上特征線方法計(jì)算的出口參數(shù)就是均勻的設(shè)計(jì)出口馬赫數(shù),只要數(shù)值過程足夠精確,理論上特征線方法得到的型面應(yīng)該就是能夠產(chǎn)生均勻出口馬赫數(shù)的型面。但是實(shí)際計(jì)算中,布點(diǎn)數(shù)、網(wǎng)格和數(shù)值過程等都帶來誤差,所以實(shí)際得到的噴管出口必然有一定的偏差,而這個(gè)偏差需要用數(shù)值模擬方法獲得。圖9~12的結(jié)果表明特征線方法實(shí)際計(jì)算中的偏差并不明顯,至少在云圖上無法明確得出,但通過分析數(shù)值模擬的出口馬赫數(shù)可以比較在相同設(shè)計(jì)參數(shù)和計(jì)算條件下,不同分布所帶來的差異。圖13給出了Case1~4噴管的出口馬赫數(shù)比較,其中面積比五次多項(xiàng)式分布和Sivells分布的最大偏差接近0.1%,B樣條曲線分布的結(jié)果靠近軸線局部區(qū)域略大于0.1%,而面積比三次多項(xiàng)式分布的結(jié)果最大偏差為0.35%。圖14給出了Case5~8噴管的結(jié)果,其中因?yàn)槊娣e比的三次多項(xiàng)式分布導(dǎo)致起點(diǎn)J附近的馬赫數(shù)增長(zhǎng)過快,所以取idx=0.8。除靠近壁面的區(qū)域外4種分布出口馬赫數(shù)偏差均小于0.1%,壁面附近偏差最大的是面積比的五次分布,約為0.3%,最小的是B樣條分布,但4種分布之間的差異小于0.03%。上述結(jié)果說明,在噴管長(zhǎng)度、特征點(diǎn)分布和數(shù)量等設(shè)計(jì)參數(shù)一致時(shí),采用不同的軸線馬赫數(shù)分布函數(shù)對(duì)噴管出口馬赫數(shù)存在影響,但只要采用合適的分布使得軸線馬赫數(shù)正確合理,可以在除壁面附近外的區(qū)域得到出口馬赫數(shù)偏差小于0.1%的結(jié)果。

    圖9 基于Sivells軸線馬赫數(shù)分布預(yù)設(shè)的噴管馬赫數(shù)等值線Fig.9 Mach number contour of the nozzle designed with the Sivells’axial Mach number distribution

    圖10 基于面積比三次多項(xiàng)式分布預(yù)設(shè)的噴管馬赫數(shù)等值線Fig.10 Mach number contour of the nozzle designed with the cubic distribution of the ratio of area

    圖11 基于面積比五次多項(xiàng)式分布預(yù)設(shè)的噴管馬赫數(shù)等值線Fig.11 Mach number contour of the nozzle designed with the quintic distribution of the ratio of area

    圖12 基于B樣條曲線馬赫數(shù)預(yù)設(shè)的噴管馬赫數(shù)等值線Fig.12 Mach number contour of the nozzle designed with the B-spline Mach number distribution

    2.2軸線特征點(diǎn)的位置分布及其數(shù)量的影響

    軸線上特征點(diǎn)的位置分布會(huì)明顯影響特征線網(wǎng),比如影響網(wǎng)格的過渡光滑性和正交性。圖15給出了馬赫12噴管,軸線上馬赫數(shù)三次多項(xiàng)式分布,軸線上點(diǎn)數(shù)101,出口馬赫數(shù)線上點(diǎn)數(shù)101,軸線上的點(diǎn)采用雙曲加密,兩端加密參數(shù)分別為(a)0.999,0.5;(b)0.99,0.5;(c)0.9,0.5的特征線網(wǎng),其中圖5(a)和5(b)圖在特征線網(wǎng)起始區(qū)域(J點(diǎn)附近)過渡明顯好于

    圖13 不同軸線馬赫數(shù)預(yù)設(shè)方法的噴管出口馬赫數(shù)Fig.13 Mach number at exit of the nozzle designed with different axial Mach number presetting methods

    圖14 不同軸線馬赫數(shù)預(yù)設(shè)方法的噴管出口馬赫數(shù)cF=3.0Fig.14 Mach number at exit of the nozzle designed with different axial Mach number presetting methods with cF=3.0

    圖15 軸線上不同布點(diǎn)的特征線網(wǎng)和噴管型面比較Fig.15 Comparision of characteristic networks and contours of different point distributions

    圖5 (c)。圖5(d)圖中噴管無粘型面的差異反映了圖5(a)、5(b)和5(c)不同特征線網(wǎng)所導(dǎo)致的差異,表明特征線網(wǎng)是否過渡光滑是其質(zhì)量的重要影響因素。

    特征點(diǎn)的數(shù)量既影響特征線網(wǎng)疏密也影響其是否光滑過渡。圖16給出了采用4種不同布點(diǎn)數(shù)時(shí)在F點(diǎn)附近的特征線網(wǎng),軸線上兩端加密參數(shù)為0.999和0.5,其中軸線和出口馬赫線分別為(a)101,101;(b)101,166;(c)151,251;(d)201,333。在圖16(a)中F點(diǎn)附近軸線和出口馬赫線上的右行特征線過渡不夠光滑,而圖16(b)、16(c)和16(d)相對(duì)較好。隨著點(diǎn)數(shù)的增加特征線網(wǎng)也隨之變密。圖17給出了噴管某位置的型面比較,可以看到采用不同點(diǎn)數(shù)時(shí)的細(xì)微差異。

    圖16 不同點(diǎn)數(shù)的特征線網(wǎng)(F點(diǎn)附近)Fig.16 Comparision of characteristic networks of different point numbers(near Fpoint)

    圖17 不同點(diǎn)數(shù)的型面比較Fig.17 Comparision of contour of different point numbers

    采用數(shù)值模擬方法考察軸線和出口馬赫線上點(diǎn)的數(shù)量和分布對(duì)于噴管流場(chǎng)的影響,考慮馬赫12噴管,三次多項(xiàng)式軸線馬赫數(shù)分布cF=3.0,初始線特征點(diǎn)數(shù)81,喉道曲率半徑比12,轉(zhuǎn)折角12°,泉流半徑為1.0。計(jì)算如下狀態(tài):(Case1)軸線101雙曲函數(shù)分布b1=0.999,b2=0.5,出口馬赫線166,平均分布;(Case2)軸線151雙曲函數(shù)分布b1=0.999,b2=0.5,出口馬赫線251,平均分布;(Case3)軸線201雙曲函數(shù)分布b1=0.999,b2=0.5,出口馬赫線333,平均分布;(Case4)軸線151雙曲函數(shù)分布b1=0.99,b2=0.5,出口馬赫線251,平均分布;(Case5)軸線151多點(diǎn)控制分布(見1.2節(jié)),出口馬赫線251,平均分布。

    圖18的噴管出口馬赫數(shù)表明,5個(gè)狀態(tài)的結(jié)果在除靠近上壁面的區(qū)域外,偏差均在0.1%以內(nèi),不同狀態(tài)在不同的y位置偏差各異。從圖19噴管上壁面的馬赫數(shù)看,狀態(tài)3最接近出口馬赫數(shù)12,偏差約0.19%,而狀態(tài)1的結(jié)果偏離最大達(dá)0.45%,說明當(dāng)網(wǎng)格點(diǎn)數(shù)偏少時(shí),計(jì)算結(jié)果偏離明顯增大,隨著網(wǎng)格點(diǎn)數(shù)增加,特征線網(wǎng)的精度明顯增加。而狀態(tài)2、4和5在相同點(diǎn)數(shù)情況下,軸線上特征點(diǎn)的不同分布會(huì)

    圖18 噴管出口馬赫數(shù)Fig.18 Mach number at exit of the nozzle

    圖19 噴管上壁面馬赫數(shù)Fig.19 Mach number at wall of the nozzle

    產(chǎn)生一定的差異,其中狀態(tài)4的結(jié)果最接近出口馬赫數(shù),但3種結(jié)果偏差小于0.02%,結(jié)合圖15分析,首端加密參數(shù)b1大于0.99,使得特征線網(wǎng)過渡光滑,是特征線網(wǎng)質(zhì)量的重要保證。

    3 結(jié) 論

    對(duì)于高超聲速風(fēng)洞噴管的無粘型面設(shè)計(jì),基于預(yù)設(shè)軸線馬赫數(shù)分布的直接設(shè)計(jì)方法是目前最普遍采用的方法。通過本文的研究得出如下結(jié)論:

    (1)對(duì)基于面積比的軸線馬赫數(shù)分布設(shè)定的方法進(jìn)行改進(jìn),使得在原方法保證軸線馬赫數(shù)不超過出口馬赫數(shù)的同時(shí),提高了軸線上起點(diǎn)附近馬赫數(shù)增長(zhǎng)的合理性。采用多點(diǎn)控制的軸線特征點(diǎn)位置分布方法,提供了軸線特征點(diǎn)分布更大的自由度和可控度。

    (2)軸線上馬赫數(shù)分布函數(shù)、特征點(diǎn)的分布和邊界上特征點(diǎn)的數(shù)量影響特征線網(wǎng)的質(zhì)量,進(jìn)而影響噴管無粘型面。當(dāng)采用某種分布函數(shù)保證軸線上馬赫數(shù)分布合理,則可以保證出口馬赫數(shù)除壁面附近外的絕大部分區(qū)域偏差小于0.1%。特征點(diǎn)數(shù)增加或者采用一定的特征點(diǎn)分布,使得特征線網(wǎng)過渡更光滑,能進(jìn)一步減小壁面附近的出口馬赫數(shù)偏差。

    [1]Shope F L.Contour design techniques for super/hypersonic wind tunnel nozzles[C].24th AIAA Applied Aerodynamics Conference,San Francisco,CA,2006.

    [2]Korte J J,Kumar A,Singh D J,et al.CAN-DO,CFD-based aerodynamic nozzle design optimization program for supersonic/hypersonic wind tunnels[C].AIAA 17thAerospace Groud Testing Conference,Nashville,TN,1992.

    [3]Shope F L.Design optimization of hypersonic test facility nozzle contours using splined corrections[R].AEDC-TR-04-2,2005.

    [4]Zucrow M J,Hoffman J D.Gas dynamics volume I~I(xiàn)I[M].New York:John Wiley &Sons Inc,1976.

    [5]Foelsch K.The analytical design of an axially symmetriclaval nozzle for a parallel and uniform jet[J].Journal of the Aeronautical Sciences,1949,16:161-166.

    [6]Cresci R J.Tabulation of coordinates for hypersonic axisymmetric nozzle Part I-Analysis and coordinates for test section Mach number of 8,12and 20[R].TN58-300,1958.

    [7]Sivells J C.Aerodynamic design of axisymmetric hypersonic wind-tunnel nozzles[J].Journal of Spacecraft,1970,7(11):1292-1299.

    [8]Sivells J C.A computer program for the aerodynamic design of axisymmetric and planar nozzles for supersonic and hypersonic wind tunnels[R].AEDC-TR-78-63,1978.

    [9]易仕和,趙玉新,何霖,等.超聲速和高超聲速噴管設(shè)計(jì)[M].北京:國(guó)防工業(yè)出版社,2013.Yi S H,Zhao Y X,He L,et al.Supersonic and hypersonic nozzle design[M].Beijing:National Defense Industry Press,2013.

    [10]Tajfar A H,Hall I M.Design of a nozzle for a hypersonic wind tunnel[M].AERO-REPT-9113,ETN-92-92779,1991.

    [11]張敏莉,易仕和,趙玉新.超聲速短化噴管的設(shè)計(jì)和試驗(yàn)研究[J].空氣動(dòng)力學(xué)學(xué)報(bào),2007,25(4):500-503.Zhang M L,Yi S H,Zhao Y X.The design and experimental investigations of supersonic length shorted nozzle[J].Acta Aerodynamica Sinica,2007,25(4):500-503.

    [12]趙一龍,趙玉新,王振國(guó),等.超聲速型面可控噴管設(shè)計(jì)方法[J].國(guó)防科技大學(xué)學(xué)報(bào),2012,34(5):1-4.Zhao Y L,Zhao Y X,Wang Z G,et al.Designing method of supersonic nozzle with controllable contour[J].Journal of National University of Defense Technology,2012,34(5):1-4.

    [13]Shapiro A H.The dynamics and thermodynamics of compressible fluid flow[M].New York:The Ronald Press Company,1953.

    [14]Moger W C,Ramsay D B.Supersonic axisymmetric nozzle design by mass flow techniques utilizing a digital computer[R].AEDC-TDR-64-110,U.S.Air Force,1964.

    [15]王永.非結(jié)構(gòu)化網(wǎng)格生成技術(shù)及在SIMPLE算法中的應(yīng)用研究[D].天津:天津大學(xué),2005.Wang Y.The study on generation techniques of unstructured grid and application in SIMPIE algorithm[D].Tianjin:Tianjin University,2005.

    [16]Hall I M.Transonic flow in two-dimensional and axially-symmetric nozzles[J].Quarterly Journal of Mechanics and Applied Mathematics.XV,1962,15(4):487-508.

    [17]Kligel J R,Levine J N.Transonic flow in small throat radius of curvature nozzles[J].AIAA,1969,7(7):1375-1378.

    [18]施法中.計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理B樣條(CAGD&NURBS)[M].北京:北京航空航天大學(xué)出版社,1994.

    [19]Piegl L,Tiller W.The NURBS book[M].New York:Springer,1997.

    Design of nozzle inviscid contour based on axial Mach number distribution

    Hu Zhenzhen*,Li Zhenqian,Shi Yilei,Chen Aiguo(Hypersonic Aerodynamics Research Institute,China Aerodynamics Research and Development Center,Mianyang Sichuan 621000,China)

    The supersonic inviscid contour design methods were studied for the hypersonic wind tunnel axial-symmetric nozzle.The Direct-design technique based on the axial Mach number distribution presetting was introduced,the axial Mach number distribution presetting method based on the ratio of area was improved,and a new multipoint-control characteristic point distribution method was proposed.The designed nozzle characteristic network was triangulated and compared with the numerical simulation result,and influencing factors were analyzed for the nozzle inviscid contour.Results indicate that:the feasibility of the axial Mach number distribution presetting method based on the ratio of area is guaranteed by the present improvement;the nozzle inviscid contour is significantly affected by the axial Mach number distribution,axial characteristic point distribution and the amount of characteristic points on the boundary.

    hypersonic wind tunnel;supersonic inviscid contour;direct-design;axial Mach number distribution

    V211.74

    :A

    (編輯:楊 娟)

    1672-9897(2016)04-0097-08

    10.11729/syltlx20150115

    2015-09-04;

    2015-12-04

    *通信作者E-mail:hzzmail@163.com

    Hu Z Z,Li Z Q,Shi Y L,et al.Design of nozzle inviscid contour based on axial Mach number distribution.Journal of Experiments in Fluid Mechanics,2016,30(4):97-104.胡振震,李震乾,石義雷,等.基于軸線馬赫數(shù)分布的噴管擴(kuò)張段無粘型面設(shè)計(jì).實(shí)驗(yàn)流體力學(xué),2016,30(4):97-104.

    胡振震(1984-),男,浙江武義人,助理研究員。研究方向:高超聲速試驗(yàn)技術(shù)。通信地址:四川省綿陽市二環(huán)路南段6號(hào)15信箱506分箱(621000)。E-mail:hzzmail@163.com

    猜你喜歡
    線網(wǎng)馬赫數(shù)軸線
    一維非等熵可壓縮微極流體的低馬赫數(shù)極限
    曲軸線工件劃傷問題改進(jìn)研究
    載荷分布對(duì)可控?cái)U(kuò)散葉型性能的影響
    新型線網(wǎng)城軌乘客信息系統(tǒng)的研究與分析
    軌道交通COCC線網(wǎng)信號(hào)系統(tǒng)設(shè)計(jì)
    基于回歸分析的水電機(jī)組軸線曲折預(yù)判斷分析
    行書章法淺析(十五)書寫應(yīng)把握行軸線
    凸輪軸孔軸線與止推面垂直度超差問題研究
    河南科技(2014年16期)2014-02-27 14:13:21
    緊湊型大都市區(qū)軌道線網(wǎng)形態(tài)配置研究
    自動(dòng)售檢票線網(wǎng)化維修管理系統(tǒng)的構(gòu)建
    国产亚洲一区二区精品| 久久久a久久爽久久v久久| 亚洲av中文字字幕乱码综合| 国产精品一区二区在线观看99| 自拍偷自拍亚洲精品老妇| 久久久国产一区二区| 97热精品久久久久久| 久久精品久久久久久久性| 久热这里只有精品99| 午夜福利在线观看免费完整高清在| 少妇被粗大猛烈的视频| 一二三四中文在线观看免费高清| 久久久久久久大尺度免费视频| 熟女av电影| 国产爱豆传媒在线观看| 久久久精品94久久精品| 狂野欧美激情性bbbbbb| 免费看不卡的av| 身体一侧抽搐| 国产亚洲欧美精品永久| 亚洲欧美中文字幕日韩二区| 国产欧美日韩精品一区二区| 久久精品人妻少妇| 国产精品爽爽va在线观看网站| a 毛片基地| 成人美女网站在线观看视频| 纵有疾风起免费观看全集完整版| 国产国拍精品亚洲av在线观看| 春色校园在线视频观看| 亚洲人成网站高清观看| 国产av精品麻豆| 少妇人妻 视频| 久久青草综合色| 日本黄色片子视频| 亚洲国产日韩一区二区| 亚洲欧美成人综合另类久久久| 久久久久久久精品精品| 国产亚洲一区二区精品| 看免费成人av毛片| 99久久精品热视频| 大片免费播放器 马上看| 亚洲精品国产av蜜桃| 日韩,欧美,国产一区二区三区| 亚洲欧美日韩无卡精品| 天堂8中文在线网| 久久99热6这里只有精品| 少妇的逼水好多| 国产乱人视频| 超碰av人人做人人爽久久| 建设人人有责人人尽责人人享有的 | 国产精品国产三级国产专区5o| av专区在线播放| av在线蜜桃| 亚洲国产欧美在线一区| 狂野欧美白嫩少妇大欣赏| 久久婷婷青草| 亚洲精品日韩在线中文字幕| 97在线视频观看| 久久久久精品性色| 人人妻人人看人人澡| 男女边摸边吃奶| 欧美国产精品一级二级三级 | 精品一品国产午夜福利视频| av在线观看视频网站免费| 精品一区在线观看国产| 成人二区视频| 久久久久久久久久久免费av| 日韩大片免费观看网站| 国产亚洲精品久久久com| tube8黄色片| 成年美女黄网站色视频大全免费 | 26uuu在线亚洲综合色| av在线播放精品| 国产精品蜜桃在线观看| 99视频精品全部免费 在线| 99久久中文字幕三级久久日本| 国产免费一区二区三区四区乱码| 精品久久久噜噜| 在线看a的网站| 大码成人一级视频| 99热全是精品| 国产伦理片在线播放av一区| 国产高清不卡午夜福利| 99热国产这里只有精品6| 久久99热这里只频精品6学生| 熟女人妻精品中文字幕| 亚洲性久久影院| 国产精品无大码| 日韩成人av中文字幕在线观看| 91精品国产九色| 亚洲最大成人中文| 亚洲精品色激情综合| 99久国产av精品国产电影| 日本欧美视频一区| 亚洲精品乱码久久久久久按摩| 91久久精品国产一区二区三区| 久久久久精品性色| 国产在线视频一区二区| 久久久久精品久久久久真实原创| 王馨瑶露胸无遮挡在线观看| 涩涩av久久男人的天堂| 日韩一本色道免费dvd| 国产精品爽爽va在线观看网站| 性色av一级| 日韩国内少妇激情av| 久久久久精品久久久久真实原创| 女性被躁到高潮视频| 秋霞伦理黄片| 日本色播在线视频| 亚洲婷婷狠狠爱综合网| 日韩视频在线欧美| 少妇人妻一区二区三区视频| 天天躁日日操中文字幕| 丝袜喷水一区| 久久热精品热| 日日摸夜夜添夜夜爱| av线在线观看网站| 边亲边吃奶的免费视频| 免费大片18禁| 极品少妇高潮喷水抽搐| 日韩av在线免费看完整版不卡| 亚洲国产精品一区三区| 亚洲国产成人一精品久久久| 欧美区成人在线视频| 日韩三级伦理在线观看| a级毛色黄片| 久久鲁丝午夜福利片| 日韩不卡一区二区三区视频在线| 国产精品99久久99久久久不卡 | 99热这里只有是精品50| 一区二区三区精品91| 久久99热这里只频精品6学生| 丝瓜视频免费看黄片| 免费观看av网站的网址| 丰满迷人的少妇在线观看| 亚洲精品国产av蜜桃| 国产av码专区亚洲av| 少妇的逼好多水| 性色av一级| 亚洲成色77777| 久久久久视频综合| 日本黄色片子视频| 成年美女黄网站色视频大全免费 | 国产欧美另类精品又又久久亚洲欧美| 人妻系列 视频| www.av在线官网国产| 亚洲国产精品专区欧美| 在线免费观看不下载黄p国产| 久久久久久久久久成人| 国产精品国产三级国产av玫瑰| 下体分泌物呈黄色| 人妻制服诱惑在线中文字幕| 高清av免费在线| 妹子高潮喷水视频| 黑丝袜美女国产一区| 一级片'在线观看视频| 久久精品国产亚洲网站| 亚洲成色77777| 久久人人爽人人爽人人片va| 国产精品久久久久久久久免| 大话2 男鬼变身卡| 久久久久人妻精品一区果冻| 高清毛片免费看| 亚洲欧美一区二区三区黑人 | 久久精品国产亚洲网站| 国产熟女欧美一区二区| 国产精品一区二区在线不卡| av黄色大香蕉| 自拍偷自拍亚洲精品老妇| 亚洲成人手机| 国产成人午夜福利电影在线观看| 亚洲美女视频黄频| 免费黄色在线免费观看| 国产精品一区二区三区四区免费观看| 国产国拍精品亚洲av在线观看| 国产无遮挡羞羞视频在线观看| 人妻制服诱惑在线中文字幕| 免费av中文字幕在线| 在线亚洲精品国产二区图片欧美 | 成人亚洲精品一区在线观看 | 97在线人人人人妻| 亚洲精品,欧美精品| av国产久精品久网站免费入址| 国产伦精品一区二区三区四那| 观看av在线不卡| 成人无遮挡网站| 久久这里有精品视频免费| av免费观看日本| 国内精品宾馆在线| 内射极品少妇av片p| 一级毛片 在线播放| 欧美xxxx性猛交bbbb| 免费大片黄手机在线观看| 国产精品人妻久久久影院| 一级黄片播放器| 久久久精品免费免费高清| 免费在线观看成人毛片| 亚洲国产色片| 中文在线观看免费www的网站| 色综合色国产| 亚洲一区二区三区欧美精品| 少妇人妻久久综合中文| 偷拍熟女少妇极品色| 大香蕉97超碰在线| 欧美日韩亚洲高清精品| 尾随美女入室| 国产免费福利视频在线观看| 少妇人妻一区二区三区视频| 国产精品欧美亚洲77777| 嫩草影院入口| 亚洲国产色片| 老熟女久久久| 波野结衣二区三区在线| 国产精品av视频在线免费观看| 国产成人aa在线观看| 国产黄片美女视频| 91久久精品国产一区二区三区| 国产男人的电影天堂91| 欧美 日韩 精品 国产| 亚洲欧美精品自产自拍| 久久午夜福利片| av视频免费观看在线观看| 成人国产麻豆网| 久久精品夜色国产| 欧美xxⅹ黑人| 卡戴珊不雅视频在线播放| 日韩一区二区视频免费看| 少妇被粗大猛烈的视频| 性色av一级| 身体一侧抽搐| 一级毛片我不卡| 在线播放无遮挡| 一本色道久久久久久精品综合| 国产精品一二三区在线看| 99视频精品全部免费 在线| 超碰av人人做人人爽久久| 亚洲第一av免费看| 观看av在线不卡| 黄色视频在线播放观看不卡| 熟女电影av网| 91在线精品国自产拍蜜月| 日韩视频在线欧美| 丰满人妻一区二区三区视频av| 一级av片app| 交换朋友夫妻互换小说| 我要看黄色一级片免费的| 免费看光身美女| 亚洲欧洲日产国产| 成人美女网站在线观看视频| 亚洲av.av天堂| 国产一级毛片在线| 免费高清在线观看视频在线观看| 在线免费十八禁| 性高湖久久久久久久久免费观看| 99九九线精品视频在线观看视频| 夜夜爽夜夜爽视频| 中国三级夫妇交换| 国产深夜福利视频在线观看| 最黄视频免费看| 水蜜桃什么品种好| 欧美三级亚洲精品| 能在线免费看毛片的网站| 国产毛片在线视频| 激情五月婷婷亚洲| 联通29元200g的流量卡| 久久久久久久久大av| 亚洲成人一二三区av| 午夜福利高清视频| 丰满少妇做爰视频| 五月开心婷婷网| 中国国产av一级| 欧美一级a爱片免费观看看| 精品国产三级普通话版| 国产无遮挡羞羞视频在线观看| 成人毛片a级毛片在线播放| 99久久综合免费| 中文字幕久久专区| 中文字幕av成人在线电影| 舔av片在线| 美女高潮的动态| 国产又色又爽无遮挡免| 国产精品久久久久久久电影| 哪个播放器可以免费观看大片| 噜噜噜噜噜久久久久久91| 精品人妻熟女av久视频| 夜夜爽夜夜爽视频| 熟女电影av网| 亚洲一级一片aⅴ在线观看| 好男人视频免费观看在线| 国产精品国产三级国产专区5o| 内射极品少妇av片p| 国产精品无大码| 91久久精品国产一区二区成人| 最近最新中文字幕大全电影3| 婷婷色综合大香蕉| 国产无遮挡羞羞视频在线观看| 日韩人妻高清精品专区| 国产淫片久久久久久久久| 人体艺术视频欧美日本| 欧美成人a在线观看| 久久久久久久久久成人| 国产精品女同一区二区软件| 日韩在线高清观看一区二区三区| 多毛熟女@视频| 国产高清三级在线| 成人国产av品久久久| 最新中文字幕久久久久| 久久久久性生活片| 91午夜精品亚洲一区二区三区| 亚洲av免费高清在线观看| 亚洲av二区三区四区| 观看美女的网站| 午夜视频国产福利| 高清日韩中文字幕在线| 成年av动漫网址| 国产毛片在线视频| 日本欧美国产在线视频| 国产亚洲午夜精品一区二区久久| 菩萨蛮人人尽说江南好唐韦庄| 日韩人妻高清精品专区| 国产精品精品国产色婷婷| 国产欧美日韩一区二区三区在线 | 成人无遮挡网站| 久久久久国产精品人妻一区二区| 嘟嘟电影网在线观看| tube8黄色片| 国产在线男女| 色网站视频免费| 在线观看免费日韩欧美大片 | a级一级毛片免费在线观看| 亚洲成人av在线免费| 久久国产乱子免费精品| 黄色怎么调成土黄色| 女的被弄到高潮叫床怎么办| 日日啪夜夜撸| 国产精品嫩草影院av在线观看| av天堂中文字幕网| 嫩草影院入口| 成人免费观看视频高清| 久久久精品94久久精品| av专区在线播放| 爱豆传媒免费全集在线观看| 亚洲av成人精品一区久久| 亚洲怡红院男人天堂| 青春草国产在线视频| 在线观看av片永久免费下载| 丝袜脚勾引网站| 亚洲电影在线观看av| 国产免费视频播放在线视频| 国产成人精品福利久久| 亚洲aⅴ乱码一区二区在线播放| 亚洲性久久影院| 欧美少妇被猛烈插入视频| 国产精品麻豆人妻色哟哟久久| 99久久精品热视频| 亚洲av在线观看美女高潮| 少妇人妻一区二区三区视频| 婷婷色av中文字幕| 国产 一区精品| 国产视频内射| av播播在线观看一区| 少妇人妻精品综合一区二区| 国产爱豆传媒在线观看| 国产美女午夜福利| 日韩av不卡免费在线播放| 女性被躁到高潮视频| av.在线天堂| 能在线免费看毛片的网站| 三级经典国产精品| 嫩草影院入口| 一级毛片黄色毛片免费观看视频| 国产精品熟女久久久久浪| 女人久久www免费人成看片| 日本爱情动作片www.在线观看| 国产伦理片在线播放av一区| 一级毛片久久久久久久久女| 精品少妇黑人巨大在线播放| 日本爱情动作片www.在线观看| 十八禁网站网址无遮挡 | 天美传媒精品一区二区| 亚洲精华国产精华液的使用体验| 久久久精品94久久精品| 观看免费一级毛片| 美女主播在线视频| 99久久人妻综合| 精品国产三级普通话版| 免费黄网站久久成人精品| 看非洲黑人一级黄片| 嘟嘟电影网在线观看| av免费观看日本| 免费看日本二区| 日本av免费视频播放| 国产视频内射| 最黄视频免费看| 国产精品蜜桃在线观看| 免费播放大片免费观看视频在线观看| 国产深夜福利视频在线观看| 色婷婷av一区二区三区视频| 一个人免费看片子| 自拍欧美九色日韩亚洲蝌蚪91 | 天堂中文最新版在线下载| 一本久久精品| 亚洲精品aⅴ在线观看| 精品一区二区三区视频在线| 国产永久视频网站| 一级黄片播放器| 日日啪夜夜撸| 国产v大片淫在线免费观看| 一区在线观看完整版| 美女视频免费永久观看网站| 26uuu在线亚洲综合色| 纯流量卡能插随身wifi吗| 看非洲黑人一级黄片| 内地一区二区视频在线| 亚洲综合色惰| 日韩电影二区| 伦理电影免费视频| 伦精品一区二区三区| 观看美女的网站| 国产成人freesex在线| av国产精品久久久久影院| 国产欧美亚洲国产| 18禁裸乳无遮挡免费网站照片| 高清午夜精品一区二区三区| 简卡轻食公司| 汤姆久久久久久久影院中文字幕| 欧美 日韩 精品 国产| 精品国产乱码久久久久久小说| 18禁动态无遮挡网站| 最新中文字幕久久久久| 新久久久久国产一级毛片| 国产精品一及| 中文字幕久久专区| 91久久精品国产一区二区三区| 亚洲精华国产精华液的使用体验| 免费观看性生交大片5| 国产在线免费精品| 成人无遮挡网站| 男女免费视频国产| 97精品久久久久久久久久精品| freevideosex欧美| 久热久热在线精品观看| 亚洲国产av新网站| 国产日韩欧美在线精品| 丰满乱子伦码专区| 亚洲久久久国产精品| 久久97久久精品| 日本黄色日本黄色录像| 18禁在线无遮挡免费观看视频| 国产av精品麻豆| 国内精品宾馆在线| 亚洲精品视频女| 成年免费大片在线观看| 亚洲怡红院男人天堂| 国产成人a区在线观看| 精品久久久久久久久av| 91狼人影院| 久久6这里有精品| 国产高清不卡午夜福利| 在线 av 中文字幕| 2022亚洲国产成人精品| 欧美日韩亚洲高清精品| 97超视频在线观看视频| 日韩欧美一区视频在线观看 | 99热这里只有是精品在线观看| 少妇人妻一区二区三区视频| 中文字幕人妻熟人妻熟丝袜美| 亚洲国产欧美在线一区| 另类亚洲欧美激情| av国产精品久久久久影院| 最近中文字幕2019免费版| 国产精品嫩草影院av在线观看| 2022亚洲国产成人精品| 九九爱精品视频在线观看| 十分钟在线观看高清视频www | 内地一区二区视频在线| 深夜a级毛片| 国产爽快片一区二区三区| 久久久久久伊人网av| 精品一品国产午夜福利视频| 欧美一区二区亚洲| 一个人看的www免费观看视频| 亚洲精品国产av成人精品| 国产av码专区亚洲av| 亚洲国产精品成人久久小说| av在线观看视频网站免费| 在线亚洲精品国产二区图片欧美 | 国产一区二区在线观看日韩| 又黄又爽又刺激的免费视频.| 夫妻午夜视频| 伊人久久精品亚洲午夜| 国产精品久久久久成人av| 精品人妻偷拍中文字幕| 美女视频免费永久观看网站| 一本色道久久久久久精品综合| 秋霞伦理黄片| 丰满迷人的少妇在线观看| 女人久久www免费人成看片| 国产精品熟女久久久久浪| 亚洲精品,欧美精品| 国产爽快片一区二区三区| 777米奇影视久久| 色5月婷婷丁香| 91精品国产国语对白视频| 久久人妻熟女aⅴ| 人妻系列 视频| 最新中文字幕久久久久| 一级毛片aaaaaa免费看小| 青青草视频在线视频观看| 夫妻性生交免费视频一级片| 丝袜脚勾引网站| 婷婷色综合大香蕉| 精品久久久噜噜| 丰满迷人的少妇在线观看| 久久久久久久久久人人人人人人| 国产精品一区www在线观看| 国产精品无大码| 最近最新中文字幕大全电影3| 一级毛片电影观看| 人人妻人人澡人人爽人人夜夜| 亚洲精品日本国产第一区| 国产精品偷伦视频观看了| 99久久人妻综合| 免费av中文字幕在线| 看免费成人av毛片| 色婷婷久久久亚洲欧美| 91aial.com中文字幕在线观看| 久久精品国产鲁丝片午夜精品| 一边亲一边摸免费视频| 日韩国内少妇激情av| 亚洲精品国产av成人精品| 亚洲真实伦在线观看| 老熟女久久久| 丰满乱子伦码专区| 老司机影院毛片| 国产精品99久久99久久久不卡 | 亚洲精品亚洲一区二区| 精品99又大又爽又粗少妇毛片| xxx大片免费视频| 伊人久久精品亚洲午夜| 日韩亚洲欧美综合| 国产 一区精品| 夜夜爽夜夜爽视频| 精品少妇久久久久久888优播| 一级毛片黄色毛片免费观看视频| 欧美丝袜亚洲另类| 精品国产乱码久久久久久小说| 国产淫语在线视频| 91精品国产国语对白视频| 在线观看av片永久免费下载| 老熟女久久久| 青春草亚洲视频在线观看| 黑丝袜美女国产一区| 免费黄色在线免费观看| 亚洲av.av天堂| 亚洲欧美成人精品一区二区| 欧美日韩国产mv在线观看视频 | 国产爽快片一区二区三区| 国国产精品蜜臀av免费| 国产男女超爽视频在线观看| 亚洲精品成人av观看孕妇| 国产成人免费观看mmmm| 欧美老熟妇乱子伦牲交| 国产亚洲精品久久久com| 一区在线观看完整版| 成年女人在线观看亚洲视频| av卡一久久| 亚洲精品自拍成人| 国产乱人偷精品视频| 国产高清国产精品国产三级 | 麻豆乱淫一区二区| 国产精品av视频在线免费观看| 赤兔流量卡办理| 欧美日韩国产mv在线观看视频 | videos熟女内射| 亚洲av电影在线观看一区二区三区| 最近中文字幕2019免费版| 国产69精品久久久久777片| 亚洲精品乱久久久久久| 免费大片黄手机在线观看| 狂野欧美白嫩少妇大欣赏| 国产成人免费观看mmmm| 久久精品国产亚洲av涩爱| 丝袜脚勾引网站| 亚洲伊人久久精品综合| 久久久久久久久久成人| 热99国产精品久久久久久7| 在线观看一区二区三区| 亚洲丝袜综合中文字幕| 男男h啪啪无遮挡| 精品人妻偷拍中文字幕| av黄色大香蕉| 一级毛片久久久久久久久女| 国产免费一级a男人的天堂| 欧美极品一区二区三区四区| 欧美国产精品一级二级三级 | 日韩三级伦理在线观看| 国产探花极品一区二区| 国产精品久久久久久精品电影小说 | 一级毛片 在线播放| 伦精品一区二区三区| 简卡轻食公司| av免费在线看不卡| 亚洲在久久综合| 最近2019中文字幕mv第一页| 国产成人91sexporn| 国产伦在线观看视频一区| 国内精品宾馆在线| 国产精品av视频在线免费观看| 人人妻人人爽人人添夜夜欢视频 | 久久99热6这里只有精品| 又粗又硬又长又爽又黄的视频| 性高湖久久久久久久久免费观看| 青春草国产在线视频| 高清av免费在线|