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

高超聲速氣流中頭錐逆噴防熱流熱耦合分析 *

2022-01-26 12:56:36劉洪鵬劉偉強
國防科技大學(xué)學(xué)報 2022年1期
關(guān)鍵詞:結(jié)構(gòu)

尹 亮,劉洪鵬,劉偉強

(1. 湖南文理學(xué)院 機械工程學(xué)院, 湖南 常德 415000; 2. 國防科技大學(xué) 空天科學(xué)學(xué)院, 湖南 長沙 410073)

對于高超聲速飛行器,特別是其頭錐和前緣等尖銳部件,由于高速空氣滯止和摩擦作用,其表面面臨嚴(yán)重的氣動加熱環(huán)境,這給飛行器的安全設(shè)計帶來巨大挑戰(zhàn)[1-4]。對于長時間飛行,傳統(tǒng)的燒蝕結(jié)構(gòu)不再適于承擔(dān)防熱任務(wù):一方面,燒蝕由于質(zhì)量的減損,飛行器表面結(jié)構(gòu)的改變帶來了氣動環(huán)境的改變,表面粗糙度的增加還增大了摩擦阻力;另一方面,燒蝕結(jié)構(gòu)不滿足可重復(fù)性使用,增加了飛行器的使用成本。而作為一種主動冷卻技術(shù)[1,5-6],逆向噴流技術(shù)在減阻和防熱方面具有顯著效果[7-8]。

逆向噴流技術(shù)的減阻和防熱前景,最初由Lopatoff[9]和 Warren[10]提出,他們采用先進的試驗和數(shù)值計算工具,進行了相關(guān)機理性研究。Hayashi等[11-12]采用試驗和數(shù)值方法對單孔逆向噴流的減阻和防熱特性進行了研究,其中自由來流馬赫數(shù)為3.98,數(shù)值計算采用k-ω湍流模型,其研究結(jié)果表明逆向噴流總壓和自由來流總壓的比值對逆向噴流的流動模式起著關(guān)鍵作用。Li等[13]研究了乘波體前緣多孔逆噴結(jié)構(gòu)的減阻和防熱特性,采用k-ω湍流模型進行數(shù)值計算,結(jié)果表明多孔噴流相對于縫隙噴流的防熱效果更加優(yōu)越;當(dāng)噴流孔徑一定時,存在最佳噴孔間隔,使得減阻和防熱效果達到最優(yōu)。Huang等[14-15]還研究了逆向噴流和其他防熱措施的組合熱防護性能,包括激波針和凹腔等結(jié)構(gòu),其研究結(jié)果表明組合結(jié)果的防熱效果更佳,并且由逆噴所帶來的流場振蕩和其他結(jié)構(gòu)之間的非穩(wěn)態(tài)模型之間存在相互耦合效應(yīng)。

目前,大部分的數(shù)值計算研究是基于固體結(jié)構(gòu)表面的均勻冷壁溫假設(shè)而開展的,忽略了流體和固體結(jié)構(gòu)之間的熱交互作用,這主要會帶來氣動熱環(huán)境的過考核問題。此外,這種與物理事實不相符合的假設(shè),不利于對防熱機理的進一步認(rèn)識。事實上,固體結(jié)構(gòu)氣動加熱的降低是和壁面溫度相關(guān)的。

為了對逆向噴流中的流動和傳熱特性進行更好的理解,本文采用共軛傳熱方法,對6馬赫飛行條件下的頭錐逆噴結(jié)構(gòu)中的流熱交互過程進行數(shù)值分析。首先介紹所采用的數(shù)值計算模型和流熱耦合方法,并分別與相關(guān)試驗結(jié)果對比,驗證其準(zhǔn)確性。然后對不同噴流壓力情況下的減阻和防熱效果進行數(shù)值計算,并進一步討論不同攻角對防熱效果的影響。為了研究固體結(jié)構(gòu)的熱響應(yīng)情況,分別采用IN718和C-103高溫難熔合金材料,對其傳熱特性進行對比分析。

1 物理模型和數(shù)值方法

1.1 物理模型

