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

    豎直平面激波與水平熱層作用后流場的理論計算方法研究*

    2020-01-02 06:20:22賈雷明
    爆炸與沖擊 2019年12期
    關(guān)鍵詞:激波計算結(jié)果流場

    賈雷明,田 宙

    (1. 清華大學(xué)工程物理系,北京 100084;2. 西北核技術(shù)研究院,陜西 西安 710024)

    在爆轟激波管[1]或強爆炸[2]等環(huán)境中,激波陣面或爆炸火球釋放出的熱輻射,先于激波到達固壁表面并發(fā)生能量沉積,致使壁面附近形成具有較高溫度的氣體層,稱為“熱層”(thermal layer)。由于熱層內(nèi)氣體聲速大于周圍空氣,熱層內(nèi)激波陣面會超越空氣中激波陣面?zhèn)鞑?,?dǎo)致流場中出現(xiàn)前驅(qū)波和渦旋等結(jié)構(gòu)。與不考慮熱層時的情形相比,固壁附近的流場參量也會發(fā)生改變。因此,激波與熱層作用常被應(yīng)用于流場診斷和控制[3-5]。

    為了認識激波與熱層作用機制,人們針對不同形式的激波與熱層作用,開展了大量的實驗和數(shù)值模擬工作[6-18]。在理論分析方面,已有的工作則主要關(guān)注豎直平面激波與水平熱層作用,在這類作用中,熱層高度是唯一的長度尺寸,隨著激波傳播距離遠大于熱層高度,流動會進入準(zhǔn)自相似階段,流場結(jié)構(gòu)近似呈比例發(fā)展,此時可假定流動定常,從而使理論分析簡化。

    最早,Griffith[6]針對熱層與周圍空氣物性差異較小的情形,在與入射激波固連的坐標(biāo)系中,假定流動定常,采用線性小擾動理論求解二維定常Euler 方程組,得到了激波陣面形狀和波后流場參量分布。針對熱層與空氣物性差異較大的情形,小擾動方法不再適用,因此人們又建立了新的理論方法,主要用于計算固壁附近的流場參量。Shreffler 等[7]將熱層高度近似為零,認為入射激波后的流體在一維簡單稀疏波的作用下膨脹加速,形成沿固壁運動的楔形“活塞”,驅(qū)動物質(zhì)界面偏折和前驅(qū)波傳播;在與楔形“活塞”固連的坐標(biāo)系內(nèi),假定流動定常,給定前驅(qū)波傳播方向,就可以求得楔形“活塞”內(nèi)流體的壓力和速度。但是,該理論方法不能計算熱層內(nèi)的激波強度,且需利用實驗或數(shù)值模擬確定前驅(qū)波傳播方向。Hess[19]和Mirels[20]等則認為,隨著流動進入準(zhǔn)自相似階段,熱層內(nèi)激波傳播速度應(yīng)與入射激波傳播速度相等,由此計算熱層內(nèi)激波的強度。為了計算固壁附近的流場參量,Mirels[20]進一步假定:(1)在與入射激波陣面固連的坐標(biāo)系中,入射激波后流體在定常等熵波作用下膨脹加速,形成沿固壁運動的“活塞”;(2)在與熱層內(nèi)激波陣面固連的坐標(biāo)系中,穿越激波后的熱層氣體在定常等熵波作用下減速增壓,直至與“活塞”內(nèi)流體壓力相等,而速度則滿足一定的線性比例關(guān)系。但是,速度比例系數(shù)依賴于具體的工況條件,需借助于實驗或數(shù)值模擬確定;同時,在流動準(zhǔn)自相似階段,熱層內(nèi)激波速度會趨近于某一值,而該值通常大于入射激波速度。

    綜上,Shreffler[7]、Hess[19]和Mirels[20]等的理論方法均基于某種流動定常假定,且由于未對熱層內(nèi)激波傳播進行詳細求解,這些方法大都依賴于特定的經(jīng)驗參數(shù),而經(jīng)驗參數(shù)又需要借助于實驗或數(shù)值模擬確定。本文中圍繞豎直平面激波與水平熱層作用,對Mirels[20]理論方法進行改進:(1)分析熱層內(nèi)激波傳播過程,然后基于激波動力學(xué)理論計算熱層內(nèi)激波傳播,從而舍棄“熱層內(nèi)激波速度與入射激波速度相等”的假定;(2)考慮到整個流場是在入射激波后流體而非入射激波的推動下發(fā)展,因此在與入射激波后流體固連的坐標(biāo)系中,假定入射激波后流體在定常等熵波作用下形成沿固壁運動的“活塞”;(3)“活塞”內(nèi)流體與其毗鄰的熱層氣體,應(yīng)滿足壓力和速度連續(xù),從而不再引入速度比例系數(shù)。本文中將通過與數(shù)值模擬結(jié)果、實驗數(shù)據(jù)和已有理論方法的結(jié)果進行對比,驗證上述改進的合理性。

    1 理論計算方法

    初始時刻,固壁附近布有性質(zhì)均勻的水平熱層,其高度為h,并與周圍空氣處于等壓靜止?fàn)顟B(tài),見圖1(a),圖中實線表示激波,短劃線表示物質(zhì)界面。在t=0 時刻,豎直平面激波I 到達熱層左側(cè)界面處,并與之發(fā)生作用,產(chǎn)生透射激波T 和反射波R。波T 和R 分別向物質(zhì)界面左右兩側(cè)傳播,見圖1(b),圖中點P1為波T 陣面與熱層上側(cè)界面的交點。記波I 和T 的傳播速度分別為DI、DT,由于DT>DI,當(dāng)波T 超越上方空氣中的波I 后,位于波T 后的高壓流體會向上方空氣中膨脹,驅(qū)動形成透射激波P,該透射激波P 被稱為“前驅(qū)波”,見圖1(c)。同時,由于激波與熱層界面作用存在斜壓效應(yīng),即 ?p×?ρ ≠0,流場中會有渦量產(chǎn)生并沿物質(zhì)界面沉積,最終在流場不穩(wěn)定性的影響下,形成渦旋。與不考慮熱層時的情形相比,流場演化更加復(fù)雜。

    1.1 熱層內(nèi)透射激波傳播

    在t=0 時刻,波T 和R 后的流場參量可借助于一維Riemann 問題的精確求解方案[21]給出。

    在t>0 時刻,波T 超越波I 沿?zé)釋由蟼?cè)界面?zhèn)鞑?。由于初始時刻波T 陣面與熱層上界面垂直,因此波T 沿?zé)釋由蟼?cè)界面的折射為非正規(guī)折射。本文中利用Whitham 提出的幾何激波動力學(xué)(geometrical shock dynamics,后文簡稱為GSD)理論[22-23],計算波T 沿?zé)釋由蟼?cè)界面的非正規(guī)折射。GSD 理論常用于近似求解二維/三維激波陣面的傳播。

    非正規(guī)折射的發(fā)生,是因為波后流場中的擾動追趕上了波T。記波T 獨自沿?zé)釋由蟼?cè)界面折射(即不考慮波I 的存在)時,點P1沿界面的移動速度為DP1。如果DP1>DI,則波T 沿界面的折射僅與波T 強度和界面兩側(cè)的物性參數(shù)有關(guān);如果DP1<DI,則波T 沿界面的折射還會受到波I 的影響。接下來,分別對這兩種情形進行討論。

    1.1.1 折射不受波I 影響

    此時,在點P1附近,熱層氣體經(jīng)過波T 后壓力增加,然后向上方空氣中膨脹,形成前驅(qū)波P 和反射稀疏波R1;稀疏波R1沿波T 陣面?zhèn)鞑ィ沟貌═ 強度減小、波陣面發(fā)生彎曲,記彎曲后的波陣面為波T′;在熱層上方,前驅(qū)波P 與波I 作用,形成激波R2、R3、T1和新的物質(zhì)界面,見圖2,圖中實線表示激波,短劃線表示物質(zhì)界面,點線表示稀疏波。

    圖2 波T 沿界面的非正規(guī)折射(DP1>DI)Fig.2 Irregular refraction of wave T at the horizontal material interface (DP1>DI)

    依據(jù)GSD 理論,建立波T 與T′之間的近似關(guān)系,

    式中:Ma為激波馬赫數(shù),θ 為激波傳播方向與水平方向的夾角,A=A(Ma)為射線管面積函數(shù),下標(biāo)T、T′分別表示波T 和T′的相關(guān)參數(shù)。本文中約定,所有的角度均以逆時針方向為正。

    在點P1附近,波T′后流場并不均勻,因此不能利用雙波或三波理論來求解點P1附近的流場??紤]到擾動是從界面下方的流場中產(chǎn)生的,然后跨越界面向上方流場傳播。當(dāng)穩(wěn)定的折射結(jié)構(gòu)形成時,在點P1附近,波P 與波T′陣面應(yīng)保持連續(xù),且流場不再發(fā)生變化。此時,界面下方流場中的擾動不會再跨越界面向上方傳播,否則波P 后流場將再次發(fā)生改變。因此,界面下方流場中擾動與波T′陣面交點的運動軌跡應(yīng)與波后物質(zhì)界面重合,即

    式中:αT′為擾動運動軌跡與水平方向的夾角,δT′為波T′后物質(zhì)界面與水平方向的夾角。根據(jù)GSD 理論,有

    式中:

    γ 為比熱比。建立與點P1固連的坐標(biāo)系,記為 ?1,根據(jù)斜激波關(guān)系式,有

    式中:ξ 為波T′后壓力與波前壓力之比,Ma0為坐標(biāo)系 ?1中波T′前流場馬赫數(shù)。波T′與P 陣面在點P1處保持連續(xù),且波后流場壓力相等,再聯(lián)立式(1)~(4),即可求得點P1附近波T′、P 后流場以及點P1沿界面的傳播速度DP1。

    在熱層上方,可利用激波關(guān)系式求解波P 與波I 作用后的流場分布。

    1.1.2 折射受波I 影響

    當(dāng)DP1<DI時,波T 沿界面的折射會受到熱層上方波I 的影響。此時,除了在界面下方流場有稀疏波產(chǎn)生外,在界面上方的流場中還產(chǎn)生了壓縮波,記為R4,它會沿著波I 陣面向上傳播,使得波I 陣面向前彎曲,形成前驅(qū)波P,見圖3。流動本身具有極強的非線性,隨著時間發(fā)展壓縮波R4會逐漸演化成激波,并與波I 和P 組成三波結(jié)構(gòu),類似于斜激波在固壁的馬赫反射。三波點附近的流場近似采用三波理論描述。

    圖3 波T 沿界面的非正規(guī)折射(DP1<DI)Fig.3 Irregular refraction of wave T at the horizontal material interface (DP1<DI)

    波T 與T′之間仍近似滿足式(1)。在點P1附近,由于折射受到波I 的影響,因此式(2)不再成立。本文中采用Catherasoo 等[24]建立的激波在非均勻介質(zhì)中的傳播理論,描述界面兩側(cè)波T′和波P 參數(shù)所滿足的關(guān)系,有

    式中:MaP為波P 的馬赫數(shù),c1、c2分別為波前未擾動空氣和熱層氣體的聲速,

    式中:DP為波P 傳播速度,θP為波P 傳播方向與水平方向夾角,其余物理量意義同式(3)。假定波P 陣面保持平直,聯(lián)立式(1)、(5)和三波理論,即可求得點P1附近波T′和P 后的流場。

    上述理論方法適用于波R1到達固壁之前的時刻。記波T 后流場聲速為cT,在t=h/cT時,波R1到達固壁并發(fā)生反射,形成新的反射稀疏波,稀疏波后激波陣面仍保持與固壁垂直,可利用式(1)求解激波參數(shù)。在t>h/cT時,新的反射稀疏波會沿著波T′陣面向上傳播。當(dāng)反射稀疏波到達熱層上側(cè)界面后,又會發(fā)生反射,流場可利用上述理論方法進行求解。新生成的反射稀疏波又會向著固壁運動。如果忽略流場中其他擾動的影響,這個過程會不斷持續(xù)下去,熱層內(nèi)激波陣面的構(gòu)型也會發(fā)生周期性變化,激波強度逐漸減小,見圖4,圖中t1~t5表示不同的時刻,點劃線表示稀疏波與波T 陣面交點軌跡,為了記述方便,將任意時刻熱層內(nèi)激波陣面統(tǒng)一記為波T。真實情況中,流場左側(cè)的擾動勢必會趕上波T,圖4 所示的過程并不會一直持續(xù)下去。

    圖4 熱層內(nèi)激波傳播Fig.4 Shock propagation with time in the thermal layer

    1.2 固壁附近流場

    當(dāng)激波走過的距離遠大于熱層高度h時,流動進入準(zhǔn)自相似階段。此時,流場參量可以近似看作是變量x/t和y/t的函數(shù)。將固壁附近流場劃分成五個區(qū)域,見圖5,其中區(qū)域1 是驅(qū)動渦旋發(fā)展的流體區(qū)域,區(qū)域2、3 是渦旋區(qū)域,渦旋與固壁之間的流體即為沿固壁運動的“活塞”,區(qū)域4 是經(jīng)過波T 壓縮后的熱層氣體,區(qū)域5 是處于未擾動狀態(tài)的熱層氣體。圖5 中,點P1為波T 與熱層上側(cè)物質(zhì)界面的交點,點P2~P5為區(qū)域1~5 的分界點。

    圖5 固壁附近流場區(qū)域劃分示意圖Fig.5 Illustration of flow field division near the wall

    在區(qū)域1 中,初始時刻波I 與界面作用形成反射波R,波R 沿固壁向界面左側(cè)傳播。由于pR≠pI,界面附近又會產(chǎn)生擾動分別向熱層上方和固壁傳播。這些擾動傳播的結(jié)果,是使得流場恢復(fù)到均衡狀態(tài),本文中近似選取波I 后流體狀態(tài)作為均衡狀態(tài)。

    在區(qū)域2、3 中,物質(zhì)界面與固壁組成了收縮-擴張管道,點P3處為收縮與擴張管道的連接處。區(qū)域1中的流體先經(jīng)過收縮管道進行加速減壓,再經(jīng)過擴張管道減速升壓后,到達點P4附近。建立與區(qū)域1中流體固連的坐標(biāo)系 ?2,在坐標(biāo)系 ?2中觀察區(qū)域2、3 處的流動,假定流動滿足準(zhǔn)一維定常絕熱條件,有

    式中:ρ 為密度,u為速度,下標(biāo)2、4L 分別表示點P2和點P4附近界面左側(cè)的流場參量。點P2附近流場參量即為區(qū)域1 中流場參量。

    在區(qū)域4 中,熱層氣體穿過波T 后,壓力和速度增加。波T 與固壁相交于點P5,建立與點P5固連的坐標(biāo)系 ?5。 在坐標(biāo)系 ?5中觀察,物質(zhì)界面、固壁與波T 組成擴張管道,假定管道內(nèi)流動滿足準(zhǔn)一維定常絕熱條件,則有

    式中:下標(biāo)4R、5 分別表示點P4附近界面右側(cè)、點P5附近的流場參量,D5為波T 沿固壁的運動速度。在點P4附近,界面兩側(cè)流場滿足連續(xù)條件,即

    Hess[19]認為流動從初始時刻到進入準(zhǔn)自相似階段的時間尺度與h和流體聲速之比具有相同的量級。同時考慮到熱層左側(cè)流場也應(yīng)達到均衡狀態(tài),本文中選取流動進入準(zhǔn)自相似階段的時間為

    式中:cR為波R 后流場聲速,N為常數(shù)。至此,根據(jù)上一小節(jié)的理論方法計算固壁附近波T 強度,聯(lián)立式(6)~(10)即可求解點P4附近的流場參量。

    2 結(jié)果與討論

    2.1 MaI=2.00,0.1≤ρtl/ρair≤0.9 時的計算結(jié)果

    給定波I 馬赫數(shù)MaI=2.00,波前流場壓力p0=0.100 MPa,空氣密度ρair=1.00 kg/m3,熱層高度h=0.50 m,空氣和熱層氣體比熱比均取γ=1.40。計算區(qū)域0<x<15.00 m 和0<y<10.00 m,熱層左側(cè)界面位于x=2.00 m 處。設(shè)定九種不同的工況,即熱層氣體密度ρtl分別取0.10、0.20、0.30、0.40、0.50、0.60、0.70、0.80、0.90 kg/m3,利用本文理論方法計算固壁附近流場參量。

    同時,利用動力學(xué)軟件ANSYS Autodyn Euler-FCT 求解器進行數(shù)值模擬。圖6 所示為三種工況條件下相應(yīng)時刻壓力p和水平速度u沿固壁分布的數(shù)值模擬結(jié)果,Grid 1 和Grid 2 分別表示不同的網(wǎng)格劃分方案,即dx=dy=0.01 m 和dx=dy=0.02 m,dx和dy分別為沿x和y方向的網(wǎng)格尺寸。兩種網(wǎng)格方案所得到的結(jié)果基本吻合,驗證了網(wǎng)格收斂性,后續(xù)的數(shù)值模擬均選用Grid 1。

    圖6 不同工況條件下壓力p 和水平速度u 沿固壁的分布Fig.6 Profiles of pressure and horizontal velocity along the wall in different cases

    2.1.1 流場中的激波結(jié)構(gòu)

    將DP1>DI和DP1<DI時流場中的激波結(jié)構(gòu)(見圖2~3)分別記為W1 和W2。兩者的區(qū)別主要體現(xiàn)為,界面上方流場中的激波結(jié)構(gòu)分別為五波和三波結(jié)構(gòu)。表1 為不同工況條件下,t<h/cT時刻熱層內(nèi)波T 相關(guān)參數(shù)的理論計算結(jié)果,pT、pT′分別為波T 以及點P1附近波T′后的壓力。對于MaI=2.00,當(dāng)ρtl=0.87 kg/m3時,DP1=DI=748.33 m/s。圖7 為t=0.5 ms 時刻壓力等值線圖的數(shù)值模擬結(jié)果,等值線間隔取為Δp=(pmax?pmin)/20,pmax、pmin為對應(yīng)工況和時刻的壓力最大值和最小值。此時,波R1未到達固壁。

    表1 流場中激波結(jié)構(gòu)類型Table 1 Wave structure types above material interface

    圖7 t=0.5 ms,不同工況下的壓力等值線圖Fig.7 Pressure contour lines at t =0.5 ms in different cases

    當(dāng)ρtl≤0.60 kg/m3時,界面上方為五波結(jié)構(gòu),理論計算結(jié)果與數(shù)值模擬結(jié)果基本吻合。受到波后流場中擾動的影響,波R2、R3和T1后的流場區(qū)域并不均勻。當(dāng)ρtl=0.70 kg/m3和0.80 kg/m3時,理論計算結(jié)果顯示界面上方為五波結(jié)構(gòu),而在數(shù)值模擬結(jié)果中并未明顯觀測到波R2。這是因為兩種工況處于波系類型由W1 向W2 的過渡區(qū)附近,根據(jù)理論方法得到的兩種工況中波R2馬赫數(shù)MaR2分別為1.20 和1.16,強度較弱。也可能是由于理論方法采用激波動力學(xué)理論近似描述波T 陣面彎曲,所得到的DP1略大于真實值,此時波系類型實為W2,理論計算結(jié)果與數(shù)值模擬結(jié)果存在偏差。

    當(dāng)ρtl=0.90 kg/m3時,界面上方的波系為三波結(jié)構(gòu),理論計算結(jié)果與數(shù)值模擬結(jié)果基本一致。波P 陣面略向前彎曲,理論計算得到波R4馬赫數(shù)MaR4=1.02,接近聲波。

    2.1.2 固壁附近流場參量

    選取工況ρtl=0.50 kg/m3,闡釋固壁附近流場特征及參量的演化。圖8 為t=13.0 ms 時刻流場密度云圖的數(shù)值模擬結(jié)果。空氣和熱層氣體經(jīng)過激波壓縮后,被卷入流場后方的渦旋中。隨著時間發(fā)展,波T 傳播距離大于14h,流動進入準(zhǔn)自相似階段。圖9 為t=10.0、12.0、14.0 和16.0 ms 時刻固壁附近壓力p和密度ρ 分布的數(shù)值模擬結(jié)果。熱層左側(cè)的流場,近似恢復(fù)到波I 后的狀態(tài),這與本文理論方法關(guān)于圖5 中區(qū)域1 的假定是一致的。點P4、P5附近的壓力和密度隨時間基本不變。在時間間隔Δt=2.0 ms 內(nèi),波T 的位移分別為1.61、1.65 和1.62 m,也隨時間基本不變。

    圖8 t=13.0 ms,流場密度云圖分布(ρtl=0.50 kg/m3)Fig.8 Density contour at t=13.0 ms with ρtl=0.50 kg/m3

    圖9 不同時刻,物理量沿固壁分布(ρtl=0.50 kg/m3)Fig.9 Parameters vs. x along the wall at different time instants with ρtl=0.50 kg/m3

    接下來,利用本文的理論方法,對不同工況下流場進入準(zhǔn)相似階段后點P4、P5附近的流場參量進行計算。其中,流場進入自相似階段的時刻按式(11)計算,并取N=2.5,假定波T 后稀疏波以波T 后流場聲速在熱層內(nèi)往返傳播。表2 為不同工況下物理量p4、u4、ρ4L、ρ4R和p5的理論計算結(jié)果與數(shù)值模擬結(jié)果,其中,TA 表示理論計算結(jié)果,NS 表示數(shù)值模擬結(jié)果,ε=|TA-NS|/NS 為理論結(jié)果與數(shù)值模擬結(jié)果的偏差,t*為數(shù)值模擬結(jié)果的取值時刻。t*的確定應(yīng)保證物理量趨近于漸近值,本文采用“在時間間隔Δt=0.5 ms 內(nèi)物理量ρ4L數(shù)值變化不超過2%”作為判斷標(biāo)準(zhǔn)。

    表2 不同工況條件下固壁附近流場參量Table 2 Parameter values near the rigid wall for different cases

    當(dāng)熱層不存在時,波I 后固壁附近流場的壓力為0.450 MPa。熱層的出現(xiàn),使得固壁附近激波峰值壓力降低,且峰值壓力隨ρtl減小而減小。對于不同的工況條件,流場參量的理論計算結(jié)果與數(shù)值模擬結(jié)果的偏差均小于10%,驗證了理論方法的合理性。

    2.2 與Shreffler 和Mirels 理論方法的比較

    選取與上一小節(jié)相同的工況參數(shù),分別利用Shreffler[7]理論方法和Mirels[20]理論方法對固壁附近流場參量p4、u4和p5進行計算,并與本文理論計算結(jié)果進行比較,見圖10。在Shreffler 理論方法中,前驅(qū)波陣面與水平方向的夾角α 根據(jù)sin2α=ρtl/ρair確定,該式的計算結(jié)果已被證實與實驗數(shù)據(jù)較為吻合[3,11]。在Mirels 理論方法中,速度比例系數(shù)k(k=u4L/u4R)描述的是點P4附近界面兩側(cè)速度之比,且0<k≤1.0。k值依賴于具體的工況條件,Mirels 根據(jù)其理論計算結(jié)果與實驗數(shù)據(jù)或數(shù)值模擬結(jié)果是否吻合來確定k值,但并沒有給出k與工況條件的具體關(guān)系。此處,初步選取k=1.0,0.3,0.01 進行計算,以涵蓋k可能的取值范圍。

    圖10 不同理論方法的計算結(jié)果(MaI=2.00)Fig.10 Results from various theoretical methods with MaI=2.00

    從圖10 可以發(fā)現(xiàn),Shreffler 理論方法僅能給出p4和u4,其中p4的結(jié)果與數(shù)值模擬結(jié)果偏差較大,u4的結(jié)果則與數(shù)值模擬結(jié)果較為接近。Mirels 理論方法給出的p4與數(shù)值模擬結(jié)果的偏差大小與k值有關(guān),當(dāng)k=0.3 時,ρth≈0.28 kg/m3的結(jié)果與數(shù)值模擬結(jié)果最接近,當(dāng)k=0.01 時,ρth≈0.70 kg/m3的結(jié)果與數(shù)值模擬結(jié)果最接近;對于任意的k值,u4的結(jié)果均與數(shù)值模擬結(jié)果存在較大偏差;隨著ρtl/ρair趨近于1,p5的結(jié)果與數(shù)值模擬結(jié)果偏差逐漸減小。本文理論方法所給出的p4、u4和p5,與數(shù)值模擬結(jié)果偏差均小于10%,表明本文理論方法在此類工況條件下優(yōu)于Shreffler 理論方法和Mirels 理論方法。

    2.3 Zaslavskii 實驗

    Zaslavskii[14]利用氮氣與氫氣的混合氣體構(gòu)建熱層,在激波管中開展了MaI=1.36 的豎直平面激波與水平熱層作用實驗,并測量了不同 ρtl/ρair(ρtl/ρair<0.5)條件下波T 沿固壁走過約12h后的強度ξT=pT/p0,按照2.1 節(jié)所采用的準(zhǔn)則判斷,流動進入準(zhǔn)自相似階段。Zaslavskii 并未詳細給出流場初始狀態(tài)參數(shù)p0、ρair等,但由于ξT是一個無量綱量且不依賴于p0、ρair等物理量的具體數(shù)值,因此本文在開展理論計算和數(shù)值模擬時,不失一般性,取p0=0.101 MPa、ρair=1.225 kg/m3、γ=1.4,且N=2.5。圖11 給出了ρtl/ρair=0.1、0.2、0.3、0.4 和0.5 時波T 強度ξT和固壁附近p4、u4的理論計算結(jié)果和數(shù)值模擬結(jié)果,并與Zaslavskii 測量得到的ξT進行比較。

    對于MaI=1.36 和ρtl/ρair≤0.5,由于入射激波速度小于熱層內(nèi)氣體聲速,Shreffler 理論方法和Mirels 理論方法均不再適用,因此圖11 中僅列出了本文理論方法的計算結(jié)果。圖11(a)為波T 強度ξT,對于不同的ρtl/ρair,數(shù)值模擬結(jié)果與實驗結(jié)果基本吻合,理論計算結(jié)果與數(shù)值模擬結(jié)果最大偏差約8.60%。圖11(b)、(c)分別為點P4附近的流場壓力p4和速度u4,理論計算結(jié)果與數(shù)值模擬結(jié)果較為吻合,其中p4對應(yīng)的最大偏差為13.8%,u4對應(yīng)的最大偏差為20.6%,隨著ρtl/ρair增大,偏差逐漸減小。

    圖11 不同理論方法的計算結(jié)果(MaI=1.36)Fig.11 Results from various theoretical methods with MaI=1.36

    3 結(jié) 論

    當(dāng)豎直平面激波與固壁附近水平熱層作用時,流場中會出現(xiàn)渦旋和前驅(qū)波等結(jié)構(gòu)。為了更為準(zhǔn)確地計算流動處于準(zhǔn)自相似階段時固壁附近的流場參量,本文對已有的Mirels[20]理論方法進行了以下三個方面的改進:

    (1)舍棄“熱層內(nèi)激波速度與入射激波速度相等”的假定,基于激波動力學(xué)理論計算熱層內(nèi)激波與波后稀疏波作用;

    (2)在與入射激波后流體而非入射激波陣面固連的坐標(biāo)系中,入射激波后流體在定常等熵波作用下,形成沿固壁運動的“活塞”;

    (3)不再引入速度比例系數(shù),而是假定“活塞”內(nèi)流體與其毗鄰的熱層氣體,滿足壓力和速度連續(xù)。

    設(shè)定不同的工況,利用本文中改進后的理論方法計算固壁附近流場參量,并與數(shù)值模擬結(jié)果、已有的實驗數(shù)據(jù)以及Shreffler 理論方法和Mirels 理論方法的計算結(jié)果進行對比。給定入射激波馬赫數(shù)2.00 和熱層與空氣密度之比0.10≤ρtl/ρair≤0.90,對于熱層上界面處擾動未到達固壁之前的流場激波類型,本文理論方法與數(shù)值模擬所給出的結(jié)果基本一致,但由于激波動力學(xué)理論自身存在一定近似,因此在激波類型發(fā)生轉(zhuǎn)變的區(qū)域附近,理論計算結(jié)果與數(shù)值模擬結(jié)果還存在一定的偏差;當(dāng)流動進入準(zhǔn)自相似階段后,本文理論方法得到的熱層內(nèi)激波強度以及物質(zhì)界面兩側(cè)流場壓力、速度和密度等流場參量,與數(shù)值模擬結(jié)果偏差均小于10%,明顯優(yōu)于Shreffler 理論方法和Mirels 理論方法。

    給定入射激波馬赫數(shù)1.36 和熱層與空氣密度之比0.10≤ρtl/ρair≤0.50,由于此時入射激波速度小于熱層內(nèi)氣體聲速,Shreffler 理論方法和Mirels 理論方法不再適用,而本文的理論方法依然可以應(yīng)用;本文理論計算得到的熱層內(nèi)激波強度,與數(shù)值模擬結(jié)果和已有的實驗數(shù)據(jù)較為吻合,與數(shù)值模擬結(jié)果最大偏差約8.60%,物質(zhì)界面處流場壓力、速度的理論計算結(jié)果與數(shù)值模擬結(jié)果最大偏差約20%。

    通過以上對比可以發(fā)現(xiàn),與現(xiàn)有的Shreffler 和Mirels 理論方法相比,本文的計算結(jié)果與數(shù)值模擬結(jié)果和已有的實驗數(shù)據(jù)更為吻合,且其適用范圍也更為廣泛。本文的研究結(jié)果可加深對激波與熱層作用機理的認識,為建立斜面或球面激波與熱層作用流場的理論計算方法提供參考。

    猜你喜歡
    激波計算結(jié)果流場
    大型空冷汽輪發(fā)電機轉(zhuǎn)子三維流場計算
    一種基于聚類分析的二維激波模式識別算法
    基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
    不等高軟橫跨橫向承力索計算及計算結(jié)果判斷研究
    甘肅科技(2020年20期)2020-04-13 00:30:40
    斜激波入射V形鈍前緣溢流口激波干擾研究
    轉(zhuǎn)杯紡排雜區(qū)流場與排雜性能
    適于可壓縮多尺度流動的緊致型激波捕捉格式
    基于HYCOM的斯里蘭卡南部海域溫、鹽、流場統(tǒng)計分析
    基于瞬態(tài)流場計算的滑動軸承靜平衡位置求解
    超壓測試方法對炸藥TNT當(dāng)量計算結(jié)果的影響
    嫩草影院新地址| 国产精品久久久久久人妻精品电影| 12—13女人毛片做爰片一| a级一级毛片免费在线观看| 成人亚洲精品av一区二区| 51午夜福利影视在线观看| 国产主播在线观看一区二区| 亚洲中文字幕一区二区三区有码在线看| 国产久久久一区二区三区| 高清毛片免费观看视频网站| 色哟哟·www| 欧美一区二区精品小视频在线| 亚洲国产日韩欧美精品在线观看| 九色国产91popny在线| 午夜精品久久久久久毛片777| 国产伦人伦偷精品视频| 波多野结衣高清作品| 亚洲精品成人久久久久久| 久久香蕉精品热| 亚洲国产日韩欧美精品在线观看| 尤物成人国产欧美一区二区三区| 久久久国产成人精品二区| 在现免费观看毛片| 亚洲精华国产精华精| 国产精品爽爽va在线观看网站| 一区二区三区激情视频| 伊人久久精品亚洲午夜| 国产av麻豆久久久久久久| 尤物成人国产欧美一区二区三区| 在线免费观看不下载黄p国产 | 国产白丝娇喘喷水9色精品| 精品一区二区三区av网在线观看| 精品人妻熟女av久视频| www.www免费av| 一区福利在线观看| 成年女人看的毛片在线观看| 亚洲国产日韩欧美精品在线观看| 看黄色毛片网站| netflix在线观看网站| 国产久久久一区二区三区| 免费观看的影片在线观看| 极品教师在线免费播放| 久久精品影院6| 男女下面进入的视频免费午夜| 午夜两性在线视频| 看十八女毛片水多多多| 欧美最新免费一区二区三区 | 亚州av有码| 欧美性猛交╳xxx乱大交人| 欧美xxxx黑人xx丫x性爽| 成人性生交大片免费视频hd| 蜜桃亚洲精品一区二区三区| av在线老鸭窝| 国产伦一二天堂av在线观看| 久久香蕉精品热| 欧美日本亚洲视频在线播放| 久久精品国产清高在天天线| 国产精品精品国产色婷婷| 一区二区三区免费毛片| 午夜激情欧美在线| 蜜桃久久精品国产亚洲av| 真实男女啪啪啪动态图| 在线观看舔阴道视频| 国产精品1区2区在线观看.| 午夜福利在线观看免费完整高清在 | 午夜免费男女啪啪视频观看 | 久久亚洲真实| 俺也久久电影网| 成年免费大片在线观看| 久久国产精品人妻蜜桃| 一个人看视频在线观看www免费| 亚洲 欧美 日韩 在线 免费| 村上凉子中文字幕在线| 97超视频在线观看视频| 我要搜黄色片| 久久九九热精品免费| 国产单亲对白刺激| 怎么达到女性高潮| 黄色配什么色好看| 中文字幕av成人在线电影| 国产精品av视频在线免费观看| 国产一区二区亚洲精品在线观看| 最近最新免费中文字幕在线| 老司机午夜十八禁免费视频| 国产又黄又爽又无遮挡在线| 免费看光身美女| 国产一区二区三区视频了| 18禁黄网站禁片免费观看直播| 亚洲熟妇中文字幕五十中出| 毛片一级片免费看久久久久 | 久久午夜福利片| 国产亚洲欧美98| 怎么达到女性高潮| 婷婷精品国产亚洲av| 欧美日本亚洲视频在线播放| 十八禁网站免费在线| 最近中文字幕高清免费大全6 | 欧美在线一区亚洲| 精品无人区乱码1区二区| 日韩欧美国产在线观看| 欧美绝顶高潮抽搐喷水| 1000部很黄的大片| 热99在线观看视频| 深爱激情五月婷婷| 精品欧美国产一区二区三| 午夜福利成人在线免费观看| 在线国产一区二区在线| 成年女人看的毛片在线观看| 最新中文字幕久久久久| 深爱激情五月婷婷| 精品人妻熟女av久视频| 观看免费一级毛片| 激情在线观看视频在线高清| 国产精品影院久久| 国产精品一区二区三区四区久久| 男女视频在线观看网站免费| xxxwww97欧美| 脱女人内裤的视频| 国产精品98久久久久久宅男小说| 免费看a级黄色片| www.999成人在线观看| 免费搜索国产男女视频| xxxwww97欧美| a级一级毛片免费在线观看| 日韩中文字幕欧美一区二区| 少妇被粗大猛烈的视频| 国产一区二区三区在线臀色熟女| 久久午夜亚洲精品久久| 国产av一区在线观看免费| 搞女人的毛片| 此物有八面人人有两片| 亚洲av免费在线观看| 日韩av在线大香蕉| 久久99热6这里只有精品| 国产伦精品一区二区三区四那| 全区人妻精品视频| 黄色女人牲交| 亚洲国产精品成人综合色| 日本撒尿小便嘘嘘汇集6| 国产精品电影一区二区三区| 亚洲人成网站在线播放欧美日韩| 五月伊人婷婷丁香| 黄色视频,在线免费观看| 国产淫片久久久久久久久 | 日韩欧美 国产精品| 免费电影在线观看免费观看| 欧美激情在线99| 国产视频内射| 精品久久久久久久人妻蜜臀av| 国产一级毛片七仙女欲春2| 草草在线视频免费看| 亚洲av熟女| 桃色一区二区三区在线观看| 少妇的逼水好多| 久久久久久九九精品二区国产| 亚洲av美国av| 中亚洲国语对白在线视频| 亚洲av不卡在线观看| 国产午夜精品论理片| 久久热精品热| 亚洲自偷自拍三级| 国产精华一区二区三区| 黄色丝袜av网址大全| 欧美日韩瑟瑟在线播放| 欧美日本亚洲视频在线播放| 淫秽高清视频在线观看| 18禁黄网站禁片午夜丰满| 看免费av毛片| av在线蜜桃| 亚洲熟妇中文字幕五十中出| 国产亚洲精品久久久com| www.熟女人妻精品国产| 欧美不卡视频在线免费观看| 深夜a级毛片| 国内毛片毛片毛片毛片毛片| 免费在线观看成人毛片| 免费一级毛片在线播放高清视频| 欧美日韩黄片免| 色吧在线观看| 午夜福利在线观看吧| 亚洲国产精品久久男人天堂| 757午夜福利合集在线观看| 亚洲美女黄片视频| 内射极品少妇av片p| 大型黄色视频在线免费观看| 国产精品三级大全| 熟女人妻精品中文字幕| 国产高清三级在线| 国产亚洲精品久久久久久毛片| 男女做爰动态图高潮gif福利片| 国产高清有码在线观看视频| 女人被狂操c到高潮| 久久久久久久久久成人| 日本五十路高清| 黄色视频,在线免费观看| 99久久精品国产亚洲精品| 在线国产一区二区在线| 国产精品久久久久久亚洲av鲁大| 国产一区二区亚洲精品在线观看| 精品人妻熟女av久视频| 精品久久久久久久人妻蜜臀av| 在线观看av片永久免费下载| 一二三四社区在线视频社区8| 亚洲五月婷婷丁香| 日韩欧美国产在线观看| 全区人妻精品视频| 深夜精品福利| av女优亚洲男人天堂| 怎么达到女性高潮| 国产高清视频在线播放一区| 国产国拍精品亚洲av在线观看| 99精品久久久久人妻精品| 丁香六月欧美| 国产三级中文精品| 国产精品女同一区二区软件 | xxxwww97欧美| 小蜜桃在线观看免费完整版高清| ponron亚洲| 天堂√8在线中文| 真人一进一出gif抽搐免费| 最近最新中文字幕大全电影3| 国产乱人伦免费视频| 久久这里只有精品中国| 欧美成人a在线观看| 在线播放国产精品三级| 国产精品国产高清国产av| 午夜影院日韩av| 国产精品精品国产色婷婷| 日本黄色片子视频| 国产精品99久久久久久久久| 亚洲国产精品999在线| 嫩草影院新地址| 精品日产1卡2卡| 亚洲国产高清在线一区二区三| 我要搜黄色片| 日日夜夜操网爽| 在线观看一区二区三区| 18禁黄网站禁片午夜丰满| 在线a可以看的网站| 噜噜噜噜噜久久久久久91| 午夜久久久久精精品| 亚洲 国产 在线| 免费搜索国产男女视频| 中文字幕av在线有码专区| www.色视频.com| 最新在线观看一区二区三区| 亚洲精品在线美女| 亚洲七黄色美女视频| 麻豆一二三区av精品| 他把我摸到了高潮在线观看| 亚洲精品456在线播放app | 波多野结衣高清作品| 亚洲av不卡在线观看| 好男人电影高清在线观看| 日韩av在线大香蕉| 好看av亚洲va欧美ⅴa在| 国产主播在线观看一区二区| 97热精品久久久久久| 69人妻影院| 国产一区二区亚洲精品在线观看| 最后的刺客免费高清国语| 欧美另类亚洲清纯唯美| 国产精品国产高清国产av| 久久精品夜夜夜夜夜久久蜜豆| 欧美性感艳星| 99久国产av精品| 欧美激情国产日韩精品一区| 此物有八面人人有两片| 久久久精品欧美日韩精品| 两人在一起打扑克的视频| 国产午夜精品论理片| 免费看a级黄色片| 51午夜福利影视在线观看| 91麻豆av在线| 国产亚洲精品久久久久久毛片| 男女床上黄色一级片免费看| 国产69精品久久久久777片| 亚洲一区二区三区色噜噜| 69av精品久久久久久| 18+在线观看网站| 日本精品一区二区三区蜜桃| 日韩欧美在线二视频| 精品人妻1区二区| 精品一区二区三区视频在线| 中文亚洲av片在线观看爽| 欧美在线一区亚洲| 国产探花极品一区二区| 欧美性猛交黑人性爽| 九色国产91popny在线| 久久精品久久久久久噜噜老黄 | 两人在一起打扑克的视频| 国产高潮美女av| a在线观看视频网站| 少妇熟女aⅴ在线视频| 久久久精品欧美日韩精品| 看免费av毛片| 深夜精品福利| 午夜久久久久精精品| 久久精品久久久久久噜噜老黄 | 老熟妇乱子伦视频在线观看| 免费观看人在逋| a级一级毛片免费在线观看| 国产成年人精品一区二区| 免费黄网站久久成人精品 | 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 久久久色成人| 亚洲国产高清在线一区二区三| 男女视频在线观看网站免费| 757午夜福利合集在线观看| 亚洲性夜色夜夜综合| 精品国内亚洲2022精品成人| 日本精品一区二区三区蜜桃| av在线老鸭窝| 国产私拍福利视频在线观看| 内地一区二区视频在线| 搡老岳熟女国产| 日本黄色片子视频| 国产av在哪里看| 国产真实乱freesex| 小说图片视频综合网站| 51午夜福利影视在线观看| 天堂影院成人在线观看| 精品久久久久久久久亚洲 | netflix在线观看网站| 美女高潮喷水抽搐中文字幕| 日本三级黄在线观看| 国产不卡一卡二| 老司机午夜福利在线观看视频| 欧美又色又爽又黄视频| 麻豆成人午夜福利视频| 一区二区三区激情视频| 一本久久中文字幕| 中文亚洲av片在线观看爽| 一区二区三区免费毛片| 变态另类成人亚洲欧美熟女| 夜夜爽天天搞| 亚洲经典国产精华液单 | 亚洲电影在线观看av| 久久久久国内视频| 小说图片视频综合网站| 亚洲精品久久国产高清桃花| 亚洲经典国产精华液单 | 国产高潮美女av| 成年人黄色毛片网站| 成人特级黄色片久久久久久久| 亚洲精品久久国产高清桃花| 九九在线视频观看精品| 午夜免费激情av| 欧美区成人在线视频| 久久午夜福利片| 亚洲精品一卡2卡三卡4卡5卡| 亚洲黑人精品在线| 日韩精品中文字幕看吧| 91在线精品国自产拍蜜月| 两性午夜刺激爽爽歪歪视频在线观看| 一级毛片久久久久久久久女| 国产探花极品一区二区| av视频在线观看入口| 亚洲欧美日韩东京热| 最近中文字幕高清免费大全6 | 精品一区二区三区视频在线观看免费| 嫁个100分男人电影在线观看| 亚洲精品一卡2卡三卡4卡5卡| 欧美黑人巨大hd| 国产乱人伦免费视频| 老鸭窝网址在线观看| 欧美不卡视频在线免费观看| 中亚洲国语对白在线视频| 黄色日韩在线| 国产精品久久视频播放| 网址你懂的国产日韩在线| 嫩草影院新地址| 人妻制服诱惑在线中文字幕| 免费搜索国产男女视频| 免费av不卡在线播放| 久久亚洲真实| 在线观看66精品国产| 欧美xxxx黑人xx丫x性爽| 一级黄色大片毛片| 最后的刺客免费高清国语| 熟女电影av网| 国产欧美日韩一区二区精品| 成人毛片a级毛片在线播放| 精品人妻熟女av久视频| 神马国产精品三级电影在线观看| 久99久视频精品免费| 亚洲成人久久性| 99riav亚洲国产免费| 亚洲电影在线观看av| 久久午夜亚洲精品久久| 在线观看av片永久免费下载| 两人在一起打扑克的视频| 亚洲自拍偷在线| 91麻豆av在线| 国产免费av片在线观看野外av| 中文字幕人妻熟人妻熟丝袜美| 久久中文看片网| 亚洲欧美日韩高清在线视频| 久久久久久久午夜电影| 亚洲午夜理论影院| 精品免费久久久久久久清纯| 亚洲经典国产精华液单 | 神马国产精品三级电影在线观看| 亚洲一区高清亚洲精品| 18禁黄网站禁片免费观看直播| 精品99又大又爽又粗少妇毛片 | 桃红色精品国产亚洲av| 日本黄大片高清| 久99久视频精品免费| 国产精品久久久久久亚洲av鲁大| 国产成人欧美在线观看| 久久这里只有精品中国| 久久久久久久亚洲中文字幕 | 日韩欧美在线乱码| 一级毛片久久久久久久久女| 亚洲av五月六月丁香网| 日韩欧美精品免费久久 | 男插女下体视频免费在线播放| 中文亚洲av片在线观看爽| 国产欧美日韩一区二区精品| av天堂在线播放| 亚洲精品456在线播放app | 亚洲av第一区精品v没综合| 真人做人爱边吃奶动态| 色吧在线观看| 中出人妻视频一区二区| 国产免费一级a男人的天堂| 制服丝袜大香蕉在线| 国产激情偷乱视频一区二区| 国产免费av片在线观看野外av| 亚洲人与动物交配视频| 激情在线观看视频在线高清| 女生性感内裤真人,穿戴方法视频| 99久国产av精品| 国产一区二区三区视频了| 日本a在线网址| 色视频www国产| 村上凉子中文字幕在线| 日韩欧美在线二视频| 在线免费观看的www视频| 欧美最黄视频在线播放免费| 舔av片在线| 亚洲最大成人中文| 97热精品久久久久久| 国产精品一区二区性色av| 热99re8久久精品国产| 色在线成人网| 性欧美人与动物交配| 亚洲片人在线观看| 国产在线精品亚洲第一网站| 噜噜噜噜噜久久久久久91| 国产探花极品一区二区| 国产成人影院久久av| 麻豆久久精品国产亚洲av| 欧美xxxx性猛交bbbb| 成人三级黄色视频| 一本精品99久久精品77| a级毛片免费高清观看在线播放| 内射极品少妇av片p| 午夜福利在线观看吧| 精品一区二区三区人妻视频| 很黄的视频免费| 欧美bdsm另类| 国产精品一区二区免费欧美| 我要看日韩黄色一级片| 国产精品98久久久久久宅男小说| 国产淫片久久久久久久久 | 国产精华一区二区三区| 国产精品影院久久| 蜜桃久久精品国产亚洲av| 亚洲18禁久久av| 91在线观看av| 一区二区三区免费毛片| 日本免费一区二区三区高清不卡| 91麻豆av在线| 亚洲18禁久久av| 在线观看舔阴道视频| 免费搜索国产男女视频| 99久久久亚洲精品蜜臀av| 亚洲成a人片在线一区二区| 亚洲片人在线观看| 免费人成在线观看视频色| 观看美女的网站| 在线观看免费视频日本深夜| 午夜激情福利司机影院| 久久这里只有精品中国| 国产伦精品一区二区三区四那| 一夜夜www| 三级男女做爰猛烈吃奶摸视频| 麻豆成人午夜福利视频| 精华霜和精华液先用哪个| 久久国产精品人妻蜜桃| av在线天堂中文字幕| 免费av不卡在线播放| 成人av在线播放网站| 亚洲精品日韩av片在线观看| 国产色爽女视频免费观看| 国产一区二区三区在线臀色熟女| av福利片在线观看| 51国产日韩欧美| avwww免费| 夜夜躁狠狠躁天天躁| 一级毛片久久久久久久久女| 成人特级黄色片久久久久久久| 一级黄片播放器| 99riav亚洲国产免费| 成人毛片a级毛片在线播放| 蜜桃亚洲精品一区二区三区| 日韩欧美三级三区| 9191精品国产免费久久| 日本 欧美在线| 国产精品综合久久久久久久免费| 成熟少妇高潮喷水视频| www.999成人在线观看| 欧美绝顶高潮抽搐喷水| 男人和女人高潮做爰伦理| 日韩成人在线观看一区二区三区| 99热这里只有是精品50| 国产黄色小视频在线观看| 欧美日韩综合久久久久久 | 中出人妻视频一区二区| www.999成人在线观看| 特级一级黄色大片| 毛片女人毛片| 91av网一区二区| 午夜久久久久精精品| 最近视频中文字幕2019在线8| 日韩精品中文字幕看吧| 精品人妻熟女av久视频| 久久这里只有精品中国| 午夜免费男女啪啪视频观看 | 51午夜福利影视在线观看| 动漫黄色视频在线观看| 一级av片app| 亚洲国产欧洲综合997久久,| 757午夜福利合集在线观看| 精品一区二区三区视频在线| 国产一区二区在线观看日韩| 免费在线观看日本一区| 免费人成视频x8x8入口观看| 成人毛片a级毛片在线播放| 校园春色视频在线观看| 久久精品综合一区二区三区| 看十八女毛片水多多多| 美女 人体艺术 gogo| 他把我摸到了高潮在线观看| 国产熟女xx| 美女cb高潮喷水在线观看| 亚洲自拍偷在线| 午夜精品一区二区三区免费看| 日日干狠狠操夜夜爽| 成人毛片a级毛片在线播放| 精品久久久久久成人av| 久久久久久大精品| 三级男女做爰猛烈吃奶摸视频| aaaaa片日本免费| 在线免费观看的www视频| 成人永久免费在线观看视频| 精品欧美国产一区二区三| 午夜福利欧美成人| 赤兔流量卡办理| www.www免费av| 村上凉子中文字幕在线| 高潮久久久久久久久久久不卡| 女人十人毛片免费观看3o分钟| 悠悠久久av| 亚洲精品乱码久久久v下载方式| a级毛片a级免费在线| 国产亚洲精品久久久久久毛片| 首页视频小说图片口味搜索| 亚洲自偷自拍三级| 少妇高潮的动态图| 亚洲人成伊人成综合网2020| 99热只有精品国产| 久久久久国产精品人妻aⅴ院| 亚洲人成电影免费在线| 丰满的人妻完整版| 国产成+人综合+亚洲专区| 亚洲最大成人手机在线| 一级a爱片免费观看的视频| 美女被艹到高潮喷水动态| 中国美女看黄片| 亚洲自拍偷在线| 啦啦啦韩国在线观看视频| 最近最新免费中文字幕在线| 村上凉子中文字幕在线| 国产精品电影一区二区三区| 18美女黄网站色大片免费观看| 亚洲自拍偷在线| 成人欧美大片| 精品久久久久久久末码| 人人妻人人澡欧美一区二区| 国内精品久久久久精免费| 夜夜爽天天搞| 日本与韩国留学比较| 怎么达到女性高潮| 中亚洲国语对白在线视频| 村上凉子中文字幕在线| 99精品在免费线老司机午夜| 国模一区二区三区四区视频| 少妇人妻一区二区三区视频| 亚洲av中文字字幕乱码综合| 久久久久国内视频| 亚洲精品成人久久久久久| 人人妻人人看人人澡| 免费观看的影片在线观看| 午夜福利欧美成人| 亚洲午夜理论影院| 真人一进一出gif抽搐免费| 黄色女人牲交| 免费看光身美女| 一级黄片播放器| 欧美精品国产亚洲| 看十八女毛片水多多多| 美女高潮的动态|