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

基于混合RANS/LES 方法的空腔流動仿真分析

2024-12-13 00:00:00謝露雷武濤張昌九廉佳趙博文
航空科學技術 2024年11期
關鍵詞:模態方法

摘要:飛機在飛行過程中打開起落架艙或武器艙會形成空腔流動,并產生高強度的壓力脈動和氣動噪聲,因此對空腔流動開展研究并分析其非定常流動特點具有重要意義。本文基于混合RANS/LES數值仿真方法,以空腔標模M219為研究對象,進行了空腔非定常流動的研究,并與風洞試驗結果進行了對比分析,驗證了方法的可靠性,詳細研究了搭接網格劃分策略,并進行網格敏感性分析,在保證計算精度的同時,提高了計算效率;接著從多個方面對計算得到的非定常流動進行后處理,得到了艙內流場的總聲壓、頻譜特性、瞬時和時均流場結構,為充分理解空腔流動的特點和艙內高強度噪聲的產生機理提供了依據;最后,分析了時間步長對空腔流動仿真結果的影響。通過研究,對空腔流動特點進行詳細分析,充分驗證了混合RANS/LES仿真方法在計算空腔流動特點方面的優越性與實用性,提高了對于混合RANS/LES仿真方法在工程應用中的認識。

關鍵詞:空腔流動;混合RANS/LES方法;搭接網格;網格敏感性;總聲壓級;頻譜特性;POD

中圖分類號:V249文獻標識碼:ADOI:10.19452/j.issn1007-5453.2024.11.006

空腔流動現象廣泛地存在于航空工程實際中,如飛機的起落架艙和武器艙、發動機燃燒室的火焰穩定器[1]??涨涣鲃哟嬖诿黠@的流動不穩定現象,在某些情況下,流動剪切層從空腔前緣脫離、破裂,并與空腔后壁面相互作用,產生反射聲波[2-4],進一步影響剪切層的穩定性,從而形成反饋回路。該反饋回路傾向于產生包括寬帶噪聲(broadbandnoise)和單音調噪聲(narrow-bandnoise)的強聲場,其中單音調噪聲又稱為Rossiter模式[5]。這些單音調噪聲根據空腔的幾何形狀、所受的流動狀態以及內埋物的存在而變化,對空腔結構、內埋武器及其他敏感設備會產生嚴重的危害。因此有必要使用計算流體力學(CFD)進行空腔非定常流動數值模擬研究,在滿足工程應用的前提下,盡可能準確地模擬空腔非定常流動現象。

對于空腔流動,工程上經常使用的CFD方法有雷諾平均N-S法(RANS)和混合雷諾平均N-S/大渦模擬(RANS/LES)方法。秦浩等[6]使用RANS方法進行帶彈彈艙三維流動數值計算研究。對于混合RANS/LES方法,兼顧了RANS方法和LES方法的優點,計算結果表明,與RANS方法相比,即使使用較粗的計算網格,也可以得到更好的計算結果,并能捕捉到更精細的流動分離和更高頻的噪聲特性。隨著計算網格的加密,該方法能夠得到更精細的空腔流動細節和更準確的聲學特性[7-8]。目前越來越多的研究人員開始使用RANS/LES混合類方法進行各種復雜形式的空腔流動的數值仿真。Loupy等[9]基于RANS/LES混合類方法進行了武器艙艙門開閉對艙內流動和噪聲的影響研究。國內對空腔流動的研究起步較晚,而且2007年以前主要以RANS和URANS數值仿真進行簡單的空腔流動模擬為主[6,10-13]。2008年以后,逐漸出現了使用LES或者RANS/LES混合類方法進行空腔非定常流動及噪聲的數值模擬,且RANS/LES混合類方法逐漸成為主要的數值模擬方法[14-19]。

本文針對空腔非定常流動,開展基于混合RANS/LES方法的仿真分析。以M219空腔標模為研究對象,使用混合RANS/LES仿真方法計算非定??涨涣鲃樱鶕鹘yRANS方法和混合RANS/LES仿真方法計算結果的對比驗證后者在非定常流動特征捕捉方面的可靠性,并通過網格敏感性、時間補償影響、空間流動特征等方面分析,進一步驗證混合RANS/LES仿真方法應用于空腔流動計算時的優越性。本文的研究工作對于混合RANS/LES方法的空腔流動特點分析與工程應用具有重要意義。