物理模型如圖1所示,由于對稱性,圖中僅給出了縱向截面的1/2結(jié)構(gòu)。其中,頭部半徑為6 mm,半楔角φ為10°,頭錐長度L為100 mm,逆噴噴口直徑d為1.2 mm。飛行高度假定為25 km,飛行馬赫數(shù)為6。逆噴相比于自由來流的總壓,即總壓比(Pressure Ratio, PR)為0.092。逆噴氣流溫度為300 K,逆噴馬赫數(shù)為1。頭錐材料采用高溫難熔合金IN718和C-103,頭錐結(jié)構(gòu)厚度設(shè)定為3 mm。

圖1 物理模型簡圖Fig.1 Diagram of physical model

1.2 數(shù)值方法

對于氣動環(huán)境的預(yù)測,采用三維Navier-Stokes (N-S)方程作為控制方程。為了提高計算準(zhǔn)確性,對流項采用AUSM-DV格式和MUSCL方法近似處理,黏性項采用中心差分格式。由于逆噴和自由來流之間,以及邊界層附近強烈的流動剪切作用,數(shù)值計算采用k-ω(SST)模型。固體域采用熱傳導(dǎo)控制方程,并在固體表面增加向周圍環(huán)境的熱輻射項,金屬固體表面輻射系數(shù)假定為0.8。

流熱耦合采用如下方法:

步驟1:初始時刻t0,為固體結(jié)構(gòu)指定均溫溫度300 K,并將流固交界面的溫度傳輸至流體域,并作為流動計算的邊界條件。

步驟2:采用傳輸至流體域的溫度邊界條件,對流場進行穩(wěn)態(tài)計算,得到流固交界面上的氣動熱環(huán)境。

步驟3:將流固交界面上的氣動熱分布傳輸至固體域,并基于該氣動熱邊界條件計算t0至t0+Δtc時刻固體域的瞬態(tài)傳熱。

步驟4:將計算所得的t0+Δtc時刻流固交界面溫度傳輸至流體域,并作為新的邊界條件進行流場計算。

步驟5:重復(fù)以上步驟2~4,直至完成所有時間步計算。

2 數(shù)值方法驗證

2.1 湍流模型驗證

流體域的計算采用FLUENT商業(yè)軟件,結(jié)構(gòu)傳熱采用Abaqus軟件,通過編程實現(xiàn)數(shù)據(jù)傳輸。為了驗證湍流模型的準(zhǔn)確性,將計算結(jié)果與文獻[11]中的實驗結(jié)果進行了對比。相關(guān)實驗設(shè)置參數(shù)如表1所示。圖2給出了數(shù)值計算和實驗測量的流場密度分布對比,可以看出流場計算中,馬赫盤、回流區(qū)和再壓縮激波均得到了較好的捕捉,并且和實驗結(jié)果較為一致。

表1 實驗參數(shù)

圖2 密度分布的計算和實驗結(jié)果對比Fig.2 Comparison of density distribution between numerical and experimental results

圖3 St分布對比Fig.3 Comparison of St distribution

圖3為固體表面的斯坦頓數(shù)(St)分布,其中θ為球頭外表面和自由來流方向間的夾角。斯坦頓數(shù)是一個無量綱數(shù),用于測量傳遞到流體中的熱量與流體的熱容量之比。由圖3可以看出,數(shù)值計算和實驗結(jié)果的分布趨勢較為一致。除了左邊第一個測點,大部分的結(jié)果偏差在15%以內(nèi)。該偏差可能是由數(shù)值計算模型的假設(shè)以及實驗測量誤差所引起的。總的來說,數(shù)值計算和實驗測量的流場和傳熱分布結(jié)果吻合較好,證明了逆向噴流所采用的數(shù)值模型的適用性。

2.2 流熱耦合方法驗證

為了驗證流熱耦合方法的準(zhǔn)確性,將計算結(jié)果與實驗測量結(jié)果進行對比。其中,實驗結(jié)果采用1987年于NASA蘭利中心開展的高速氣流中的圓管氣動加熱實驗,來流馬赫數(shù)為6.47,管材為不銹鋼,其他實驗參數(shù)參照文獻[16]。

為了節(jié)省計算資源,僅對氣動加熱嚴(yán)重的圓管前半部分進行數(shù)值計算。為了保證氣動熱計算的準(zhǔn)確性,流體域近壁面第一層網(wǎng)格厚度取1×10-5m,固體域第一層網(wǎng)格厚度取1×10-4m,流固交界面附近的固體域和流體域均進行了網(wǎng)格加密,如圖4所示。

