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

    加權(quán)緊致非線性格式在熱化學(xué)非平衡流數(shù)值模擬中的應(yīng)用*

    2016-11-28 01:18:22葛明明趙小宇
    關(guān)鍵詞:熱化學(xué)駐點(diǎn)激波

    葛明明,曾 明,趙小宇

    (國防科技大學(xué) 航天科學(xué)與工程學(xué)院, 湖南 長沙 410073)

    ?

    加權(quán)緊致非線性格式在熱化學(xué)非平衡流數(shù)值模擬中的應(yīng)用*

    葛明明,曾 明,趙小宇

    (國防科技大學(xué) 航天科學(xué)與工程學(xué)院, 湖南 長沙 410073)

    高精度格式因其在較低的計(jì)算資源消耗下仍能保持足夠的精度而在精細(xì)流動(dòng)結(jié)構(gòu)的模擬中具有獨(dú)特優(yōu)勢(shì)。將顯式5階加權(quán)緊致非線性格式(WCNS-E-5)引入二維高溫非平衡流計(jì)算,進(jìn)行自由流速度5.0~5.7 km/s的高焓風(fēng)洞和高空條件下的圓柱非平衡流場(chǎng)數(shù)值模擬,得到正確的流場(chǎng)參數(shù)分布、壁面壓力和壁面熱流等結(jié)果與實(shí)驗(yàn)值吻合較好。計(jì)算結(jié)果表明,WCNS-E-5格式較2階的MUSCL格式具有更好的熱流網(wǎng)格收斂特性。分析并針對(duì)算例條件初步解決了高階格式在剛性化學(xué)反應(yīng)流模擬中的魯棒性問題,為進(jìn)一步采用高階格式開展流場(chǎng)結(jié)構(gòu)更加復(fù)雜的熱化學(xué)非平衡流數(shù)值研究奠定了良好基礎(chǔ)。

    高精度格式;熱化學(xué)非平衡流;數(shù)值模擬

    高溫條件下的熱化學(xué)非平衡現(xiàn)象是高超聲速流場(chǎng)的重要特征之一,這方面的數(shù)值研究始于20世紀(jì)六七十年代。經(jīng)過半個(gè)多世紀(jì)的發(fā)展,各方面都得到了較為成熟并被普遍接受的理論方法。解算器方面,國內(nèi)外也發(fā)展出了較為完善、功能全面的解算器,如NASA Langley中心發(fā)展的蘭利氣動(dòng)熱力學(xué)迎風(fēng)松弛算法(Langley Aerothermodynamic Upwind Relaxation Algorithm, LAURA)程序[1],Gaitonde等發(fā)展的通用氣動(dòng)熱力學(xué)模擬程序(General Aerothermodynamic Simulation Program, GASP)[2-3]等。不過目前高超聲速熱化學(xué)非平衡流模擬所采用的數(shù)值算法以2階格式為主,這對(duì)外形較為簡單的流場(chǎng)模擬和流動(dòng)機(jī)理研究,尚可以達(dá)到精度需求。而21世紀(jì)發(fā)展的新型高超聲速飛行器,外形相對(duì)復(fù)雜,飛行器流場(chǎng)中激波/邊界層、激波/激波相互作用等流動(dòng)現(xiàn)象普遍存在。此時(shí)不僅存在宏觀流動(dòng)和微觀物理化學(xué)過程的緊密耦合,還存在復(fù)雜的流動(dòng)相互作用,多種效應(yīng)并存且彼此制約,互相影響。激波位置、形狀和強(qiáng)度,邊界層厚度和形狀都直接影響到相互作用的性質(zhì)和強(qiáng)度,影響到流動(dòng)分離的特征和分離區(qū)大小,這些又都進(jìn)一步影響到流場(chǎng)各點(diǎn)的熱化學(xué)狀態(tài)。為深入理解存在復(fù)雜波系和相互作用的非平衡流動(dòng)的內(nèi)在機(jī)理,準(zhǔn)確預(yù)測(cè)這類情況下飛行器氣動(dòng)力、氣動(dòng)熱和氣動(dòng)物理特性,要求非平衡流的數(shù)值模擬具有更高的精度和效率。

    高階格式因其兼具高空間分辨率和計(jì)算效率的特點(diǎn)而在精細(xì)流動(dòng)結(jié)構(gòu)的模擬中具有獨(dú)特優(yōu)勢(shì)[4],將其應(yīng)用于包含復(fù)雜流場(chǎng)結(jié)構(gòu)的熱化學(xué)非平衡流的數(shù)值模擬很有意義。近期國外應(yīng)用高階格式模擬化學(xué)反應(yīng)流的工作主要是采用Shu[5]的加權(quán)本質(zhì)無振蕩(Weighted Essentially Non Oscillatory, WENO)格式。Wang[6]發(fā)展了一個(gè)基于5階WENO格式的高速化學(xué)非平衡流動(dòng)求解器,通過對(duì)圓柱非平衡流的模擬及與2階總變差減小(Total Variation Diminishing, TVD)格式結(jié)果的對(duì)比驗(yàn)證了WENO求解高超聲速化學(xué)非平衡流動(dòng)的可行性。Lani[7]將5階、7階、9階的WENO格式應(yīng)用于非定常三維可壓縮流求解器ADPDIS 3D,進(jìn)行了高超聲速熱化學(xué)非平衡湍流數(shù)值模擬。Prakash[8]采用一種空間5階、時(shí)間3階精度的激波裝配法進(jìn)行化學(xué)非平衡流計(jì)算,流場(chǎng)捕捉到了不同頻率的擾動(dòng)波疊加效果,結(jié)論指出足夠高的精度對(duì)捕捉擾動(dòng)的必要性。

    加權(quán)緊致非線性格式[9-11](Weighted Compact Nonlinear Scheme, WCNS)是鄧小剛等在20世紀(jì)90年代提出的具有激波捕捉能力的高精度有限差分格式,已成功應(yīng)用于低速、超聲速和高超聲速(量熱完全氣體)流場(chǎng)的數(shù)值模擬,在大量復(fù)雜流動(dòng)問題中展示出優(yōu)良性能[12-13]。國外學(xué)者Nonomura[14]指出,采用對(duì)稱守恒網(wǎng)格導(dǎo)數(shù)算法(Symmetrical Conservative Metric Method, SCMM)的WCNS能滿足幾何守恒,因而在復(fù)雜網(wǎng)格問題中更具優(yōu)勢(shì);Matsukawa[15]評(píng)價(jià):WCNS既具有高的精度,又具有激波捕捉能力,形式簡便,計(jì)算耗費(fèi)也較少。

    將顯式5階加權(quán)緊致非線性格式(WCNS-E-5)引入高超聲速熱化學(xué)非平衡流的數(shù)值求解,探討高階格式在非平衡流中應(yīng)用的數(shù)值困難并采取相應(yīng)措施,編制采用WCNS-E-5的二維熱化學(xué)非平衡流計(jì)算程序。通過圓柱非平衡流場(chǎng)的模擬初步考察程序的有效性,并研究高階和2階格式在網(wǎng)格收斂性上表現(xiàn)的差異。為進(jìn)一步采用高階格式開展具有復(fù)雜流場(chǎng)結(jié)構(gòu)的熱化學(xué)非平衡流動(dòng)數(shù)值研究奠定基礎(chǔ)。

    1 熱化學(xué)模型與控制方程

    對(duì)高溫空氣考慮7組元(O2, N2, N, O, NO, NO+, e-)。化學(xué)反應(yīng)模型采用Gupta的7組元6反應(yīng)模型[16],采用Park的雙溫度模型[17]實(shí)現(xiàn)化學(xué)與熱力學(xué)非平衡的耦合?;瘜W(xué)反應(yīng)方程、振動(dòng)松弛模型及化學(xué)反應(yīng)源項(xiàng)和振動(dòng)松弛源項(xiàng)的計(jì)算方法詳見文獻(xiàn)[18-19]。

    采用時(shí)間相關(guān)的二維熱化學(xué)非平衡納維斯托克斯(Navier-Stokes, NS)方程,求解流場(chǎng)的定常解。計(jì)算坐標(biāo)系下的無量綱控制方程為:

    (1)

    (2)

    式中,ρi和ρ分別代表組元i和混合氣體的密度,E為單位體積總能,ev為單位質(zhì)量振動(dòng)能,J為坐標(biāo)變換雅可比行列式。

    (3)

    (4)

    (5)

    (6)

    bx=uτxx+vτxy+qx

    (7)

    (8)

    (9)

    式中,ns為組元個(gè)數(shù),nm為分子組元個(gè)數(shù)。by,qy,qyv的表達(dá)式類似。

    2 數(shù)值方法

    2.1 時(shí)間格式

    控制方程的時(shí)間離散格式采用Shu[20]提出的三階精度具有TVD 性質(zhì)的Runge-Kutta方法。

    (10)

    式中,U代表原始量,RHS(U)為控制方程右端項(xiàng),則:

    (11)

    2.2 對(duì)流通量

    根據(jù)WCNS-E-5對(duì)無粘通量項(xiàng)進(jìn)行空間離散,通量導(dǎo)數(shù)的計(jì)算采用六階精度的中心差分格式。

    (12)

    (13)

    根據(jù)左右原始量得到通量:

    (14)

    本文嘗試了多種通量格式。AUSMPW+通量格式因其較高的計(jì)算效率與在間斷和邊界層內(nèi)良好的分辨率被廣泛運(yùn)用,對(duì)于非平衡流場(chǎng)熱流的模擬表現(xiàn)較好[21]。但在計(jì)算實(shí)踐中發(fā)現(xiàn),在含有強(qiáng)激波的算例中,WCNS-E-5與AUSMPW+通量格式的組合可能導(dǎo)致激波附近流場(chǎng)出現(xiàn)振蕩,甚至?xí)褂?jì)算發(fā)散。這可能是因?yàn)楦叻直媛矢袷讲蹲郊げ〞r(shí)耗散不夠,引起了數(shù)值計(jì)算不穩(wěn)定。而Shu研究[22]指出,隨著格式精度的提高,通量格式差別對(duì)于計(jì)算結(jié)果的影響逐漸減小。故本文最后選用了耗散較大的Steger-Warming格式。

    左右原始量利用5階加權(quán)非線性插值得到。這里以左值為代表給出具體插值方法,右值類似。

    (15)

    其中

    (16)

    (17)

    其中

    (18)

    式中,CLk為最優(yōu)權(quán)重系數(shù),ISk代表模板k的光滑度量函數(shù),ε為小量,一般取10-6。

    一般情況下,這種非線性加權(quán)方法能夠很好地捕捉激波,并在光滑區(qū)域恢復(fù)精度。但在熱化學(xué)非平衡流動(dòng)中,激波附近微量組元的質(zhì)量分?jǐn)?shù)量級(jí)可能在10-10以下,如果ε仍然取10-6,則光滑度量函數(shù)失效,可能導(dǎo)致微量組元密度在激波附近劇烈振蕩。因此需要減小ε的取值。但有學(xué)者[23]針對(duì)WENO研究指出,ε取值過小也會(huì)降低光滑區(qū)域求解精度。因此,本文采用

    ε=max[10-26,10-6×min(IS1,IS2,IS3)]

    (19)

    這樣針對(duì)間斷處不同物理量的量級(jí)設(shè)定ε,避免在間斷區(qū)間影響光滑度量函數(shù)的效果,又不至于在光滑區(qū)域過小而影響求解精度。

    2.3 黏性通量

    黏性項(xiàng)的差分格式與對(duì)流項(xiàng)的處理一致,采用6階中心差分。半節(jié)點(diǎn)的通量可以用6階中心插值得到的原始量計(jì)算。

    (20)

    式中

    (21)

    在非平衡流動(dòng)中,跨過激波后的化學(xué)反應(yīng)使部分組元質(zhì)量分?jǐn)?shù)發(fā)生量級(jí)變化,此時(shí)高階線性插值會(huì)在激波附近產(chǎn)生非物理振蕩,甚至插值出負(fù)值,導(dǎo)致計(jì)算失敗。文獻(xiàn)[13]提出了一種應(yīng)用于完全氣體條件下的非線性插值格式,它利用3個(gè)線性模板,通過構(gòu)造合理的權(quán)重,來消除間斷導(dǎo)致的振蕩,并且在光滑區(qū)域理論上可達(dá)到6階精度。本文則采用了類似2.2節(jié)中給出的對(duì)流項(xiàng)插值的方法,首先得到半節(jié)點(diǎn)的左右值,對(duì)二者取平均得到半節(jié)點(diǎn)的原始量,進(jìn)而計(jì)算出黏性通量。算例計(jì)算實(shí)踐中發(fā)現(xiàn),對(duì)于非平衡流動(dòng),本文做法的魯棒性更強(qiáng)。

    3 算例與結(jié)果分析

    根據(jù)上述數(shù)值方法編制了應(yīng)用WCNS-E-5的二維熱化學(xué)非平衡流計(jì)算程序,開展了高焓條件下的圓柱非平衡流場(chǎng)數(shù)值模擬,通過與實(shí)驗(yàn)和文獻(xiàn)計(jì)算結(jié)果的比較,初步考察程序的有效性。

    這里給出三個(gè)代表性算例。算例1為德國航空宇航中心(Deutsches zentrum für Luftund Raumfahrt, DLR)進(jìn)行的高焓風(fēng)洞(High Enthalpy shock tunnel Gōttingen, HEG)圓柱繞流實(shí)驗(yàn)[24],重點(diǎn)對(duì)比分析實(shí)驗(yàn)測(cè)得的和本文計(jì)算得到的壁面壓力、熱流分布結(jié)果。算例2為中科院力學(xué)所高焓風(fēng)洞自由流條件[18]下的圓柱繞流,主要比較2階(采用MUSCL格式,記作MUSCL2)和高階格式的計(jì)算結(jié)果及二者關(guān)于駐點(diǎn)熱流的網(wǎng)格收斂性。算例3采用與算例2自由流總焓相同的高空自由流條件,與算例1、算例2(均為高焓風(fēng)洞非平衡自由流條件)不同的是,算例3的自由流中僅含N2和O2組元,跨過激波后的化學(xué)反應(yīng)引起原子和離子組元的質(zhì)量分?jǐn)?shù)發(fā)生很大的量級(jí)變化。這對(duì)應(yīng)用高階格式計(jì)算程序的魯棒性提出了更高要求。

    3.1 算例條件與計(jì)算網(wǎng)格

    三個(gè)算例的來流參數(shù)見表1。

    對(duì)于算例1,文獻(xiàn)[25]采用了5組元(不含NO+和e-)的高溫空氣模型,本文將自由流的NO+質(zhì)量分?jǐn)?shù)設(shè)定為10-20。算例3為高空自由流條件,只含有N2和O2組元,本文將其他組元的質(zhì)量分?jǐn)?shù)人為設(shè)定為10-10。

    算例1的圓柱半徑為45 mm,算例2和算例3的均為35 mm。計(jì)算域的外邊界在駐點(diǎn)區(qū)離開壁面0.5個(gè)半徑長度。計(jì)算網(wǎng)格量為121×121或121×61,在壁面附近采用指數(shù)拉伸的方法加密,采用了多個(gè)不同的網(wǎng)格拉伸因子生成多套網(wǎng)格(對(duì)應(yīng)壁面處網(wǎng)格雷諾數(shù)在1~1000范圍內(nèi)),以驗(yàn)證計(jì)算結(jié)果的網(wǎng)格收斂性。壁面網(wǎng)格雷諾數(shù)ReΔη采用式(22)計(jì)算:

    (22)

    式中:Δη1為駐點(diǎn)無量綱法向第一層網(wǎng)格高度;Rn為圓柱半徑;ρw,Vs,w,μw分別為壁面處混合氣體密度、聲速以及黏性系數(shù)。下面對(duì)三個(gè)算例分別給出ReΔη=8,10,8的121×121網(wǎng)格對(duì)應(yīng)的計(jì)算結(jié)果。對(duì)于算例2,還以駐點(diǎn)熱流為代表,詳細(xì)考察了計(jì)算結(jié)果的網(wǎng)格收斂性。

    表1 來流參數(shù)

    計(jì)算中除5階的WCNS-E-5外,還采用了2階的MUSCL格式,對(duì)比分析二者結(jié)果。通量格式為Steger-Warming格式,熵修正系數(shù)取0.05,MUSCL2格式中采用minmod限制器。三個(gè)算例均設(shè)為等壁溫(300 K)、全催化壁條件。

    3.2 算例1

    將本文計(jì)算結(jié)果分別與文獻(xiàn)[24]實(shí)驗(yàn)測(cè)量結(jié)果及文獻(xiàn)[25]計(jì)算進(jìn)行對(duì)比,文獻(xiàn)[25]是應(yīng)用GASP程序[2],采用2階的MUSCL格式,通量格式為AUSMPW+格式,使用minmod限制器,網(wǎng)格量為161×151。從表面壓力分布看,本文結(jié)果、文獻(xiàn)[25]的數(shù)值結(jié)果均與實(shí)驗(yàn)[24]測(cè)量值吻合很好。從表面熱流分布(如圖1所示)看,本文的WCNS-E-5與2階MUSCL格式結(jié)果基本一致,較文獻(xiàn)[25]的數(shù)值結(jié)果更加接近實(shí)驗(yàn)[24]值。另外,本文的5階和2階格式熱流結(jié)果在駐點(diǎn)區(qū)存在一定差異,2階格式得到的熱流分布曲線的光滑性不如5階格式。

    圖1 算例1壁面熱流分布Fig.1 Surface heat flux distribution for case 1

    (a) p (b) T (c) Tv

    (d) CO (e) CNO+圖2 算例2流場(chǎng)參數(shù)分布Fig. 2 Flow field parameters distribution for case 2

    3.3 算例2

    圖2分別給出了本文采用2階和5階格式得到的流場(chǎng)壓力、溫度和O,NO+組元質(zhì)量分?jǐn)?shù)分布,圖中上半部分為2階格式結(jié)果,下半部分為5階格式結(jié)果。在此算例條件下,兩種格式的結(jié)果基本一致。注意到5階格式得到的流場(chǎng)壓力等值線在圖中標(biāo)出的位置存在抖動(dòng),這一現(xiàn)象在文獻(xiàn)[26]的量熱完全氣體流場(chǎng)結(jié)果中也存在,是由于激波與網(wǎng)格走向不一致造成的,文獻(xiàn)[26]通過優(yōu)化網(wǎng)格來消除。

    2階和5階格式得到的駐點(diǎn)線上組元質(zhì)量分?jǐn)?shù)分布總體吻合,只有CNO+在激波附近和過激波后的峰值存在微小差別,如圖3所示。

    圖3 算例2駐點(diǎn)線上組元NO+質(zhì)量分?jǐn)?shù)分布Fig.3 NO+mass fraction distribution along the stagnation line for case 2

    下面對(duì)比分析5階WCNS-E-5和2階MUSCL格式在熱流計(jì)算的網(wǎng)格收斂性上的表現(xiàn)。圖4給出了二者得到的駐點(diǎn)熱流隨ReΔη的變化曲線。在5%的誤差范圍內(nèi),5階格式在ReΔη達(dá)到300左右基本收斂,而2階格式則要求ReΔη低至60。

    圖4 算例2駐點(diǎn)熱流網(wǎng)格收斂性曲線Fig.4 Grid convergence curve of the heat flux at stagnation point for case 2

    可見同樣精度要求下高階格式對(duì)網(wǎng)格密度要求明顯降低,因此可提高計(jì)算效率。圖5給出了5階格式采用121×61(網(wǎng)格雷諾數(shù)ReΔη=90),121×61(ReΔη=60)和2階格式采用121×121(ReΔη=30)網(wǎng)格的殘差收斂曲線。5階格式可比2階格式減少網(wǎng)格量,從而減少每個(gè)時(shí)間步的計(jì)算時(shí)間;壁面附近網(wǎng)格也可更疏,這樣能顯著加大每個(gè)時(shí)間步推進(jìn)的Δt,因此得到定常結(jié)果需要的時(shí)間推進(jìn)步數(shù)大大減少。雖然網(wǎng)格量相同時(shí)5階格式的單步計(jì)算時(shí)間比2階格式長(該算例條件下高22%),但殘差達(dá)到10-12量級(jí)時(shí)5階格式采用ReΔη=90網(wǎng)格的總計(jì)算時(shí)間僅為2階格式采用ReΔη=30網(wǎng)格的22.5%。

    圖5 算例2密度殘差收斂曲線Fig.5 Convergence histories of density residuals for case 2

    3.4 算例3

    (a) CO,CN,CNO

    (b) CNO+圖6 算例3駐點(diǎn)線組元質(zhì)量分?jǐn)?shù)分布Fig.6 Species mass fraction distribution along the stagnation line for case 3

    與算例2類似,5階和2階格式得到的壓力和溫度分布沒有差異。圖6給出了兩種格式得到的駐點(diǎn)線組元質(zhì)量分?jǐn)?shù)分布。與算例2不同的是,兩種格式得到的微量組元分布在激波附近有細(xì)小差別,峰值也略有不同。這是因?yàn)楦呖諄砹鳁l件下微量組元跨過激波變化更大,而高階格式在同樣網(wǎng)格下對(duì)梯度的分辨精度更高,能夠更好地捕捉梯度。以圖6(b)給出的CNO+分布為例,5階格式曲線上升的起始位置略滯后于2階格式,而峰值則比2階格式稍高。其實(shí)兩種格式結(jié)果差異的這個(gè)特點(diǎn)在其他參數(shù)分布上也有體現(xiàn),但由于相對(duì)差值很小,表現(xiàn)不明顯。

    5階和2階格式在熱流計(jì)算的網(wǎng)格收斂性上表現(xiàn)與算例2類似。采用ReΔη等于320,160,80,40,20,8的121×121網(wǎng)格的計(jì)算結(jié)果表明,5階格式在ReΔη低于160后駐點(diǎn)熱流結(jié)果的變化就在2%以內(nèi),而2階格式需要ReΔη低至40。

    4 結(jié)論

    將顯式五階加權(quán)緊致非線性格式WCNS-E-5引入非平衡流的數(shù)值求解,針對(duì)高階格式在化學(xué)反應(yīng)流應(yīng)用中的魯棒性問題采取措施,編制了采用WCNS-E-5格式的二維熱化學(xué)非平衡流計(jì)算程序,開展了高焓條件下的圓柱流場(chǎng)模擬。得到以下結(jié)論:

    1)關(guān)于非線性加權(quán)插值方法的改進(jìn)(即針對(duì)間斷處不同物理量的量級(jí)設(shè)定權(quán)系數(shù)計(jì)算中的小量),確保了微量化學(xué)組元在間斷區(qū)的光滑度量函數(shù)效果;黏性項(xiàng)插值計(jì)算中借鑒對(duì)流項(xiàng)非線性插值方法,解決了激波附近組元質(zhì)量分?jǐn)?shù)量級(jí)變化引起的插值振蕩問題。

    2)在高焓風(fēng)洞條件下的圓柱繞流計(jì)算中,5階WCNS-E-5格式結(jié)果優(yōu)于2階的MUSCL格式。HEG風(fēng)洞算例中本文結(jié)果比文獻(xiàn)計(jì)算結(jié)果更接近實(shí)驗(yàn)值,5階格式得到的熱流分布曲線比2階格式光滑。根據(jù)駐點(diǎn)熱流考察的網(wǎng)格收斂性,5階格式也明顯優(yōu)于2階格式。

    3)高空自由流條件下,跨過激波后的化學(xué)反應(yīng)使來流中微量組元質(zhì)量分?jǐn)?shù)發(fā)生量級(jí)變化。算例結(jié)果表明,5階格式對(duì)組元梯度的捕捉比2階格式精細(xì),在熱流計(jì)算的網(wǎng)格收斂性上也有明顯優(yōu)勢(shì)。

    4)在同樣精度要求下,采用5階格式可比2階格式減少網(wǎng)格量,并且顯著放寬對(duì)壁面附近網(wǎng)格密度的要求,從而可加大每個(gè)時(shí)間步推進(jìn)的Δt。雖然同樣網(wǎng)格量時(shí)5階格式單步計(jì)算時(shí)間更長,但得到定常結(jié)果需要的總計(jì)算時(shí)間減少。在本文算例2中可降至2階格式的1/4。

    5)未來還需在多方面努力,包括較寬范圍速度高度條件下流場(chǎng)激波附近網(wǎng)格的調(diào)試,組元連續(xù)方程和振動(dòng)能方程的特殊處理,進(jìn)一步將高階格式應(yīng)用于存在復(fù)雜干擾現(xiàn)象的熱化學(xué)非平衡流場(chǎng)數(shù)值模擬,發(fā)揮高階格式在捕捉精細(xì)流場(chǎng)結(jié)構(gòu)方面的優(yōu)勢(shì)。

    References)

    [1] Cheatwood F M, Gnoffo P A. User′s manual for the langley aerothermodynamic upwind relaxation algorithm (LAURA)[R]. NASA TM-4674, 1996.

    [2] Knight D, Longo L, Drikakis D, et al. Assessment of CFD capability for prediction of hypersonic shock in interactions[J]. Progress in Aerospace Sciences, 2012, S48/49(2): 8-26.

    [3] Olynick D, Tam T. Trajectory based validation of the shuttle heating environment [C]//Proceedings of 31st Thermophysics Conference, AIAA 96-1891, 1996.

    [4] Wang Z J, Fidkowski K, Abgrall R, et al. High-order CFD methods:current status and perspective [J]. International Journal for Numerical Methods in Fluids, 2013, 72(8): 811-845.

    [5] Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes [J]. Journal of Computational Physics, 1995, 126(1): 202-228.

    [6] Wang X W, Zhong X L. Nonequilibrium and reactive high-speed flow simulations with a fifth-order WENO scheme[C]//Proceedings of 39th AIAA Fluid Dynamics Conference, AIAA 2009-4041, 2009.

    [7] Lani A, Sjogreen B, Yee H C, et al. High-order simulation of hypersonic nonequilibrium flows on overset grids [M]//Center for Turbulence Research Annual Research Briefs, 2010.

    [8] Prakash A, Parsons N, Wang X, et al. High-order shock-fitting methods for direct numerical simulation of hypersonic flow with chemical and thermal nonequilibrium [J]. Journal of Computational Physics, 2011, 230(23): 8474-8507.

    [9] Deng X, Maekawa H. Compact high-order accurate nonlinear schemes [J]. Journal of Computational Physics, 1997, 130(1): 77-91.

    [10] Deng X G, Mao M L. Weighted compact high-order nonlinear schemes for the Euler equations[C]//Proceedings of 13th Computational Fluid Dynamics Conference, AIAA 97-1941, 1997.

    [11] Deng X G, Zhang H X. Developing high-order weighted compact nonlinear schemes [J]. Journal of Computational Physics, 2000, 165(1): 22-44.

    [12] Deng X G, Wang G X, Tu G H, et al. Applications of high-order weighted compact nonlinear scheme for complex transonic flows[C]//Proceedings of 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace, AIAA 2011-364, 2011.

    [13] Liu H Y, Deng X G, Mao M L, et al. High order nonlinear schemes for viscous terms and the application to complex conguration problems[C]//Proceedings of 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace, AIAA 2011-369, 2011.

    [14] Nonomura T, Iizuka N, Fujii K.Free-stream and vortex preservation properties of high-order WENO and WCNS on curvilinear grids [J]. Computers and Fluids, 2010, 39(2): 197-214.

    [15] Matsukawa Y. Implicit large eddy simulation of a supersonic flat-plate boundary layer flow by weighted compact nonlinear scheme [J]. International Journal of Computational Fluid Dynamics, 2011, 25(2): 47-57.

    [16] Gupta R N, Yos J M, Thompson R A. Review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K[R]. NASA RP-1232, 1990.

    [17] Park C. Problems of rate chemistry in the flight regimes of aeroassisted orbital transfer vehicles [C]//Proceedings of 19th Thermophysics Conference, AIAA 1984-1730, 1984.

    [18] 曾明. 高焓風(fēng)洞流場(chǎng)測(cè)量的數(shù)值重建和非平衡效應(yīng)的數(shù)值分析[D].北京:中國科學(xué)院力學(xué)研究, 2007.

    ZENG Ming. Numerical rebuilding of free-stream measurement and analysis of nonequilibrium effects in high-enthalpy tunnel [D]. Beijing: Institute of Mechanics, Chinese Academy of Sciences, 2007. (in Chinese)

    [19] Zeng M, Xu D, Liu J. Novel method to calculate vibrational thermal conduction in hypersonic nonequilibrium flow [J]. Journal of Thermophysics and Heat Transfer, 2016, 30(1): 12-24.

    [20] Shu C W, Osher S. Efficient implementation of essentially non-oscillatory shock capturing schemes[J]. Journal of Computational Physics, 1988, 77(2): 493-471.

    [21] 柳軍. 熱化學(xué)非平衡流及其輻射現(xiàn)象的實(shí)驗(yàn)和數(shù)值計(jì)算研究[D].長沙:國防科學(xué)技術(shù)大學(xué), 2004.LIU Jun. Experimental and numerical research on thermo-chemical nonequilibrium flow with radiation phenomenon[D]. Changsha: National University of Defense Technology, 2004. (in Chinese)

    [22] Shu C W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws [R]. NASA/CR-97-206253, 1997: 19-26

    [23] 侯中喜. 超聲速復(fù)雜流場(chǎng)并行數(shù)值分析及高階格式研究[D].長沙: 國防科學(xué)技術(shù)大學(xué), 2000.

    HOU Zhongxi. Parallel numerical analysis of hypersonic complex flow and research of high-oder schemes[D]. Changsha: National University of Defense Technology, 2000. (in Chinese)

    [24] Karl S, Martinez-Schramm J, Hannemann K. High enthalpy cylinder flow in HEG: a basis for CFD validation [C]//Proceedings of 33rd AIAA Fluid Dynamics Conference and Exhibit,AIAA 2003-4252, 2003.

    [25] Knight D, Longo J, Drikakis D, et al.Assessment of CFD capability for prediction of hypersonic shock interactions [J]. Progress in Aerospace Sciences, 2012, 48/49: 8-26.

    [26] 董義道. WCNS高階格式在典型流動(dòng)狀態(tài)下的應(yīng)用研究[D]. 長沙:國防科學(xué)技術(shù)大學(xué), 2014.

    DONG Yidao. Application research of weighted compact high-oder nonlinear schemes in some typical flow conditions[D]. Changsha: National University of Defense Technology, 2014.(in Chinese)

    Applications of high-order weighted compact nonlinear scheme for thermo-chemical nonequilibrium flow

    GE Mingming, ZENG Ming, ZHAO Xiaoyu

    (College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China)

    High-order scheme has advantages in the simulations of complex flow for its good behavior of preserving high spatial resolution with lower cost than the second order methods. WCNS-E-5(weighted compact nonlinear scheme) was applied in the simulations of thermo-chemical nonequilibrium flow around a cylinder with the freestream velocity from 5 km/s to 5.7 km/s. The distribution of correct flow field parameter was acquired, the pressure and heat flux on the surface agreed well with the experimental data. WCNS-E-5 also exhibited better performance on the grid convergence of heat flux than the second order MUSCL scheme. Results show that the present attempt to combine high order scheme with the calculations of stiff chemically reacting flow is basically successful, and it has established a good foundation for the further study.

    high-order scheme; thermo-chemical nonequilibrium flow; numerical simulation

    10.11887/j.cn.201605026

    http://journal.nudt.edu.cn

    2016-03-21

    國家自然科學(xué)基金資助項(xiàng)目(11572348);國防科技大學(xué)重大應(yīng)用基礎(chǔ)研究資助項(xiàng)目(ZDYYJCYJ20140101)

    葛明明(1991—),男,江蘇南通人,博士研究生,E-mail:owen2024@sina.cn;曾明(通信作者),女,副教授,博士,碩士生導(dǎo)師,E-mail:ming_z@163.com

    V211.3

    A

    1001-2486(2016)05-163-07

    猜你喜歡
    熱化學(xué)駐點(diǎn)激波
    “熱化學(xué)方程式”知識(shí)揭秘
    一種基于聚類分析的二維激波模式識(shí)別算法
    基于HIFiRE-2超燃發(fā)動(dòng)機(jī)內(nèi)流道的激波邊界層干擾分析
    稠油硫酸鹽熱化學(xué)還原生成H2S實(shí)驗(yàn)研究
    斜激波入射V形鈍前緣溢流口激波干擾研究
    基于游人游賞行為的留園駐點(diǎn)分布規(guī)律研究
    中國園林(2018年7期)2018-08-07 07:07:48
    適于可壓縮多尺度流動(dòng)的緊致型激波捕捉格式
    高超聲速熱化學(xué)非平衡對(duì)氣動(dòng)熱環(huán)境影響
    例析熱化學(xué)中焓變計(jì)算的常見題型
    利用遠(yuǎn)教站點(diǎn),落實(shí)駐點(diǎn)干部帶學(xué)
    日韩精品青青久久久久久| av.在线天堂| 国产又黄又爽又无遮挡在线| 最近视频中文字幕2019在线8| 欧美又色又爽又黄视频| 欧美+日韩+精品| av专区在线播放| 国产一级毛片在线| 国产人妻一区二区三区在| 69人妻影院| 婷婷色综合大香蕉| 我的老师免费观看完整版| 成人永久免费在线观看视频| 级片在线观看| 亚洲欧美日韩高清专用| 久久久精品欧美日韩精品| 久久亚洲精品不卡| 人体艺术视频欧美日本| 国产午夜精品一二区理论片| 中文字幕精品亚洲无线码一区| 免费大片18禁| 国产精品,欧美在线| 国产伦在线观看视频一区| 久久人人爽人人片av| 乱人视频在线观看| 五月玫瑰六月丁香| 久久久久久久久大av| 两个人的视频大全免费| 亚洲欧美成人精品一区二区| 亚洲精品456在线播放app| 两个人视频免费观看高清| 亚洲精品乱码久久久久久按摩| 在线a可以看的网站| 九九爱精品视频在线观看| 亚洲成a人片在线一区二区| 好男人在线观看高清免费视频| 亚洲精品影视一区二区三区av| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 午夜精品国产一区二区电影 | 午夜精品国产一区二区电影 | 尾随美女入室| 日韩视频在线欧美| 别揉我奶头 嗯啊视频| 夜夜夜夜夜久久久久| 一级毛片久久久久久久久女| 色哟哟·www| 国产高潮美女av| 高清午夜精品一区二区三区 | 成年女人永久免费观看视频| 国内揄拍国产精品人妻在线| 国产不卡一卡二| 国产午夜福利久久久久久| 亚洲,欧美,日韩| 国产日韩欧美在线精品| 中文字幕人妻熟人妻熟丝袜美| 一进一出抽搐动态| 一本久久精品| 久久久久九九精品影院| 美女黄网站色视频| 中文字幕av成人在线电影| 免费观看在线日韩| .国产精品久久| 18+在线观看网站| 青春草视频在线免费观看| 午夜福利成人在线免费观看| 热99re8久久精品国产| 久久精品夜色国产| 精品无人区乱码1区二区| 亚洲成a人片在线一区二区| 热99re8久久精品国产| 国产精品久久久久久精品电影小说 | 2022亚洲国产成人精品| 国产老妇女一区| 变态另类成人亚洲欧美熟女| 女同久久另类99精品国产91| 免费无遮挡裸体视频| 可以在线观看毛片的网站| 国产av麻豆久久久久久久| 久久久久久久亚洲中文字幕| 国产在线精品亚洲第一网站| 国产又黄又爽又无遮挡在线| 丝袜喷水一区| 色尼玛亚洲综合影院| 亚洲av电影不卡..在线观看| 亚洲第一电影网av| 午夜精品一区二区三区免费看| 欧美一级a爱片免费观看看| 真实男女啪啪啪动态图| 欧美高清性xxxxhd video| 麻豆av噜噜一区二区三区| 国产精品一二三区在线看| 99热这里只有是精品在线观看| 午夜激情福利司机影院| 亚洲精品456在线播放app| 99热全是精品| 国产亚洲欧美98| 久久精品久久久久久久性| 成人鲁丝片一二三区免费| 99热全是精品| 热99re8久久精品国产| 干丝袜人妻中文字幕| 最近2019中文字幕mv第一页| 精品免费久久久久久久清纯| 久久久久久久午夜电影| 青青草视频在线视频观看| 免费人成视频x8x8入口观看| 夜夜爽天天搞| 人妻少妇偷人精品九色| 日韩 亚洲 欧美在线| 秋霞在线观看毛片| 嘟嘟电影网在线观看| 极品教师在线视频| 午夜福利成人在线免费观看| 蜜桃久久精品国产亚洲av| 亚洲美女搞黄在线观看| 男女那种视频在线观看| 啦啦啦韩国在线观看视频| 在线国产一区二区在线| 欧美日本视频| 黄色一级大片看看| 免费无遮挡裸体视频| 人体艺术视频欧美日本| 亚洲成人久久性| 午夜福利在线在线| 精品人妻视频免费看| 久久精品国产亚洲av香蕉五月| 国产精品av视频在线免费观看| 亚洲综合色惰| 丰满人妻一区二区三区视频av| 麻豆国产97在线/欧美| 麻豆国产97在线/欧美| 我要搜黄色片| av在线播放精品| 国产精品精品国产色婷婷| 国国产精品蜜臀av免费| 日日干狠狠操夜夜爽| 97热精品久久久久久| 禁无遮挡网站| 亚洲人与动物交配视频| 日本在线视频免费播放| 欧美潮喷喷水| 在线播放无遮挡| 国内精品久久久久精免费| av在线亚洲专区| 国产午夜精品论理片| 成人欧美大片| 小蜜桃在线观看免费完整版高清| 亚洲国产日韩欧美精品在线观看| 久久久久国产网址| 全区人妻精品视频| 日本免费一区二区三区高清不卡| 亚洲自拍偷在线| 亚洲美女视频黄频| 一区福利在线观看| 国产精品人妻久久久影院| 免费看美女性在线毛片视频| 亚洲五月天丁香| 在线免费十八禁| 免费看日本二区| 日本在线视频免费播放| 亚洲第一电影网av| 午夜视频国产福利| 熟妇人妻久久中文字幕3abv| h日本视频在线播放| 国内少妇人妻偷人精品xxx网站| 成人三级黄色视频| www日本黄色视频网| 99热只有精品国产| 亚洲国产欧美人成| 青青草视频在线视频观看| 亚洲aⅴ乱码一区二区在线播放| 给我免费播放毛片高清在线观看| 亚洲国产色片| 国产单亲对白刺激| 99热这里只有是精品50| 免费av观看视频| 亚洲,欧美,日韩| 国产免费男女视频| 在线播放无遮挡| 国产精品99久久久久久久久| 日本黄色视频三级网站网址| 亚洲中文字幕一区二区三区有码在线看| 成人毛片a级毛片在线播放| 亚洲av中文av极速乱| 高清日韩中文字幕在线| 少妇裸体淫交视频免费看高清| 久久午夜亚洲精品久久| 男人狂女人下面高潮的视频| 久久精品国产鲁丝片午夜精品| 国产精品电影一区二区三区| 欧美激情久久久久久爽电影| 午夜福利在线观看免费完整高清在 | 免费观看在线日韩| 免费观看a级毛片全部| 色哟哟·www| 日韩av在线大香蕉| 久久精品国产亚洲网站| 国产又黄又爽又无遮挡在线| 国产精品.久久久| 欧美日韩综合久久久久久| av天堂在线播放| 在线免费观看的www视频| 亚洲欧美日韩卡通动漫| 日韩高清综合在线| 久久久久性生活片| 成人亚洲精品av一区二区| 亚洲丝袜综合中文字幕| 久久久久久久久中文| 国产精品99久久久久久久久| 夫妻性生交免费视频一级片| 国内少妇人妻偷人精品xxx网站| 97人妻精品一区二区三区麻豆| 精品午夜福利在线看| 国国产精品蜜臀av免费| 99riav亚洲国产免费| 亚洲av不卡在线观看| 小蜜桃在线观看免费完整版高清| 99视频精品全部免费 在线| 亚洲自拍偷在线| 99在线人妻在线中文字幕| 欧美xxxx黑人xx丫x性爽| 美女 人体艺术 gogo| 亚洲美女搞黄在线观看| 久久人妻av系列| 在线播放无遮挡| 寂寞人妻少妇视频99o| 久久99热这里只有精品18| a级毛片a级免费在线| 久久午夜福利片| 亚洲国产精品成人久久小说 | 22中文网久久字幕| 免费av毛片视频| 69人妻影院| 久久人人爽人人片av| 黄色欧美视频在线观看| 日本熟妇午夜| 最好的美女福利视频网| 国产一区二区在线av高清观看| 给我免费播放毛片高清在线观看| 午夜精品国产一区二区电影 | 精品人妻熟女av久视频| 成人鲁丝片一二三区免费| 精品久久久久久久久久久久久| 亚洲在久久综合| 亚洲真实伦在线观看| 亚洲自偷自拍三级| 成人无遮挡网站| 亚洲国产精品久久男人天堂| 亚洲最大成人中文| 神马国产精品三级电影在线观看| 国产精品一及| 18禁在线播放成人免费| 久99久视频精品免费| 男女那种视频在线观看| 亚洲av二区三区四区| 精品人妻视频免费看| 乱人视频在线观看| 麻豆国产av国片精品| 中文在线观看免费www的网站| 九草在线视频观看| 91午夜精品亚洲一区二区三区| 亚洲av免费在线观看| 亚洲精品国产av成人精品| 中文字幕制服av| 成人高潮视频无遮挡免费网站| 日韩精品青青久久久久久| 欧美一区二区国产精品久久精品| 一进一出抽搐动态| 一级毛片aaaaaa免费看小| 你懂的网址亚洲精品在线观看 | 亚洲成人中文字幕在线播放| 一级二级三级毛片免费看| 国产视频内射| 少妇的逼水好多| 亚洲激情五月婷婷啪啪| 天堂√8在线中文| 99riav亚洲国产免费| 国产乱人偷精品视频| 国产精品一区二区在线观看99 | 日韩av在线大香蕉| 婷婷色av中文字幕| 国产精品不卡视频一区二区| 神马国产精品三级电影在线观看| 欧美性猛交黑人性爽| 国产精品美女特级片免费视频播放器| 大又大粗又爽又黄少妇毛片口| 全区人妻精品视频| 久久久成人免费电影| 久久韩国三级中文字幕| 欧美成人免费av一区二区三区| 在线播放无遮挡| 午夜福利在线在线| 观看美女的网站| 国产精品久久久久久久久免| 熟妇人妻久久中文字幕3abv| 深夜a级毛片| 一个人免费在线观看电影| 麻豆av噜噜一区二区三区| 搞女人的毛片| 3wmmmm亚洲av在线观看| 99国产精品一区二区蜜桃av| 美女高潮的动态| 久久精品国产亚洲av涩爱 | 99热全是精品| 国产午夜精品论理片| 黄片无遮挡物在线观看| 久久久精品欧美日韩精品| 桃色一区二区三区在线观看| 国产午夜精品论理片| 哪个播放器可以免费观看大片| 91在线精品国自产拍蜜月| 男人狂女人下面高潮的视频| 国产伦在线观看视频一区| 日本与韩国留学比较| 亚洲欧美日韩卡通动漫| 高清在线视频一区二区三区 | 女人被狂操c到高潮| 国产探花极品一区二区| 欧美日韩精品成人综合77777| 久久6这里有精品| 国产乱人偷精品视频| 午夜福利在线观看免费完整高清在 | av在线观看视频网站免费| 综合色丁香网| 精品99又大又爽又粗少妇毛片| 久久久久久久久久久免费av| av国产免费在线观看| 日本av手机在线免费观看| 日本撒尿小便嘘嘘汇集6| 日韩国内少妇激情av| 国产一区二区三区av在线 | 国内精品久久久久精免费| 国产午夜福利久久久久久| 国产av不卡久久| 亚洲精品国产av成人精品| 日本一本二区三区精品| 不卡一级毛片| 成年女人永久免费观看视频| 国产老妇伦熟女老妇高清| 国产在线男女| 精华霜和精华液先用哪个| 免费人成视频x8x8入口观看| 久久精品国产鲁丝片午夜精品| 高清午夜精品一区二区三区 | 国产精品永久免费网站| 久久久久久久久久久免费av| 久久这里只有精品中国| 男人和女人高潮做爰伦理| 欧美潮喷喷水| 国产免费一级a男人的天堂| 中文精品一卡2卡3卡4更新| 非洲黑人性xxxx精品又粗又长| 国产成人freesex在线| 精品一区二区免费观看| 久久久色成人| 日韩欧美在线乱码| 人人妻人人澡欧美一区二区| 青青草视频在线视频观看| 蜜桃亚洲精品一区二区三区| 婷婷精品国产亚洲av| www.色视频.com| 日韩亚洲欧美综合| 亚洲精华国产精华液的使用体验 | 成人av在线播放网站| 可以在线观看毛片的网站| 亚洲精品亚洲一区二区| 国产黄色小视频在线观看| 国产av一区在线观看免费| 悠悠久久av| 日韩强制内射视频| 非洲黑人性xxxx精品又粗又长| 中文资源天堂在线| 草草在线视频免费看| 99久久精品国产国产毛片| 精品国内亚洲2022精品成人| 午夜福利在线在线| 女人十人毛片免费观看3o分钟| 蜜桃久久精品国产亚洲av| 亚洲国产高清在线一区二区三| 国产黄a三级三级三级人| 久久精品夜色国产| www.av在线官网国产| 美女国产视频在线观看| 最近2019中文字幕mv第一页| 熟妇人妻久久中文字幕3abv| 国产精品,欧美在线| 成人性生交大片免费视频hd| 尾随美女入室| 男人舔女人下体高潮全视频| 18禁裸乳无遮挡免费网站照片| 中文字幕久久专区| 中文欧美无线码| 97热精品久久久久久| 亚洲三级黄色毛片| 黄片无遮挡物在线观看| 99久久精品国产国产毛片| 日韩一区二区视频免费看| 91av网一区二区| 一本久久精品| 18+在线观看网站| 男人舔女人下体高潮全视频| 99久久无色码亚洲精品果冻| 一级毛片aaaaaa免费看小| 精华霜和精华液先用哪个| 一边摸一边抽搐一进一小说| 精品久久久久久久人妻蜜臀av| 欧美日本亚洲视频在线播放| 男插女下体视频免费在线播放| 国产伦精品一区二区三区四那| 国产真实伦视频高清在线观看| 男人和女人高潮做爰伦理| 91狼人影院| 日韩三级伦理在线观看| a级毛色黄片| 精品熟女少妇av免费看| 欧美性感艳星| 级片在线观看| 一区二区三区高清视频在线| 日日干狠狠操夜夜爽| 免费观看精品视频网站| www日本黄色视频网| 国产精品久久视频播放| 18+在线观看网站| 久久久成人免费电影| 精品一区二区三区视频在线| 国产在线男女| 在线国产一区二区在线| 日本在线视频免费播放| 能在线免费观看的黄片| 午夜激情欧美在线| 国产精品爽爽va在线观看网站| 男女啪啪激烈高潮av片| 久久精品国产亚洲av香蕉五月| 99久久无色码亚洲精品果冻| 国产欧美日韩精品一区二区| 99久久精品热视频| 狂野欧美白嫩少妇大欣赏| 国产免费男女视频| 亚洲av男天堂| 51国产日韩欧美| 久久久久久久久久久免费av| 97在线视频观看| 高清在线视频一区二区三区 | 成熟少妇高潮喷水视频| 九九在线视频观看精品| 我要看日韩黄色一级片| 内地一区二区视频在线| 天堂中文最新版在线下载 | 久久亚洲精品不卡| 简卡轻食公司| 国语自产精品视频在线第100页| 日韩强制内射视频| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 亚洲国产欧美在线一区| 婷婷色av中文字幕| 欧美精品一区二区大全| 18禁在线无遮挡免费观看视频| 亚洲丝袜综合中文字幕| 搡女人真爽免费视频火全软件| 午夜激情欧美在线| 在线观看免费视频日本深夜| 免费大片18禁| 在线播放国产精品三级| 蜜臀久久99精品久久宅男| 99久国产av精品国产电影| 看免费成人av毛片| 久久久久久久午夜电影| 久久久a久久爽久久v久久| www日本黄色视频网| 国产片特级美女逼逼视频| 免费无遮挡裸体视频| 嫩草影院精品99| 晚上一个人看的免费电影| 国产精品一区二区性色av| 亚洲婷婷狠狠爱综合网| 午夜免费男女啪啪视频观看| 99热只有精品国产| 免费黄网站久久成人精品| 久久人人精品亚洲av| 熟女人妻精品中文字幕| 免费无遮挡裸体视频| 欧美一区二区亚洲| 啦啦啦韩国在线观看视频| 91麻豆精品激情在线观看国产| 成年av动漫网址| 成人三级黄色视频| 久久精品国产亚洲网站| 日本免费一区二区三区高清不卡| 最后的刺客免费高清国语| 国产精品一区二区在线观看99 | 久久草成人影院| 1000部很黄的大片| 97热精品久久久久久| 国产精品.久久久| 综合色丁香网| 日日摸夜夜添夜夜添av毛片| 我的女老师完整版在线观看| 国产91av在线免费观看| 亚洲激情五月婷婷啪啪| 久久午夜亚洲精品久久| 欧美最新免费一区二区三区| 偷拍熟女少妇极品色| 美女被艹到高潮喷水动态| 日韩欧美三级三区| 99热这里只有精品一区| 久久99精品国语久久久| 国产麻豆成人av免费视频| 成人毛片60女人毛片免费| 久久国内精品自在自线图片| 婷婷精品国产亚洲av| 嫩草影院新地址| 91狼人影院| 丰满乱子伦码专区| 天堂网av新在线| 久久中文看片网| 精品久久久久久久久久久久久| 热99在线观看视频| 女的被弄到高潮叫床怎么办| 最好的美女福利视频网| 精品欧美国产一区二区三| 在线观看免费视频日本深夜| 嫩草影院入口| 老熟妇乱子伦视频在线观看| 国产黄片视频在线免费观看| 啦啦啦观看免费观看视频高清| 丰满人妻一区二区三区视频av| 精品久久久久久久久久久久久| 亚洲美女视频黄频| 女的被弄到高潮叫床怎么办| 最好的美女福利视频网| 精品无人区乱码1区二区| 午夜福利高清视频| 国产片特级美女逼逼视频| 亚洲中文字幕一区二区三区有码在线看| 插阴视频在线观看视频| 看十八女毛片水多多多| 丝袜美腿在线中文| 国产毛片a区久久久久| 狠狠狠狠99中文字幕| 日韩国内少妇激情av| 午夜亚洲福利在线播放| 日韩欧美在线乱码| 91久久精品国产一区二区三区| 大又大粗又爽又黄少妇毛片口| 亚洲国产欧洲综合997久久,| 精品国内亚洲2022精品成人| 99在线人妻在线中文字幕| 国产精品野战在线观看| 国产一区二区亚洲精品在线观看| 精品国内亚洲2022精品成人| 国产国拍精品亚洲av在线观看| 青春草视频在线免费观看| 亚洲国产欧美在线一区| 国产欧美日韩精品一区二区| 变态另类丝袜制服| 欧美最黄视频在线播放免费| 又粗又硬又长又爽又黄的视频 | 国产白丝娇喘喷水9色精品| 黄色视频,在线免费观看| 国产亚洲av片在线观看秒播厂 | 一本久久中文字幕| 欧美日韩一区二区视频在线观看视频在线 | 日韩一区二区三区影片| 又粗又爽又猛毛片免费看| 日本撒尿小便嘘嘘汇集6| 久久99蜜桃精品久久| 嘟嘟电影网在线观看| 26uuu在线亚洲综合色| 高清在线视频一区二区三区 | 两个人的视频大全免费| 在现免费观看毛片| 国内精品美女久久久久久| 日韩亚洲欧美综合| 国产乱人视频| 秋霞在线观看毛片| 日韩欧美精品v在线| 成熟少妇高潮喷水视频| 在线观看一区二区三区| 亚洲第一电影网av| 免费电影在线观看免费观看| 国产成人a区在线观看| 亚洲精品粉嫩美女一区| 18+在线观看网站| 午夜亚洲福利在线播放| 伦理电影大哥的女人| 在线观看午夜福利视频| 欧美三级亚洲精品| 能在线免费观看的黄片| 最好的美女福利视频网| 欧美成人精品欧美一级黄| 国产精品av视频在线免费观看| 在线播放国产精品三级| 夜夜爽天天搞| 黄色一级大片看看| 国产精品久久视频播放| av.在线天堂| 日韩 亚洲 欧美在线| 人人妻人人澡欧美一区二区| 欧美成人精品欧美一级黄| 国产成人91sexporn| 国产色爽女视频免费观看| 中文精品一卡2卡3卡4更新| a级一级毛片免费在线观看| 嫩草影院精品99| 亚洲最大成人中文| 欧美xxxx黑人xx丫x性爽| 你懂的网址亚洲精品在线观看 | 欧美色视频一区免费| 边亲边吃奶的免费视频| 精品欧美国产一区二区三| 日韩欧美 国产精品| 久久99热6这里只有精品| 亚洲精品成人久久久久久| 亚洲图色成人| 国产亚洲91精品色在线|