1M219空腔模型

仿真計算采用M219空腔試驗標模(見圖1),該試驗模型全長1.8288m,寬0.4318m,高0.14605m;模型內矩形空腔寬度W為0.1016m,長度L為0.508m,深度D為0.1016m,空腔前壁距離模型前緣0.7874m;該模型整體為具有尖劈形頭部的平板,且內部矩形空腔中心線與平板中心線平行不共線,偏移量為0.0254m。

空腔模型內脈動壓力測量點位置如圖2所示,各個測點等距分布,以K20~K29命名;測點之間相距0.0508m(空腔長度的10%),兩端測點距離空腔的前后壁長度均為0.0254m,且各測點分布在偏離空腔中心線0.0254m的平行線上。

風洞的自由來流Ma0.85,自由來流靜壓p=6.21×104Pa,自由來流靜溫T=263K?;诳涨婚L度的試驗雷諾數為Re=6.84×106。

2網格劃分

為提高計算資源利用效率,使用更高效的網格劃分策略十分重要,尤其對于非定常數值模擬計算而言,其本身對于計算資源的要求較高,更需要慎重選擇網格劃分方法。相較于非結構網格,傳統結構化網格在計算精度方面勝出一籌;但由于其本身具備傳導效應,對于需要局部加密的重點區域進行調整時將會同等程度地影響遠場網格分布,顯著提高網格數量,降低網格的計算效率,不利于計算資源的高效利用。相比之下,搭接網格[20-23]能夠有效規避這一缺陷,在保證重點區域的網格量和網格精度的同時控制總網格數量和合理提高網格效率。

圖3所示為搭接網格劃分區域示意圖,圖中紅色線條標注區域即為通過搭接網格局部加密的重點區域。在加密區域與外部區域的交界面兩側網格的位置并不一一對應,因此也稱搭接網格為非共形網格??紤]到空腔流動計算的重點關注區域應為空腔內部受剪切層影響區域,因此在內部網格劃分時盡量保持網格為立方體網格,各向尺度保持相近。外部區域同樣采用結構化網格,但網格尺度較大,且網格形狀不必保持標準立方體外形。

根據網格量和網格尺寸不同,將搭接的內部網格劃分為三類,其中粗網格(coarse)的最大網格尺寸為5mm,中等網格(medium)的最大網格尺寸為2.5mm,細網格(fine)的最大網格尺寸為1.25mm。表1為三類搭接網格的具體網格分布情況,根據表1中的數據可知,內部搭接網格的局部加密并不會影響到外部區域的網格劃分,恰好印證了搭接網格具有克服傳統結構網格傳導效應缺陷的能力。

3計算方法

常用數值模擬方法主要包括雷諾平均(RANS)、直接數值模擬(DNS)和大渦模擬(LES)等。在這些方法中,RANS方法因為對流場中部分時空流動細節進行了平均運算,其中的一些流動特征會因平均化處理而無法顯現,因而不適用于捕捉空腔流動中的小尺度渦;DNS方法能夠捕捉到流場中的脈動細節,但也正因如此,對網格和計算量有著極大要求;LES方法能夠較好地捕捉流場中的脈動特性,從而反映出流場中分離流動的真實特性,但是直接采用LES方法進行計算,其網格量和計算量仍然有著難以實現的需求,對于真實情況下的航空器其飛行雷諾數往往較高并有著較薄的邊界層,其中的小型渦結構的尺度則會更小,無法通過LES方法應用于工程實踐。

綜合上述各類方法優缺點,為較好模擬航空器邊界層中較小尺度的渦結構,可以采用混合RANS/LES方法進行模擬計算。該方法既涵蓋了RANS方法計算量小的優點,又能在遠離物面的位置處將湍流模型中的耗散項中湍流尺度用網格尺度參數和常數的乘積來代替,通過對大尺度的渦結構進行數值模擬來捕捉流場中的分離流特征。

混合RANS/LES方法經過近20年來的研究,目前已經發展出了眾多計算方法分支,主要包括應力混合渦流模擬(SBES)、尺度自適用模擬(SAS)、脫體渦模擬(DES)等;本文采用SBES方法對空腔流動進行數值模擬計算。由于RANS方程和空間過濾的LES方程在形式上具有一定相似性,在使用SBES湍流量混合模型時可以通過在不同流動區域對二者進行權重分配來更好地實現數值模擬計算,具體加權方法為