圖4 流體域和固體域網(wǎng)格分布Fig.4 Grid distribution of fluid domain and solid domain

圖5為圓柱表面熱流分布隨時間的變化,實驗數(shù)據(jù)為冷壁溫時的測量結(jié)果[16]。其中q為任意時刻圓柱駐點熱流,q0為0 s時圓柱駐點熱流。由圖5可以看出,0 s時數(shù)值計算和實驗結(jié)果較為一致。數(shù)值結(jié)果表明,在加熱時間1 s內(nèi),駐點位置的氣動熱流值迅速下降,而后部區(qū)域的熱流降低幅度較小。這表明,采用流熱耦合方法對準(zhǔn)確計算氣動熱分布是非常重要的,而采用均勻的冷壁溫邊界條件假設(shè)將使氣動熱計算過考核。

圖5 熱流分布對比Fig.5 Comparison of heat flux distribution

圖6 表面溫度分布對比Fig.6 Comparison of surface temperature distribution

圖6為不同時間圓管表面的溫度分布,其中實驗值為約5 s時的數(shù)據(jù)。由圖6可以看出,大部分的實驗數(shù)據(jù)位于計算所得的4 s和5 s曲線之間,其中5 s時的駐點溫度比實驗值高約3%,證明了數(shù)值計算和實驗結(jié)果吻合良好。從圖6還可以看出,隨著氣動加熱時長增加,駐點溫度上升速度遠(yuǎn)遠(yuǎn)高于圓管尾部溫度上升速度,即駐點和尾部溫差隨之增大,這意味圓管要承受的熱應(yīng)力也將隨之增加。

此外,駐點溫度隨時間的變化如圖7所示。在氣動加熱的前1 s內(nèi)駐點溫度快速上升,1 s后溫度上升速度較為趨緩。當(dāng)駐點固體壁面溫度升高時,空氣滯止溫度和固體壁面溫度間的溫差隨之減小,對固體域的氣動加熱熱流也隨之減小,這可與圖5中的熱流分布加以相互印證。 5 s時固體結(jié)構(gòu)中的溫度云圖如圖8所示。駐點區(qū)域固體壁面由于氣動加熱最為嚴(yán)重,其溫度最高。溫度云圖的溫度梯度分布清晰表明固體內(nèi)由圓管外壁面向內(nèi)壁面的導(dǎo)熱以及沿著圓管周向上由駐點向尾部的導(dǎo)熱。

圖7 駐點溫度隨時間的變化Fig.7 Stagnation temperature changes with time

圖8 固體溫度分布(t=5 s)Fig.8 Solid domain temperature distribution(t=5 s)

3 計算結(jié)果

3.1 流動與傳熱性能

采用IN718作為結(jié)構(gòu)材料,對逆噴總壓比PR=0.092的工況進行數(shù)值計算,溫度、壓力和流體速度計算結(jié)果如圖9所示。由圖9可以看出,氣體自噴口逆向噴出后迅速膨脹,形成馬赫盤,該馬赫盤主要是激波后壓力和逆噴壓力的平衡結(jié)果。逆噴和自由來流相互作用結(jié)果對減低氣動加熱起到了直接作用:一是逆向噴流將激波推離固體壁面,并且將高溫氣體封裝于遠(yuǎn)離固體壁面處;二是環(huán)噴口附近形成了低溫低壓回流區(qū),這對頭部區(qū)域冷卻具有積極作用。此外還可以看到,在氣體再附層的下游形成了再壓縮激波,此處高溫高壓的氣流流向近壁面區(qū)域,會造成此處氣動熱的升高。

(a) 速度場分布(a) Velocity distributions

圖9 速度場、溫度場及壓力場分布(PR=0.092)Fig.9 Velocity,temperature and pressure distributions (PR=0.092)

圖10 不同時刻熱流分布Fig.10 Heat flux distribution at different times

