999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

高超聲速飛行器跨流域氣動力/熱預(yù)測技術(shù)研究

2016-04-06 03:02:38楊彥廣李中華李緒國戴金雯
空氣動力學(xué)學(xué)報 2016年1期
關(guān)鍵詞:測量方法

楊彥廣,李 明,李中華,李緒國,戴金雯

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

高超聲速飛行器跨流域氣動力/熱預(yù)測技術(shù)研究

楊彥廣*,李 明,李中華,李緒國,戴金雯

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

分析了跨流域流動特點、工程背景需求、國內(nèi)外研究現(xiàn)狀,針對高超聲速飛行器跨流域氣動力/熱預(yù)測存在的主要問題,從粘性干擾參數(shù)與克努森數(shù)搭接的跨流域模擬準(zhǔn)則、微量天平結(jié)構(gòu)優(yōu)化及測力技術(shù)、紅外大面積中低量值熱流測量技術(shù)、流場顯示與測量技術(shù)和N-S/DSMC自適應(yīng)緊耦合數(shù)值模擬方法等方面,總結(jié)了近年來取得的研究進展,并討論了下一步研究方向。

跨流域;氣動力;氣動熱;流動顯示;耦合算法

0 引 言

跨流域一般是指介于連續(xù)流和自由分子流之間的流動區(qū)域,從氣體動力學(xué)角度看,主要包括黏性主導(dǎo)的近連續(xù)滑移流和稀薄過渡流,從空域范圍看,與臨近空間基本重疊。跨流域流動的雷諾數(shù)很低而馬赫數(shù)很高,流動往往是層流,而且粘性層很厚,脫體激波后的極大熵梯度明顯影響粘性流的結(jié)構(gòu)。物面對流動的影響可以傳播到遠(yuǎn)離物面的流動中,使粘性影響向外擴張并與外圍無粘流動相互作用,并導(dǎo)致二者之間的界限難以嚴(yán)格劃分。同時又因馬赫數(shù)很高,激波后的氣體分子狀態(tài)遠(yuǎn)離平衡態(tài),不僅氣體分子的平動溫度、轉(zhuǎn)動溫度和振動溫度存在著顯著差別,有些情況下還存在氣體分子的化學(xué)反應(yīng)。這種稀薄氣體效應(yīng)、強粘性效應(yīng)和高溫真實氣體效應(yīng)高度耦合的流動,已成為稀薄氣體動力學(xué)的前沿研究課題。

雖然飛行器在該區(qū)域的氣動力、熱載荷比稠密大氣層內(nèi)低,但大氣的作用依然顯著,且飛行器周圍氣體的流動特性與連續(xù)流相比有了顯著差異,目前尚缺少系統(tǒng)可靠的跨流域氣動預(yù)測方法。傳統(tǒng)的彈道再入由于時間短、載荷小而將這類氣動問題忽略或進行保守處理,從而回避了預(yù)測的困難。而新型高超聲速飛行器具有非彈道、高可控、長時間飛行的特點,對跨流域氣動特性的精細(xì)預(yù)測提出迫切需求。美國HTV2等新型飛行器試飛失敗進一步表明,跨流域飛行的稀薄氣體與高溫真實氣體效應(yīng)對于飛行器的發(fā)展而言仍然有很多科學(xué)問題亟待深入研究。另外,美國X-37B、X-51A的試飛成功極大地激發(fā)了各國對于開發(fā)臨近空間飛行器的熱情,使跨流域問題成為當(dāng)今熱門研究課題之一。

國內(nèi)外對跨流域的研究重點主要集中在以下幾個方面:地球大氣環(huán)境稀薄氣體效應(yīng)對飛行器氣動力/熱特性的影響、過渡區(qū)RCS噴流干擾對飛行器氣動特性及控制效率的影響、行星(火星)大氣環(huán)境稀薄氣動特性研究、稀薄氣體高溫?zé)峄瘜W(xué)非平衡效應(yīng)及分子內(nèi)部能態(tài)的激發(fā)與傳遞建模研究等。

本文在分析國內(nèi)外研究現(xiàn)狀的基礎(chǔ)上,給出了風(fēng)洞試驗、數(shù)值模擬兩個方面近十多年的研究進展。主要包括:跨流域氣動力試驗?zāi)M準(zhǔn)則、微量天平結(jié)構(gòu)優(yōu)化和測力技術(shù)、紅外大面積中低熱流測量技術(shù)、流場顯示測量技術(shù)、過渡流區(qū)復(fù)雜流動的N-S/DSMC耦合算法等。最后討論了下一步研究方向。

1 國內(nèi)外研究現(xiàn)狀

開展高聲速飛行器跨流域氣動力/熱特性研究主要有地面試驗、數(shù)值模擬和工程估算三類手段。其中,工程估算主要采用當(dāng)?shù)鼗椒▽橛谶B續(xù)流與自由分子流之間的過渡區(qū)域進行搭橋,由于其計算結(jié)果置信度不高、橋函數(shù)依賴于風(fēng)洞試驗結(jié)果修正、對不同外形的適應(yīng)性差,一般只用于簡單外形或飛行器概念設(shè)計的氣動估算,本文不再贅述。由于地面試驗存在靜溫低、雷諾數(shù)低、有洞壁干擾、支架干擾等,而數(shù)值模擬又需要風(fēng)洞試驗數(shù)據(jù)驗證,因此高超聲速飛行器跨流域氣動力/熱預(yù)測研究有賴于多種手段的結(jié)合。

1.1 地面試驗研究方面

地面試驗方面,主要利用低密度風(fēng)洞、真空艙設(shè)備等開展試驗研究。德國DLR高超聲速真空風(fēng)洞中的v1g、v2g試驗段主要進行飛行器70~120km高度范圍的氣動力/氣動熱試驗、反推力控制裝置射流試驗,及高馬赫數(shù)稀薄氣體流動的基礎(chǔ)研究。所發(fā)展的測試技術(shù)包括利用三分量、五分量應(yīng)變天平及電磁阻力天平進行稀薄流區(qū)氣動力測量;采用瞬態(tài)薄膜技術(shù)及熱電偶技術(shù)進行稀薄流區(qū)氣動熱測量;采用微壓傳感器和真空計進行模型表面壓力測量;配置電子束流場診斷系統(tǒng),可開展流場顯示、流場密度測量、轉(zhuǎn)動溫度和振動溫度測量等試驗。法國的SR3低密度風(fēng)洞配有電子槍進行定量的密度測量和流動顯示,用外支架天平進行廣泛的氣動載荷測量,還配置有壁面壓力和熱流計、紅外測熱系統(tǒng)、干涉儀和紋影儀等。SR3低密度風(fēng)洞的雷諾數(shù)覆蓋范圍廣,覆蓋了從近自由分子流到連續(xù)流,擴展到強干擾和混合層區(qū)。