?hybrid=f?RANS+(1-f)?LES,0≤f≤1

圖4、圖5通過給出RANS和混合RANS/LES兩種方法的空間分離渦計算結果對比,更加明晰地體現了混合RANS/LES方法在計算空腔流動方面的優越性。

4計算結果分析

4.1網格敏感性分析

為分析不同尺度網格對空腔流動計算結果的敏感性,本文通過粗網格(coarse)、中等網格(medium)、細網格(fine)三種尺度的網格計算結果進行網格敏感性分析。在非定常計算中控制時間步長統一為2e-5s,并通過空腔底部K29監測點處的脈動壓力結果進行分析;根據圖6所示結果,空腔底部脈動壓力隨時間推進逐步收斂,在0.06s后基本達到完全收斂狀態。

圖7所示結果為空腔底部K20~K29監測點的總聲壓級監測結果,各監測點平行于空腔中心線等距分布。根據計算結果,空腔底部的總聲壓級分布沿流向呈遞增趨勢,從155dB逐漸增大到165dB。不同網格尺度下的計算結果差異并不大,誤差小于1.6%,且整體發展規律保持一致。

與此同時,針對不同尺度網格的計算結果進行快速傅里葉變換(FFT)分析,以監測點K29的聲壓頻譜特性為例,圖8所示為計算結果與試驗結果(EXP_1024)的對比情況。受限于計算資源等因素,實際數值計算模擬時長有限,對于第一階Rossiter模態(mode1)的捕捉不夠精確,在計算結果中沒有明顯體現。除了一階模態外,其余的第2~4階空腔流動Rossiter模態(mode2、mode3和mode4)的模擬都較明顯,從結果中能夠較清晰地識別出來,各階模態對應的頻率與聲壓級見表2。

根據計算結果,三種尺度網格均能夠捕捉到空腔內部自激振蕩模態,但與試驗結果仍然存在一定差異。受限于計算資源等原因,數值模擬計算的時間歷程遠小于試驗中的實際采樣時間歷程,通過增大計算時長,能夠有效改善這種誤差。

圖9所示為三種不同尺度網格下計算得到的瞬時空間渦結構;通過Q等值面體現出某時刻下空間渦的形態與位置特征。根據圖示結果可以發現,網格越精細,空間渦結構細節越精細。

4.2時間步長的影響

以中等規模網格為研究對象,在選取物理時間步長t=2e-5s基礎上,將物理時間步長變為2t、10t,研究時間步長對計算結果的影響。

不同時間步長對應的瞬時空間渦結構如圖10~圖12所示。可以看出,時間步長t和2t對應的空間渦結構差別不大;時間步長10t對應的空間渦結構明顯偏大,尤其是位于空腔后部的渦,明顯沒有得到很好的刻畫。

不同時間步長對應的總聲壓級分布如圖13所示。從圖13可以看出,時間步長t和2t對應的空腔底部的總聲壓級(OASPL)分布差別不大,總體趨勢與試驗結果符合得很好。但是時間步長10t對應的空腔底部總聲壓級OASPL分布與試驗結果的差別較大。典型監測點K29的頻譜特性分別如圖14所示。可以看出,時間步長過大,無法準確獲取高頻部分的噪聲特性。

因此,時間步長對計算結果有較大的影響,要與網格的尺度匹配,不能過大。

4.3艙內流動特征分析

圖15為中等規模網格計算得到的中間剖面的時均流場示意圖。符合典型的開式空腔流動特點,即從空腔前緣分離的剪切層直接撞擊到空腔并撞擊后緣,后壁面附近的壓力增加。由于剪切層與空腔底面之間有足夠的距離,剪切層得到充分發展,讓回流在空腔內可以自由發生,從而其脈動壓力的模態特性表現出全局性的特點。

在計算的過程中,對不同時刻的流場進行統計處理,可以得到任意位置的脈動壓力大小的分布。本部分使用總聲壓級OASPL表示腔內脈動壓力的大小。

圖16給出了空腔附近壁面上的總聲壓級分布,圖17為中間剖面上的總聲壓級分布,可以看出總聲壓級最大的區域位于空腔的后壁面及其附近的區域,其次是剪切層跨過空腔空向下游運動的區域。這也進一步說明了腔內產生高強度噪聲的噪聲源位于剪切層與后壁面的撞擊區域。