圖10為熱流分布隨時間變化的計算結(jié)果。可以看出,在2 s時,固體壁面的最大熱流由1.29 MW/m2降低至0.78 MW/m2,降低幅度約為40%;60 s時,最大熱流值降低至0.29 MW/m2,降低幅度約為78%。此外,60 s時尾端的熱流降低幅度約為68%。顯然,采用均勻的冷壁溫邊界條件假設(shè),將帶來氣動熱環(huán)境的嚴(yán)重高估,這對結(jié)構(gòu)設(shè)計提出了更為苛刻的要求。固體表面的溫度變化如圖11所示。可以看出,固體壁面的最高溫度位于氣流再壓縮區(qū)域,同時可以看到固體壁面受氣動加熱后的升溫過程。60 s時固體表面的溫度分布云圖如圖12所示。可以看出,60 s時固體表面最高溫度位于氣流的再壓縮區(qū)。

圖11 不同時刻結(jié)構(gòu)表面溫度分布Fig.11 Surface temperature distribution at different times

圖12 60 s時固體表面溫度分布云圖Fig.12 Solid surface temperature distribution at t=60 s

3.2 總壓比的影響

圖13比較了總壓比PR=0.092 (圖中上半部分)和PR=0.184(圖中下半部分)的流場計算結(jié)果。可以看出,當(dāng)總壓比增大到0.184時,弓形激波和高溫氣體均被進一步推離頭錐的頭部區(qū)域。從圖中還可以看出,總壓比增大時,逆向噴流的出流膨脹角有所增大,環(huán)出流的低溫回流區(qū)域也有所增大,這將有利于頭部區(qū)域的進一步冷卻。此外,再壓縮區(qū)域也受到總壓比增大的影響,總壓比增大時,再壓縮激波同樣被進一步推離固體壁面,這將使得肩部附近的最大熱流值有所降低。

圖13 不同總壓比時溫度場和速度場分布Fig.13 Temperature and velocity distributions at different total pressure ratio

圖14比較了沒有逆噴、逆噴總壓比PR=0.092和PR=0.184時壁面熱流分布的對比。可以看出,在數(shù)值模擬時間段的起始和最終時刻,采用逆噴結(jié)構(gòu)后壁面的最高熱流均有了大幅度的下降;當(dāng)總壓比PR由0.092增大到0.184時,壁面最高熱流有了更進一步的降低,這主要是由于流場結(jié)構(gòu)改變的結(jié)果。

圖14 不同總壓比時熱流分布Fig.14 Heat flux distribution at different total pressure ratio

此外,圖14中可以看到,逆噴總壓比變化時,環(huán)出流的氣流回流區(qū)和再壓縮區(qū)域是受影響最顯著的區(qū)域。由于逆向噴流冷空氣的注入,固體表面得到冷卻,其溫度也隨之得到了大幅度的降低,并且固體結(jié)構(gòu)的溫度分布整體上更加均勻:其中,PR=0.092時固體表面溫差降低了55.5%,PR=0.184時固體表面溫差降低了83.1%,如圖15所示。

圖15 不同總壓比時固體壁面溫度分布Fig.15 Solid surface temperature distribution at different total pressure ratio

3.3 攻角的影響

本小節(jié)通過數(shù)值計算比較不同攻角對氣動加熱環(huán)境的影響,其中設(shè)定PR=0.138,攻角分別取0°、2°和5°。其中5°攻角時,對稱面上的流場密度和溫度分布云圖如圖16所示。可以看出,迎風(fēng)面上的再壓縮激波明顯薄于背風(fēng)面,這意味著在肩部區(qū)域,迎風(fēng)面處的高溫氣流更加流向固體表面,這也必然帶來此處氣動加熱環(huán)境的惡化。

圖16 5°攻角時流場云圖Fig.16 Flow field contour when the attack angle is 5°