俄羅斯中央流體動力研究院的T-117高超聲速風(fēng)洞、中央通用機械研究院的U-16高超聲速低密度風(fēng)洞、ITAM的T-327高超聲速低密度風(fēng)洞等配置有天平、熱敏涂料、紅外測熱系統(tǒng)和電子束熒光流場診斷系統(tǒng)等測試設(shè)備。美國AEDC的L、M兩座低密度風(fēng)洞M-Re數(shù)覆蓋范圍較寬,并與其他高超聲速風(fēng)洞配合,滿足了其航天飛行器的試驗需求。R.J.Vidal對高超聲速楔形流中低密度效應(yīng)進行了試驗研究,將試驗獲得的熱流結(jié)果與粘性邊界層理論數(shù)據(jù)進行了比較,試驗獲得的壓力結(jié)果與邊界層位移理論、楔形流理論數(shù)據(jù)進行了對比,研究表明低密度效應(yīng)在壓力數(shù)據(jù)中是明顯的[1]。

中國空氣動力研究與發(fā)展中心(CARDC)擁有國內(nèi)唯一可以進行高空稀薄氣動力/熱試驗的高超聲速低密度風(fēng)洞。經(jīng)過三十多年的不斷發(fā)展,逐步形成了高空氣動力/氣動熱稀薄氣動特性試驗技術(shù)、高空羽流及噴流干擾試驗技術(shù)。唐志共、楊彥廣、戴金雯、李緒國等研制發(fā)展了多種結(jié)構(gòu)形式和量程、從二分量到六分量的系列微量天平,研制了國際上首套體軸系自動復(fù)位的微量天平校準(zhǔn)系統(tǒng),先后建立了錐-桿組合體、大鈍頭體、多體干擾、復(fù)雜升力體等的稀薄過渡流區(qū)氣動力試驗方法[2-8]。張金濤、華威等針對鈍錐外形開展了薄壁量熱技術(shù)與液晶熱圖技術(shù),李明、祝智偉等人在國內(nèi)率先建立了紅外熱圖測熱試驗技術(shù),并成功應(yīng)用于大升阻比復(fù)雜外形的稀薄氣動熱測量試驗[9-11]。楊彥廣、陳愛國等建立了平板橫向噴流、姿控發(fā)動機羽流干擾試驗技術(shù);李明、廖俊必、李震乾等開展了羽流撞擊平板試驗技術(shù)研究,并利用輝光放電裝置進行了流場顯示,發(fā)展了低壓力測量技術(shù)[]。

盡管如此,地面試驗還需要進一步完善試驗?zāi)M準(zhǔn)則,以實現(xiàn)從近連續(xù)流到稀薄過渡流跨流動區(qū)域的氣動力/熱試驗?zāi)M和數(shù)據(jù)關(guān)聯(lián);針對小尺寸模型內(nèi)腔狹小和氣動載荷小的特點,解決小尺寸、小量程內(nèi)式六分量應(yīng)變天平結(jié)構(gòu)布局難和靈敏度低的技術(shù)難題;針對中低熱流測量問題,解決局部干擾區(qū)與駐點熱流的精細(xì)測量問題;針對跨流域模擬高度寬的特點,需要建立多種流場顯示與測量手段。

1.2 數(shù)值模擬研究方面

數(shù)值模擬方面,主要采用三種方法:DSMC、N-S/DSMC耦合、基于Boltzmann及其模型方程的數(shù)值模擬方法。

DSMC方法不僅能夠模擬三維稀薄氣體流動的復(fù)雜流場,而且能夠真實地模擬包括熱化學(xué)非平衡效應(yīng)以及熱輻射等物理化學(xué)過程在內(nèi)的稀薄氣體流動問題,其結(jié)果得到微觀細(xì)節(jié)和宏觀氣動力試驗的驗證,是公認(rèn)可用于復(fù)雜工程問題研究預(yù)測的稀薄流數(shù)值模擬方法。這種方法發(fā)展較早,可查閱的文獻較多[1344]。直接求解Boltzmann方程一直是學(xué)者們追求的目標(biāo),但其解析求解極其困難,目前數(shù)值求解方法主要有矩方程方法、離散速度坐標(biāo)法、模型方程方法、有限差分方法、BGK格式方法等,可參閱香港科技大學(xué)的徐昆[45]、CARDC的李志輝等開展的研究工作[46-49]。DSMC方法與基于Boltzmann及其模型方程的數(shù)值模擬方法不是本文討論的重點,下面僅對N-S/DSMC耦合方法國內(nèi)外研究現(xiàn)狀作一回顧。

為開展較低高度過渡流區(qū)(通常為90km以下)流動的數(shù)值模擬,研究者們進行了各種嘗試。其中N-S/DSMC耦合算法得到快速發(fā)展,國內(nèi)外學(xué)者開展了多種耦合算法的研究。該算法基于N-S方程與DSMC方法的分區(qū)耦合模擬,擴展這兩種數(shù)值方法的應(yīng)用范圍,并使計算效率在工程上可接受,是近期有望率先達到工程實用化的跨流域數(shù)值模擬方法。Hash和Hassan采用松耦合方式對高超聲速鈍錐繞流進行了模擬,這種耦合方式僅適合分界面一次性確定并且界面上參數(shù)不變的情況[50]。另一種耦合方法是緊耦合方式,CFD和DSMC計算區(qū)域時刻進行流場信息交換。Wadsworth和Erwin最先利用參數(shù)插值方法開展這種緊耦合方法的研究[51]。Hash和Hassan采用Marshak邊界條件[52],建立了一種耦合方法,應(yīng)用于Couette流和高超聲速繞流。Marshak邊界在耦合邊界上采用通量守恒邊界條件,通量一半來自CFD,一半來自DSMC。為了發(fā)展耦合算法,后加濾鏡方法來降低DSMC統(tǒng)計波動的影響。Oh和Oran研究表明通量修正傳輸(FCT)方法能夠濾掉高頻波動,信息通過雙線性插值方法從DSMC向N-S區(qū)域傳遞[53]。Thomas提出了一種MPC(modular particle-continuum)耦合技術(shù)[54],N-S方程和DSMC方法幾乎不用修改就可以耦合在一個程序中。在信息交換過程中,采用亞松弛技術(shù),可以有效抑制DSMC統(tǒng)計波動對CFD計算的影響,隨后,該耦合方法被推廣到熱、化學(xué)非平衡流動中。Timothy對MPC耦合方法實現(xiàn)并行化,進一步提高了計算效率[55]。李中華根據(jù)MPC耦合思想,建立了能夠自動分區(qū)的耦合算法,并在三維氣動特性預(yù)測中進行了驗證[56-58]。N-S/DSMC耦合算法需要解決CFD方法與DSMC方法自動分區(qū)與信息交換問題。