4.4空腔主要流場結構提取

開式空腔流動在一定條件下表現出強烈的自持振蕩,會引起高度的氣聲噪聲和結構疲勞。為了識別開式空腔非定常流動的主導模態和流動結構,采用本征正交分解(POD)方法[24-28]對非定常空腔流動進行了分析。

在去掉平均流場之后,圖18顯示了空腔流動前5階能量最大的POD模態。左列顯示了它們的三維渦結構,用Q準則的等值面表示,并用X方向的速度著色。中間列顯示了在中間剖面上X方向的速度云圖。右列顯示了在中間剖面上,Y方向的速度云圖。模態1在空腔的后部包含一個大型結構。模態2和模態3形成模態對,從Y方向的速度云圖可以清晰地看到,在空腔口顯示兩個波長,并有一定的相位差。這種相位差是穿過空腔的對流運動的標志。高階的POD模態顯示更精細的空間尺度,它包含少量的能量。

5結論

通過研究,可以得出如下結論:

(1)湍流量混合模型SBES作為混合RANS/LES方法的一種,可以很好地預測空腔非定常流動。艙內的總聲壓級OASPL分布、脈動壓力的模態特征與試驗結果吻合得很好。

(2)使用搭接網格可以很好地降低結構網格的數量,從而提高空腔非定常數值模擬的技術效率;通過網格敏感性的研究,證實了SBES方法在計算空腔非定常流動方面具有很好的魯棒性,即使在非常粗的網格上,也可以得到很好的結果。

(3)時間步長對計算結果有較大的影響,要與網格的尺度匹配,可以根據空腔內部網格的平均尺寸和來流速度進行確定。

(4)后續將借助混合RANS/LES類方法進行復雜構型空腔流動的仿真分析,并進行流動控制方法的研究,以降低艙內的高強度脈動壓力。

參考文獻

[1]ThomasCF,MingT,JackRE.ValidationofahybridReynolds-Averaged/Large-Eddysimulationmethodforsimulatingcavityflameholderconfigurations[R].AIAA2001-2929,2001.

[2]SchmitR,GroveJ,SemmelmayerF,etal.Nonlinearfeedbackmechanismsinsidearectangularcavity[J].AIAAJournal,2014,52(10):2127-2142.

[3]LoupyG,BarakosG.Modellingoftransonicshallowcavityflowsandstroereleasesimulationsfromweaponbays[R].AIAA2017-3252,2017.

[4]JeffreyMM.Effectsofcavitydimensions,boundarylayerandtemperatureoncavitynoisegenerationandcontrol[D].Georgia:GeorgiaInstituteofTechnology,1997.

[5]RossiterJE.Windtunnelexperimentsontheflowoverrectangularcavitiesatsubsonicandtransonicspeeds[R].TechnicalReport64037,1964.

[6]秦浩,周長悅.帶彈彈艙三維流動數值計算研究[J].航空科學技術,2017,28(7):25-29.

QinHao,ZhouChangyue.Studyonnumericalcalculationof3Dflowinmissilebay[J].AeronauticalScienceamp;Technology,2017,28(7):25-29.(inChinese)

[7]NayyarP,BarakosGN,BadcockJN,etal.Numericalstudyoftransoniccavityflowsusinglarge-eddyanddetached-eddysimulation[J].AeronauticalJournal,2007,111:153-164.

[8]RichardA,FredM.RANSandDESturbulencemodelpredictions

[9]LoupyGJM,BarakosGN,TaylorNJ.Cavityflowoveratransonicweaponbayduringdooroperation[J].JournalofAircraft,2018,55(1):339-354.

[10]侯中喜,易仕和,王承堯.超聲速開式空腔流動的數值模擬[J].推進技術,2001,22(5):400-403.

HouZhongxi,YiShihe,WangChengyao.Numericalanalysisofsupersonicopencavity[J].JournalofPropulsionTechnology,2001,22(5):400-403.(inChinese)

[11]司海青,王同光.邊界條件對三維空腔流動振蕩的影響[J].南京航空航天大學學報,2006,38(5):595-599.

SiHaiqing,WangTongguang.Influenceofboundaryconditionon3Dcavityflow-inducedoscillations[J].JournalofNanjingUniversityofAeronauticsamp;Astronautics,2006,38(5):595-599.(inChinese)