圖17為0 s和60 s時迎風(fēng)面和背風(fēng)面中心線上的熱流分布。由圖17(a)可以看出,隨著攻角的增加,迎風(fēng)面肩部區(qū)域的熱流值增大幅度較大,當(dāng)攻角增大到5°時,最大熱流值相比于無攻角時增大超過60%,而背風(fēng)面的熱流值相應(yīng)得到大幅度降低。這主要是攻角的增大使得迎風(fēng)面和背風(fēng)面的激波厚度均得以改變所致,如圖16所示。相應(yīng)地,在頭錐尾部區(qū)域的迎風(fēng)面和背風(fēng)面上,氣動熱流也分別有不同程度的增大和減小。60 s時固體表面的熱流分布如圖17(b)所示,在球頭的肩部區(qū)域,迎風(fēng)面中心線上的熱流隨攻角的增大而增大,背風(fēng)面中心線上的熱流隨攻角的增大而減小。當(dāng)攻角為2°時,在頭錐的尾部區(qū)域,60 s時迎風(fēng)面中心線上的熱流分布相較于0°時變化不大,這主要是由于攻角增大幅度較小時,尾部區(qū)域激波厚度變化相對不顯著,并且受流體和固體之間的耦合傳熱以及固體內(nèi)部的導(dǎo)熱所影響。此外還可以看出,60 s時噴口附近的熱流為負(fù)值,并且在背風(fēng)面處負(fù)值區(qū)域相對更大,這是由于攻角導(dǎo)致逆向出流更容易流向流場壓力較低的背風(fēng)面,從而使背風(fēng)面得到冷氣流更好的冷卻效果。

(a) 0 s時的熱流分布(a) Heat flux distribution at 0 s

(b) 60 s時的熱流分布(b) Heat flux distribution at 60 s圖17 不同攻角下的熱流分布Fig.17 Heat flux distribution at different attack angle

圖18為60 s時固體表面的溫度分布,可以看出攻角的增大使得固體表面的溫差得以大幅度增大,這對材料的耐高溫性能和力學(xué)性能提出了更嚴(yán)苛的要求。

圖18 不同攻角下溫度分布Fig.18 Temperature distribution at different attack angle

3.4 材料的影響

圖19 不同材料熱流分布Fig.19 Heat flux distribution of different materials

如前所述,逆向噴流中的傳熱實際上是流體和固體結(jié)構(gòu)間的動態(tài)交互過程,所以本小節(jié)對材料的影響進行對比研究,其中頭錐飛行攻角為0°,PR=0.184。圖19為10 s和60 s時刻固體壁面的熱流分布,固體材料分別為IN7118和C-103。可以看出,不同材料下頭部區(qū)域的熱流分布比較一致,尾部區(qū)域差異相對明顯,這是因為:頭部區(qū)域受逆向出流的影響較為顯著,回流區(qū)中的對流冷卻對頭部區(qū)域的防熱占主導(dǎo)因素;而尾部區(qū)域則以固體結(jié)構(gòu)內(nèi)部的熱傳導(dǎo)為主,不同材料的熱物性導(dǎo)致固體內(nèi)部的傳熱速度不同,從而進一步影響了流體和固體表面的耦合傳熱情況。圖20為固體結(jié)構(gòu)內(nèi)表面和外表面的溫度分布,可以看出選用C-103為固體結(jié)構(gòu)材料時結(jié)構(gòu)溫度分布更為均勻,這主要是因為C-103有著更大的導(dǎo)熱系數(shù),從而使加熱過程中固體結(jié)構(gòu)內(nèi)部的溫差更小。圖21為固體結(jié)構(gòu)最高溫度隨時間的變化,可以看出不同材料下固體最高溫度的變化較為一致,這與圖19中不同時刻下不同材料表面的最高熱流較為一致是相吻合的。

圖20 結(jié)構(gòu)內(nèi)外表面溫度分布(t=60 s)Fig.20 Inner surface and outer face temperature distribution (t=60 s)

圖21 不同材料結(jié)構(gòu)表面最高溫度隨時間的變化Fig.21 Maximum surface temperature for different materials with time

4 結(jié)論

采用流熱耦合方法,對6馬赫下的頭錐逆噴結(jié)構(gòu)的流動與傳熱特性進行了數(shù)值研究。基于對湍流模型和流熱耦合方法的驗證,針對不同逆噴總壓比、攻角和材料對頭錐防熱效果的影響進行了對比研究。主要結(jié)論如下:

1)對于逆向噴流結(jié)構(gòu),為了獲得準(zhǔn)確的氣動加熱環(huán)境,采用流熱耦合進行相關(guān)計算非常必要。與流熱耦合數(shù)值計算方法相比,基于均勻冷壁溫邊界條件假設(shè)的熱流計算值被嚴(yán)重高估。

2)冷氣流回流區(qū)和再壓縮激波對逆噴結(jié)構(gòu)的防熱起著關(guān)鍵作用。隨著總壓比的增大,冷氣流回流區(qū)域增大,再壓縮激波被進一步推離固體表面,致使固體結(jié)構(gòu)得到更大幅度的冷卻。采用逆向噴流結(jié)構(gòu)以后,固體結(jié)構(gòu)溫度下降顯著,結(jié)構(gòu)整體溫度分布更加均勻。