2 主要難點與研究進展

2.1 跨流域氣動力試驗?zāi)M準(zhǔn)則

按照錢學(xué)森提出的流態(tài)劃分準(zhǔn)則,流動分區(qū)如圖1所示。跨流域氣動力/熱試驗?zāi)M高度范圍跨度較大,試驗流場涵蓋近連續(xù)流、滑移流和過渡流。

圖1 流動分區(qū)圖Fig.1 Map of different flow region

在稀薄流區(qū),表征流動稀薄程度的模擬參數(shù)是克努森數(shù)Kn,其定義為平均分子自由程與流動的特征尺度之比。一般而言,120km以上為自由分子流區(qū)。隨著飛行高度的降低,逐步進入過渡流區(qū)(90km左右),激波和邊界層開始出現(xiàn),此時,表征流動的特征尺度由飛行器特征長度變?yōu)檫吔鐚雍穸龋琄n的計算公式變?yōu)椋?/p>

式中,Re∞為雷諾數(shù),L為飛行器特征長度,δ為邊界層厚度。

在連續(xù)流區(qū),飛行高度較低(高Re數(shù))時,飛行器繞流的邊界層很薄,粘性對流動的影響小。隨著飛行高度增加,邊界層增厚,與外部無粘流的相互作用不斷增強,粘性干擾逐漸成為流動的主導(dǎo)因素,流動狀態(tài)向近連續(xù)的滑移流和過渡流(60km左右)演化,此時,一般采用粘性干擾參數(shù)進行表征。該參數(shù)有多種表達形式,其中考慮邊界層內(nèi)參考溫度影響的第三粘性干擾參數(shù)V*在航天飛機氣動力系數(shù)的關(guān)聯(lián)中取得較好的一致性,可作為近連續(xù)流氣動模擬的相似參數(shù)。V*由下式確定:

式中,M∞為馬赫數(shù),C*為Chapman-Rubesin因子,ReL∞為雷諾數(shù)。

過渡流與連續(xù)流的分界,因飛行器的外形及特征尺度、馬赫數(shù)等因素而異,分界點高度不盡相同,一般在50~70km范圍。在不同流區(qū)分別采用克努森數(shù)、粘性干擾參數(shù)進行表征,可以對地面試驗的模擬狀態(tài)進行很好地搭接,實現(xiàn)從稀薄流到連續(xù)流的統(tǒng)一模擬。圖2給出了某乘波體飛行器在模型縮比J=0.05、馬赫數(shù)M∞=12條件下,兩種模擬參數(shù)的關(guān)聯(lián)結(jié)果,搭接高度為52km。

圖2 某乘波體不同模擬參數(shù)的模擬高度關(guān)聯(lián)圖Fig.2 Simulation altitude map using different simulation parameters

2.2 微量天平測力技術(shù)

測力天平是氣動力試驗的關(guān)鍵測試儀器,由機械結(jié)構(gòu)本體、應(yīng)變計和測量線路組成。研制用于飛行器跨流域氣動力測量的微量天平,其結(jié)構(gòu)主要存在以下難點:1)小尺寸的結(jié)構(gòu)限制;2)氣動載荷小與模型自重大的矛盾;3)滾轉(zhuǎn)力矩、軸向力單元結(jié)構(gòu)小與天平縱向剛度/強度要求矛盾突出;4)縱橫向載荷不匹配問題。此外,在高超聲速流動中,由于氣動加熱,試驗?zāi)P捅砻鏈囟壬撸€容易產(chǎn)生溫度效應(yīng)問題。

天平要完成氣動力的6分量分解,其結(jié)構(gòu)非常復(fù)雜,其測力元件結(jié)構(gòu)形式、元件間布局等方面都必須綜合考慮。天平元件采用對稱布局方式:將法向/側(cè)向力、俯仰/偏航力矩元件在兩端分段布置,使這4個分量既縮短了元件長度,又提高了靈敏度;將軸向力、滾轉(zhuǎn)力矩元件居中垂直布置,并采用多片梁結(jié)構(gòu),使這兩個分量既有較高的靈敏度,又降低了其他分量帶來的干擾。上述設(shè)計使天平的整體強度、剛度,各分量的應(yīng)力、應(yīng)變及靈敏度等性能比較均衡,同時有效減小了分量間干擾和單元熱輸出,降低了小尺寸天平加工難度。通過優(yōu)化,大大降低了天平的設(shè)計裕度,解決了天平靈敏度提高和設(shè)計風(fēng)險增加之間的矛盾,有效提高了天平測量的精準(zhǔn)度。

采取在天平與模型、支桿之間采取隔熱套、水冷、選擇溫度自補償?shù)闹袦貞?yīng)變計等方法,有效降低了溫度效應(yīng)對微量天平測量精度的影響。

通過上述優(yōu)化設(shè)計和反復(fù)實踐,研制出可實現(xiàn)復(fù)雜外形跨流域氣動力測量的尾支撐內(nèi)式天平。測量分量數(shù)從三分量擴展到六個分量,A、N、Mz的測量精度從3%提高到1%以內(nèi),Z、My、Mx達3%水平。圖3為典型的小量程六分量天平結(jié)構(gòu),圖4給出了某臨近空間飛行器(CAV)在不同俯仰舵偏角下軸向力系數(shù)CA隨迎角的變化數(shù)據(jù)。

圖3 天平結(jié)構(gòu)與網(wǎng)格劃分Fig.3 Framework of microbalance and partitioned grids

圖4 CAV外形CA隨迎角的變化Fig.4 CAvaried with attack of angle of CAV model

2.3 中低量值熱流紅外測量技術(shù)

