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

    振動(dòng)激發(fā)對(duì)高超聲速氣動(dòng)力/熱影響1)

    2017-07-03 15:00:44張子健劉云峰姜宗林
    力學(xué)學(xué)報(bào) 2017年3期
    關(guān)鍵詞:摩阻氣動(dòng)力邊界層

    張子健劉云峰 姜宗林

    (中國(guó)科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京100190)(中國(guó)科學(xué)院大學(xué),北京100049)

    -流體力學(xué)

    振動(dòng)激發(fā)對(duì)高超聲速氣動(dòng)力/熱影響1)

    張子健2)劉云峰 姜宗林

    (中國(guó)科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京100190)(中國(guó)科學(xué)院大學(xué),北京100049)

    隨著飛行馬赫數(shù)的不斷提高,空氣的高溫氣體效應(yīng)越來(lái)越明顯,對(duì)高超聲速飛行器的氣動(dòng)力/熱特性產(chǎn)生重要影響.高溫氣體效應(yīng)對(duì)氣動(dòng)力/熱的影響機(jī)理復(fù)雜,影響參數(shù)眾多,迄今為止國(guó)內(nèi)外尚未完全研究清楚.發(fā)生高溫氣體效應(yīng)時(shí),多個(gè)非線性物理過(guò)程耦合在一起,地面試驗(yàn)和數(shù)值模擬無(wú)法將這些過(guò)程解耦,無(wú)法給出關(guān)鍵物理機(jī)理.為了解決這一問(wèn)題,文章提出一種理論分析與數(shù)值模擬相結(jié)合的兩步漸進(jìn)新方法:先通過(guò)牛頓迭代法得到發(fā)生振動(dòng)激發(fā)過(guò)程的斜激波無(wú)黏解;再將該無(wú)黏解的結(jié)果作為邊界條件,求解邊界層的黏性解.利用該方法研究了振動(dòng)激發(fā)過(guò)程對(duì)二維斜劈的氣動(dòng)力/熱特性的影響規(guī)律.研究結(jié)果表明,振動(dòng)激發(fā)過(guò)程對(duì)斜激波后的溫度、密度、馬赫數(shù)、雷諾數(shù)和斜激波角影響較大,而對(duì)壓力和速度影響較小.斜激波波后的無(wú)黏流動(dòng)與邊界層流動(dòng)是耦合在一起的.發(fā)生振動(dòng)激發(fā)后,斜激波波后雷諾數(shù)的增大會(huì)導(dǎo)致邊界層厚度減小,結(jié)合多個(gè)物理量的變化,如速度增大和溫度減小,共同對(duì)邊界層內(nèi)的摩擦阻力和氣動(dòng)熱產(chǎn)生影響.對(duì)比完全氣體的結(jié)果發(fā)現(xiàn),振動(dòng)激發(fā)使壁面摩阻升高,而使壁面熱流降低.分別通過(guò)影響激波層和邊界層,振動(dòng)激發(fā)對(duì)摩阻的影響是弱耦合的,而對(duì)熱流的影響則是強(qiáng)耦合的.

    振動(dòng)激發(fā),氣動(dòng)力/熱,邊界層,高溫氣體效應(yīng)

    引言

    近十幾年來(lái),高超聲速飛行器由于具有巨大的戰(zhàn)略價(jià)值,受到了許多國(guó)家的高度重視,成為世界航空航天領(lǐng)域的重要發(fā)展方向和研究熱點(diǎn).高超聲速飛行器是指在大氣層內(nèi)實(shí)現(xiàn)高超聲速機(jī)動(dòng)飛行的飛行器,其飛行速度一般超過(guò)5倍聲速,主要包括航天運(yùn)載器、太空作戰(zhàn)飛行器、高超聲速巡航導(dǎo)彈、高超聲速飛機(jī)等.

    在高超聲速飛行器的飛行過(guò)程中,高速空氣經(jīng)過(guò)激波壓縮和黏性阻滯后減速,空氣巨大的動(dòng)能轉(zhuǎn)變?yōu)闊崮?,?dǎo)致飛行器表面激波層和邊界層內(nèi)空氣的溫度急劇升高.空氣分子在如此高的溫度下將發(fā)生振動(dòng)激發(fā)、解離和電離等熱化學(xué)反應(yīng),這就是高溫氣體效應(yīng)[14].在高溫氣體效應(yīng)的影響下,空氣介質(zhì)性質(zhì)發(fā)生變化,進(jìn)而影響飛行器的外部流動(dòng),最終導(dǎo)致高超聲速飛行器的氣動(dòng)力/熱特性偏離完全氣體理論的預(yù)測(cè)結(jié)果.

    高溫氣體效應(yīng)對(duì)高超聲速飛行器氣動(dòng)力特性的影響最早是在航天飛機(jī)再入時(shí)所發(fā)生的“上仰異?!爆F(xiàn)象中發(fā)現(xiàn)的[511].與航天飛機(jī)大鈍頭體、大攻角再入不同的是,新一代高超聲速飛行器的飛行空域?yàn)?0~100km高度的臨近空間,多采用乘波體或升力體構(gòu)型,飛行器頭部不會(huì)形成很強(qiáng)的弓形激波,而主要是斜激波.飛行馬赫數(shù)為6~15時(shí),斜激波后空氣溫度為600~2500K,發(fā)生的熱化學(xué)反應(yīng)主要是空氣分子振動(dòng)激發(fā)[34].另一方面,與航天飛機(jī)相比,新一代高超聲速飛行器飛行時(shí)間長(zhǎng)、飛行距離遠(yuǎn),對(duì)其所受的氣動(dòng)力/熱的預(yù)測(cè)精度要求非常高,希望實(shí)現(xiàn)精細(xì)控制.因此,非常有必要深入細(xì)致地研究空氣分子振動(dòng)激發(fā)對(duì)高超聲速飛行器氣動(dòng)力/熱特性的影響.

    中國(guó)科學(xué)院力學(xué)研究所的JF12激波風(fēng)洞[1213],是國(guó)際上首座可復(fù)現(xiàn)25~40km高空、馬赫數(shù)5~9飛行條件、噴管出口直徑Φ2.5/Φ1.5m、試驗(yàn)氣體為潔凈空氣、試驗(yàn)時(shí)間超過(guò)100ms的超大型高超聲速風(fēng)洞.目前,在JF12復(fù)現(xiàn)激波風(fēng)洞上正大量開(kāi)展高溫氣體效應(yīng)對(duì)氣動(dòng)力/熱特性影響的研究,并取得了一些成果[14-15].

    高溫氣體效應(yīng)對(duì)高超聲速飛行器氣動(dòng)力/熱特性的影響主要體現(xiàn)在以下兩個(gè)方面[34]:一是對(duì)激波形狀和激波層內(nèi)壓力、溫度等的影響,這將會(huì)影響飛行器表面的壓力分布以及激波與激波相互作用、激波與邊界層相互作用的位置和強(qiáng)度;二是對(duì)邊界層內(nèi)流動(dòng)的影響,這將會(huì)直接影響飛行器的摩阻和熱流.更重要的是,高溫氣體效應(yīng)對(duì)激波和邊界層的影響會(huì)相互耦合疊加,使其對(duì)氣動(dòng)力/熱的影響機(jī)理變得極其復(fù)雜.

    本文主要研究高溫下空氣分子振動(dòng)激發(fā)對(duì)高超聲速氣動(dòng)力/熱的影響,暫時(shí)不考慮解離和電離等過(guò)程.先給出考慮振動(dòng)激發(fā)的空氣熱力學(xué)模型,以及斜激波與邊界層的求解方法;然后分別單獨(dú)討論振動(dòng)激發(fā)對(duì)斜激波波后流場(chǎng)和邊界層內(nèi)壁面摩阻、熱流的影響;最后采用兩種方法分析斜激波和邊界層共同考慮振動(dòng)激發(fā)時(shí),對(duì)氣動(dòng)力/熱特性的影響,并討論其影響機(jī)制.

    1 熱力學(xué)模型與研究方法

    1.1 高溫空氣的熱力學(xué)模型

    常溫下,空氣一般不考慮熱化學(xué)反應(yīng),按完全氣體來(lái)處理,其比內(nèi)能、比焓、比定壓熱容和比熱比等滿足以下熱力學(xué)關(guān)系

    式中,R是空氣的氣體常數(shù),取值288.28J/(kg·K).但考慮振動(dòng)激發(fā)后,空氣的比定壓熱容和比熱比不再是常數(shù),而是溫度的非線性函數(shù).為了簡(jiǎn)化分析,假設(shè)空氣是由79%的N2和21%的O2組成的,則考慮振動(dòng)激發(fā)后,比內(nèi)能和比焓中將多出振動(dòng)能項(xiàng)ev.

    式中,χN2和χO2分別是N2和 O2的摩爾分?jǐn)?shù),而Tve,N2和Tve,O2是振動(dòng)特征溫度,取值為[16]

    則比定壓熱容和比熱比分別為

    根據(jù)以上熱力學(xué)關(guān)系,靜止空氣比振動(dòng)能占比焓的比例以及空氣比熱比隨溫度的變化關(guān)系分別如圖1所示.可見(jiàn),當(dāng)溫度上升到800K時(shí),空氣振動(dòng)能所占比例已達(dá)2.5%,比熱比從1.4降低到1.353,此時(shí)振動(dòng)激發(fā)對(duì)空氣熱力學(xué)性質(zhì)的影響已不可忽略[2].當(dāng)溫度達(dá)到2500K時(shí),空氣振動(dòng)能所占比例達(dá)12.6%,比熱比降低到1.296,此時(shí)振動(dòng)激發(fā)是相當(dāng)可觀的,空氣的熱力學(xué)性質(zhì)將會(huì)發(fā)生較大變化,進(jìn)而影響飛行器的外部流動(dòng),最終影響高超聲速飛行器的氣動(dòng)力/熱特性.

    圖1 靜止空氣比振動(dòng)能占比焓的比例和空氣比熱比隨溫度的變化Fig.1 Variation of proportion of specifi vibration energy in specifi enthapy of stationary air and the specifi heat ratio with temperature

    1.2 激波---邊界層組合結(jié)構(gòu)的兩步漸進(jìn)方法

    本文提出一種理論求解斜激波與數(shù)值求解邊界層相結(jié)合的兩步漸進(jìn)方法,將振動(dòng)激發(fā)對(duì)激波和邊界層的影響分開(kāi),來(lái)分析振動(dòng)激發(fā)效應(yīng)在激波層與邊界層內(nèi)的傳遞和干擾機(jī)理.高超聲速飛行器的氣動(dòng)力/熱特性是由其周圍流場(chǎng)決定的,基本的流場(chǎng)結(jié)構(gòu)是來(lái)流先經(jīng)過(guò)一道激波壓縮,然后通過(guò)邊界層作用在飛行器壁面上.本文采用斜劈作為基本幾何構(gòu)型進(jìn)行分析,將振動(dòng)激發(fā)當(dāng)成是一種擾動(dòng)施加在流場(chǎng)上,如圖2.

    圖2 流場(chǎng)示意圖Fig.2 Schematic of the fl w fiel

    根據(jù)邊界層理論,壁面壓力主要由斜激波波后壓力決定,而摩阻和熱流則主要由斜激波波后流動(dòng)和邊界層共同決定.于是,振動(dòng)激發(fā)先對(duì)斜激波波后流場(chǎng)產(chǎn)生影響,使其壓力、溫度、速度等發(fā)生變化,從而對(duì)壁面壓力產(chǎn)生影響;接著斜激波波后流動(dòng)作為邊界層的外流,其速度、溫度的變化結(jié)合振動(dòng)激發(fā)共同影響邊界層內(nèi)流動(dòng),最終影響壁面摩阻、熱流.作為近似,忽略邊界層對(duì)斜激波的反饋?zhàn)饔?這種近似方法在振動(dòng)激發(fā)效應(yīng)分析中的誤差將在后續(xù)2.3.2節(jié)中進(jìn)行簡(jiǎn)要說(shuō)明.

    因此,在漸進(jìn)求解時(shí),先求解斜激波關(guān)系,然后將斜激波波后流場(chǎng)參數(shù)作為平板邊界層問(wèn)題的來(lái)流參數(shù)進(jìn)一步求解,最后根據(jù)邊界層內(nèi)流場(chǎng)得到壁面摩阻和熱流.由于邊界層內(nèi)壓力與外流壓力非常接近,因此壁面壓力可以直接取為斜激波后壓力.另外,本文還直接求解黏性斜劈問(wèn)題來(lái)驗(yàn)證以上兩步漸進(jìn)分析方法的準(zhǔn)確性.

    1.3 精確求解斜激波關(guān)系

    考慮振動(dòng)激發(fā)后,空氣的焓是溫度的非線性函數(shù),難以得到解析的斜激波關(guān)系.本文從基本方程出發(fā),推導(dǎo)出精確求解斜激波關(guān)系的牛頓迭代公式

    其中,β是斜激波角,θ是楔面角,其余流場(chǎng)參數(shù)的下標(biāo)1表示激波前,下標(biāo)2表示激波后,下標(biāo)n表示垂直于激波的分量.對(duì)于完全氣體,h(T)采用式(2)進(jìn)行計(jì)算,而對(duì)于振動(dòng)激發(fā),h(T)則采用式(6)進(jìn)行計(jì)算.根據(jù)以上牛頓迭代公式可以求出斜激波角β,進(jìn)而求出所有斜激波波后流場(chǎng)參數(shù).

    1.4 數(shù)值求解邊界層問(wèn)題

    控制方程采用二維Navier-Stokes方程[17-18]

    其中,U為守恒變量,F(xiàn),G分別為x,y方向的無(wú)黏守恒通量,F(xiàn)v,Gv分別為x,y方向的黏性守恒通量,具體形式為

    其中對(duì)于內(nèi)能e,完全氣體采用式(1),而振動(dòng)激發(fā)采用式(5).黏性系數(shù)μ隨溫度T變化的近似關(guān)系由Sutherland公式給出,而熱傳導(dǎo)系數(shù)k則由普朗特?cái)?shù)Pr給出[19]

    式中,C1=1.458×10-6,C2=110.4,Pr=0.72.

    計(jì)算中采用有限差分方法、LUSGS隱式格式[20]進(jìn)行求解.其中對(duì)無(wú)黏通量 F和 G的求解采用AUSMPW+格式[2122],將流動(dòng)通量分為對(duì)流項(xiàng)和壓力項(xiàng)兩部分,分別進(jìn)行近似求解.該數(shù)值格式在高超聲速流動(dòng)的求解中,對(duì)激波和邊界層的捕捉均表現(xiàn)出較好的性能[2324].迎風(fēng)格式基本上是一階的,為了取得更好的空間精度,采用MUSCL格式[25]對(duì)原始變量進(jìn)行重構(gòu),并引入Van Albada限制器來(lái)限制重構(gòu)時(shí)產(chǎn)生過(guò)大或過(guò)小的梯度

    其中,W是原始變量.當(dāng)κ=1/3時(shí),該方法具有三階空間精度[26].

    對(duì)于平板邊界層問(wèn)題,計(jì)算域和網(wǎng)格示意圖如圖3.其中總網(wǎng)格數(shù)為858×812≈7.0×105,邊界層第一層網(wǎng)格間距為5μm.壁面邊界條件為等溫壁,Tw=1500K.

    圖3 計(jì)算域和網(wǎng)格示意圖Fig.3 Computational domain and grids

    2 結(jié)果與討論

    2.1 振動(dòng)激發(fā)對(duì)斜激波的影響

    這里采用1.3節(jié)提出的牛頓迭代公式分別對(duì)完全氣體和考慮振動(dòng)激發(fā)的斜激波關(guān)系進(jìn)行精確求解,以分析振動(dòng)激發(fā)對(duì)斜激波波后流場(chǎng)的定量影響.計(jì)算時(shí),來(lái)流采用40km高空空氣參數(shù)(T=250K,p=287Pa),楔面角為 30?.

    圖4給出了來(lái)流馬赫數(shù)從6增大到12時(shí),斜激波后的溫度、壓力、密度、速度、馬赫數(shù)和激波角的變化.從圖中可以看出,與完全氣體相比,考慮振動(dòng)激發(fā)后斜激波波后溫度降低、壓力變小、密度變大、速度和馬赫數(shù)變大、激波角變小.且隨著來(lái)流馬赫數(shù)的增大,激波壓縮增強(qiáng),波后溫度增大,振動(dòng)激發(fā)變得越來(lái)越顯著,斜激波波后流場(chǎng)參數(shù)的偏差也越來(lái)越大.事實(shí)上,其他使斜激波波后溫度增大的過(guò)程,如升高來(lái)流溫度或增大楔面角,都會(huì)使振動(dòng)激發(fā)變得顯著,進(jìn)而導(dǎo)致振動(dòng)激發(fā)對(duì)斜激波的影響變大.

    圖4 斜激波前后流場(chǎng)參數(shù)Fig.4 Flow fiel parameters across oblique shock

    取Ma=10的斜激波關(guān)系如表1所示.可見(jiàn),振動(dòng)激發(fā)對(duì)斜激波波后溫度、密度、馬赫數(shù)和激波角的影響較大,超過(guò)12%;而對(duì)波后壓力和速度的影響較小,在4%以內(nèi).事實(shí)上,考慮振動(dòng)激發(fā)后,斜激波波后氣體能量有較大一部分轉(zhuǎn)化成氣體分子的振動(dòng)能,使得氣體分子平動(dòng)能和轉(zhuǎn)動(dòng)能顯著變小,于是作為表征分子平動(dòng)能的溫度比完全氣體的相應(yīng)值要小得多.這是振動(dòng)激發(fā)對(duì)流場(chǎng)最直接的影響.溫度顯著變小,則聲速也將顯著變小,由于速度基本保持不變,導(dǎo)致波后馬赫數(shù)顯著變大.另外,由于壓力變化不大,從狀態(tài)方程可得密度將顯著變大.

    表1 Ma=10斜激波關(guān)系Table 1 Oblique shock properties at Ma=10

    注意到,考慮振動(dòng)激發(fā)后,斜激波波后流動(dòng)的單位雷諾數(shù)Rex變化明顯,達(dá)到23.65%.雷諾數(shù)的增大會(huì)導(dǎo)致邊界層厚度減小,結(jié)合多個(gè)物理量的變化,如速度和溫度,將共同對(duì)邊界層內(nèi)的摩擦阻力和氣動(dòng)熱產(chǎn)生一定影響.由式(26)可知,溫度顯著降低將導(dǎo)致黏性系數(shù)μ顯著變小,且由于密度顯著變大,則式(31)給出的單位雷諾數(shù)Rex顯著變大.

    振動(dòng)激發(fā)對(duì)斜激波波后流場(chǎng)參數(shù)產(chǎn)生顯著影響,該流場(chǎng)又作為邊界層的外流,將與振動(dòng)激發(fā)對(duì)邊界層流動(dòng)的影響進(jìn)行耦合,共同對(duì)邊界層內(nèi)流場(chǎng)產(chǎn)生影響,最終影響壁面處的氣動(dòng)特性.

    2.2 振動(dòng)激發(fā)對(duì)邊界層的影響

    這里分析邊界層外流取斜激波波后流場(chǎng)參數(shù)時(shí),振動(dòng)激發(fā)對(duì)邊界層流動(dòng)的影響.采用1.4節(jié)中的數(shù)值方法求解平板邊界層問(wèn)題,即可得到邊界層內(nèi)流場(chǎng)以及壁面熱流和摩阻.

    當(dāng)將斜激波波后流場(chǎng)取為邊界層外流時(shí),由表1可知,完全氣體與考慮振動(dòng)激發(fā)的流場(chǎng)溫度差異很大,從而聲速差異很大,馬赫數(shù)差異也很大.即在高溫氣體效應(yīng)顯著時(shí),采用速度描述流動(dòng)比馬赫數(shù)更準(zhǔn)確.另外,與邊界層內(nèi)壁面摩阻直接相關(guān)的是速度而不是馬赫數(shù).因此,在本小節(jié)采用速度替代馬赫數(shù)來(lái)描述流動(dòng).斜激波前自由來(lái)流為U=3176m/s(Ma=10),T=250K,p=287Pa,θ=30?時(shí),完全氣體或考慮振動(dòng)激發(fā),斜激波波后流場(chǎng)參數(shù)約為U=2530m/s,T=2000K,p=12750Pa.

    2.2.1 對(duì)邊界層流場(chǎng)參數(shù)的影響

    取U=2530m/s,T=2000K,p=12750Pa為邊界層外流,壁面溫度為Tw=1500K,此時(shí)平板邊界層流場(chǎng)的溫度云圖如圖5所示.距離平板前緣x=150mm處,完全氣體或考慮振動(dòng)激發(fā)的溫度邊界層及速度邊界層的對(duì)比如圖6.該位置的邊界層厚度及與壁面熱流、摩阻相關(guān)的參數(shù)匯總?cè)绫?和表3所示.其中速度邊界層厚度δ、壁面熱流qw和壁面摩阻τw分別按式(32)和式(33)定義為

    圖5 平板邊界層(完全氣體)流場(chǎng)溫度云圖Fig.5 Temperature contour of flat-plat boundary layer of perfect gas

    從圖6和表2可見(jiàn),與完全氣體對(duì)比,考慮振動(dòng)激發(fā)后,邊界層內(nèi)最高溫度 Tmax下降明顯,從2390K降低到2278K,降低了4.69%.而邊界層內(nèi)最高溫度所在位置y|Tmax變化較小.由此可得,振動(dòng)激發(fā)使壁面溫度梯度顯著降低,計(jì)算結(jié)果給出變化為-13.11%.然而壁面是等溫壁,Tw=1500K,由式(9)可得,壁面處的比定壓熱容cp|w考慮振動(dòng)激發(fā)后多出與分子振動(dòng)能相關(guān)的項(xiàng),從而顯著增大.而黏性系數(shù)僅與溫度有關(guān),等溫壁處的黏性系數(shù)μw保持不變.因此,由式(26)可得壁面處的熱傳導(dǎo)系數(shù)kw顯著增大,變化為20.18%.根據(jù)壁面熱流的定義式(33),壁面熱流qw的變化為

    從圖6和表3可知,與完全氣體對(duì)比,考慮振動(dòng)激發(fā)后,速度邊界層厚度稍微減小,減小了2.82%.在外流速度相同時(shí),壁面速度梯度也稍微增大,變化為0.95%.另外,等溫壁的黏性系數(shù)μw保持不變.因此,根據(jù)壁面摩阻的定義式(33),壁面摩阻τw與壁面速度梯度的變化趨勢(shì)相同,也有0.95%的小變化.

    綜上,在邊界層外流條件相同時(shí),振動(dòng)激發(fā)對(duì)邊界層內(nèi)流動(dòng)以及壁面熱流、摩阻都產(chǎn)生一定影響.

    圖6 邊界層的溫度和速度剖面(x=150mm)Fig.6 Profile of temperature and velocity in boundary layer at x=150mm

    表2 壁面溫度梯度與熱流(x=150mm)Table 2 Temperature gradient and heat flu at wall of x=150mm

    表3 壁面速度梯度與摩阻(x=150mm)Table 3 Velocity gradient and friction at wall of x=150mm

    2.2.2 隨外流參數(shù)的變化規(guī)律

    邊界層內(nèi)流動(dòng)和壁面熱流、摩阻主要受外流速度 U、溫度 T以及壓力 p的共同影響.在 U =2530m/s,T=2000K,p=12750Pa條件下,固定其中兩個(gè),改變另一個(gè),可以得到完全氣體或考慮振動(dòng)激發(fā)的壁面熱流qw、摩阻τw以及振動(dòng)激發(fā)對(duì)其影響大小隨外流參數(shù)的變化規(guī)律,如圖7.其中,分別按式(35)和式(36)將振動(dòng)激發(fā)對(duì)壁面熱流和摩阻的影響定義成兩個(gè)物理量:δqw和δτw.

    由圖7可見(jiàn),壁面熱流qw對(duì)外流速度U、溫度T、壓力p都比較敏感,而壁面摩阻τw對(duì)外流速度U和壓力p比較敏感,外流溫度T對(duì)τw影響較小.壁面熱流qw、摩阻τw隨外流速度U、溫度T和壓力p而變化說(shuō)明,斜激波波后流動(dòng)作為邊界層外流,在振動(dòng)激發(fā)使斜激波波后流場(chǎng)參數(shù)發(fā)生變化時(shí),最終通過(guò)邊界層的傳遞作用,會(huì)對(duì)壁面熱流和摩阻產(chǎn)生影響.

    從圖7還可得到,振動(dòng)激發(fā)的影響大小δqw和δτw也會(huì)隨外流速度U和外流溫度T而變化.這又說(shuō)明,邊界層的存在不僅傳遞振動(dòng)激發(fā)對(duì)外流產(chǎn)生的影響,還將其與振動(dòng)激發(fā)對(duì)自身的影響進(jìn)行耦合,使振動(dòng)激發(fā)對(duì)氣動(dòng)力/熱特性的影響機(jī)理變得復(fù)雜.其中從圖7(a)和圖7(b)可得,δqw隨U,T變化較大,即振動(dòng)激發(fā)通過(guò)斜激波使邊界層外流發(fā)生變化時(shí),會(huì)極大地干擾振動(dòng)激發(fā)通過(guò)改變邊界層內(nèi)流動(dòng)而對(duì)壁面熱流qw產(chǎn)生的影響大小.這說(shuō)明,振動(dòng)激發(fā)通過(guò)在斜激波中改變邊界層外流而對(duì)qw產(chǎn)生的影響與通過(guò)改變邊界層內(nèi)流動(dòng)而對(duì)qw產(chǎn)生的影響是強(qiáng)烈耦合在一起的,不是簡(jiǎn)單的疊加.

    從圖7(d)和圖7(e)可得,δτw隨U,T變化相對(duì)較小,即在振動(dòng)激發(fā)通過(guò)斜激波使邊界層外流發(fā)生變化時(shí),振動(dòng)激發(fā)通過(guò)改變邊界層內(nèi)流動(dòng)而對(duì)壁面摩阻τw產(chǎn)生的影響大小相對(duì)較穩(wěn)定,所受干擾較小.這說(shuō)明,振動(dòng)激發(fā)通過(guò)在斜激波中改變邊界層外流而對(duì)τw產(chǎn)生的影響,與通過(guò)改變邊界層內(nèi)流動(dòng)而對(duì)τw產(chǎn)生的影響的耦合作用較弱,基本是解耦的.

    圖7 qw,|δqw|,τw,|δτw|隨 U,T,p 的變化Fig.7 Variation of qw,|δqw|,τw,|δτw|as functions of U,T,p

    2.3 綜合分析

    2.3.1 漸進(jìn)求解的結(jié)果

    為了分析斜激波和邊界層同時(shí)考慮振動(dòng)激發(fā)時(shí)對(duì)壁面壓力、摩阻和熱流的影響,這里選擇速度U=3176m/s(Ma=10)、溫度T=250K、壓力p=287Pa的空氣流過(guò)角度θ=30?,長(zhǎng)度L=20cm的斜劈作為實(shí)例,采用1.2節(jié)中提出的漸進(jìn)分析方法分兩步求解.分析結(jié)果如圖8,圖中點(diǎn)標(biāo)志X-X(X=P,R)中,第1個(gè)字母表示斜激波的計(jì)算模型,第2個(gè)字母表示邊界層的計(jì)算模型,P表示完全氣體模型,R表示考慮振動(dòng)激發(fā)的熱完全氣體模型,兩點(diǎn)之間連線上的百分比代表的是兩者的差別大小.

    由圖8(a)可見(jiàn),振動(dòng)激發(fā)對(duì)壁面壓力系數(shù)的影響基本上是通過(guò)改變斜激波波后壓力實(shí)現(xiàn)的,邊界層的作用較小,符合邊界層理論.

    從圖8(b)可見(jiàn),斜激波或邊界層單獨(dú)考慮振動(dòng)激發(fā)時(shí),都會(huì)使壁面摩阻τw增大,分別為1.85%和0.92%,斜激波的影響較大.注意到,P-P,R-P,P-R,R-R四個(gè)點(diǎn)組成的圖形非常接近平行四邊形.這說(shuō)明斜激波考慮振動(dòng)激發(fā)時(shí)對(duì)壁面摩阻τw產(chǎn)生的影響與邊界層考慮振動(dòng)激發(fā)產(chǎn)生的影響近似可以疊加,疊加結(jié)果為:1.85%+0.92%=2.77%,與斜激波和邊界層同時(shí)考慮振動(dòng)激發(fā)時(shí)的結(jié)果(2.85%)非常接近.即振動(dòng)激發(fā)通過(guò)在斜激波中改變邊界層外流而對(duì)τw產(chǎn)生的影響,與通過(guò)改變邊界層內(nèi)流動(dòng)而對(duì)τw產(chǎn)生的影響的耦合作用較弱,基本是解耦的.2.2.2節(jié)中指出,這是因?yàn)檎駝?dòng)激發(fā)通過(guò)改變邊界層內(nèi)流動(dòng)而對(duì)壁面摩阻τw產(chǎn)生的影響大小相對(duì)較穩(wěn)定,受邊界層外流參數(shù)變化的干擾較小,即受振動(dòng)激發(fā)對(duì)斜激波影響的干擾較小.

    由圖8(c)可見(jiàn),斜激波考慮振動(dòng)激發(fā)時(shí)會(huì)使壁面熱流qw減小,影響為-6.51%.而邊界層考慮振動(dòng)激發(fā)時(shí)會(huì)使壁面熱流qw增大,影響為5.18%.即在斜激波或邊界層這兩種流場(chǎng)結(jié)構(gòu)中,振動(dòng)激發(fā)對(duì)壁面熱流qw的影響方向是相反的,當(dāng)斜激波和邊界層同時(shí)考慮振動(dòng)激發(fā)時(shí),兩者的作用會(huì)發(fā)生一定的抵消,最終削弱了振動(dòng)激發(fā)對(duì)壁面熱流qw產(chǎn)生的影響(-3.46%).注意到,與壁面摩阻的結(jié)果不同的是,P-P,R-P,P-R,R-R四個(gè)點(diǎn)組成的圖形與平行四邊形(圖中虛線)差異較大.當(dāng)將斜激波考慮振動(dòng)激發(fā)時(shí)對(duì)壁面熱流qw產(chǎn)生的影響與邊界層考慮振動(dòng)激發(fā)產(chǎn)生的影響進(jìn)行疊加,結(jié)果為:-6.51%+5.18%=-1.33%,圖中(+)點(diǎn).對(duì)比斜激波和邊界層同時(shí)考慮振動(dòng)激發(fā)時(shí)的結(jié)果(-3.46%),差異較大,相差-1.33%-(-3.46%)=2.13%.這說(shuō)明振動(dòng)激發(fā)通過(guò)在斜激波中改變邊界層外流而對(duì)qw產(chǎn)生的影響與通過(guò)改變邊界層內(nèi)流動(dòng)而對(duì)qw產(chǎn)生的影響是強(qiáng)烈耦合在一起的,不能簡(jiǎn)單疊加.2.2.2節(jié)中指出,這是因?yàn)樵谡駝?dòng)激發(fā)通過(guò)斜激波使邊界層外流發(fā)生變化時(shí),會(huì)極大地干擾振動(dòng)激發(fā)通過(guò)改變邊界層內(nèi)流動(dòng)而對(duì)壁面熱流qw產(chǎn)生的影響大小.

    圖8 振動(dòng)激發(fā)對(duì)氣動(dòng)力/熱特性影響的路徑圖Fig.8 Paths of vibration e ff ects on aerodynamic and aerothermodynamic characteristics

    2.3.2 直接求解的結(jié)果

    為了驗(yàn)證以上兩步漸進(jìn)分析方法的合理性,對(duì)同樣的高超聲速黏性斜劈問(wèn)題采用1.4節(jié)中給出的數(shù)值方法進(jìn)行直接求解.圖9是完全氣體和考慮振動(dòng)激發(fā)的斜劈流場(chǎng)溫度云圖的對(duì)比圖.從圖中可見(jiàn),考慮振動(dòng)激發(fā)斜激波波后溫度從完全氣體的2143K降低到1844K,變化顯著.另外,與完全氣體進(jìn)行對(duì)比,考慮振動(dòng)激發(fā)后激波層厚度明顯變薄,激波角明顯變小.這些變化規(guī)律與2.1節(jié)的結(jié)果吻合.

    距離斜劈前緣x=150mm處,完全氣體和考慮振動(dòng)激發(fā)兩種氣體模型的壁面壓力系數(shù)Cp、摩阻τw和熱流 qw的計(jì)算結(jié)果如表 4.可見(jiàn),對(duì)比完全氣體的結(jié)果,考慮振動(dòng)激發(fā)后壁面壓力系數(shù)Cp變化-3.08%,摩阻變化3.45%,熱流變化-3.09%.將該結(jié)果與圖8給出的結(jié)果對(duì)比發(fā)現(xiàn),漸進(jìn)求解得到的振動(dòng)激發(fā)對(duì)高超聲速氣動(dòng)力/熱特性的影響大小與該直接求解得到的結(jié)果非常接近,誤差較小.這說(shuō)明以上提出的兩步漸進(jìn)分析方法及其結(jié)果是合理的,且該方法對(duì)在斜激波與邊界層組合結(jié)構(gòu)中分析振動(dòng)激發(fā)對(duì)壁面摩阻、熱流的影響規(guī)律以及估計(jì)其影響大小具有很大的工程應(yīng)用價(jià)值.

    圖9 完全氣體(a)和考慮振動(dòng)激發(fā)(b)的斜劈流場(chǎng)溫度云圖的對(duì)比圖Fig.9 Comparison of temperature contour of wedge of perfect gas(a)and vibration excitation(b)

    表4 斜劈的氣動(dòng)力/熱特性(x=150mm)Table 4 Aerodynamic and aerothermodynamic characteristics of wedge at x=150mm

    3 結(jié)論

    高超聲速氣流通過(guò)斜激波與邊界層的加熱誘導(dǎo)空氣分子振動(dòng)激發(fā),從而影響飛行器表面的氣動(dòng)力/熱特性.本文新提出一種理論求解斜激波與數(shù)值求解邊界層相結(jié)合的兩步漸進(jìn)方法,研究了振動(dòng)激發(fā)過(guò)程對(duì)二維斜劈的氣動(dòng)力/熱特性的影響規(guī)律,并分析了振動(dòng)激發(fā)效應(yīng)在激波層與邊界層內(nèi)的傳遞和干擾機(jī)理.

    研究結(jié)果表明,對(duì)比完全氣體,振動(dòng)激發(fā)使壁面摩阻升高,使壁面熱流降低.在斜激波中,振動(dòng)激發(fā)使波后溫度、壓力降低,密度、速度、馬赫數(shù)升高,激波角減小,單位雷諾數(shù)增大.其中,振動(dòng)激發(fā)對(duì)斜激波波后溫度、密度、馬赫數(shù)、激波角和單位雷諾數(shù)的影響較大.發(fā)生振動(dòng)激發(fā)后,斜激波波后雷諾數(shù)的增大會(huì)導(dǎo)致邊界層厚度減小,結(jié)合多個(gè)物理量的變化,如速度和溫度,共同對(duì)邊界層內(nèi)的摩擦阻力和氣動(dòng)熱產(chǎn)生影響.在邊界層流動(dòng)中,壁面熱流、摩阻隨外流速度、溫度和壓力而變化,表明邊界層的存在會(huì)傳遞斜激波的振動(dòng)激發(fā)效應(yīng)而對(duì)壁面熱流和摩阻產(chǎn)生影響.另外,振動(dòng)激發(fā)對(duì)壁面熱流、摩阻的影響大小也隨外流速度和溫度而變化,說(shuō)明邊界層內(nèi)的振動(dòng)激發(fā)效應(yīng)會(huì)與斜激波的振動(dòng)激發(fā)效應(yīng)發(fā)生相互干擾.采用漸進(jìn)方法分析發(fā)現(xiàn),振動(dòng)激發(fā)通過(guò)在斜激波中改變邊界層外流而對(duì)壁面熱流產(chǎn)生的影響與通過(guò)改變邊界層內(nèi)流動(dòng)而產(chǎn)生的影響是強(qiáng)烈耦合的,不能簡(jiǎn)單疊加.而在對(duì)壁面摩阻的影響中,耦合作用較弱,基本是解耦的.對(duì)高超聲速黏性斜劈問(wèn)題的直接求解結(jié)果驗(yàn)證了提出的兩步漸進(jìn)分析方法的合理性.

    1姜宗林.觸摸高溫氣體動(dòng)力學(xué).力學(xué)與實(shí)踐,2006,28(5):1-7(Jiang Zonglin.Feeling high temperature gas dynamics.Mechanics in Engineering,2006,28(5):1-7(in Chinese))

    2樊菁.高超聲速高溫氣體效應(yīng)判據(jù).力學(xué)學(xué)報(bào),2010,42(4):591-596(Fan Jing.Criteria on high-temperature gas e ff ects around hypersonic vehicles.Chinese Journal of Theoretical and Applied Mechanics,2010,42(4):591-596(in Chinese))

    3 Anderson JD.Hypersonic and High-Temperature Gas Dynamics,2nd Edition.AIAA,2006

    4童秉綱,孔祥言,鄧國(guó)華.氣體動(dòng)力學(xué).第2版.北京:高等教育出版社,2012(Tong Binggang,Kong Xiangyan,Deng Guohua.Gas Dynamics.2nd Edition.Beijing:Higher Education Press,2012(in Chinese))

    5 Hirschel EH,Weiland C.Selected Aerothermodynamic Design Problems of Hypersonic Flight Vehicles.Berlin:Springer,2009

    6 Compton HR,Schiess JR,Suit WT,et al.Stability and control over the supersonic and hypersonic speed range.Conference Paper,NASA Langley,1983

    7 Brauckmann GJ,Paulson JW.Experimental and computational analysis of shuttle orbiter hypersonic trim anomaly.Journal of Spacecraft and Rockets,1995,32(5):758-764

    8 Muylaert J,Walpot L,Rostand P,et al.Extrapolation from wind tunnel to flight shuttle orbiter aerodynamics.Technical Report,NASA Langley,1998

    9 Griffith BJ,Maus JR,Best JT.Explanation of the hypersonic longitudinal stability problem:lessons learned.NASA,1983

    10 Romere PO,Young JC.Space shuttle entry longitudinal aerodynamic comparisons of fligh 2 with prefligh predictions.Journal of Spacecraft and Rockets,1983,20(6):518-523

    11 Calloway RL.Real-gas simulation for the shuttle orbiter and planetary entry configuration including fligh results.AIAA,1984

    12 Jiang ZL,Yu HR.Experimental and development of the longtest-duration hypervelocity detonation-driven shock tunnel(LHDst).AIAA,2014

    13姜宗林,李進(jìn)平,趙偉等.長(zhǎng)試驗(yàn)時(shí)間爆轟驅(qū)動(dòng)激波風(fēng)洞技術(shù)研究.力學(xué)學(xué)報(bào),2012,44(5):824-831(Jiang Zonglin,Li Jinping,Zhao Wei,et al.Investigating into techniques for extending the test-duration of detonation-driven shock tunnels.Chinese Journal of Theoretical and Applied Mechanics,2012,44(5):824-831(in Chinese))

    14劉云峰,汪運(yùn)鵬,苑朝凱等.復(fù)現(xiàn)高超聲速飛行條件下10?尖錐標(biāo)模氣動(dòng)力特性試驗(yàn)研究.中國(guó)力學(xué)大會(huì),上海,2015(Liu Yunfeng,Wang Yunpeng,Yuan Chaokai,et al.Aerodynamic force and moment characteristics of 10 degree sharp cone at Mach 7.0 in JF12 shock tunnel.The Chinese Congress of Theoretical and Applied Mechanics,Shanghai,2015(in Chinese))

    15汪運(yùn)鵬,劉云峰,苑朝凱等.長(zhǎng)試驗(yàn)時(shí)間激波風(fēng)洞測(cè)力技術(shù)研究.力學(xué)學(xué)報(bào),2016,48(3):545-556(Wang Yunpeng,Liu Yunfeng,Yuan Chaokai,et al.Study on force measurement in long-test duration shock tunnel.Chinese Journal of Theoretical and Applied Mechanics,2016,48(3):545-556(in Chinese))

    16汪志誠(chéng).熱力學(xué)·統(tǒng)計(jì)物理.第4版.北京:高等教育出版社,2008(Wang Zhicheng.Thermodynamics&Statistical Physics.4th Edition.Beijing:Higher Education Press,2008(in Chinese))

    17張德良.計(jì)算流體力學(xué)教程.北京:高等教育出版社,2010(Zhang Deliang.A Course in Computational Fluid Dynamics.Beijing:Higher Education Press,2010(in Chinese))

    18任玉新,陳海昕.計(jì)算流體力學(xué)基礎(chǔ).北京:清華大學(xué)出版社,2006(Ren Yuxin,Chen Haixin.Fundamentals of Computational Fluid Dynamics.Beijing:Tsinghua University Press,2006(in Chinese))

    19卞蔭貴,徐立功.氣動(dòng)熱力學(xué).第2版.合肥:中國(guó)科學(xué)技術(shù)大學(xué)出版社,2010(Bian Yingui,Xu Ligong.Aerothermodynamics.2nd Edition.Hefei:University of Science and Technology of China Press,2010(in Chinese))

    20 Seokkwan Y,Antony J.Lower-upper symmetric-gauss-seidel method for the Euler and Navier-Stokes equations.AIAA Journal,1988,26(9):1025-1026

    21 KimKH,KimC,RhoOH.Methodsfortheaccuratecomputationsof hypersonic fl ws I.AUSMPW+scheme.Journal of ComputationalPhysics,2001,174(1):38-80

    22 Kim KH,Lee JH,Rho OH.An improvement of AUSM schemes by introducing the pressure-based weight functions.Computers&Fluids,1988,27(3):311-346

    23 Keiichi K,Eiji S,Yoshiaki N,et al.Evaluation of Euler flu for hypersonic heating computations.AIAA Journal,2010,48(4):763-776

    24 Ren Yuxin.A robust shock-capturing scheme based on rotated Riemann solvers.Computers&Fluid,2003,32(10):1379-1403

    25 Van Leer B.Towards to the ultimate conservation di ff erence schemes:V.a second-order sequel to Godunov method.Journal of Computational Physics,1979,32(1):101-136

    26 Collela P.A direct Eulerian MUSCL scheme for gas dynamics.SIAM Journal on Scientifi and Statistical Computing,1985,6(1):104-117

    EFFECT OF VIBRATION EXCITATION ON HYPERSONIC AERODYNAMIC AND AEROTHERMODYNAMIC1)

    Zhang Zijian2)Liu Yunfeng Jiang Zonglin
    (State Key Laboratory of High-Temperature Gas Dynamics,Institute of Mechanics,CAS,Beijing 100190,China)(University of Chinese Academy of Sciences,Beijing 100049,China)

    With the increasing of fligh Mach number,the high-temperature gas e ff ect of air has becoming remarkable,which has significan impacts on the aerodynamics and aerothermodynamic characteristics of hypersonic vehicles.Because of the complex mechanism and numerous key parameters of high-temperature gas e ff ect,it has not been fully studied at home and abroad.When the high-temperature gas e ff ect occurs,multiple nonlinear physical processes are coupled together.However,ground tests and numerical simulations can not decouple these processes and can not explain the key physical mechanisms.To solve this problem,a new two-step asymptotic approximation method combining theoretic analysis and numerical simulation is proposed.In this method,the oblique shock relation with vibration excitation e ff ect is obtained by Newton iterative method,then the results are used as the boundary conditions of the boundary layer andit is solved numerically.By using this method,the e ff ect of vibration excitation on the aerodynamics and aerothermodynamic characteristics of a two dimensional wedge is studied.The results show that,the vibration excitation process has great e ff ect on the shock angle,the temperature,density,Mach number,and Reynolds number behind the oblique shock,but little influenc on the pressure and velocity.The inviscid fl w behind the oblique shock is coupled together with the boundary layer fl w.The changes of multiple physical quantities,including the increase of velocity and the decrease of the temperature behind the oblique shock,and the decrease of the boundary layer thickness due to the increase of the Reynolds number,have an e ff ect on the friction and aerodynamic heating in the boundary layer.Comparing with perfect gas model,vibration excitation increases the wall friction and decreases the wall heat flu of the wedge.By influencin the shock layer and the boundary layer respectively,the e ff ects of vibration excitation on heat flu are strong coupled,while they are weak coupled on friction.

    vibration excitation,aerodynamic and aerothermodynamic,boundary layer,high-temperature gas e ff ect

    V211.22

    :A

    10.6052/0459-1879-16-307

    2016–11–02 收稿,2017–01–24 錄用,2017–01–24 網(wǎng)絡(luò)版發(fā)表.

    1)國(guó)家自然科學(xué)基金資助項(xiàng)目(11672312,11532014).

    2)張子健,博士研究生,主要研究方向:高溫氣體動(dòng)力學(xué).E-mail:zhangzijian@imech.ac.cn

    張子健,劉云峰,姜宗林,振動(dòng)激發(fā)對(duì)高超聲速氣動(dòng)力/熱影響.力學(xué)學(xué)報(bào),2017,49(3):616-626

    Zhang Zijian,Liu Yunfeng,Jiang Zonglin.E ff ect of vibration excitation on hypersonic aerodynamic and aerothermodynamic.Chinese Journal of Theoretical and Applied Mechanics,2017,49(3):616-626

    猜你喜歡
    摩阻氣動(dòng)力邊界層
    飛行載荷外部氣動(dòng)力的二次規(guī)劃等效映射方法
    基于HIFiRE-2超燃發(fā)動(dòng)機(jī)內(nèi)流道的激波邊界層干擾分析
    市政橋梁預(yù)應(yīng)力管道摩阻系數(shù)測(cè)試研究
    江西建材(2018年4期)2018-04-10 12:37:20
    側(cè)風(fēng)對(duì)拍動(dòng)翅氣動(dòng)力的影響
    一類具有邊界層性質(zhì)的二次奇攝動(dòng)邊值問(wèn)題
    非特征邊界的MHD方程的邊界層
    高速鐵路接觸線覆冰后氣動(dòng)力特性的風(fēng)洞試驗(yàn)研究
    計(jì)算隱式摩阻系數(shù)方程數(shù)值解的簡(jiǎn)便方法
    考慮扶正器影響的套管摩阻計(jì)算方法研究
    降低壓裂施工摩阻技術(shù)研究
    自拍偷自拍亚洲精品老妇| 亚洲三级黄色毛片| 国产亚洲一区二区精品| 午夜老司机福利剧场| 国产成人a∨麻豆精品| 69人妻影院| 亚洲国产欧美人成| 一级片'在线观看视频| 男女边吃奶边做爰视频| 又粗又硬又长又爽又黄的视频| 亚洲av成人精品一二三区| 久久精品国产自在天天线| 亚洲精品成人久久久久久| 国产一区二区三区av在线| 午夜福利在线在线| 国产极品天堂在线| 成人毛片60女人毛片免费| 亚洲电影在线观看av| 看免费成人av毛片| 看黄色毛片网站| 成人一区二区视频在线观看| 中文在线观看免费www的网站| 高清在线视频一区二区三区| 国产爽快片一区二区三区| 亚洲成人一二三区av| 精品久久久精品久久久| 国产免费一区二区三区四区乱码| 男男h啪啪无遮挡| 日本爱情动作片www.在线观看| 嫩草影院新地址| 国产黄片视频在线免费观看| 午夜亚洲福利在线播放| 男女国产视频网站| 亚洲精品国产色婷婷电影| 男女国产视频网站| 午夜免费鲁丝| 超碰av人人做人人爽久久| 久久久久久国产a免费观看| 日韩三级伦理在线观看| 久久久久久伊人网av| 久久久久久久久大av| 国产精品一及| 日日撸夜夜添| 亚洲av福利一区| 亚洲欧洲日产国产| 色视频在线一区二区三区| 久久精品久久精品一区二区三区| 国产免费视频播放在线视频| 卡戴珊不雅视频在线播放| 亚洲精品国产av成人精品| 亚洲欧美清纯卡通| 国产人妻一区二区三区在| 热re99久久精品国产66热6| 国产欧美日韩一区二区三区在线 | 午夜激情久久久久久久| 啦啦啦在线观看免费高清www| 亚洲精品乱久久久久久| 亚洲自偷自拍三级| 国产精品三级大全| 国产精品爽爽va在线观看网站| 黄色怎么调成土黄色| 少妇猛男粗大的猛烈进出视频 | 91狼人影院| 看免费成人av毛片| 寂寞人妻少妇视频99o| 国产乱人视频| 国产黄色免费在线视频| 日韩中字成人| 99热网站在线观看| 高清日韩中文字幕在线| 久久人人爽人人爽人人片va| 日本黄色片子视频| 成年人午夜在线观看视频| 亚洲人与动物交配视频| 又爽又黄无遮挡网站| 国产爽快片一区二区三区| 免费黄网站久久成人精品| 天堂网av新在线| 色5月婷婷丁香| 精品酒店卫生间| 精品熟女少妇av免费看| 我的老师免费观看完整版| 国产午夜精品久久久久久一区二区三区| 五月伊人婷婷丁香| 亚洲美女视频黄频| 各种免费的搞黄视频| 自拍欧美九色日韩亚洲蝌蚪91 | kizo精华| 免费看光身美女| h日本视频在线播放| 全区人妻精品视频| 国产成人精品久久久久久| 2022亚洲国产成人精品| 日韩精品有码人妻一区| 亚洲精品色激情综合| av一本久久久久| 99九九线精品视频在线观看视频| 91在线精品国自产拍蜜月| 男插女下体视频免费在线播放| 狂野欧美白嫩少妇大欣赏| 成年免费大片在线观看| 亚洲天堂国产精品一区在线| 男女无遮挡免费网站观看| 久久这里有精品视频免费| 不卡视频在线观看欧美| 国产午夜福利久久久久久| 国产精品偷伦视频观看了| 成人国产av品久久久| 99久国产av精品国产电影| 一本一本综合久久| 男人狂女人下面高潮的视频| 又爽又黄a免费视频| 中文资源天堂在线| 夜夜看夜夜爽夜夜摸| 亚洲精品,欧美精品| 亚洲伊人久久精品综合| 久久影院123| 好男人视频免费观看在线| 丝袜脚勾引网站| 日本三级黄在线观看| 国产伦精品一区二区三区四那| 成人免费观看视频高清| 欧美亚洲 丝袜 人妻 在线| 99久久精品热视频| 一级毛片电影观看| 中国国产av一级| 亚洲高清免费不卡视频| a级毛色黄片| 一级爰片在线观看| 日韩,欧美,国产一区二区三区| 高清午夜精品一区二区三区| 国产午夜精品久久久久久一区二区三区| 亚洲国产欧美在线一区| 国精品久久久久久国模美| 水蜜桃什么品种好| 熟妇人妻不卡中文字幕| 成人午夜精彩视频在线观看| 亚洲av二区三区四区| 日本与韩国留学比较| 亚洲av男天堂| 免费大片18禁| 亚洲四区av| 成年av动漫网址| 人人妻人人爽人人添夜夜欢视频 | 69人妻影院| 亚洲aⅴ乱码一区二区在线播放| 国产人妻一区二区三区在| 久久精品国产亚洲av天美| 国产精品人妻久久久久久| 欧美日韩一区二区视频在线观看视频在线 | 中文资源天堂在线| 人人妻人人澡人人爽人人夜夜| 永久免费av网站大全| 亚洲精品成人久久久久久| 丰满乱子伦码专区| 欧美日韩国产mv在线观看视频 | 又爽又黄无遮挡网站| 国产淫语在线视频| 久久精品国产亚洲av涩爱| 成人毛片60女人毛片免费| 久久精品国产自在天天线| 男女国产视频网站| 大香蕉97超碰在线| 亚洲人成网站在线播| kizo精华| 永久网站在线| 精品久久久久久电影网| 国产精品女同一区二区软件| 男人和女人高潮做爰伦理| 亚洲精品乱码久久久v下载方式| 97超碰精品成人国产| 精品熟女少妇av免费看| 亚洲av电影在线观看一区二区三区 | 在线观看美女被高潮喷水网站| 日本一二三区视频观看| 亚洲精品456在线播放app| 成人美女网站在线观看视频| 搞女人的毛片| 日本免费在线观看一区| 欧美bdsm另类| 我要看日韩黄色一级片| 一本一本综合久久| 日本三级黄在线观看| 欧美人与善性xxx| 午夜福利高清视频| 最近最新中文字幕免费大全7| 亚洲色图综合在线观看| 久久精品久久久久久久性| 另类亚洲欧美激情| 免费观看在线日韩| 亚洲国产精品专区欧美| 亚洲av电影在线观看一区二区三区 | 欧美3d第一页| 直男gayav资源| 久久精品国产亚洲av涩爱| 天天躁夜夜躁狠狠久久av| xxx大片免费视频| 大片免费播放器 马上看| 国产精品久久久久久精品古装| 亚洲丝袜综合中文字幕| 五月开心婷婷网| freevideosex欧美| 亚洲av福利一区| 插阴视频在线观看视频| 女人被狂操c到高潮| 亚洲天堂国产精品一区在线| 免费看光身美女| 午夜福利在线在线| 国产亚洲91精品色在线| 国产高清有码在线观看视频| 国产视频内射| 免费看不卡的av| 91午夜精品亚洲一区二区三区| 国产精品人妻久久久影院| 久久午夜福利片| 女人被狂操c到高潮| 国产成人一区二区在线| 成人毛片a级毛片在线播放| 草草在线视频免费看| 久久97久久精品| 亚洲成人一二三区av| 赤兔流量卡办理| 国产成人精品久久久久久| 亚州av有码| 久久午夜福利片| 精品99又大又爽又粗少妇毛片| 尾随美女入室| 一区二区三区精品91| 久久久久久久午夜电影| 国产在视频线精品| videos熟女内射| 在线天堂最新版资源| 亚洲欧洲国产日韩| 毛片一级片免费看久久久久| 欧美成人精品欧美一级黄| 尾随美女入室| 国产精品国产三级专区第一集| 国产国拍精品亚洲av在线观看| 国产亚洲精品久久久com| 亚洲精品一区蜜桃| 色网站视频免费| 九九久久精品国产亚洲av麻豆| 91久久精品国产一区二区成人| 18禁动态无遮挡网站| 精品酒店卫生间| 精品人妻视频免费看| 国产av国产精品国产| 欧美日韩亚洲高清精品| 国产精品国产av在线观看| 国内少妇人妻偷人精品xxx网站| 大片电影免费在线观看免费| 男人和女人高潮做爰伦理| 涩涩av久久男人的天堂| 国产老妇女一区| 亚洲怡红院男人天堂| 97热精品久久久久久| 亚洲精品国产av蜜桃| 成人高潮视频无遮挡免费网站| 街头女战士在线观看网站| 嫩草影院精品99| 秋霞在线观看毛片| 欧美高清成人免费视频www| 欧美xxⅹ黑人| 国产男女内射视频| 2021少妇久久久久久久久久久| 国产高清三级在线| 熟女电影av网| 国产亚洲91精品色在线| 国产老妇伦熟女老妇高清| 国产精品偷伦视频观看了| av.在线天堂| 成人亚洲精品av一区二区| 老女人水多毛片| 国产永久视频网站| 日韩欧美精品免费久久| 熟女av电影| 日韩大片免费观看网站| 亚洲国产成人一精品久久久| 中文字幕av成人在线电影| 久久6这里有精品| 男女无遮挡免费网站观看| 亚洲精品成人久久久久久| 老司机影院成人| 成年版毛片免费区| 久久99热这里只频精品6学生| 国产综合懂色| 国产精品秋霞免费鲁丝片| 伊人久久精品亚洲午夜| 午夜视频国产福利| 少妇猛男粗大的猛烈进出视频 | 丝瓜视频免费看黄片| 欧美成人午夜免费资源| 国产有黄有色有爽视频| 一个人观看的视频www高清免费观看| 午夜视频国产福利| 国产男女内射视频| 一级爰片在线观看| 我的女老师完整版在线观看| 九九爱精品视频在线观看| 亚洲真实伦在线观看| 男人和女人高潮做爰伦理| 久久久久久久久久久免费av| 亚洲精品久久午夜乱码| 日韩成人伦理影院| 狂野欧美激情性xxxx在线观看| 中文在线观看免费www的网站| 深夜a级毛片| 水蜜桃什么品种好| 欧美三级亚洲精品| 一本久久精品| 欧美最新免费一区二区三区| 色播亚洲综合网| 色5月婷婷丁香| 亚洲内射少妇av| 高清在线视频一区二区三区| 少妇丰满av| a级毛色黄片| 91久久精品电影网| 精品一区二区三区视频在线| 日本-黄色视频高清免费观看| 91精品伊人久久大香线蕉| 亚洲无线观看免费| 舔av片在线| 中文字幕久久专区| 亚洲国产欧美在线一区| 在线观看人妻少妇| 亚洲精品,欧美精品| 国产一区亚洲一区在线观看| 91午夜精品亚洲一区二区三区| av福利片在线观看| 国产永久视频网站| freevideosex欧美| 大码成人一级视频| av国产久精品久网站免费入址| 亚洲不卡免费看| 精品酒店卫生间| 亚洲天堂av无毛| 亚洲婷婷狠狠爱综合网| 看十八女毛片水多多多| 一区二区三区四区激情视频| 内射极品少妇av片p| 一级毛片电影观看| 18+在线观看网站| 麻豆国产97在线/欧美| tube8黄色片| 少妇人妻久久综合中文| 伊人久久国产一区二区| 日本wwww免费看| 中文字幕亚洲精品专区| av国产免费在线观看| 国产69精品久久久久777片| 午夜日本视频在线| 免费在线观看成人毛片| 久久99精品国语久久久| 性色av一级| 国产极品天堂在线| 麻豆乱淫一区二区| av在线观看视频网站免费| 狂野欧美激情性bbbbbb| 亚洲精品亚洲一区二区| 秋霞在线观看毛片| 汤姆久久久久久久影院中文字幕| 一级黄片播放器| 国产永久视频网站| 亚洲国产色片| 亚洲一区二区三区欧美精品 | 一级毛片黄色毛片免费观看视频| 国产亚洲av嫩草精品影院| 欧美区成人在线视频| 天天躁日日操中文字幕| 久久女婷五月综合色啪小说 | 久久久欧美国产精品| 亚洲欧洲日产国产| 18+在线观看网站| 综合色av麻豆| 91精品一卡2卡3卡4卡| 91久久精品国产一区二区三区| 国产探花在线观看一区二区| 亚洲内射少妇av| 日韩一区二区视频免费看| 日韩三级伦理在线观看| 免费少妇av软件| 97超视频在线观看视频| 免费少妇av软件| 国产成人精品婷婷| 人妻 亚洲 视频| 97在线视频观看| 80岁老熟妇乱子伦牲交| 久久ye,这里只有精品| 国产又色又爽无遮挡免| av在线天堂中文字幕| 看黄色毛片网站| 狂野欧美激情性xxxx在线观看| 免费av观看视频| 男女国产视频网站| 中国美白少妇内射xxxbb| 欧美性感艳星| 日韩国内少妇激情av| 一区二区三区四区激情视频| .国产精品久久| 一级毛片我不卡| 精品国产一区二区三区久久久樱花 | av免费在线看不卡| 亚洲人成网站高清观看| 少妇的逼好多水| 日韩欧美一区视频在线观看 | 18禁裸乳无遮挡动漫免费视频 | 嫩草影院精品99| 免费看不卡的av| 搡老乐熟女国产| 91精品一卡2卡3卡4卡| 偷拍熟女少妇极品色| 特级一级黄色大片| 亚洲欧美日韩卡通动漫| av福利片在线观看| 老司机影院成人| 久久热精品热| 女人被狂操c到高潮| 欧美日韩视频精品一区| 女的被弄到高潮叫床怎么办| 国产伦精品一区二区三区视频9| 99久久人妻综合| 日日撸夜夜添| 黄色欧美视频在线观看| av一本久久久久| 日韩在线高清观看一区二区三区| 在线观看美女被高潮喷水网站| 国产精品99久久久久久久久| 亚洲,一卡二卡三卡| 国产黄片美女视频| 春色校园在线视频观看| 国内少妇人妻偷人精品xxx网站| 51国产日韩欧美| 免费看a级黄色片| 嘟嘟电影网在线观看| 最新中文字幕久久久久| 特大巨黑吊av在线直播| 欧美日韩精品成人综合77777| 国产爽快片一区二区三区| 国产精品久久久久久精品电影| 国产高清三级在线| 夫妻午夜视频| 欧美人与善性xxx| 欧美zozozo另类| 亚洲经典国产精华液单| 高清欧美精品videossex| 男的添女的下面高潮视频| 麻豆乱淫一区二区| 日韩精品有码人妻一区| 男女那种视频在线观看| 欧美另类一区| 国产探花极品一区二区| 一区二区三区免费毛片| 三级国产精品欧美在线观看| 少妇 在线观看| 国产精品偷伦视频观看了| 亚洲av不卡在线观看| 菩萨蛮人人尽说江南好唐韦庄| 国产色爽女视频免费观看| 又粗又硬又长又爽又黄的视频| a级一级毛片免费在线观看| 亚洲av电影在线观看一区二区三区 | 国产成人精品婷婷| 尾随美女入室| 少妇猛男粗大的猛烈进出视频 | 天堂俺去俺来也www色官网| 日本色播在线视频| 国产日韩欧美亚洲二区| 有码 亚洲区| 熟女av电影| 在线观看免费高清a一片| 亚洲电影在线观看av| 简卡轻食公司| 中文字幕av成人在线电影| 亚洲成人一二三区av| 免费播放大片免费观看视频在线观看| 国产高清三级在线| 中国国产av一级| 国产精品一区二区三区四区免费观看| 午夜福利视频1000在线观看| 午夜福利在线在线| 国产91av在线免费观看| xxx大片免费视频| 亚洲欧美成人精品一区二区| 97在线人人人人妻| 精品午夜福利在线看| 另类亚洲欧美激情| 在线观看免费高清a一片| 国产精品一区二区性色av| 国产老妇伦熟女老妇高清| 免费高清在线观看视频在线观看| 国产91av在线免费观看| 免费不卡的大黄色大毛片视频在线观看| 国产在线一区二区三区精| 菩萨蛮人人尽说江南好唐韦庄| 中文字幕免费在线视频6| 亚洲一区二区三区欧美精品 | 天堂中文最新版在线下载 | www.av在线官网国产| 能在线免费看毛片的网站| 国产永久视频网站| 成年av动漫网址| 久久久国产一区二区| 国产毛片在线视频| 大片免费播放器 马上看| 国产免费一级a男人的天堂| 欧美3d第一页| 日韩大片免费观看网站| 91狼人影院| 18禁在线无遮挡免费观看视频| 1000部很黄的大片| 热re99久久精品国产66热6| 国产精品久久久久久精品电影| 精品久久久久久久人妻蜜臀av| 青青草视频在线视频观看| 一级毛片黄色毛片免费观看视频| 国产视频内射| 日日撸夜夜添| 色5月婷婷丁香| 最近手机中文字幕大全| 黄片无遮挡物在线观看| 免费黄频网站在线观看国产| 丰满人妻一区二区三区视频av| 丝袜喷水一区| 波野结衣二区三区在线| 国产一区亚洲一区在线观看| 国产乱人偷精品视频| 亚洲欧洲国产日韩| 欧美极品一区二区三区四区| 日韩人妻高清精品专区| 精品一区二区三区视频在线| 边亲边吃奶的免费视频| 国产淫片久久久久久久久| 免费av不卡在线播放| videossex国产| 高清午夜精品一区二区三区| 亚洲成人久久爱视频| 3wmmmm亚洲av在线观看| 久久这里有精品视频免费| 亚洲色图综合在线观看| 三级国产精品欧美在线观看| 18禁动态无遮挡网站| 少妇的逼好多水| 黄色配什么色好看| 亚洲欧美中文字幕日韩二区| 最新中文字幕久久久久| 国产精品国产av在线观看| 大片免费播放器 马上看| 亚洲精品日韩在线中文字幕| 又爽又黄a免费视频| 99九九线精品视频在线观看视频| 精品熟女少妇av免费看| 一级av片app| 美女国产视频在线观看| 国产成人一区二区在线| 毛片一级片免费看久久久久| 22中文网久久字幕| 亚洲av在线观看美女高潮| 成人特级av手机在线观看| 人人妻人人澡人人爽人人夜夜| 水蜜桃什么品种好| www.色视频.com| 男女无遮挡免费网站观看| 国产高清国产精品国产三级 | 国产探花在线观看一区二区| 欧美另类一区| 中国美白少妇内射xxxbb| 深爱激情五月婷婷| 欧美xxⅹ黑人| 中文字幕制服av| 一级a做视频免费观看| 亚洲欧美日韩东京热| av国产久精品久网站免费入址| 免费av不卡在线播放| 国产亚洲5aaaaa淫片| 狂野欧美激情性bbbbbb| 亚洲精品久久久久久婷婷小说| 亚洲精品成人av观看孕妇| 2021天堂中文幕一二区在线观| 又爽又黄无遮挡网站| 3wmmmm亚洲av在线观看| 精品亚洲乱码少妇综合久久| 亚洲国产欧美在线一区| 人人妻人人看人人澡| 日本爱情动作片www.在线观看| 三级经典国产精品| 热99国产精品久久久久久7| 亚洲一区二区三区欧美精品 | 特大巨黑吊av在线直播| 中文天堂在线官网| 97精品久久久久久久久久精品| 亚洲一区二区三区欧美精品 | 国内少妇人妻偷人精品xxx网站| 三级国产精品片| 国产一区二区三区av在线| 蜜桃久久精品国产亚洲av| 色吧在线观看| 久热这里只有精品99| 精品少妇久久久久久888优播| 国产免费一级a男人的天堂| av线在线观看网站| 久久人人爽av亚洲精品天堂 | 久久久精品94久久精品| 色网站视频免费| 国产精品久久久久久精品古装| 大片电影免费在线观看免费| 看十八女毛片水多多多| 国产69精品久久久久777片| 一级a做视频免费观看| 精品国产露脸久久av麻豆| 亚洲国产成人一精品久久久| 激情 狠狠 欧美| 欧美xxⅹ黑人| 国产在线一区二区三区精| 精品少妇久久久久久888优播|