[12]司海青,王同光.數值模擬有外掛物的空腔流動[J].空氣動力學學報,2007,25(3):404-409.

SiHaiqing,WangTongguang.Numericalsimulationsofthecavitywithastore[J].ActaAerodynamicaSinica,2007,25(3):404-409.(inChinese)

[13]李曉東,劉靖東,高軍輝.空腔流激振蕩發聲的數值模擬研究[J].力學學報,2006,38(5):599-604.

LiXiaodong,LiuJingdong,GaoJunhui.Numericalsimulationofflow-inducedoscillationandsoundgenerationinacavity[J].ChineseJournalofTheoreticalandAppliedMechanics,2006,38(5):599-604.(inChinese)

[14]賴煥新,周邵萍,蘇永升,等.空腔流動的大渦模擬及氣動噪聲控制[J].工程熱物理學報,2008,29(2):228-232.

LaiHuanxin,ZhouShaoping,SuYongsheng,etal.Large-eddysimulationandcontrollingofnoiseincavity[J].JournalofEngineeringThermophysics,2008,29(2):228-232.(inChinese)

[15]譚玉婷.空腔非定常流動特性的數值模擬研究[D].南京:南京航空航天大學,2009.

TanYuting.Numericalsimulationoftheunsteadycavityflowfeatures[D].Nanjing:NanjingUniversityofAeronauticsandAstronautics,2009.(inChinese)

[16]譚玉婷,伍貽兆,田書玲.基于DES的二維和三維空腔流動特性研究[J].航空計算技術,2010,40(1):67-70.

TanYuting,WuYizhao,TianShuling.Numericalsimulationof2D/3DcavityflowsusingDES[J].AeronauticalComputingTechnique,2010,40(1):67-70.(inChinese)

[17]肖志祥,羅堃宇,劉健.寬速域RANS-LES混合方法的發展及應用[J].空氣動力學學報,2017,35(3):338-353.

XiaoZhixiang,LuoKunyu,LiuJian.DevelopmentsandapplicationsofhybridRANS-LESmethodsforwide-speedrangeflows[J].ActaAerodynamicaSinica,2017,35(3):338-353.(inChinese)

[18]余培汛,白俊強,郭博智,等.射流對空腔噪聲抑制效果研究[J].計算力學學報,2014(5):663-669.

YuPeixun,BaiJunqiang,GuoBozhi,etal.Suppressioneffectofjetflowonaerodynamicnoiseofcavity[J].ChineseJournalofComputationalMechanics,2014(5):663-669.(inChinese)

[19]劉俊,楊黨國,王顯圣,等.湍流邊界層厚度對三維空腔流動的影響[J].航空學報,2015,37(2):475-483.

LiuJun,YangDangguo,WangXiansheng,etal.Effectofturbulentboundarylayerthicknessonathree-dimensionalcavityflow[J].ActaAeronauticaetAstronauticaSinica,2015,37(2):475-483.(inChinese)

[20]趙軻,高正紅,黃江濤,等.拼接網格技術在復雜流場數值模擬中的應用研究[J].應用力學學報,2011,28(1):69-74.

ZhaoKe,GaoZhenghong,HuangJiangtao,etal.Applicationsofthepatched-gridtechnologyinnumericalsimulationofflowfield[J].ChineseJournalofAppliedMechanics,2011,28(1):69-74.(inChinese)

[21]崔英俊,潘若癡,孟德君,等.不同類型搭接網格對周向槽處理機匣模擬結果影響研究[J].推進技術,2021,42(6):1265-1275.

CuiYingjun,PanRuochi,MengDejun,etal.Effectsofdifferenttypesofpatched-gridonsimulatedresultsrelatedtocircumferentialgrooves[J].JournalofPropulsionTechnology,2021,42(6):1265-1275.(inChinese)

[22]艾俊強,謝露.基于混合RANS/LES方法的亞聲速空腔流動主要影響因素的數值研究[J].南京航空航天大學學報,2022,54(5):927-936.

AiJunqiang,XieLu.NumericalstudyonmaininfluencingfactorsofsubsoniccavityflowbasedonhybridRANS/LESmethod[J].JournalofNanjingUniversityofAeronauticsamp;Astronautics,2022,54(5):927-936.(inChinese)

[23]謝露,艾俊強,李權,等.長深比對空腔流動與聲學特性的影響分析[J].航空工程進展,2014,5(1):18-24.