中低量值熱流測量具有以下難點:1)模型表面熱流密度低。由于熱流載荷小,各種因素的干擾影響相對變大,這種小參數(shù)、高精度測量是一個棘手的問題。2)模型小。為了盡量模擬高空環(huán)境,通常需要盡量減小模型縮尺。采用小模型,使熱電偶布點少、安裝困難、測量精度難以保證。3)采用非接觸方法進行模型表面熱流測量,需要解決大極角發(fā)射率的修正、駐點區(qū)域熱流測量、模型物面與熱圖像素坐標(biāo)的對應(yīng)關(guān)系等問題。針對以上難點,采取了以下措施:采用紅外熱圖非接觸方法,一次試驗可以獲得連續(xù)的熱流分布圖譜,解決熱電偶測點不連續(xù)、數(shù)據(jù)量小且不便安裝的問題;通過降低模型背景噪聲、控制模型表面溫度范圍,提高低熱流測量精度;通過密封旋轉(zhuǎn)裝置及內(nèi)外氣流循環(huán)方式,確保紅外窗口冷卻與紅外熱像儀內(nèi)外氣壓平衡,提出主動適應(yīng)方法,解決大極角情況下模型表面發(fā)射率的修正問題;采用紅外反射光路對測量區(qū)域進行拓展和放大,解決試驗?zāi)P途植扛蓴_區(qū)和駐點區(qū)域熱流測量問題[10];采用參考網(wǎng)格法,確定模型物面與熱圖像素坐標(biāo)的對應(yīng)關(guān)系。分析紅外測熱各誤差項權(quán)重,建立基于泰勒級數(shù)誤差傳遞的測熱不確定度評估方法,指導(dǎo)改進重點誤差環(huán)節(jié)。通過紅外熱圖與熱電偶同時測量的平板斜坡薄壁模型熱流二者符合得較好,如圖5。圖6給出了測量方法改進前后的熱流測量區(qū)域,過去只能獲得半球圓柱模型表面θ角30°以后區(qū)域熱流,改進后可以得到整個區(qū)域熱流。半球圓柱模型大面積區(qū)域測量精度從20%提高到12%、駐點/局部區(qū)域從30%提高到20%。

圖5 平板斜坡模型紅外與熱電偶結(jié)果比較Fig.5 Comparison of heat transfer rates between with thermocouple and with infrared mapping

圖6 半球圓柱球頭測量區(qū)域過去與現(xiàn)在比較Fig.6 Comparison of h/h0distributions between in the past and at present on semi-cylinder testing model

2.4 流場顯示與測量技術(shù)

盡管流場顯示與測量技術(shù)發(fā)展很快,但對臨近空間跨流域模擬,高度從30km變化到100km左右,流場密度及靜壓相差幾個數(shù)量級,至今還沒有任何一種技術(shù)能完全滿足其流動顯示與測量的需要。為此,根據(jù)不同用途開發(fā)了適用于低密度流場的三種流動顯示與測量技術(shù):納米粒子平面激光散射(NPLS)測量技術(shù)、輝光放電流場顯示技術(shù)、低雷諾數(shù)油流顯示技術(shù)。

NPLS技術(shù)由國防科技大學(xué)易仕和教授研發(fā)并成功應(yīng)用于超聲速流動的顯示與測量,該技術(shù)與PIV相比具有粒子跟隨性好的優(yōu)勢,與PLIF相比又有流動成像信噪比高的優(yōu)點[59-60]。將該技術(shù)拓展運用于高超聲速流動,需要解決高超聲速試驗流場激光散射光信號微弱、流場示蹤粒子與試驗主氣流均勻混合等問題。當(dāng)流場密度較低時,可通過輝光放電流場、低雷諾數(shù)油流顯示技術(shù)對流場進行定性顯示。輝光放電技術(shù)結(jié)構(gòu)簡單,操作方便,便于推廣運用。

2.4.1 高超聲速流中的NPLS測量技術(shù)

高超聲速流場密度比駐室密度低2~3個數(shù)量級,導(dǎo)致隨氣流運動的TiO2納米示蹤粒子濃度也很低,散射光非常弱。例如在馬赫數(shù)M5試驗條件下,流場密度大約為駐室密度的1/88,而在馬赫數(shù)M10時,流場密度大約為駐室密度的1/2000。為了解決這個問題,在氣流中注入CO2氣體,CO2氣體通過噴管膨脹后迅速冷凝成納米量級的粒子,與TiO2納米粒子共同作為流場示蹤粒子,提高散射光強度。此外,通過加長管路、加大管徑、多孔注入等方式,解決流場示蹤粒子與試驗主氣流均勻混合等問題。通過上述措施,將NPLS技術(shù)拓展運用于高超聲速激波/邊界層干擾研究[61]。

圖7、圖8分別給出了入射斜激波與平板激波/邊界層干擾的紋影和納米粒子激光散射流場。由于紋影技術(shù)檢測的是光學(xué)擾動沿整個光路的積分值,加之模型邊界的衍射效應(yīng),很難觀察到邊界層內(nèi)部流場結(jié)構(gòu)。在相同試驗條件下,NPLS技術(shù)獲得了模型邊界層內(nèi)湍流的演變過程,這是紋影方法等常規(guī)手段所難以得到的(圖8中兩圖間隔時間為5μs)。

圖7 馬赫數(shù)M5條件下平板模型紋影流場圖Fig.7 Flowflield schlieren picture of the flat plate model under M5experimental condition

圖8 馬赫數(shù)M5條件下平板模型邊界層流場圖Fig.8 Flowfield picture of the flat plate model in boundary layer under M5experimental condition

2.4.2 輝光放電流場顯示技術(shù)

輝光顯示的原理為:低壓氣體中的自由電子和離子在外加電場的作用下加速,并與氣體分子發(fā)生碰撞產(chǎn)生二次電子和離子,使氣體分子激發(fā)到高能態(tài);不穩(wěn)定的受激分子通過躍遷發(fā)射出特定波長的輝光,輝光強弱是試驗流場氣體壓力的函數(shù),通過輝光強弱可以顯示試驗?zāi)P偷牧鲌鼋Y(jié)構(gòu)。通過輝光放電流場顯示技術(shù),解決了模擬高度從60km到100km的流場顯示問題,可提供高超聲速飛行器繞流波系、翼身激波干擾等直觀圖像。輝光放電裝置由交流高壓電源、放電電極及記錄裝置組成,通過其獲得的高空羽流流場圖像如圖9。

2.4.3 低雷諾數(shù)油流顯示技術(shù)

油流顯示是通過模型表面油膜在氣流作用下的變化,反映近壁極限流線的一種流動顯示技術(shù),其關(guān)鍵是油膜粘度及流動性與邊界層粘流狀態(tài)的匹配。針對低雷諾數(shù)、低剪切力的特殊流動條件,配制相應(yīng)的示蹤粒子混合油劑,其示蹤粒子濃度、油劑粘稠度等指標(biāo)都較好地適應(yīng)了稀薄流場邊界層流動及風(fēng)洞運行時的試驗條件,可以獲得滿意的壁面流動圖譜。圖10給出了平板模型側(cè)向噴流與低密度主流相互干擾的流動圖譜,及依據(jù)該圖譜分析得到的噴流干擾區(qū)流動多層次分離結(jié)構(gòu)。

