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

    細(xì)長旋成體亞聲速超大攻角非定常流動特性研究

    2022-03-20 15:52:42王方劍王宏偉李曉輝
    力學(xué)學(xué)報(bào) 2022年2期
    關(guān)鍵詞:物面油流細(xì)長

    王方劍 王宏偉 李曉輝 董 磊 黃 湛 陳 蘭

    (中國航天空氣動力技術(shù)研究院,北京 100074)

    引言

    空空導(dǎo)彈作為現(xiàn)代空戰(zhàn)的主要攻擊手段,要求比目標(biāo)飛機(jī)更高的機(jī)動性和敏捷性.新型空空導(dǎo)彈在面對新一代飛機(jī)時必須具備全方位攻擊能力,尤其對來自后方目標(biāo)的威脅,則需要更高轉(zhuǎn)彎率和更大機(jī)動包絡(luò)線的航向反轉(zhuǎn)機(jī)動等先進(jìn)高效機(jī)動方法.為了保證高效機(jī)動的順利完成,要求導(dǎo)彈在超大攻角(α=0°~ 180°)范圍內(nèi)具有飛行和機(jī)動控制能力[1].以往對超大攻角流動的觀測和研究大多集中在α=40°~ 60°范圍內(nèi),最大角度不超過90°[2-14].在超大攻角條件下,復(fù)雜的非定常分離流問題廣泛存在,包括非定常渦的產(chǎn)生、脫落、干擾和非定常非對稱等氣動現(xiàn)象[2-10].因此,提高對細(xì)長體復(fù)雜非定常氣動特性的認(rèn)識尤為重要.

    早在20 世紀(jì)50 年代[3],對細(xì)長體(錐柱體)分離流動特性的探索就已經(jīng)開始.在小攻角處,細(xì)長體背風(fēng)側(cè)形成對稱的軸向集中渦[4].隨著攻角的增大(達(dá)到30°左右),流型由對稱型轉(zhuǎn)變?yōu)榉菍ΨQ型.隨著旋渦沿軸向向后發(fā)展,旋渦一個接一個離開物面,形成旋渦脫落流動結(jié)構(gòu),并以一定角度向下游流動.此時如果觀察其橫截面,會發(fā)現(xiàn)形成類卡門渦街的流動形式[5].

    在無側(cè)滑來流條件下的側(cè)向力是非對稱渦的重要標(biāo)志[9-10].Lamont 和Hunt[11]的研究表明,在大部分攻角范圍內(nèi)(0°~ 90°),沿旋成體軸向的側(cè)向力呈振蕩分布.Dexter[12]發(fā)現(xiàn)側(cè)向力強(qiáng)烈依賴于旋成體滾轉(zhuǎn)角.在一個可獨(dú)立滾轉(zhuǎn)的頭尖部上進(jìn)行的風(fēng)洞試驗(yàn)表明,滾轉(zhuǎn)角所造成的對側(cè)向力的影響主要頭尖部帶來.Luo 等[13]在一個橢圓形截面頭尖部的旋成體上進(jìn)行了壓力分布測量,結(jié)果表明,橢圓形截面頭尖部側(cè)向力隨滾轉(zhuǎn)角的變化規(guī)律比圓形更容易預(yù)測,且兩種頭尖部的最大側(cè)向力基本相同.通過在頭尖部附近安裝一根可伸縮的鋼絲[14]或沙礫[15],可以通過改變頭尖部的周向角來有效地控制側(cè)向力的方向.大多數(shù)文獻(xiàn)采用的試驗(yàn)?zāi)P褪穷^部尖銳的旋成體外形,然而頭部形狀也會對流動產(chǎn)生一定影響.Hsieh 等[16-17]對半球形頭部旋成體在攻角從0°~50°的工況下的分離流進(jìn)行了試驗(yàn)和數(shù)值研究.令人驚訝的是,分離流包含了許多復(fù)雜的流動結(jié)構(gòu)特征,這是三維流所獨(dú)有的.結(jié)果表明,半球形頭部旋成體同時存在兩個分離區(qū)(頭部分離泡和橫流分離區(qū)).

    20 世紀(jì)80 年代以前,由于缺乏強(qiáng)大的計(jì)算軟件和先進(jìn)的風(fēng)洞試驗(yàn)方法,學(xué)者們關(guān)注點(diǎn)主要集中在附著流和穩(wěn)定分離區(qū)上[18].隨著技術(shù)的發(fā)展,一些非定常流動特性被揭示出來.Zeiger 等[19]對細(xì)長旋成體的流動特性進(jìn)行流動顯示水洞試驗(yàn),認(rèn)為在大攻角時,細(xì)長體旋成體上的流動可主要分為3 種流動形態(tài),第1 種流動形態(tài)是頭尖部附近的產(chǎn)生的較為集中的渦旋,隨著流動沿軸向向后發(fā)展,流動結(jié)構(gòu)逐漸體現(xiàn)為旋渦脫落流動形態(tài)(第2 種流動形態(tài)),渦脫落方向與旋成體呈一定角度,當(dāng)流動發(fā)展至接近底部時(橫截面位置X/D 從10 到32),脫落的旋渦渦軸與旋成體平行,形成第3 流動形態(tài).管小榮和徐誠[20]、楊云軍和周偉江[21]、張贏和劉超峰[22]、劉仙名和符松[23],采用數(shù)值模擬方法細(xì)致的研究了非對稱旋渦的演化與發(fā)展,研究表明不同的Re、湍流模型均會對旋渦非對稱流動的模擬產(chǎn)生一定影響.Degani 和Zilliac[24]的研究表明,其非定常流動主要包括類卡門渦街的低頻流動、剪切層失穩(wěn)引起的高頻流動,以及介于以上兩種頻率之間的旋渦干擾流動.Ayoub 和Karamcheti[25]同時測量了細(xì)長旋成體表面壓力和尾流速度波動.結(jié)果表明,在頭部區(qū)域存在一定的旋渦脫落現(xiàn)象.在離旋成體頭部不同距離處,渦脫落產(chǎn)生了一個接近10.5 Hz 的頻率.

    除此之外,馬赫數(shù)對旋成體分離流動細(xì)節(jié)也有一定影響,Keneer[26]對旋成體前體在不同馬赫數(shù)及攻角細(xì)長體流動進(jìn)行研究,在風(fēng)洞試驗(yàn)紋影照片中可以觀察到,第一脫落渦在高亞聲速下(Ma=0.6)的軌跡比低速(Ma=0.25)時更接近表面,側(cè)力系數(shù)從Ma=0.25 時的1.9 降至Ma=0.6 時的0.8,同時也可以看到,在Ma=0.6 時,攻角55°時后半部分出現(xiàn)不對稱,旋渦較為分散,而Ma=0.25 在此攻角下攻角還較為集中,這體現(xiàn)出旋成體高亞音速流動與低速旋渦流動相比,分離流動的旋渦位置、旋渦強(qiáng)度、渦破裂位置、分離位置均有一定的不同.

    在傳統(tǒng)的復(fù)雜分離流動數(shù)值模擬工作中,基于雷諾平均Navier-Stokes (RANS)方程的流場預(yù)測方法對附著流以及小分離流動有較好的模擬能力,計(jì)算成本較低,是工程運(yùn)用中最廣泛的一種方法,但由于RANS 方法平均了小尺度的旋渦脈動,難以精細(xì)的模擬較大范圍的分離與復(fù)雜的旋渦流動.大渦模擬方法(LES)能夠較好的模擬大范圍分離旋渦流動現(xiàn)象,但對邊界層的網(wǎng)格密度有非常高的要求,帶來很大的計(jì)算成本.所以為了結(jié)合RANS 方法與LES方法的優(yōu)勢,Spalart[27]提出了基于S-A 湍流模型的DES 方法,引入特征長度參數(shù),根據(jù)特征長度選擇不同的流動模擬方法,在近壁面區(qū)域,采用RANS 方法,在遠(yuǎn)離物面的區(qū)域,選擇LES 方法進(jìn)行模擬,能夠兼顧數(shù)值計(jì)算精度與效率.深入研究發(fā)現(xiàn),DES 方法會發(fā)生模型應(yīng)力耗散現(xiàn)象(modeled stress depletion,MSD),導(dǎo)致在流動中出現(xiàn)網(wǎng)格誘導(dǎo)分離問題(grid induced separa-tion,GIS).2006 年,Spalart 等[28-37]提出了延遲脫體渦模擬方法(DDES)有效解決MSD與GIS 問題,該方法被廣泛應(yīng)用于圓柱擾流、空腔流動和起落架分離流動能計(jì)算研究中.

    從細(xì)長體大攻角的研究情況來看,之前對細(xì)長體大攻角流動研究速域主要集中在低速,并且研究攻角范圍主要在α=0°~ 90°之間,本文針對細(xì)長旋成體在高速(Ma=0.6)超大攻角下(α=0°~ 180°)存在的旋渦非定常性、旋渦脫落和旋渦非對稱等復(fù)雜流動現(xiàn)象,采用基于結(jié)構(gòu)網(wǎng)格的DDES 數(shù)值方法,開展細(xì)長旋成體復(fù)雜氣動非定常特性研究,為新型機(jī)動導(dǎo)彈超大攻角氣動設(shè)計(jì)提供理論支撐.

    1 研究方法

    本文采用數(shù)值模擬與油流顯示風(fēng)洞試驗(yàn)相結(jié)合的方法來研究細(xì)長體超大攻角流動特性問題.其中數(shù)值模擬采用Roe 方法求解Navier-Stokes 方程,細(xì)致分析細(xì)長體超大攻角在不同攻角下的非定常流動特性及壓力脈動特性,油流顯示風(fēng)洞試驗(yàn)用于獲得細(xì)長體模型表面的流動拓?fù)?通過對比試驗(yàn)結(jié)果與數(shù)值結(jié)果,驗(yàn)證數(shù)值模擬結(jié)果的準(zhǔn)確性.

    1.1 數(shù)值方法與算例驗(yàn)證

    控制方程為三維可壓縮Navier-Stokes 方程,廣義坐標(biāo)表示為

    空間離散采用基于MUSCL 方法插值方法的FDS-Roe 格式,黏性項(xiàng)采用二階中心差分離散,對流和壓力項(xiàng)使用3 階迎風(fēng)格式.時間推進(jìn)格式為隱式LUSGS 方法,保證了較高的時間計(jì)算精度.非定常計(jì)算采用雙時間步長方法,同時采用多重網(wǎng)格算法加快子迭代的收斂速度.

    在用湍流模型進(jìn)行流動模擬時,控制方程中應(yīng)加入雷諾應(yīng)力相關(guān)項(xiàng).本文采用的湍流模型是基于Menter 的兩方程剪切應(yīng)力輸運(yùn)湍流模型(SST),表達(dá)式為

    DES 方法的擴(kuò)展形式是基于SST 湍流模型,其中湍流尺度為

    DES 長度尺度dDES表示為

    式中,Δ 為當(dāng)前網(wǎng)格單元與周邊單元的最大距離,CDES是SST 湍流模型中附加的經(jīng)驗(yàn)?zāi)P统?shù),CDES表示為

    為了解決MSD 與GIS 問題問題,Spalart[27]提出了延遲函數(shù)fd,將SST 湍流模型的湍流能量輸運(yùn)方程的耗散項(xiàng)重新定義為與網(wǎng)格尺度和湍流黏度有關(guān)的DES 長度尺度,結(jié)合延遲函數(shù)fd的長度尺度dDES被重新定義為

    因此,dDES不僅與網(wǎng)格尺度有關(guān),而且與湍流黏度有關(guān),可以防止湍流模擬過早地切換到LES 模式.

    算例驗(yàn)證采用高雷諾數(shù)圓柱擾流算例進(jìn)行計(jì)算,圓柱直徑為D=457 mm,圓柱長度為5.66D,圓柱兩端采用對稱邊界條件,其他邊界采用自由來流條件.圖1 為圓柱擾流計(jì)算域,為保證滿足自由來流邊界條件,在流向計(jì)算域?yàn)?0D,圖2 是圓柱橫截面網(wǎng)格.計(jì)算參數(shù)為:Ma=0.25,Re=3 × 106,物理時間步長為Δt=4 × 10-4s,一共計(jì)算5000 個物理時間步,計(jì)算物理時間為2 s.

    圖1 圓柱擾流計(jì)算域Fig.1 Domain of calculation

    圖2 圓柱橫截面網(wǎng)格Fig.2 Zoomed-in-view of mesh around the cylinder

    圖3 是圓柱擾流壓力分布試驗(yàn)值[38-40]與計(jì)算值壓力分布對比,其中縱軸是壓力系數(shù),橫軸是方位角θ,計(jì)算值是截取了圓柱中間截面(y=2.83D)的壓力分布時均值,通過計(jì)算獲得的圓柱升力斯特勞哈爾數(shù)St=0.32,并將取一個周期內(nèi)的60 個瞬時值進(jìn)行平均,得到時均值.圖中顯示,計(jì)算值與試驗(yàn)值吻合較好,說明本文所采用的DDES 方法能夠較好的捕捉圓柱尾跡的旋渦結(jié)構(gòu),見圖4.

    圖3 壓力分布試驗(yàn)值[38-40]與計(jì)算值壓力分布對比Fig.3 Comparisons of pressure coefficients distribution[38-40]

    圖4 圓柱擾流尾跡流動結(jié)構(gòu)Fig.4 Q-criterion iso-surface of wake

    1.2 油流流動顯示技術(shù)

    為了便于油流的光學(xué)測量技術(shù)的應(yīng)用,試驗(yàn)選定在中國航天空氣動力技術(shù)研究院FD-12 風(fēng)洞常規(guī)亞跨聲速試驗(yàn)段作為流場顯示試驗(yàn)的主要試驗(yàn)環(huán)境,該試驗(yàn)段長度為3.8 m,試驗(yàn)段橫截面尺寸為1.2 m ×1.2 m,兩側(cè)沿來流方向?qū)ΨQ布置4 個窗口,內(nèi)芯頂部也具有光學(xué)窗口,是目前能夠?qū)崿F(xiàn)常規(guī)攻角機(jī)構(gòu)試驗(yàn)光學(xué)條件最好的試驗(yàn)段.由于模型需要安裝在常規(guī)攻角機(jī)構(gòu)上,攻角調(diào)節(jié)范圍是-15°~ 25°,無法實(shí)現(xiàn)連續(xù)的大攻角變化,需要通過預(yù)置攻角模塊來實(shí)現(xiàn)不同階段的攻角調(diào)節(jié).

    如圖5,模型為了實(shí)現(xiàn)180°可翻轉(zhuǎn)設(shè)計(jì)安裝,頭部和尾部均可與支撐段連接,因此,頭部需要設(shè)計(jì)可拆卸段,適合尾部的支撐連接.

    圖5 模型的支撐設(shè)計(jì)Fig.5 Support design of experimental model

    由前述風(fēng)洞試驗(yàn)條件可知,常規(guī)試驗(yàn)段和攻角機(jī)構(gòu)的攻角調(diào)節(jié)范圍是-15°~ 25°,流動顯示試驗(yàn)?zāi)P蜑榱嗽陲L(fēng)洞中實(shí)現(xiàn)0°~ 180°的攻角變化,需要對模型的攻角進(jìn)行預(yù)置.在0°~ 30°范圍內(nèi),預(yù)置-15°的攻角,在60°~ 90°范圍內(nèi),預(yù)置-75°的攻角,配合模型頭尾更換,可以實(shí)現(xiàn)模型0°~ 180°的攻角變化.

    傳統(tǒng)油流試驗(yàn),一般在風(fēng)洞運(yùn)行結(jié)束后,對表面油流圖像進(jìn)行拍攝,獲取表面流動結(jié)構(gòu).在高速、大攻角條件下,較強(qiáng)的流動分離、風(fēng)洞臺階波干擾、油流自身重力等時刻在影響著表面流動圖譜,而暫沖式風(fēng)洞關(guān)車時氣流沖擊對表面流動影響非常大,風(fēng)洞運(yùn)行結(jié)束后的表面流動圖譜可能與風(fēng)洞運(yùn)行時的表面流動圖譜純在明顯差別.因此,本試驗(yàn)采用實(shí)時化表面油流圖像拍攝方法,在風(fēng)洞運(yùn)行時,按照下圖布局進(jìn)行拍攝.同時,選取特定工況進(jìn)行精細(xì)化的熒光油膜摩阻測量試驗(yàn),采用紫外光源的照明系統(tǒng)、混有熒光指示劑的油膜、帶有高通濾波片的圖像采集系統(tǒng),針對模型局部區(qū)域進(jìn)行試驗(yàn)拍攝,見圖6.

    圖6 油流流動顯示試驗(yàn)測量布局Fig.6 Oil-flow visualization test measurement layout

    2 細(xì)長旋成體外形及網(wǎng)格

    現(xiàn)代高機(jī)動導(dǎo)彈除了細(xì)長彈身外,還具有鴨舵、尾翼等控制部件,Wang 等[41]的研究結(jié)果顯示,對于細(xì)長導(dǎo)彈大攻角流動,其非定常特性、頻率特性主要由細(xì)長彈身誘導(dǎo),彈翼貢獻(xiàn)較小,所以為了能夠更清晰的研究超大攻角下的流動,本文采用的細(xì)長旋成體外形.如圖7,彈身直徑為60 mm,彈身全長為707.7 mm,彈身頭部采用的是尖型頭部設(shè)計(jì),坐標(biāo)設(shè)置為X 軸向后,Z 軸向上.

    圖7 細(xì)長旋成體外形 (單位:mm)Fig.7 Slender revolutionary body (unit:mm)

    計(jì)算網(wǎng)格采用結(jié)構(gòu)網(wǎng)格進(jìn)行求解,如圖8,采用O 型網(wǎng)格拓?fù)?外邊界距離的設(shè)置以能夠滿足自由來流條件為判斷準(zhǔn)則.貼近物面的第一層網(wǎng)格厚度保持Y+~ 1,以確保邊界層的準(zhǔn)確模擬.為了能夠精確的捕捉細(xì)長旋成體背風(fēng)側(cè)的分離渦流動,在模型的法向進(jìn)行了整體加密,在模型背風(fēng)側(cè)展向網(wǎng)格進(jìn)行適當(dāng)加密,同時當(dāng)攻角超過90°時,底部超前,這時底部引起的回流區(qū)模擬也非常重要,此時會在尾部附近進(jìn)行加密,總網(wǎng)格量為1500 萬.

    圖8 細(xì)長體網(wǎng)格Fig.8 Mesh of slender revolutionary body

    由于本文針對大攻角下非定常流動特性問題,其對旋渦尺度的分辨以及旋渦頻率特性的捕捉十分重要,所以需要對網(wǎng)格與時間步長無關(guān)性進(jìn)行檢查,即確定數(shù)值模擬結(jié)果不再隨網(wǎng)格規(guī)模與時間步長大小而變化.因此本文對典型工況(Ma=0.6,α=60°)進(jìn)行模擬,并選取升力CL與阻力Cd的平均值,以及側(cè)向力的斯特勞哈爾數(shù)St 進(jìn)行分析.首先進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證,采用了3 套計(jì)算網(wǎng)格,3 套網(wǎng)格在流向、周向和法向均進(jìn)行變化,網(wǎng)格量分別為530 萬、1500 萬和4200 萬,圖9 展示了3 套網(wǎng)格的物面網(wǎng)格分布.表1 為3 套網(wǎng)格下的升力CL、阻力Cd、側(cè)向力St 計(jì)算結(jié)果.從表1 我們可以看到,網(wǎng)格規(guī)模對CL,Cd與St 均有一定的影響,隨著網(wǎng)格量的增加,當(dāng)網(wǎng)格量達(dá)到1500 萬時,其計(jì)算結(jié)果(CL,Cd,St)與網(wǎng)格量4200 萬較為接近,證明網(wǎng)格量1500 萬對于研究該問題已滿足要求.

    表1 網(wǎng)格無關(guān)性Table 1 Sensitivity of mesh resolution

    圖9 3 套不同規(guī)模網(wǎng)格示意Fig.9 Three grids with different cell counts

    在時間步長無關(guān)性方面,分別采用物理時間步長為dt=1.2 × 10-5s,2.3 × 10-5s,5 × 10-5s 進(jìn)行計(jì)算,網(wǎng)格量為1500 萬.表2 為不同時間步長下的計(jì)算那結(jié)果,從表2 中我們可以看到,時間步長的變化對升力CL、阻力Cd影響不大,但對St 影響相對較大,隨著時間步長不斷減小,當(dāng)時間步長達(dá)到dt=2.3 × 10-5s 時,CL,Cd與St 基本不隨時間步長而變化,證明時間步長dt=2.3 × 10-5s 能夠滿足研究要求.

    表2 時間步長無關(guān)性Table 2 Sensitivity of physical time step

    3 超大攻角瞬時流動特性

    本文采用CFD 數(shù)值模擬手段,計(jì)算了細(xì)長旋成體在攻角α=0°~ 180°下的背風(fēng)側(cè)流動形態(tài),側(cè)滑角β=0°,雷諾數(shù)Re=1.29 × 106m-1,參考長度為直徑D=60 mm.由于細(xì)長體在大攻角下具有較強(qiáng)的非定常效應(yīng),所以這里采用數(shù)值模擬得到的瞬時流場與試驗(yàn)進(jìn)行對比.

    圖10 為攻角30°,CFD 模擬的表面極限流線分布,以及CFD 與油流試驗(yàn)圖像比較,從圖中我們可以看到,油流試驗(yàn)與數(shù)值模擬得到的主分離線位置較為一致.主分離線在模型頭尖部起始并較為靠近背風(fēng)側(cè),隨著流動沿軸向向后發(fā)展,分離線周向位置逐漸向迎風(fēng)側(cè)移動,當(dāng)流動發(fā)展至彈身圓柱段時,分離線周向位置沿著軸向幾乎不變.

    圖10 CFD 物面流線與油流試驗(yàn)對比(α=30°)Fig.10 Surface streamlines comparisons between CFD and oil-flow test (α=30°)

    圖11 是攻角30°是旋成體背風(fēng)側(cè)流動及物面流線,其中背風(fēng)側(cè)流動通過采用不同截面下的渦量值來體現(xiàn).從圖中可以看到,旋成體背風(fēng)側(cè)主要呈現(xiàn)出了對稱渦的流動形態(tài),旋渦從頭部起始,一直沿軸向向后發(fā)展至尾部.并且在主旋渦下方靠近物面的位置還存在二次渦流動結(jié)構(gòu).從物面流線分布也能夠看到明顯的主分離線以及二次分離線.圖12 是X/L=0.1,0.2,0.4,0.6,0.8 截面的壓力分布情況,圖中顯示5 個截面下壓力分布均為對稱分布形式,體現(xiàn)出旋成體在該攻角的對稱旋渦流動結(jié)構(gòu).

    圖11 物面流線及空間流動(α=30°)Fig.11 Surface streamlines and flow (α=30°)

    圖12 各截面壓力分布(α=30°)Fig.12 Pressure distribution of cross sections (α=30°)

    圖13 是CFD 物面流線與油流試驗(yàn)對比圖,圖中有4 張圖,從左至右分別為油流試驗(yàn)拍攝的照片(圖13(a)左)、通過試驗(yàn)照片處理出的物面流線圖(圖13(a)右)、CFD 手段獲得物面流線圖(圖13(b)左)、CFD 獲得的橫截面流線圖(圖13(b)右).從圖中試驗(yàn)獲得的物面流線圖(圖13(a)左)與CFD 獲得的物面流線圖(圖13(b)左)進(jìn)行對比發(fā)現(xiàn),存在清晰的主分離線、二次分離線、二次再附線.試驗(yàn)與CFD 獲得分離位置與再附位置對比較好.圖13(b)也展示了橫截面流線,從圖中可以看到較為清晰的主分離渦流動結(jié)構(gòu),同時在物面附近存在較小的二次渦流動結(jié)構(gòu).

    圖13 CFD 物面流線與油流試對比(α=60°)Fig.13 Surface streamlines comparisons between CFD and oil-flow test (α=60°)

    圖14 是攻角在α=60°下的背風(fēng)側(cè)旋渦流動形態(tài),圖中可以看到,在旋成體的前半段(約X/L=0.4 之前),其旋渦流動還較為對稱,隨著旋渦流動沿軸向向后發(fā)展,背風(fēng)側(cè)旋渦流動逐漸轉(zhuǎn)化為非對稱旋渦流動形態(tài).當(dāng)達(dá)到X/L=0.6 截面左右時,旋渦流動變的不再集中,逐漸發(fā)展為旋渦脫落的流動形態(tài),此時二次旋渦也變的不再穩(wěn)定,在物面上也觀察不到較為明顯的二次分離線.圖15 是各截面的壓力分布情況,從圖中也可以看到,在截面X/L=0.1~ 0.2,其壓力分布較為對稱.在X/L=0.4 位置,旋渦流動開始出現(xiàn)輕微的非對稱形態(tài),隨著旋渦向后發(fā)展,非對稱愈加明顯清晰.

    圖14 物面流線及空間流動(α=60°)Fig.14 Surface streamlines and flow (α=60°)

    圖15 各截面壓力分布(α=60°)Fig.15 Pressure distribution of cross sections (α=60°)

    圖16 是攻角α=90°工況下的流動圖像.圖中我們可以看到,由于此時旋成體軸線與來流成垂直關(guān)系,在背風(fēng)側(cè)的軸向流動速度很低,所以背風(fēng)側(cè)無法形成集中渦流動結(jié)構(gòu),在不同截面下基本都體現(xiàn)出了旋渦脫落的流動現(xiàn)象,并且具有較強(qiáng)的非對稱性.從壓力分布中(圖17)我們也可以看到,由于背風(fēng)側(cè)旋渦流動主要是渦脫落流動形態(tài),旋渦在物面附近形成的吸力較弱,沒有出現(xiàn)吸力峰值.

    圖16 物面流線及空間流動(α=90°)Fig.16 Surface streamlines and flow (α=90°)

    圖17 各截面壓力分布(α=90°)Fig.17 Pressure distribution of cross sections (α=90°)

    圖18 是攻角α=120°下油流試驗(yàn)物面流線與CFD 物面流線的對比.從圖中可以看到明顯的主分離線,CFD 計(jì)算結(jié)果與油流試驗(yàn)結(jié)果吻合很好,沒有看到比較明顯的二次分離線在模型底部附近均出現(xiàn)了主分離線“彎曲”的現(xiàn)象,這主要是旋成體底部誘導(dǎo)的分離流動的影響.圖18 還展示了旋成體在對稱面的流動特性,圖中我們看到,當(dāng)來流流過模型的底部時會背風(fēng)側(cè)產(chǎn)生一個回流區(qū),該回流區(qū)在軸向的影響范圍約為1.2D.

    圖18 CFD 物面流線與油流試驗(yàn)對比(α=120°)Fig.18 Surface streamlines comparisons between CFD and oil-flow test (α=120°)

    圖19 是α=120°攻角下空間流動與物面流線.圖中我們可以看到,當(dāng)自由來流在底部附近產(chǎn)生回流區(qū)后,隨著流動沿軸向向后發(fā)展,背風(fēng)側(cè)流動逐漸發(fā)展非對稱旋渦流動,并伴隨著一定的旋渦脫落現(xiàn)象.圖20 是不同截面下的壓力分布情況,可以看到在不同截面下壓力分布始終為非對稱.

    圖19 物面流線及空間流動(α=120°)Fig.19 Surface streamlines and flow (α=120°)

    圖20 各截面壓力分布(α=120°)Fig.20 Pressure distribution of cross sections (α=120°)

    圖21 是攻角150°下試驗(yàn)與數(shù)值模擬結(jié)果的對比情況.從圖中我們可以看到,隨著攻角進(jìn)一步增加,底部分離帶來的回流區(qū)在模型軸向的影響范圍將進(jìn)一步加大,影響范圍約為2.5D.隨著流動沿軸向向后發(fā)展,模型背風(fēng)側(cè)逐漸發(fā)展為旋渦流動,油流試驗(yàn)中的分離線逐漸清晰起來.圖22 是在攻角150°下的空間流動,從圖中我們可以看到,旋渦流動在回流區(qū)后開始逐漸發(fā)展起來,在流動剛流過回流區(qū)時,旋渦渦位較低,并且仍然受底部回流區(qū)的影響.隨著流動繼續(xù)向后發(fā)展,非對稱旋渦流動越來越清晰起來,受底部回流區(qū)的影響逐漸減弱.

    圖21 CFD 物面流線與油流試驗(yàn)對比(α=150°)Fig.21 Surface streamlines comparisons between CFD and oil-flow test (α=150°)

    圖22 物面流線及空間流動(α=150°)Fig.22 Surface streamlines and flow (α=150°)

    圖23 是攻角180°工況下物面流線與油流試驗(yàn)的對比情況,圖中一共分為上下兩張附圖,上面附圖為CFD 得到物面流線以及對稱面的渦量分布圖,下面附圖為油流試驗(yàn)照片,在此攻角下主要流動特性為底部分離引起的復(fù)雜回流流動.從試驗(yàn)以及數(shù)值模擬結(jié)果我們都可以看到,在底部附近會出現(xiàn)一個十分明顯的分離線,并且回流區(qū)再附位置主要在距離底部約2.5D 的位置,試驗(yàn)與計(jì)算結(jié)果對比較好.同時在對稱面渦量分布中我們可以看到,來流經(jīng)過底部時會產(chǎn)生較強(qiáng)的剪切層,隨著流動向后發(fā)展逐漸失穩(wěn),與回流區(qū)中的小尺度旋渦流動相互作用,形成較為復(fù)雜的具有較強(qiáng)非線性的局部流動.

    圖23 CFD 物面流線與油流試驗(yàn)對比(α=180°)Fig.23 Surface streamlines comparisons between CFD and oil-flow test (α=180°)

    4 超大攻角旋渦流動時序演化

    上一節(jié)對超大攻角瞬時流動進(jìn)行分析,但細(xì)長體在超大攻角下體現(xiàn)為復(fù)雜的非定常流動現(xiàn)象,其背風(fēng)側(cè)旋渦的強(qiáng)度、位置、渦脫落均會隨著時間產(chǎn)生變化,所以本節(jié)對超大攻角旋渦流動的時序演化進(jìn)行分析.本文每個計(jì)算物理時間步為Δt=2.3 ×10-5s,共計(jì)算2000 個物理時間步,總計(jì)算物理時間為0.046 s,以保證足夠的時長使得流動進(jìn)入穩(wěn)態(tài),并選取后1200 步進(jìn)行分析.

    圖24 是攻角α=60°下在X=350 mm 截面背風(fēng)側(cè)渦量圖隨時間的變化.從圖中我們可以看到,在該截面下背風(fēng)側(cè),呈現(xiàn)出了兩個集中渦以及兩個二次渦,集中渦的渦量和渦位隨著時間增加較小,呈現(xiàn)出輕微的非對稱流動形態(tài).隨著流動沿著軸向向后發(fā)展,當(dāng)流動發(fā)展至X=550 mm 截面,如圖25,集中渦體現(xiàn)出明顯的非對稱流動特性,同時隨著時間的變化,體現(xiàn)出旋渦脫落特性,產(chǎn)生類卡門渦街式的旋渦交替脫落現(xiàn)象.

    圖24 截面渦量隨時間變化(α=60°,X=350 mm)Fig.24 Vortice of cross section (α=60°,X=350 mm)

    圖25 截面渦量隨時間變化(α=60°,X=550 mm)Fig.25 Vortice of cross section (α=60°,X=350 mm)

    圖26 是攻角α=90°下在X=350 mm 截面下的渦量圖隨時間變化.從圖中我們可以看到,此時彈身垂直于來流,來流流過彈身后會產(chǎn)生較長的剪切層,隨著剪切層逐漸發(fā)展,慢慢卷起旋渦,旋渦強(qiáng)度相對較弱,并且旋渦距離物面較遠(yuǎn),旋渦對物面的誘導(dǎo)能力較低,沒有發(fā)現(xiàn)明顯的二次渦結(jié)構(gòu).彈身兩側(cè)剪切層誘導(dǎo)的旋渦呈現(xiàn)較為明顯的非對稱性,并且隨著時間增加,旋渦逐漸脫離剪切層向后脫落,左右旋渦的渦脫落現(xiàn)象交替進(jìn)行.

    圖26 截面渦量隨時間變化(α=90°,X=350 mm)Fig.26 Vortice of cross section (α=90°,X=350 mm)

    圖27 是攻角α=120°下的X=350 mm 截面渦量圖.從圖中我們可以看到,此時旋渦距離物面較近,出現(xiàn)了明顯的非對稱旋渦流動結(jié)構(gòu),同時在物面附近會誘導(dǎo)出較小尺度的二次渦結(jié)構(gòu),隨著時間發(fā)展,左右旋渦的渦位會出現(xiàn)切換現(xiàn)象,渦位較高的旋渦逐漸脫落并飄向尾跡區(qū).圖28 是Y=0 mm 的渦量分布圖,圖中可以看到,在底部附近的回流區(qū)存在較多的小尺度旋渦結(jié)構(gòu),隨著時間增加,距離物面較近小尺度渦在物面附近形成回流區(qū),距離物面較遠(yuǎn)的小尺度渦慢慢脫落飄向尾跡區(qū).彈身(圓柱等直段)附近也存在較強(qiáng)的渦量分布,這是圓柱段集中渦渦量在Y=0 mm 平面的投影,反應(yīng)出了超大攻角下旋渦具有復(fù)雜的三維效應(yīng).

    圖27 截面渦量隨時間變化(α=120°,X=350 mm)Fig.27 Vortice of cross section (α=120°,X=350 mm)

    圖28 截面渦量隨時間變化(α=120°,Y=0 mm)Fig.28 Vortice of cross section (α=120°,Y=0 mm)

    圖29 是攻角α=150°下的X=350 mm 截面渦量圖.圖中我們可以看到,旋成體背風(fēng)側(cè)產(chǎn)生了兩個集中渦,渦位較低,隨著時間增加,旋成體背風(fēng)側(cè)的渦位出現(xiàn)左右切換的現(xiàn)象,并且隨著時間的增加逐漸脫體向下游運(yùn)動.圖30 是Y=0 mm 截面渦量圖,來流流過底部后會出現(xiàn)較為明顯的剪切層,剪切層逐漸向后發(fā)展逐漸失穩(wěn)變成小尺度渦飄向下游,同時在底部回流區(qū)內(nèi)依然為很多小尺度旋渦流動.

    圖29 截面渦量隨時間變化(α=150°,X=350 mm)Fig.29 Vortice of cross section (α=150°,X=350 mm)

    圖30 截面渦量隨時間變化(α=150°,Y=0 mm)Fig.30 Vortice of cross section (α=150°,Y=0 mm)

    圖31 是攻角α=180°時在Y=0 mm 截面的渦量分布圖,從圖中我們可以看到,來流會在底部處拖出明顯的剪切層,剪切層逐漸向后發(fā)展逐漸失穩(wěn),一部分渦量會變成小尺度旋渦向下游發(fā)展,另一部分會再附至物面,在頭部形成一個回流區(qū),回流區(qū)內(nèi)的小尺度旋渦相互干擾作用,形成復(fù)雜的非定?;亓髁鲃咏Y(jié)構(gòu).

    圖31 截面渦量隨時間變化(α=180°,Y=0 mm)Fig.31 Vortice of cross section (α=180°,Y=0 mm)

    5 超大攻角物面壓力脈動特性

    為了更好的對細(xì)長體超大攻角下非定常流動頻率特性進(jìn)行分析,如圖32 所示,取兩個主要截面,一個是X=350 mm 截面,位于模型中間位置,能夠較好的描述背風(fēng)側(cè)由圓柱引起的旋渦流動;另一個是X=550 mm 截面,能夠較好的描述模型底部附近流動.在這兩個平面內(nèi),分別取了一個檢測點(diǎn)(m1 點(diǎn)和m2 點(diǎn)),通過提取檢測點(diǎn)的壓力脈動,可以更深入的分析非定常流動的頻域特性.

    圖32 監(jiān)測點(diǎn)及截面位置示意圖Fig.32 Location of monitored point and cross section

    圖33 是攻角α=60°下監(jiān)測點(diǎn)m2 的壓力脈動,在此截面流動已經(jīng)體現(xiàn)出了明顯的非定常特性,從測壓點(diǎn)的壓力脈動就可以看出來,通過對m2 點(diǎn)壓力脈動進(jìn)行頻譜分析,發(fā)現(xiàn)壓力脈動具有主頻,主頻為f=543 Hz.

    圖33 監(jiān)測點(diǎn)壓力脈動(α=60°,X=550 mm)Fig.33 Fluctuating pressure of monitored point (α=60°,X=550 mm)

    圖34 是兩個監(jiān)測點(diǎn)(m1 點(diǎn)和m2 點(diǎn))的壓力系數(shù)脈動數(shù)據(jù).從圖中可以看到,壓力脈動較為明顯,體現(xiàn)出了明顯的非定常性,通過壓力脈動進(jìn)行頻譜分析,發(fā)現(xiàn)在攻角α=90°下沒有發(fā)現(xiàn)明顯的主頻現(xiàn)象,這主要是由于該攻角下,旋渦距離物面較遠(yuǎn),對物面誘導(dǎo)能力較弱的原因.

    圖34 監(jiān)測點(diǎn)壓力脈動(α=90°)Fig.34 Fluctuating pressure of monitored point (α=90°)

    圖35 展示了監(jiān)測點(diǎn)m1 和m2 的壓力脈動情況,m1 點(diǎn)的壓力脈動體現(xiàn)出了圓柱段旋渦的切換情況,從頻譜分析中可以看到,該點(diǎn)的壓力脈動有一定的主頻,主頻為f=717 Hz 左右,但是該主頻并沒有特別集中,這里面也體現(xiàn)出了圓柱段的旋渦結(jié)構(gòu)不僅僅存在較大尺度的集中渦結(jié)構(gòu),同時也存在較小尺度的旋渦,小尺度渦與大尺度渦相互干擾影響,造成一定程度上的頻率分散.m2 點(diǎn)體現(xiàn)出了底部分離回流區(qū)的壓力脈動特性(如圖35(b)),可以看到在m2 點(diǎn)在較低頻率(1000 Hz 以內(nèi))沒有發(fā)現(xiàn)主頻,但是在高頻段出現(xiàn)了較為明顯的主頻,主頻為5391 Hz左右,這主要是在該截面的壓力脈動主要由底部分離回流區(qū)中的小尺度旋渦帶來,小尺度旋渦頻率較高,對應(yīng)壓力脈動頻率也較高.

    圖35 監(jiān)測點(diǎn)壓力脈動(α=120°)Fig.35 Fluctuating pressure of monitored point (α=120°)

    圖36 為α=150°監(jiān)測點(diǎn)m1 和m2 的壓力脈動情況,圖中可以看到,集中渦所產(chǎn)生的脈動存在主頻,主頻為543 Hz 左右.從監(jiān)測點(diǎn)m2 的頻譜分析來看,類似α=120°所發(fā)現(xiàn)的高頻脈動,在α=150°下的高頻脈動為5086 Hz 左右,但是此壓力脈動能量較低.

    圖36 監(jiān)測點(diǎn)壓力脈動(α=150°)Fig.36 Fluctuating pressure of monitored point (α=150°)

    圖37 是攻角α=180°下m2 點(diǎn)壓力脈動,用于監(jiān)測底部回流區(qū)的壓力脈動情況.從圖中可以看到,底部回流區(qū)引起了較高頻的壓力脈動,脈動頻率為f=5173 Hz.

    圖37 監(jiān)測點(diǎn)壓力脈動(α=180°)Fig.37 Fluctuating pressure of monitored point (α=180°)

    從以上分析可以看到,不同攻角下,超大攻角非定常流動體現(xiàn)出了不同的物面壓力脈動特性,主要有兩類主頻,一類為圓柱段非對稱旋渦誘導(dǎo)的較低頻率,另一類為攻角超過90°時,底部分離流動誘導(dǎo)的較高的頻率.本文將特征頻率進(jìn)行無量綱化,提取出脈動主頻的斯特勞哈爾數(shù)St

    式中f 為流動頻率,D 為直徑,V 為特征速度.對于非對稱旋渦誘導(dǎo)的頻率,類比圓柱擾流流動,特征速度取橫截面流動速度V=V∞sinα,對于底部分離流動誘導(dǎo)的流動頻率,特征速度取V=V∞.圖38 是不同攻角下的流動St,圖中我們可以看到,由非對稱渦引起的物面壓力脈動主要集中在攻角60°~ 150°,其中流動主頻St 范圍為St=0.19~ 0.33,而由底部分離誘導(dǎo)的非定常流動在120°~ 180°都有主頻,其脈動頻率較高,無量綱主頻范圍為St=1.55~ 1.64.

    圖38 不同攻角下的非定常流動頻率St 數(shù)Fig.38 St of unsteady flow versus angle of attack

    6 結(jié)論

    本文結(jié)合數(shù)值模擬DDES 方法與油流顯示試驗(yàn)方法,研究了細(xì)長旋成體在Ma=0.6 時超大攻角范圍內(nèi)(α=0°~ 180°)的流動特性,得到以下結(jié)論.

    (1)數(shù)值模擬與油流顯示試驗(yàn)獲得的物面流線吻合較好,在攻角范圍α=30°~ 60°和α=150°背風(fēng)側(cè)旋渦較為集中并且離物面較近時,均能捕捉到明顯的分離線、再附線等物面流動拓?fù)湫问?在攻角范圍α=90°~ 120°,背風(fēng)側(cè)旋渦距離物面較遠(yuǎn),物面流動拓?fù)湟砸淮畏蛛x線為主,當(dāng)攻角α=180°時,物面的分離線與再附線主要由底部分離回流區(qū)產(chǎn)生.

    (2)在攻角范圍α=0°~ 90°,細(xì)長體背風(fēng)側(cè)流動特性主要體現(xiàn)為由圓柱段產(chǎn)生的分離渦流動結(jié)構(gòu),當(dāng)攻角α=30°時,分離渦位對稱集中渦,攻角至α=60°時,細(xì)長體頭部附近依然為對稱渦,在尾部附近非定常非對稱渦結(jié)構(gòu),攻角至90°時,細(xì)長體背風(fēng)側(cè)均為非對稱旋渦結(jié)構(gòu),并伴有多尺度旋渦干擾、渦脫落等現(xiàn)象.

    (3)攻角α=120°~ 150°范圍內(nèi),細(xì)長體底部朝前,底部產(chǎn)生較為明顯的回流區(qū),回流區(qū)內(nèi)存在較多小尺度旋渦相互作用干擾,隨著流動逐漸沿軸向向后發(fā)展,流動逐漸發(fā)展為非對稱渦流動.當(dāng)攻角達(dá)到α=180°時,來流經(jīng)過底部時會產(chǎn)生較強(qiáng)的剪切層,隨著流動向后發(fā)展逐漸失穩(wěn),與回流區(qū)中的小尺度旋渦流動相互作用.

    (4)由非對稱渦引起的物面壓力脈動主要集中在攻角60°~ 150°,其中流動主頻St 范圍為St=0.19~0.33,而由底部分離誘導(dǎo)的非定常流動在120°~180°都有主頻,其脈動頻率較高,無量綱主頻范圍為St=1.55~ 1.64.

    本文此次研究旨在對超大攻角下細(xì)長體的流動特性有個較為清晰的認(rèn)識,下一步將對非定常氣動特性對機(jī)動飛行的影響及其控制策略做進(jìn)一步研究.

    猜你喜歡
    物面油流細(xì)長
    主變壓器油流繼電器指針頻繁抖動的原因分析
    寧夏電力(2022年4期)2022-11-10 04:13:30
    激波/湍流邊界層干擾壓力脈動特性數(shù)值研究1)
    帶擾動塊的細(xì)長旋成體背部繞流數(shù)值模擬
    正交車銑細(xì)長軸的切削穩(wěn)定性研究
    讓吸盤掛鉤更牢固
    脂肪流油流油 快瘦快瘦“脂肪炸彈”“炸出”財(cái)富一片片
    牽引變壓器繞組溫升與油流的關(guān)聯(lián)性
    新型單面陣自由曲面光學(xué)測量方法成像特性仿真
    彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究
    圓筒內(nèi)有接頭的細(xì)長桿穩(wěn)定性問題
    精品国产亚洲在线| 国产成人免费观看mmmm| 国产在线精品亚洲第一网站| 黑人巨大精品欧美一区二区mp4| 建设人人有责人人尽责人人享有的| 久久国产精品男人的天堂亚洲| 免费高清在线观看日韩| 中文字幕人妻熟女乱码| 丰满迷人的少妇在线观看| 一级毛片高清免费大全| 亚洲全国av大片| 久久精品aⅴ一区二区三区四区| 亚洲九九香蕉| 最新在线观看一区二区三区| 黑人操中国人逼视频| 亚洲国产欧美日韩在线播放| 一二三四社区在线视频社区8| 中国美女看黄片| 777米奇影视久久| 亚洲av成人一区二区三| 国产高清视频在线播放一区| 美女高潮到喷水免费观看| 如日韩欧美国产精品一区二区三区| 夫妻午夜视频| 国产精品免费视频内射| 国产不卡一卡二| 美女高潮到喷水免费观看| 欧美乱妇无乱码| 在线播放国产精品三级| 妹子高潮喷水视频| 捣出白浆h1v1| 99香蕉大伊视频| 精品久久久久久久久久免费视频 | 一边摸一边做爽爽视频免费| 日韩欧美一区二区三区在线观看 | 黑人操中国人逼视频| 国产精品综合久久久久久久免费 | 女人被躁到高潮嗷嗷叫费观| 国产精品九九99| 欧美成狂野欧美在线观看| 搡老乐熟女国产| 777久久人妻少妇嫩草av网站| 国产一区在线观看成人免费| 国产精品久久久久久人妻精品电影| 国产97色在线日韩免费| 欧美日韩瑟瑟在线播放| 在线观看免费视频网站a站| 亚洲一码二码三码区别大吗| 自拍欧美九色日韩亚洲蝌蚪91| 好看av亚洲va欧美ⅴa在| 美女国产高潮福利片在线看| 黄色毛片三级朝国网站| 99精国产麻豆久久婷婷| 免费高清在线观看日韩| cao死你这个sao货| 久久精品国产亚洲av香蕉五月 | 国产一区二区激情短视频| 国产一区二区激情短视频| 欧美激情高清一区二区三区| 亚洲五月天丁香| 亚洲熟女精品中文字幕| 人妻久久中文字幕网| 精品国产乱子伦一区二区三区| 看黄色毛片网站| www日本在线高清视频| 黑人巨大精品欧美一区二区mp4| 韩国精品一区二区三区| 免费少妇av软件| 美女高潮喷水抽搐中文字幕| 亚洲熟妇熟女久久| 久久人妻熟女aⅴ| 亚洲五月婷婷丁香| 精品熟女少妇八av免费久了| 一级作爱视频免费观看| 日本黄色视频三级网站网址 | 国产一区二区三区综合在线观看| 精品免费久久久久久久清纯 | 老熟妇仑乱视频hdxx| 久久国产乱子伦精品免费另类| 国产精品二区激情视频| 欧美日韩乱码在线| 欧美亚洲 丝袜 人妻 在线| 一级作爱视频免费观看| 成熟少妇高潮喷水视频| 窝窝影院91人妻| 麻豆乱淫一区二区| 飞空精品影院首页| 久久久久久人人人人人| 搡老岳熟女国产| 亚洲av美国av| 久久中文字幕一级| 久久精品成人免费网站| 欧美日韩一级在线毛片| 中文字幕色久视频| 人人澡人人妻人| 久久精品亚洲精品国产色婷小说| 国产亚洲精品久久久久久毛片 | 久久这里只有精品19| 国产片内射在线| 在线观看舔阴道视频| 在线观看www视频免费| 午夜福利视频在线观看免费| av线在线观看网站| 高清黄色对白视频在线免费看| 黑人猛操日本美女一级片| 手机成人av网站| 午夜福利乱码中文字幕| 桃红色精品国产亚洲av| 中文字幕高清在线视频| 久久狼人影院| 午夜精品在线福利| 亚洲片人在线观看| 中国美女看黄片| 夜夜躁狠狠躁天天躁| 久久久久久免费高清国产稀缺| 亚洲熟妇中文字幕五十中出 | 国产精品电影一区二区三区 | 国产av又大| 免费观看人在逋| 香蕉久久夜色| 国产免费男女视频| 亚洲国产欧美日韩在线播放| 一a级毛片在线观看| 国产精品九九99| 久久精品成人免费网站| 18禁裸乳无遮挡动漫免费视频| 亚洲午夜理论影院| 一边摸一边做爽爽视频免费| 国产野战对白在线观看| 捣出白浆h1v1| 夜夜爽天天搞| av片东京热男人的天堂| 久久人人爽av亚洲精品天堂| 又黄又粗又硬又大视频| 丝袜在线中文字幕| 国产主播在线观看一区二区| 国产精品自产拍在线观看55亚洲 | 亚洲国产欧美日韩在线播放| 午夜福利免费观看在线| 久久久国产一区二区| 下体分泌物呈黄色| 美女高潮到喷水免费观看| 久久久国产成人精品二区 | 久久亚洲精品不卡| 麻豆国产av国片精品| 久久国产亚洲av麻豆专区| av中文乱码字幕在线| 18禁美女被吸乳视频| www.精华液| 亚洲欧美日韩另类电影网站| 成年人免费黄色播放视频| 黄色片一级片一级黄色片| 无限看片的www在线观看| 一a级毛片在线观看| 日韩大码丰满熟妇| 亚洲一区二区三区不卡视频| 欧美激情 高清一区二区三区| 久久久国产成人免费| 精品国产国语对白av| 亚洲第一av免费看| 麻豆国产av国片精品| 在线观看舔阴道视频| 午夜福利欧美成人| 99精品久久久久人妻精品| 黑人操中国人逼视频| 国产主播在线观看一区二区| 无人区码免费观看不卡| 欧美国产精品va在线观看不卡| 天堂俺去俺来也www色官网| 亚洲熟妇中文字幕五十中出 | 国产97色在线日韩免费| 天堂俺去俺来也www色官网| 女人久久www免费人成看片| 最新在线观看一区二区三区| 女同久久另类99精品国产91| 亚洲性夜色夜夜综合| 老司机福利观看| 国产亚洲一区二区精品| 免费女性裸体啪啪无遮挡网站| 黄色丝袜av网址大全| 国产成人欧美| 久久人妻熟女aⅴ| 在线播放国产精品三级| 午夜免费成人在线视频| 夜夜爽天天搞| 在线观看免费高清a一片| 亚洲美女黄片视频| 久久精品亚洲av国产电影网| 久久久久久久久免费视频了| a级毛片黄视频| 身体一侧抽搐| 日本一区二区免费在线视频| 三级毛片av免费| 交换朋友夫妻互换小说| 黄网站色视频无遮挡免费观看| 女性被躁到高潮视频| 乱人伦中国视频| 国产欧美亚洲国产| 欧美色视频一区免费| 中文字幕最新亚洲高清| 曰老女人黄片| 亚洲精品久久成人aⅴ小说| 老鸭窝网址在线观看| 久久久久视频综合| 中文字幕av电影在线播放| 亚洲精品国产区一区二| 日韩欧美一区视频在线观看| 黄频高清免费视频| 极品人妻少妇av视频| 91精品国产国语对白视频| 一边摸一边做爽爽视频免费| 久热这里只有精品99| 香蕉丝袜av| 国产极品粉嫩免费观看在线| av国产精品久久久久影院| 美女视频免费永久观看网站| 国产成人系列免费观看| 精品一区二区三区视频在线观看免费 | 国产一卡二卡三卡精品| 满18在线观看网站| 黄色片一级片一级黄色片| www.熟女人妻精品国产| www日本在线高清视频| 80岁老熟妇乱子伦牲交| 一级,二级,三级黄色视频| 两个人看的免费小视频| 手机成人av网站| 一进一出好大好爽视频| 久久性视频一级片| 国产精品免费大片| 国产视频一区二区在线看| 亚洲精品美女久久久久99蜜臀| 亚洲人成电影免费在线| 国产免费男女视频| 国产免费av片在线观看野外av| 最新在线观看一区二区三区| 99久久国产精品久久久| 成人精品一区二区免费| 午夜精品国产一区二区电影| 成年人免费黄色播放视频| 欧美在线黄色| 亚洲欧美一区二区三区久久| 97人妻天天添夜夜摸| 制服人妻中文乱码| 亚洲中文字幕日韩| 美女高潮喷水抽搐中文字幕| 99热国产这里只有精品6| av电影中文网址| 久久人人97超碰香蕉20202| 在线免费观看的www视频| 国产极品粉嫩免费观看在线| 中文字幕高清在线视频| av中文乱码字幕在线| 亚洲av日韩在线播放| 法律面前人人平等表现在哪些方面| 一边摸一边做爽爽视频免费| 91精品三级在线观看| 国产欧美日韩综合在线一区二区| 丝袜美足系列| 成人av一区二区三区在线看| 精品久久久久久久毛片微露脸| 可以免费在线观看a视频的电影网站| 精品福利观看| av视频免费观看在线观看| 交换朋友夫妻互换小说| 制服人妻中文乱码| 成人av一区二区三区在线看| 极品教师在线免费播放| 高清av免费在线| 在线观看66精品国产| 国产午夜精品久久久久久| 99香蕉大伊视频| 高潮久久久久久久久久久不卡| 欧美成人午夜精品| 欧美黑人欧美精品刺激| 国产欧美日韩一区二区精品| 热re99久久国产66热| 亚洲少妇的诱惑av| 日韩欧美免费精品| av天堂久久9| 日韩中文字幕欧美一区二区| 久久中文看片网| 国产av又大| 看黄色毛片网站| 新久久久久国产一级毛片| 老汉色av国产亚洲站长工具| 欧美乱妇无乱码| 亚洲av成人一区二区三| 国产一区有黄有色的免费视频| 桃红色精品国产亚洲av| 国产1区2区3区精品| 国产一区在线观看成人免费| 国产成人影院久久av| 国产激情欧美一区二区| 欧美精品高潮呻吟av久久| av超薄肉色丝袜交足视频| 中文字幕av电影在线播放| 啦啦啦免费观看视频1| 成人手机av| 美国免费a级毛片| 国产一区有黄有色的免费视频| 国产伦人伦偷精品视频| 99久久综合精品五月天人人| 好看av亚洲va欧美ⅴa在| av中文乱码字幕在线| 国产成人精品在线电影| 国产97色在线日韩免费| 亚洲欧美一区二区三区久久| 乱人伦中国视频| 高清黄色对白视频在线免费看| 国产精品免费大片| 又紧又爽又黄一区二区| 成年人免费黄色播放视频| 免费在线观看亚洲国产| 天天躁日日躁夜夜躁夜夜| 午夜福利视频在线观看免费| 亚洲国产精品合色在线| 亚洲午夜精品一区,二区,三区| 久久久国产欧美日韩av| 18禁裸乳无遮挡免费网站照片 | 成人特级黄色片久久久久久久| a级片在线免费高清观看视频| 成年人黄色毛片网站| 色播在线永久视频| 国产精品影院久久| 久久99一区二区三区| 久久香蕉激情| 国产激情欧美一区二区| 国产男靠女视频免费网站| 欧美激情 高清一区二区三区| 一边摸一边抽搐一进一小说 | 国产不卡一卡二| 99精国产麻豆久久婷婷| av网站在线播放免费| 亚洲欧美激情在线| 日韩三级视频一区二区三区| 色婷婷av一区二区三区视频| 欧美精品啪啪一区二区三区| 1024香蕉在线观看| 久99久视频精品免费| 亚洲欧美一区二区三区久久| 久久国产精品男人的天堂亚洲| 国产成人影院久久av| av一本久久久久| 亚洲精品中文字幕在线视频| 老汉色av国产亚洲站长工具| 色婷婷av一区二区三区视频| 久久精品亚洲熟妇少妇任你| 午夜福利欧美成人| 狠狠婷婷综合久久久久久88av| 亚洲国产欧美一区二区综合| 久久精品aⅴ一区二区三区四区| 国产精品久久久久久人妻精品电影| 久久午夜综合久久蜜桃| 韩国精品一区二区三区| 99国产精品一区二区蜜桃av | 少妇粗大呻吟视频| 久久人妻熟女aⅴ| 国产高清videossex| 亚洲欧美日韩另类电影网站| 欧美日本中文国产一区发布| 欧美日韩福利视频一区二区| av不卡在线播放| 国产精品综合久久久久久久免费 | 99热国产这里只有精品6| 日韩免费高清中文字幕av| 18禁观看日本| cao死你这个sao货| 久久热在线av| 在线视频色国产色| 99香蕉大伊视频| 亚洲午夜精品一区,二区,三区| 午夜福利,免费看| 午夜久久久在线观看| 久久人人97超碰香蕉20202| 亚洲精华国产精华精| 两个人看的免费小视频| 国产野战对白在线观看| 最近最新免费中文字幕在线| 精品卡一卡二卡四卡免费| 美女午夜性视频免费| 男女下面插进去视频免费观看| 亚洲av成人一区二区三| 国产精品香港三级国产av潘金莲| 精品一区二区三区四区五区乱码| 国产精品乱码一区二三区的特点 | 女同久久另类99精品国产91| 亚洲国产毛片av蜜桃av| 欧美在线一区亚洲| 亚洲 欧美一区二区三区| 亚洲成人国产一区在线观看| 91麻豆av在线| 国产精品国产高清国产av | 咕卡用的链子| 村上凉子中文字幕在线| 成人黄色视频免费在线看| 老司机在亚洲福利影院| 国产伦人伦偷精品视频| 国产高清视频在线播放一区| a在线观看视频网站| 色婷婷久久久亚洲欧美| 午夜福利免费观看在线| 两个人看的免费小视频| 欧美精品高潮呻吟av久久| bbb黄色大片| 老汉色∧v一级毛片| 老鸭窝网址在线观看| 成人18禁高潮啪啪吃奶动态图| 最新美女视频免费是黄的| 久久久国产欧美日韩av| 黄色片一级片一级黄色片| 国产精品99久久99久久久不卡| 色94色欧美一区二区| 国产激情欧美一区二区| 一区二区日韩欧美中文字幕| 久久久精品国产亚洲av高清涩受| 国产成+人综合+亚洲专区| 视频区欧美日本亚洲| 9色porny在线观看| 他把我摸到了高潮在线观看| 少妇的丰满在线观看| 嫩草影视91久久| 久久久久国内视频| 国产精品国产高清国产av | 国产aⅴ精品一区二区三区波| 亚洲成人免费av在线播放| 最近最新免费中文字幕在线| 91在线观看av| 在线观看免费视频网站a站| 亚洲成国产人片在线观看| 麻豆乱淫一区二区| 久久久久久久久久久久大奶| 久久国产精品人妻蜜桃| 欧美激情 高清一区二区三区| 久久人妻av系列| 人妻一区二区av| 国产国语露脸激情在线看| 国产亚洲欧美精品永久| 欧美成人午夜精品| 午夜福利,免费看| 国产一区二区激情短视频| 精品一区二区三区四区五区乱码| 成熟少妇高潮喷水视频| 欧美黑人欧美精品刺激| 精品福利永久在线观看| 黄色a级毛片大全视频| 一级a爱片免费观看的视频| 99re在线观看精品视频| 一边摸一边抽搐一进一小说 | 亚洲少妇的诱惑av| 男人操女人黄网站| 精品高清国产在线一区| 久久精品国产亚洲av香蕉五月 | 精品国产超薄肉色丝袜足j| 久久香蕉国产精品| 老司机福利观看| 亚洲三区欧美一区| 国产av精品麻豆| 99精品在免费线老司机午夜| 国产在线观看jvid| 精品久久蜜臀av无| 极品少妇高潮喷水抽搐| 欧美成狂野欧美在线观看| 久久国产精品影院| 天天添夜夜摸| 成人18禁高潮啪啪吃奶动态图| 天天影视国产精品| 高清毛片免费观看视频网站 | videosex国产| 欧美久久黑人一区二区| 搡老熟女国产l中国老女人| 免费在线观看完整版高清| 精品久久蜜臀av无| 国产亚洲精品第一综合不卡| 在线观看日韩欧美| 欧美国产精品va在线观看不卡| 国内毛片毛片毛片毛片毛片| 久久99一区二区三区| 色尼玛亚洲综合影院| av中文乱码字幕在线| 日本一区二区免费在线视频| 亚洲专区字幕在线| 一进一出好大好爽视频| 精品无人区乱码1区二区| 国产成人av激情在线播放| 久久香蕉精品热| 天堂中文最新版在线下载| 丁香六月欧美| 午夜福利乱码中文字幕| 久久久久视频综合| 成人18禁在线播放| 不卡一级毛片| 电影成人av| 日韩免费高清中文字幕av| 人人妻人人澡人人爽人人夜夜| 人妻久久中文字幕网| 真人做人爱边吃奶动态| 天堂动漫精品| 免费久久久久久久精品成人欧美视频| 乱人伦中国视频| 18禁裸乳无遮挡免费网站照片 | 精品国产超薄肉色丝袜足j| 国产三级黄色录像| 女人被狂操c到高潮| 精品亚洲成a人片在线观看| 中文欧美无线码| 婷婷成人精品国产| 在线观看免费日韩欧美大片| 不卡一级毛片| 最新的欧美精品一区二区| 久久香蕉精品热| 天堂中文最新版在线下载| 欧美日韩瑟瑟在线播放| 成在线人永久免费视频| 看免费av毛片| 黄色视频不卡| 黄片播放在线免费| 一本一本久久a久久精品综合妖精| 777米奇影视久久| ponron亚洲| 日本黄色日本黄色录像| 正在播放国产对白刺激| 91麻豆精品激情在线观看国产 | 国产精品成人在线| 无限看片的www在线观看| 国内毛片毛片毛片毛片毛片| 黑人欧美特级aaaaaa片| 欧美黑人欧美精品刺激| 在线观看日韩欧美| 妹子高潮喷水视频| 麻豆av在线久日| 日日爽夜夜爽网站| 亚洲精品一二三| 校园春色视频在线观看| 老熟女久久久| 日本wwww免费看| 亚洲一区中文字幕在线| 中文字幕人妻丝袜制服| 国产99久久九九免费精品| 国产精品九九99| 黄色怎么调成土黄色| www.999成人在线观看| 中出人妻视频一区二区| 久久狼人影院| 亚洲熟女精品中文字幕| 免费在线观看日本一区| 91精品三级在线观看| 亚洲精品国产一区二区精华液| 欧美精品高潮呻吟av久久| 午夜老司机福利片| 自线自在国产av| 久久精品91无色码中文字幕| 757午夜福利合集在线观看| 99re6热这里在线精品视频| 国产男女内射视频| 亚洲欧美色中文字幕在线| 中亚洲国语对白在线视频| 人人妻人人添人人爽欧美一区卜| 一进一出抽搐gif免费好疼 | 日韩欧美一区视频在线观看| 欧美色视频一区免费| 精品少妇一区二区三区视频日本电影| 天堂中文最新版在线下载| 国产在线精品亚洲第一网站| 另类亚洲欧美激情| 久久香蕉激情| 999久久久国产精品视频| 国产色视频综合| av线在线观看网站| x7x7x7水蜜桃| 国产不卡一卡二| 国产区一区二久久| 亚洲少妇的诱惑av| 热99re8久久精品国产| 欧美av亚洲av综合av国产av| 中出人妻视频一区二区| 亚洲av成人不卡在线观看播放网| 久久久精品国产亚洲av高清涩受| 成人av一区二区三区在线看| 国产极品粉嫩免费观看在线| 精品国产超薄肉色丝袜足j| 香蕉国产在线看| 亚洲中文字幕日韩| 午夜亚洲福利在线播放| 50天的宝宝边吃奶边哭怎么回事| tube8黄色片| 国产精品美女特级片免费视频播放器 | tube8黄色片| 国产成人av教育| 村上凉子中文字幕在线| 精品国产一区二区久久| 曰老女人黄片| 中国美女看黄片| 国产高清视频在线播放一区| 青草久久国产| 人妻丰满熟妇av一区二区三区 | 亚洲精华国产精华精| netflix在线观看网站| avwww免费| 夫妻午夜视频| videos熟女内射| 国产亚洲欧美98| 亚洲,欧美精品.| 黄色毛片三级朝国网站| 18禁裸乳无遮挡免费网站照片 | 免费在线观看完整版高清| 日日爽夜夜爽网站| 99riav亚洲国产免费| 女人被躁到高潮嗷嗷叫费观| 伦理电影免费视频| 免费女性裸体啪啪无遮挡网站| 九色亚洲精品在线播放| 中文欧美无线码| 好男人电影高清在线观看| 国产精品一区二区免费欧美|