XieLu.AiJunqiang,LiQuan,etal.Theeffectoflength-todepthratioonflowandaeroacousticcharacteristicsofcavity[J].AdvancesinAeronauticalScienceandEngineering,2014,5(1):18-24.(inChinese)

[24]PhilippeRS.Detached-eddysimulation[J].AnnualReviewofFluidMechanics,2009,41:181-202.

[25]RokitaT,GreenbergJB,ArieliR,etal.Spatial-temporalpatternsofthree-dimensionalsubsonicturbulentcavityflow[J].InternationalJournalofHeatandFluidFlow,2018,71:260-274.

[26]MurrayN,S?llstr?mE,UkeileyL.Propertiesofsubsonicopencavityflowfields[J].PhysicsofFluids,2009,21(9):095103.

[27]RowleyCW,ColoniusT,MurrayRM.PODbasedmodelsofself-sustainedoscillationsintheflowpastanopencavity[R].AIAA2000-1969,2000.

[28]劉宏鵬.三維效應對超聲速湍流燃燒流動大渦模擬的影響研究[J].航空科學技術,2022,33(7):57-65.

LiuHongpeng.Onthethreedimensionaleffectsforthelargeeddysimulationofsupersonicturbulentcombustionflows[J].AeronauticalScienceamp;Technology,2022,33(7):57-65.(inChinese)

基金項目:航空科學基金(2019ZA003002)

猜你喜歡
模態方法
學習方法
車輛CAE分析中自由模態和約束模態的應用與對比
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
國內多模態教學研究回顧與展望
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
高速顫振模型設計中顫振主要模態的判斷
航空學報(2015年4期)2015-05-07 06:43:35
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: v天堂中文在线| 亚洲中文字幕久久精品无码一区| 熟妇丰满人妻av无码区| 欧美国产在线看| 国内精品小视频福利网址| 婷婷综合缴情亚洲五月伊| 无码精品福利一区二区三区| 亚洲欧洲日韩综合色天使| 91久久夜色精品| 日韩中文无码av超清| a国产精品| 国产剧情一区二区| 永久在线精品免费视频观看| 人妻丰满熟妇AV无码区| 真人免费一级毛片一区二区| 热思思久久免费视频| 免费毛片a| 国产精品久久久精品三级| 一级毛片免费不卡在线| 亚洲午夜国产精品无卡| 18禁高潮出水呻吟娇喘蜜芽| 亚洲日韩AV无码一区二区三区人 | 国产内射一区亚洲| 日本精品中文字幕在线不卡 | 米奇精品一区二区三区| 欧美一区二区精品久久久| 日本三级欧美三级| 99视频在线免费观看| 久草青青在线视频| 免费久久一级欧美特大黄| 91精品国产91欠久久久久| 日韩高清在线观看不卡一区二区| 亚洲成在线观看| 老司机久久精品视频| 成年人国产网站| 激情国产精品一区| 2024av在线无码中文最新| 青青操视频免费观看| 午夜色综合| 久久成人免费| 欧美人在线一区二区三区| a色毛片免费视频| 成人国产精品视频频| 亚洲无卡视频| 亚洲欧美极品| 国产原创自拍不卡第一页| 国产欧美精品一区二区 | 国产尤物视频网址导航| 日韩久草视频| 国产三级国产精品国产普男人| 久久久久亚洲精品成人网| 亚洲欧美国产视频| 国产第四页| 香蕉网久久| 久草网视频在线| 亚洲日韩Av中文字幕无码 | 热久久国产| 欧美综合一区二区三区| 九九热精品免费视频| 真实国产乱子伦视频| 亚洲永久视频| 老司国产精品视频| 人妻出轨无码中文一区二区| 国语少妇高潮| 色老二精品视频在线观看| 国产免费羞羞视频| 亚洲欧洲日韩久久狠狠爱| 欧美激情一区二区三区成人| 精品一区二区三区无码视频无码| 亚洲香蕉久久| 91久久偷偷做嫩草影院电| 精品国产污污免费网站| 亚洲国产精品无码久久一线| 99久久国产综合精品2020| 毛片在线看网站| 精品第一国产综合精品Aⅴ| 久草中文网| 亚洲中文精品人人永久免费| 99久久精品久久久久久婷婷| 国产成人一区免费观看| 美女内射视频WWW网站午夜| 精品久久久久无码|