圖9 高空羽流輝光放電流場圖Fig.9 Flow field picture of the plume at high altitude

圖10 橫噴干擾平板模型油流圖及流動分離結(jié)構(gòu)Fig.10 Oil flow picture of the flat plate model as well as flow separated structure in interaction of lateral jet flow

2.5 N-S/DSMC耦合算法

在過渡流區(qū),傳統(tǒng)的數(shù)值計算方法中,連續(xù)流求解N-S方程CFD技術(shù)和稀薄流的DSMC方法都會遇到各自的困難。在過渡流區(qū)由于稀薄度的增加,N-S方程往往會在局部或者全部區(qū)域失效,其結(jié)果可能會在流場和物面力/熱參數(shù)等計算方面產(chǎn)生比較大的誤差。DSMC方法受網(wǎng)格、時間步長以及仿真分子數(shù)等因素的限制,計算效率很低,這給數(shù)值求解該流區(qū)高超聲速繞流問題帶來了極大的困難。為拓展N-S方程和DSMC方法應(yīng)用范圍,在已有的CFD和DSMC方法和程序的基礎(chǔ)上,采用MPC(Modular Particle-Continuum)耦合技術(shù),對CFD和DSMC方法的計算程序進行耦合計算。

MPC的耦合方案為:在統(tǒng)一的網(wǎng)格系統(tǒng)下,把流場分為CFD計算區(qū)域和DSMC計算區(qū)域,CFD和DSMC在各自的區(qū)域進行計算。在分界面上,兩種方法進行信息交換,CFD和DSMC方法相互提供邊界條件。在兩個獨立的程序外部加入網(wǎng)格和信息交換模塊,就可以實現(xiàn)N-S/DSMC的耦合計算。

N-S/DSMC耦合算法需要解決CFD方法與DSMC方法自動分區(qū)問題、信息交換問題。

采用當(dāng)?shù)乜伺瓟?shù)Knl判斷連續(xù)流方程的失效問題:

式中,λ為當(dāng)?shù)亓鲌鰵怏w分子平均自由程;Q為當(dāng)?shù)亓鲌龊暧^參數(shù)。

當(dāng)克努森數(shù)Knl大于0.02,認(rèn)為連續(xù)流方程失效,該網(wǎng)格采用DSMC方法進行模擬。在計算過程中,根據(jù)CFD結(jié)果不斷調(diào)整DSMC方法計算區(qū)域,直到流場達到穩(wěn)定。

CFD的信息比較容易傳遞給DSMC方法的計算區(qū)域。對交界面上的每個網(wǎng)格上的當(dāng)?shù)睾暧^流動參數(shù),進行適時地更新,可以把CFD結(jié)果的變化適時地傳遞給DSMC求解區(qū)域。DSMC方法求解區(qū)域的信息向CFD求解區(qū)域的傳遞比較困難。這是因為在DSMC方法計算穩(wěn)定之前,統(tǒng)計結(jié)果要定時清零更新,統(tǒng)計樣本不可能太大,會使得宏觀參數(shù)有較大的統(tǒng)計波動。為了抑制DSMC方法的統(tǒng)計波動對CFD計算的影響,采用 “亞松弛”技術(shù)。在每個時間點j

式中,θ為網(wǎng)格中新參數(shù)的較小的權(quán)重;相應(yīng)地,(1- θ)是較大的權(quán)重。通過這種“亞松弛”技術(shù),能夠比較好地抑制DSMC方法的統(tǒng)計波動對CFD計算的影響。

在Case1(馬赫數(shù)為12.98、總壓1.309MPa、總溫595K、克努森數(shù)0.018 1)、Case2(馬赫數(shù)為12.90、總壓0.662MPa、總溫597K、克努森數(shù)0.025 9)和Case3(馬赫數(shù)為12.69、總壓0.280MPa、總溫595K、克努森數(shù)0.037 3)的試驗條件下,圖11給出了以返回艙低密度風(fēng)洞試驗條件為算例,對N-S/DSMC耦合算法進行了驗證。比較數(shù)值計算結(jié)果和風(fēng)洞試驗結(jié)果,以迎角10°為例,計算與試驗結(jié)果的偏差CA為3.2%,CN為3.1%。表明耦合算法的結(jié)果與試驗結(jié)果符合較好,驗證了所建立的方法和程序。上,一個網(wǎng)格中由DSMC方法得到的宏觀參數(shù)為Qj(包括ρ、u、v、w、T),傳遞給CFD的參數(shù)與上一步的參數(shù)進行加權(quán)平均:

圖11 不同狀態(tài)下氣動力計算與風(fēng)洞試驗比較Fig.11 Aerodynamic force comparisons between by computation and by wind tunnel under different conditions

3 下一步研究方向

樊菁(2013)根據(jù)2008年、2010年和2012年3屆國際稀薄氣體動力學(xué)會議的主題,指出了稀薄氣體動力學(xué)的研究方向[36]。本文在此基礎(chǔ)上,從學(xué)科發(fā)展、工程型號需求等方面,總結(jié)跨流域流動下一步研究。

3.1 地面試驗研究方面

發(fā)展電子束熒光流場診斷技術(shù),實現(xiàn)粒子數(shù)密度、轉(zhuǎn)動溫度、振動溫度、點速度和平均速度等基礎(chǔ)流動參數(shù)的定量測量,為跨流域流動機理研究及數(shù)值模擬驗證提供基礎(chǔ)手段;發(fā)展新型光纖天平測力技術(shù),進一步提高微量天平的抗溫度、電磁干擾能力及測試精度;發(fā)展跨流域的飛行器舵面鉸鏈力矩測量技術(shù)、熱噴流試驗技術(shù)、多體分離試驗技術(shù)、完善微壓測量、大面測熱試驗技術(shù),為臨近空間高超聲速飛行器跨流域氣動力、熱及控制、分離特性預(yù)測提供更精細(xì)、保真度更高的試驗研究手段。

針對空間環(huán)境的姿控發(fā)動機射流這一典型跨流域流動問題,發(fā)展羽流撞擊、污染的試驗?zāi)M環(huán)境及相關(guān)試驗技術(shù),開展羽流污染物輸運/沉積過程及其污染效應(yīng)、多噴管羽流干擾效應(yīng)研究,建立羽流干擾效應(yīng)及空間環(huán)境污染數(shù)據(jù)庫。