3)當(dāng)頭錐攻角大于0°時,迎風(fēng)面熱流得以增大,背風(fēng)面熱流減小;頭部環(huán)出流的冷氣流回流區(qū)和再壓縮激波的厚度受影響明顯。攻角增大時,固體結(jié)構(gòu)內(nèi)的溫差增大幅度顯著,這給固體材料的熱和力學(xué)性能提出了更嚴(yán)苛要求。

4)對于IN718和C-103兩種不同材料,其在頭部氣動加熱環(huán)境的差異可以忽略。整體上,采用C-103作為結(jié)構(gòu)材料時,結(jié)構(gòu)內(nèi)的溫度分布更為均勻,這主要得益于C-103較大的導(dǎo)熱系數(shù)。

猜你喜歡
結(jié)構(gòu)
DNA結(jié)構(gòu)的發(fā)現(xiàn)
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結(jié)構(gòu)的應(yīng)用
模具制造(2019年3期)2019-06-06 02:10:54
循環(huán)結(jié)構(gòu)謹(jǐn)防“死循環(huán)”
論《日出》的結(jié)構(gòu)
縱向結(jié)構(gòu)
縱向結(jié)構(gòu)
我國社會結(jié)構(gòu)的重建
人間(2015年21期)2015-03-11 15:23:21
創(chuàng)新治理結(jié)構(gòu)促進中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 日韩区欧美国产区在线观看| 免费国产好深啊好涨好硬视频| 亚洲最大福利视频网| 伊人久久精品无码麻豆精品| 91亚洲影院| 亚洲日韩精品无码专区97| 亚洲伦理一区二区| 婷婷综合在线观看丁香| 午夜爽爽视频| 国产成人综合久久| 乱人伦99久久| 毛片网站免费在线观看| 国产精选自拍| 亚洲国产精品美女| 国内熟女少妇一线天| 亚洲色图欧美一区| 嫩草国产在线| 精品无码国产一区二区三区AV| 亚洲av成人无码网站在线观看| 99精品伊人久久久大香线蕉| 26uuu国产精品视频| 欧美精品色视频| 亚洲熟女中文字幕男人总站| 免费亚洲成人| 国产成人久视频免费| 日韩毛片在线视频| 国产剧情一区二区| 浮力影院国产第一页| 亚洲色图综合在线| 天堂在线视频精品| 亚洲AV人人澡人人双人| 91原创视频在线| 中国精品久久| 老熟妇喷水一区二区三区| 91国内视频在线观看| 视频一区亚洲| 国产精品99久久久| 国产丝袜第一页| 无码'专区第一页| 国产乱子伦手机在线| 99热这里都是国产精品| 国产亚洲精品97在线观看| 拍国产真实乱人偷精品| 免费无码网站| 国产一级特黄aa级特黄裸毛片| 青青青视频91在线 | 国产午夜人做人免费视频| 国产丝袜啪啪| 最新亚洲av女人的天堂| 日本成人精品视频| 亚洲视频一区在线| 999精品视频在线| 亚洲国产中文欧美在线人成大黄瓜 | 亚洲高清在线天堂精品| 国产精品天干天干在线观看| 国产日韩欧美黄色片免费观看| 国产福利免费视频| 国产毛片不卡| 麻豆精品在线播放| 国产成人精品午夜视频'| 久久综合九色综合97婷婷| 乱色熟女综合一区二区| 波多野结衣在线se| 日本高清免费一本在线观看| 日韩精品无码免费专网站| 一区二区三区国产| 久久精品人人做人人| 亚洲丝袜中文字幕| 午夜电影在线观看国产1区| 亚洲第一成人在线| 57pao国产成视频免费播放| 日韩 欧美 小说 综合网 另类| 亚洲第一视频网| 免费又黄又爽又猛大片午夜| 中文字幕亚洲综久久2021| 欧美α片免费观看| 无码电影在线观看| 无码视频国产精品一区二区| 手机在线免费毛片| 91美女视频在线| 亚洲高清在线天堂精品| 久久久久久久久久国产精品|