3.2 數(shù)值模擬研究方面

需進一步發(fā)展完善氣-面相互作用模型、考慮量子效應(yīng)的振動激發(fā)與非平衡化學(xué)反應(yīng)模型、氣體電離與輻射模型,并深化上述復(fù)雜動力學(xué)、熱力學(xué)及熱化學(xué)作用機理及對飛行器氣動力/熱特性的影響研究。

N-S/DSMC耦合計算方面:發(fā)展完善化學(xué)非平衡流動的計算方法,解決流動與化學(xué)反應(yīng)兩種時間尺度引起的剛性問題、微量組元在兩種算法交界面信息傳遞的波動問題,建立N-S/DSMC耦合算法的復(fù)雜飛行器跨流域氣動特性通用分析軟件系統(tǒng)。

基于Boltzmann模型方程的統(tǒng)一算法方面:發(fā)展適于混合氣體和化學(xué)反應(yīng)的Boltzmann模型方程,進一步優(yōu)化大規(guī)模并行算法,提高計算效率,使之穩(wěn)定運行于數(shù)萬CPU的超大規(guī)模并行計算機系統(tǒng)。

4 結(jié) 論

本文回顧了近年來在跨流域氣動力/熱預(yù)測技術(shù)方面的研究進展。建立了粘性干擾參數(shù)與克努森數(shù)搭接的跨流域模擬準(zhǔn)則,實現(xiàn)了從連續(xù)流到稀薄流的統(tǒng)一模擬;發(fā)展了微量天平優(yōu)化設(shè)計方法和測力技術(shù),解決了六分量微小載荷測量的天平布局、靈敏度、分量間干擾等難題;發(fā)展了紅外大面積中低量值熱流測量技術(shù),采用主動適應(yīng)、光路反射等方法,解決了大極角發(fā)射率修正與駐點區(qū)域熱流準(zhǔn)確測量問題;針對跨流域流動特性變化大、單一手段難以覆蓋的特點,發(fā)展了多手段結(jié)合的流動顯示與測量技術(shù);建立了N-S/DSMC自適應(yīng)緊耦合的數(shù)值模擬方法,實現(xiàn)了連續(xù)流和稀薄流區(qū)域的自動分區(qū),發(fā)展了可有效抑制DSMC統(tǒng)計波動影響的亞松弛技術(shù)。地面試驗具備了模擬高度從30多km到100km、馬赫數(shù)從5到24的跨流域氣動力、氣動熱與流場顯示能力。通過風(fēng)洞試驗和數(shù)值模擬相結(jié)合,為新型高超聲速飛行器跨流域氣動特性模擬和預(yù)測提供了有效的研究手段。

當(dāng)前,跨流域氣動研究尚未成熟,亟待發(fā)展更為精細(xì)、保真的試驗和數(shù)值模擬方法,在包含各種復(fù)雜效應(yīng)的流動機理方面也有待開展更加深入細(xì)致的研究,本文對下一步的研究方向也進行了初步探討,其中部分工作已經(jīng)開展,在不久的將來可以取得實質(zhì)性突破。

[1] Vidal R J,Bartz J A.Experimental study of low-density effect in hypersonic wedge flow[R],AD616898,1965.

[2] 楊彥廣,唐志共.低密度風(fēng)洞飛船測力技術(shù)研究[C]//第七屆高超聲速氣動力(熱)學(xué)術(shù)交流會議論文集,1993.

[3] 楊彥廣,唐志共,許曉斌.高超聲速低密度風(fēng)洞氣動力測量技術(shù)研究[C]//第九屆高超聲速氣動力(熱)學(xué)術(shù)交流會議論文集,1997.

[4] 楊彥廣,唐志共,許曉斌.微量天平測力技術(shù)研究與應(yīng)用[C]//第7屆全國天平學(xué)術(shù)會議,1998.

[5] 戴金雯,楊彥廣,李緒國.尾支撐外式三分力微量天平研制[C]//第八屆全國天平學(xué)術(shù)會議,1999.

[6] Tang Zhigong,Yang Yanguang.Development and application of microbalance in hypersonic low densty wind tunnel[R].NASA/CP-1999-209101/PT1,March 1999:669-76

[7] Yang Yanguang,Tang Zhigong,Xu Xiaobin.Aerodynamic microbalance calibration system[C]//2th International Symposium on Strain-Gauge Balance.Bedford,England,UK.May 1999.

[8]Yang Yanguang,Dai Jinwen,Tang Zhigong.Two types of rear support aerodynamic microbalances[C]//3th International Symposium on Strain-Gauge Balance.Frankfurt Germany,April 2001.

[9] 李明,楊彥廣,等.利用紅外熱圖開展通用航空飛行器氣動熱特性試驗[J].紅外與激光工程,2013,42(2):285-289.

[10]李明,祝智偉,李志輝.紅外熱圖在高超聲速低密度風(fēng)洞測熱試驗中的應(yīng)用概述[J].實驗流體力學(xué),2013,27(3):108-112.

[11]曾學(xué)軍,李明,劉太奎,等.用紅外熱圖技術(shù)進行升力體模型氣動熱特性試驗研究[J].空氣動力學(xué)學(xué)報,2010,22(4):494-498.

[12]李明,廖俊必,祝智偉,等.紅外熱圖在羽流撞擊平板熱效應(yīng)試驗中的應(yīng)用[J].紅外與激光工程,2010,39(5):796-800.

[13]Seiff A.Secondary flow-fields embedded in hypersonic shock layers[R].NASATN D-1304,1962.

[14]Griffith B J,Boulan D E.Postflight[AS-202]Appol command module aerodynamic simulation test[R].AEDC-TR-67-238.1968.

[15]Arvel E.Gentry,Donglas N.Smyth,Wayne R.Oliver.The mark IV supersonic-h(huán)ypersonic arbitrarybody program[R].AD77844,1973.

[16]Bird G A.Application of the direct simulation Monte Carlo method to the full shuttle geometry[R].AIAA 90-1692,1990.

[17]Dietrich S.Efficient computation of particle movement in 3-D DSMC calculation on structured body-fitted grids[C]//17th International Symposium on Rarefied Gas Dynamics,Achen,1990:745-752.

[18]Bird G A.Definition of mean free path for real gases[J].Phys Fluids,1983,26:3222-3223.

[19]Musanov S V,Nikiforov A P,Omelik A L,et al.Experimental determination of momentum tran sfer coefficients in hypersonic free molecular flow and distribution function recovery of reflected molecules[A].Proceedings of the XIII International Symposium Rarefied Gas Dynamics[M],Vol.1,Plenum,New York,1985:669-676.

[20]Antonov S G,Ivanov M S,Kashkovsky A V,et al.Influence of atmospheric rarefied on aerodynamic characteristics of flying vehicles[A].In Proc.Int.Symp.Rarefied Gas Dynamics[M],17th,Aachen,ed.AE Beylich,1991:512-530.

[21]Ivanov M S,Antonov S G,Gimelshein S F,et al.Computational tools for rarefied aerodynamics[C]//The Eighteenth International Symposium on Rarefied Gas Dynamics,pp115-126,1992.

[22]Boyd I D,Penko P F,Meissner D L,et al.Experimental and numerical investigation of low-density nozzle and plume flows of nitrogen[J].AIAA J.,1992,30(10):2453-2461.

[23]Bird G A.Molecular gas dynamics and the direct simulation of gas flows[M].Clarendon Press,Oxford,(1994.).

[24]Rault D F G.Aerodynamics of the shuttle orbiter at high altitudes[J].J.Spacecraft and Rockets,1994,31(6):944-952.

[25]Pallegoix J F.Computation of tansitional flows around three-dimensional reentry bodies[A].In Rarefied Gas Dynamics:Theory and Simulations[M],ed.BD Shizgal,DP Weaver,1994,Prog.Astronaut.Aeronaut.160:161-170.

[26]Kannenberg K.C.,Boyd I.D.,Dietrich S.Development of an object-oriented parallel DSMC code for plume impingement studies[R].AIAA 95-2052,1995

[27]Wilmoth R G,Lebeau G J,Carlson A B.DSMC grid methodologies for computing low-density hypersonic flows about reusable launch vehicles[R].AIAA 96-1812,1996.

[28]Nance R P,Hash D B,Hassan H A.Role of boundary condition in monte carlo simulation of mems devices[R].AIAA 97-0375,1997.

[29]Cai Chunpei,Boid I D,F(xiàn)an Jing.Direct simulation methods for low-speed micro-channel flows[R].AIAA 99-3801,1999

[30]Moss J N,Prince J M.DSMC simulation of OREX entry[C]//In Proc.Int.Symp.Rarefied Gas Dynamics,20th,Beijing,1997.

[31]Blanchard R C,Wilmoth R G,LeBeau G J.Rarefied-flow transition regime on orbiter aerodynamic acceleration flight measurements[J].J.Spacecraft and Rockets,1997,34(1):8-15.

[32]Ivanov M S,Markelov G N,Gimelshein S F.High-altitude capsule aerodynamics with real gas effects[J].Journal of Spacecraft and Rockets,1998,35(1):16-22.

[33]Quade M,Hufnagel K.The development of modern manual calibration and measuring system for internal balances[R].NASA/CP-1999-209101/PT2,March 1999:7-18.

[34]Gatsonis N A,NaN-Son R A,LeBeau G J.Simulation N-S of cold-gas nozzle and plume flows and flight data comparison[J].J.Spacecraft and Rockets,2000,37(1):39-48.

[35]Fang Yichuan,liou W W.Microfluid flow computatio N-S using aparallel DSMC code[R].AIAA 2002-1057,2002.

[36]樊菁.稀薄氣體動力學(xué):進展與應(yīng)用[M].力學(xué)進展,2013,43(2):185-201.

[37]沈青.航天飛行中的稀薄氣體動力學(xué)[J].流體力學(xué)實驗與測量,1997,11(4):1-7.

[38]王保國,李耀華,錢耕.四種飛行器繞流的三維DSMC計算與傳熱分析[J].航空動力學(xué)報,2011,26(1):1-20.

[39]王智慧,鮑麟,童秉綱.高超聲速尖頭體駐點熱流從連續(xù)態(tài)過渡到稀薄態(tài)的變化特征和橋函數(shù)研究[J].中國科學(xué):物理學(xué) 力學(xué)天文學(xué),2009,39(8):1134-1140.

[40]梁杰.飛船返回艙稀薄高超音速繞流的DSMC方法研究[D].西安:西北工業(yè)大學(xué)研究生院,2003.

[41]梁杰,閻超,杜波強.基于兩級直角網(wǎng)格結(jié)構(gòu)的三維DSMC算法研究[J].空氣動力學(xué)學(xué)報,2010,28(4):466-471.

[42]梁杰,閻超,楊彥廣,杜波強.過渡區(qū)側(cè)向噴流干擾的并行DSMC數(shù)值模擬研究[J].宇航學(xué)報,2011,32(5):1012-1018.

[43]梁杰,閻超,李志輝,等.稀薄過渡流區(qū)橫向噴流干擾效應(yīng)數(shù)值模擬研究[J].空氣動力學(xué)學(xué)報,2013,31(1):27-33.

[44]梁杰,李志輝,杜波強,等.探月返回器稀薄氣體熱化學(xué)非平衡特性數(shù)值模擬[J].載人航天,2015,21(3):295-302.

[45]徐昆,李啟兵,黎作武.離散空間直接建模的計算流體力學(xué)方法[J].中國科學(xué):物理學(xué)力學(xué)天文學(xué),2014,44(5):519-530.

[46]李志輝.從稀薄流到連續(xù)流的氣體運動論統(tǒng)一數(shù)值算法研究[D].綿陽:中國空氣動力研究與發(fā)展中心,2001.

[47]李志輝,張涵信.跨流域三維復(fù)雜繞流問題的氣體運動論并行計算[J].空氣動力學(xué)學(xué)報,2010,28(1):7-16.

[48]李志輝,張涵信,符松.基于Boltzmann模型方程的氣體運動論HPF并行算法[J].計算物理,2003,20(1):1-8.

[49]李志輝,張涵信,基于Boltzmann模型方程各流域飛行器繞流問題大規(guī)模并行計算[C]//高性能計算戰(zhàn)略研討會報告文集,pp.114-130,中國科學(xué)院學(xué)部“高性能計算戰(zhàn)略研究”咨詢組,2006.

[50]Hash D B,Hassan H A.An decoupled DSMC/Navier-Stokes analysis of a transitional flow experiment[R].AIAA-96-0353,1996.

[51]Wadsworth D C,Erwin D A.Two-dimensional hybrid continu-um/particle simulation approach for rarefied hypersonic flows[R].AIAA Paper 1992-2975.

[52]Hash D B,Hassan H A.Assessment of schemes for coupling monte carlo and Navier-Stokes solution methods[J].Journal of Thermophysics and Heat Transfer,1996,10(2):242-249.

[53]Oh C K,Oran E S.A new hybrid algorithm:Navier-Stokes as a DSMC filter[R].AIAA Paper 1998-0849.

[54]Schwartzentruber T E,Scalabrin L C,Boyd I D.Hybrid particle-continuum simulations of non-equilibrium hypersonic blunt body flow fields[R].AIAA 2006-3602.

[55]Deschenes T R,Boyd I D.Extension of a modular particle-continuum method to vibrationally excited,hypersonic flows[J].AIAA JOURNAL,2011,49(9):1951-1959.

[56]李中華,李志輝,李海燕,等.過渡流區(qū)NS/DSMC耦合算法研究[J].空氣動力學(xué)學(xué)報.2013,31(3):282-287.

[57]Li Zhonghua,Li Zhihui,Li Haiyan,et al.Application of hybrid N-S/DSMC method in hypersonic transitional flow[C]//28th International Symposium on Rarefied Gas Dynamics.2012,pp.435-442.

[58]Li Zhonghua,Li Zhihui,Li Haiyan,et al.N-S/DSMC hybrid simulation of hypersonic flow over blunt body including wakes[C]//29th International Symposium on Rarefied Gas Dynamics.2014,pp.519-526.

[59]趙玉新,易仕和,何霖,等.激波與湍流相互作用的實驗研究[J].科學(xué)通訊,2007,52(2):140-143.

[60]Zhao Y X,Yi S H,Tian L F,et al.Supersonic flow imaging via nanoparticles[J].Sci.China Ser.E-Tech.Sci.,2009,52(12):3640-3648,doi:10.1007/s11431-009-0281-3.

[61]李明,易仕和,祝智偉.激光散射技術(shù)在高超聲速激波與邊界層干擾試驗中的應(yīng)用[J].紅外與激光工程,2013,42(S1):79-83.

Aerodynamic force/heating measurement on hypersonic vehicle across different flow regions

Yang Yanguang*,Li Ming,Li Zhonghua,Li Xuguo,Dai Jinwen
(Hypersonic Aerodynamic Institute of China Aerodynamics Research and Development Center,Mianyang Sichuan 621000,China)

The flow characteristics across different flow regions,engineering requires for aerodynamic force/heating predicting technique,and the research status in China and foreign countries are summarized.To solve the main problems of aerodynamic force/heating predicting techniques across different flow regions,some recent developments,such as the similarity criterion for experimental simulation based on viscous interaction and Knudsen number,framework optimization of microbalance as well as aerodynamic force measurement technique,the measurement of the middle and low heat transfer using infrared mapping,flow visualization and measurement,together with N-S/DSMC hybrid numerical method,have been presented,and further investigation is suggested.

across different flow region;aerodynamic force;aerodynamic heating;flow visualization;hybrid numerical method

V211.3;V211.7

Adoi:10.7638/kqdlxxb-2015.0149

0258-1825(2016)01-0005-09

2015-08-10;

2015-09-15

楊彥廣*(1969-),研究員,空氣動力學(xué)專業(yè)。E-mail:yyyyygggg@sina.com

楊彥廣,李明,李中華,等.高超聲速飛行器跨流域氣動力/熱預(yù)測技術(shù)研究[J].空氣動力學(xué)學(xué)報,2016,34(1):5-13.

10.7638/kqdlxxb-2015. Yang Y G,Li M,Li Z H,et al.Aerodynamic force/heating measurement on hypersonic vehicle across different flow regions[J].Acta Aerodynamica Sinica,2016,34(1):5-15.

猜你喜歡
測量方法
把握四個“三” 測量變簡單
學(xué)習(xí)方法
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
測量的樂趣
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
測量
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 中文字幕永久在线看| 91美女视频在线| 国产高清不卡视频| 亚洲av无码牛牛影视在线二区| 日韩小视频在线播放| 曰AV在线无码| 日韩精品无码一级毛片免费| 亚洲精品少妇熟女| 自慰高潮喷白浆在线观看| 亚洲精品无码久久久久苍井空| 国产日韩欧美视频| 中文无码毛片又爽又刺激| 久久人体视频| 在线欧美一区| 91麻豆久久久| 久爱午夜精品免费视频| 欧美日韩一区二区三区四区在线观看 | 国产亚洲精久久久久久久91| 色偷偷一区| 国产精品专区第1页| 欧美成人精品在线| 青青草欧美| 男女男精品视频| 亚洲精品欧美日本中文字幕| 亚洲swag精品自拍一区| 71pao成人国产永久免费视频| 少妇被粗大的猛烈进出免费视频| 丰满人妻中出白浆| 国产精品白浆无码流出在线看| 中文字幕一区二区人妻电影| 欧美午夜在线观看| 亚洲精品天堂自在久久77| 免费在线观看av| 亚洲国产成熟视频在线多多| 国产成人一区免费观看 | 熟妇人妻无乱码中文字幕真矢织江| 日韩区欧美国产区在线观看| 高清久久精品亚洲日韩Av| 日韩成人在线网站| 国产在线观看人成激情视频| 色窝窝免费一区二区三区| 91亚洲免费视频| 波多野结衣亚洲一区| 狠狠躁天天躁夜夜躁婷婷| 国产小视频在线高清播放 | 亚洲天堂日韩av电影| a亚洲视频| 欧美v在线| 91一级片| 直接黄91麻豆网站| 国产凹凸一区在线观看视频| 国产真实乱子伦精品视手机观看| av色爱 天堂网| 一级毛片a女人刺激视频免费| 亚洲精品无码久久毛片波多野吉| 亚洲国产精品不卡在线| 在线观看免费国产| 尤物亚洲最大AV无码网站| 亚洲色图欧美在线| 91年精品国产福利线观看久久| 国产精品女主播| 久久精品这里只有精99品| 91丝袜乱伦| 一本一本大道香蕉久在线播放| 午夜日b视频| 亚洲av无码成人专区| 91蝌蚪视频在线观看| 欧美一区二区精品久久久| 亚洲毛片网站| 久久人人爽人人爽人人片aV东京热| 伊人久综合| 色婷婷综合激情视频免费看 | 国产99欧美精品久久精品久久| 亚洲欧美成人综合| 在线国产毛片手机小视频| 丝袜高跟美脚国产1区| 成人小视频在线观看免费| 国产日韩欧美精品区性色| 亚洲看片网| 无码福利视频| 国产日韩欧美精品区性色| 成人一级免费视频|