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

    動(dòng)態(tài)艙門對空腔復(fù)雜流動(dòng)和噪聲特性影響分析

    2022-07-13 01:55:02李曉東姜蔚婉蔡晉生
    關(guān)鍵詞:艙門空腔脈動(dòng)

    李曉東,姜蔚婉,2,屈 崑,*,蔡晉生

    (1. 西北工業(yè)大學(xué) 航空學(xué)院,西安 710072;2. 中國航空工業(yè)集團(tuán) 成都飛機(jī)設(shè)計(jì)研究所,成都 610091)

    0 引 言

    內(nèi)埋式武器在投放分離過程中,彈艙將演變?yōu)椋◣撻T)大尺度空腔。大尺度空腔是航空航天武器裝備必不可少且普遍采用的一種典型結(jié)構(gòu)形式,主要用于作戰(zhàn)武器裝載、動(dòng)力燃油存儲和起落架存放等,具有尺度大、結(jié)構(gòu)復(fù)雜等特點(diǎn)。在動(dòng)載荷激勵(lì)下,空腔薄壁結(jié)構(gòu)容易產(chǎn)生嚴(yán)重的振動(dòng)和噪聲,不僅影響武器裝備的作戰(zhàn)效能,甚至危及武器裝備安全。因此,空腔繞流是空氣動(dòng)力學(xué)領(lǐng)域的一個(gè)重要且典型問題。

    針對空腔中流動(dòng)和噪聲的復(fù)雜現(xiàn)象與物理機(jī)理,國內(nèi)外科研工作者通過風(fēng)洞試驗(yàn)和計(jì)算手段對干凈彈艙(空腔)模型[1-4]進(jìn)行了廣泛深入的研究。具體包括分析流致振蕩現(xiàn)象的成因[5-6],對空腔類型[7]、來流馬赫數(shù)、邊界層厚度等[3,8]參數(shù)進(jìn)行詳細(xì)研究,并探討了多種不同的控制措施[2,9],對空腔設(shè)計(jì)具有重要的指導(dǎo)意義。Rossiter[10]提出了預(yù)測空腔噪聲模態(tài)頻率的半經(jīng)驗(yàn)?zāi)P停籋eller 等[11]改進(jìn)了此模型,將其推廣到高馬赫數(shù)(Ma= 2、3)的可壓縮流動(dòng)中。

    空腔是彈艙的一種理想模型,實(shí)際中艙門的存在會使得空腔內(nèi)非定常流動(dòng)更加復(fù)雜,因而研究學(xué)者們對帶艙門的空腔流動(dòng)開展了研究和分析。首先是對帶有靜態(tài)艙門的空腔流動(dòng)的研究:Blair 和Stallings[12]采用風(fēng)洞試驗(yàn)研究了不同艙門開度和艙門高度對導(dǎo)彈類物體在超聲速空腔中投放的影響,結(jié)果表明艙門開度對彈體的氣動(dòng)載荷有顯著影響,并且降低艙門高度會減少彈體的氣動(dòng)載荷;Lawson 和Barakos[13]采用脫體渦模擬(DES)方法對艙門開啟90°的空腔標(biāo)模M219 和1303 無人機(jī)模型進(jìn)行了模擬,并分析了艙門打開對彈體拋投的影響。Casper 等[14]采用試驗(yàn)的方式模擬了帶有艙門的復(fù)雜空腔中彈體拋投的流固耦合問題,發(fā)現(xiàn)增加空腔復(fù)雜度會引發(fā)彈體展向共振。其次在動(dòng)態(tài)艙門的研究上,Sheta 等[15]采用RANS/LES混合方法,模擬了超聲速(Ma= 1.44)空腔艙門從5°打開到35°過程中的流動(dòng),發(fā)現(xiàn)在艙門運(yùn)動(dòng)到30°時(shí),剪切層和艙內(nèi)聲波產(chǎn)生較強(qiáng)的相互作用,腔內(nèi)壁面出現(xiàn)強(qiáng)的壓力脈動(dòng);Loupy 等[16]對比分析了固定開度艙門和動(dòng)態(tài)艙門對馬赫數(shù)為0.85 的空腔流動(dòng)、載荷和氣動(dòng)噪聲的影響,發(fā)現(xiàn)艙門打開過程中的過渡階段,艙門氣動(dòng)載荷增大,流動(dòng)的脈動(dòng)和總聲壓級增強(qiáng),對空腔結(jié)構(gòu)產(chǎn)生顯著影響。

    國內(nèi)多位學(xué)者采用RANS 方法模擬分析靜態(tài)艙門[17]和艙門開啟過程中[18]空腔動(dòng)態(tài)氣動(dòng)特性和載荷的變化。吳繼飛等[19-20]對Ma= 0.6 空腔艙門的開閉問題開展了研究,設(shè)計(jì)了不同的運(yùn)動(dòng)方式,發(fā)現(xiàn)艙門運(yùn)動(dòng)時(shí)氣動(dòng)載荷發(fā)生明顯變化,艙門打開30°時(shí)氣動(dòng)載荷最大,并分析了腔內(nèi)和艙門的總聲壓級和聲壓譜特性。由于帶艙門的空腔流動(dòng)非定常效應(yīng)顯著,且流動(dòng)的湍流結(jié)構(gòu)復(fù)雜,所以傳統(tǒng)的RANS方法捕捉復(fù)雜流動(dòng)結(jié)構(gòu)比較困難。閆盼盼等[21]采用基于k?ω的分離渦模擬方法(IDDES)計(jì)算了艙門不同開啟姿態(tài)對內(nèi)埋武器投放的影響,發(fā)現(xiàn)艙門會改變超聲速流場中的波系結(jié)構(gòu),對彈體偏航影響較大。

    考慮運(yùn)動(dòng)艙門時(shí),空腔中氣動(dòng)與噪聲特性會變得更加復(fù)雜,預(yù)測難度增加,而國內(nèi)對帶運(yùn)動(dòng)艙門的復(fù)雜空腔流動(dòng)與噪聲的研究較少。為了研究艙門運(yùn)動(dòng)情況下的空腔復(fù)雜流動(dòng)現(xiàn)象,深入分析靜、動(dòng)態(tài)艙門對空腔復(fù)雜流動(dòng)的影響,本文采用嵌套重疊網(wǎng)格方法對艙門和空腔分別生成網(wǎng)格, 采用改進(jìn)的脫體渦模擬方法對亞聲速(Ma= 0.6)空腔流場進(jìn)行數(shù)值模擬,從流場特性、湍流結(jié)構(gòu)、脈動(dòng)特性上分析艙門對湍流流場和腔內(nèi)噪聲的影響。

    1 數(shù)值方法

    針對空腔流動(dòng)強(qiáng)脈動(dòng)、強(qiáng)湍流的特點(diǎn),采用RANS/LES 混合方法進(jìn)行計(jì)算,本文采用改進(jìn)的脫體渦湍流模擬方法(Improved Delayed DES,IDDES)[22],解析空腔內(nèi)分離區(qū)域的湍流結(jié)構(gòu)和脈動(dòng)特性?;谝陨系臄?shù)值方法,本文發(fā)展了一套基于有限體積方法的多塊結(jié)構(gòu)網(wǎng)格求解器,并通過多種復(fù)雜流動(dòng)問題驗(yàn)證其準(zhǔn)確性和可靠性[23-26]。

    1.1 高精度低耗散格式

    式中,

    1.2 空腔標(biāo)模(M219)驗(yàn)證

    為驗(yàn)證高精度數(shù)值格式在空腔問題中的有效性,采用本文提出的格式對馬赫數(shù)0.85 空腔標(biāo)模M219進(jìn)行模擬分析,并與參考文獻(xiàn)[32]的LES 計(jì)算結(jié)果進(jìn)行對比。網(wǎng)格單元總數(shù)量約800 萬,空腔內(nèi)部與附近的網(wǎng)格數(shù)量約380 萬。對于馬赫數(shù)為0.85 的流動(dòng),時(shí)間步長為 Δt=2×10?5s, 約為Tref=L/U∞的百分之一,足以捕捉第一階Rossiter 頻率。圖1 給出了來流馬赫數(shù)為0.85 時(shí),中截面上時(shí)均流向速度、縱向速度和湍動(dòng)能的分布,結(jié)果與參考文獻(xiàn)符合得較好。由此可見,本文采用的計(jì)算方法可以比較準(zhǔn)確地捕捉空腔復(fù)雜湍流的時(shí)均統(tǒng)計(jì)特性。

    圖1 空腔標(biāo)模M219 在Ma = 0.85 時(shí)均速度與湍動(dòng)能分布Fig. 1 Wall-normal profiles of the mean velocity and turbulence kinetic energy in canonical cavity M219 at Ma = 0.85

    2 動(dòng)態(tài)嵌套重疊網(wǎng)格方法

    針對運(yùn)動(dòng)物體的網(wǎng)格生成問題,發(fā)展了快速的動(dòng)態(tài)嵌套重疊網(wǎng)格方法[33]。首先通過采用分層嵌套重疊網(wǎng)格策略(如圖2)對網(wǎng)格進(jìn)行劃分,然后采用自適應(yīng)樹結(jié)構(gòu)完成網(wǎng)格切割和貢獻(xiàn)單元搜索。引入壁面距離概念,設(shè)計(jì)了統(tǒng)一的分層嵌套重疊策略,同時(shí)實(shí)現(xiàn)了重疊網(wǎng)格切割和嵌套網(wǎng)格切割兩個(gè)功能,完成重疊最小化的功能,使得重疊網(wǎng)格組裝質(zhì)量得到了改善。

    圖2 分層嵌套重疊網(wǎng)格策略Fig. 2 Hierarchically embedded strategy of overset grids

    根據(jù)運(yùn)動(dòng)規(guī)律的不同,可將所有運(yùn)動(dòng)物體分為若干動(dòng)態(tài)組,而所有靜態(tài)物體歸為唯一的靜態(tài)組。各組分別構(gòu)造獨(dú)立的二叉樹和八叉樹輔助網(wǎng)格,并且每個(gè)組內(nèi)的流場網(wǎng)格和二/八叉樹網(wǎng)格的坐標(biāo)均保持不變。即,各個(gè)動(dòng)態(tài)分組仍然用各自原有的局部坐標(biāo)系,但在每個(gè)時(shí)間步上記錄下各動(dòng)態(tài)組運(yùn)動(dòng)的位移。而在更新切割信息時(shí),將所有被切割的坐標(biāo)點(diǎn)變換到切割物體所在的局部坐標(biāo)內(nèi),從而完成運(yùn)動(dòng)問題的網(wǎng)格切割。對于同一組內(nèi)的切割,由于相對位置不變,只需進(jìn)行一次切割,后續(xù)不需再改變。此方法使得二叉樹和八叉樹在運(yùn)動(dòng)過程中不必重新構(gòu)造,只需進(jìn)行坐標(biāo)變換操作,大大節(jié)省了計(jì)算時(shí)間。表1 給出了艙門運(yùn)動(dòng)時(shí)網(wǎng)格組裝時(shí)間的對比。對于總計(jì)2 100 多萬的網(wǎng)格單元來說,兩種動(dòng)態(tài)嵌套重疊網(wǎng)格組裝方法均可滿足該類非定常流動(dòng)計(jì)算的效率要求。

    表1 考慮運(yùn)動(dòng)艙門的空腔動(dòng)態(tài)重疊網(wǎng)格組裝的耗時(shí)統(tǒng)計(jì)Table 1 Assembly time of overset grids for a cavity with moving doors

    3 空腔標(biāo)模C201 與艙門模型

    3.1 空腔標(biāo)模C201 與艙門模型

    空腔標(biāo)模C201 由中國空氣動(dòng)力研究與發(fā)展中心設(shè)計(jì),空腔底部、前后壁分別布置了靜壓測點(diǎn)和脈動(dòng)壓力測點(diǎn)(見圖3)。本文采用的計(jì)算模型如圖4(a)所示,與試驗(yàn)的模型外形一致??涨婚L度(L)為200 mm,深度(D)為33.3 mm,寬度(W)為66.7 mm,空腔長深比L/D為6,寬深比W/D為2;與空腔相連壁面的長度和寬度分別是Wa和La。

    圖3 空腔標(biāo)模C201 中截面的壓力測點(diǎn)Fig. 3 Distribution of pressure probes in the central section of canonical cavity C201

    參照文獻(xiàn)[1]的數(shù)據(jù)選取艙門,艙門厚度1.5 mm,艙門四周采用弧形物面進(jìn)行過渡。此構(gòu)型與空腔標(biāo)模M219 的艙門外形有一定相似性。在艙門與空腔、艙門與艙門之間保留1 mm(約為艙門深度的3%)縫隙,這樣可以保證縫隙處保留足夠的網(wǎng)格用于數(shù)值模擬。艙門三維模型如圖4(b)所示。

    圖4 空腔標(biāo)模C201 與艙門計(jì)算模型Fig. 4 Geometry of canonical cavity C201 and doors

    3.2 計(jì)算網(wǎng)格

    整個(gè)模型網(wǎng)格采用一個(gè)網(wǎng)格層,包含三個(gè)網(wǎng)格簇:空腔網(wǎng)格簇(其中包含背景網(wǎng)格)和兩個(gè)艙門網(wǎng)格簇。圖5 給出了空腔的網(wǎng)格,由于網(wǎng)格在展向(z方向)是對稱的,所以圖中顯示了一半網(wǎng)格。邊界條件設(shè)置如下:在圖5 給出的模型中,紅色為無滑移物面邊界條件;綠色表示的是空腔前后的自由端,采用對稱邊界條件;展向的兩個(gè)邊界均采用對稱邊界條件;最上面的外場邊界設(shè)置為入口/出口邊界;前后外場分別為入口/出口邊界??涨痪W(wǎng)格采用結(jié)構(gòu)網(wǎng)格,共計(jì)約1 600 萬網(wǎng)格單元。

    圖5 空腔標(biāo)模C201 計(jì)算模型與網(wǎng)格(每4 個(gè)點(diǎn)顯示一個(gè)網(wǎng)格點(diǎn))Fig. 5 Computational domain and grids for canonical cavity C201

    艙門采用沿著物面法向?qū)⑽锩婢W(wǎng)格向外推進(jìn)的方法生成貼體網(wǎng)格,然后在外圍采用矩形外場,使得外場網(wǎng)格的周向尺寸得到控制。單個(gè)艙門的網(wǎng)格單元共計(jì)約454 萬,圖6(a)給出了艙門的物面網(wǎng)格,艙y+=1門的前向剖視圖如圖6(b)所示。艙門壁面的第一層網(wǎng)格厚度由 給出。艙門開度的定義如圖6(c)所示。

    圖6 艙門計(jì)算網(wǎng)格與艙門打開角度(開度)定義Fig. 6 Computational grids for cavity doors and the definition of the door-opening angle

    剪切層流動(dòng)是空腔復(fù)雜流場的重要因素之一,所以需要保留高質(zhì)量的空腔剪切層網(wǎng)格。通過改進(jìn)空腔壁面距離的計(jì)算方法,使得在四種不同艙門打開角度下均能保留空腔剪切層網(wǎng)格(見圖7),從而保證剪切層流動(dòng)計(jì)算的精細(xì)和準(zhǔn)確。

    圖7 四種不同艙門開度在 x/L=0.5處yOz 視圖上網(wǎng)格的組裝結(jié)果Fig. 7 The grid slice of assembled overset grids in yOz planes for cavities with different door-opening angles

    4 空腔標(biāo)模C201 數(shù)值模擬

    對無艙門的空腔標(biāo)模C201 亞聲速流動(dòng)(Ma=0.6)進(jìn)行模擬,來流條件為T∞=286.66 K,Re=12.14×106。 時(shí) 間步 長 采用 Δt=1.0×10?6s,約為 參 考 時(shí) 間Tref=L/U∞的千分之一,足以分辨流動(dòng)的細(xì)節(jié)。采用雙時(shí)間步推進(jìn)求解非定常流場,每步偽時(shí)間推進(jìn)采用多重網(wǎng)格加速收斂,并采用MPI 并行技術(shù)提高計(jì)算效率。非定常計(jì)算采用IDDES 模擬方法,在廣州天河二號超級計(jì)算機(jī)上進(jìn)行。采用SA 湍流模型計(jì)算的穩(wěn)態(tài)流場作為初場,在非定常流場充分發(fā)展之后(約15Tref~20Tref),開始非定常流場的數(shù)據(jù)采集和時(shí)均場的累計(jì)。時(shí)均場的累計(jì)大約為 50Tref,約為0.05 s,足以分辨第一階Rossiter 頻率。

    4.1 流場時(shí)均特性分析

    圖13(a)給出了中截面時(shí)均流場的壓強(qiáng)系數(shù)分布和流線形狀,在空腔后部存在一個(gè)較小的第二環(huán)流。在較低速來流條件下,剪切層厚度增長較大,大范圍沖擊到后壁、角落和底部。沖擊后壁與角落的氣流較為靠近剪切層核心區(qū),流速較高,形成局部高壓環(huán)流;而沖擊底部的氣流較為遠(yuǎn)離剪切層核心區(qū),流速較低,被角落的高壓環(huán)流所阻礙,提前回流形成了主環(huán)流。

    圖8 對比了三種不同網(wǎng)格尺度下腔內(nèi)底部壓力系數(shù)的分布與試驗(yàn)值的比較。隨著網(wǎng)格的加密,計(jì)算結(jié)果趨于穩(wěn)定。盡管其值高于試驗(yàn)結(jié)果,但是計(jì)算結(jié)果可以較好地預(yù)測壓力系數(shù)的變化。從圖中可以看出,底部壓強(qiáng)分布表現(xiàn)為前段平緩下降、中段增長和后段快速上升,最后氣流沖擊后壁形成高壓區(qū)域。

    圖8 三種網(wǎng)格對腔內(nèi)底部壓力系數(shù)時(shí)均分布的影響Fig. 8 Grid convergence of the mean pressure coefficients at cavity floor compared to experimental distribution

    4.2 總聲壓級

    總聲壓級(overall sound pressure level,OASPL)可以通過湍流模擬得到的壓力脈動(dòng)來獲得,其定義如下:

    圖9 空腔標(biāo)模C201 腔內(nèi)底部總聲壓級Fig. 9 OASPL distribution at the floor of canonical cavity C201 without doors

    圖10 空腔標(biāo)模C201 底部聲壓級成分合成示意圖Fig. 10 Sketch of the development of OASPL in canonical cavity C201

    4.3 脈動(dòng)壓力功率譜

    通過脈動(dòng)壓力功率譜分析空腔底部的頻譜特性,對脈動(dòng)壓力p′采用快速傅里葉分析得到功率譜密度(PSD), 然后將PSD 用SPL[13]呈現(xiàn),進(jìn)行功率譜的分析。圖11 給出了空腔標(biāo)模C201 腔內(nèi)底部第3 和第11 監(jiān)測點(diǎn)脈動(dòng)壓強(qiáng)功率譜,圖中灰色條帶是根據(jù)Heller 公式計(jì)算的峰值頻率帶。計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合,主導(dǎo)峰值頻率都是第二峰值,因此在工程上需要特別注意第二峰值頻率的影響。

    圖11 空腔標(biāo)模C201 腔內(nèi)底部p3 和p11 測點(diǎn)脈動(dòng)壓強(qiáng)功率譜Fig. 11 Comparison of SPL of fluctuating pressure at probes p3 and p11 between numerical simulation and experiment

    5 靜態(tài)艙門空腔模擬

    在滿足準(zhǔn)定常假設(shè)時(shí),可以對靜態(tài)艙門直接模擬分析。采用發(fā)展的嵌套重疊網(wǎng)格方法和IDDES 湍流模擬方法,對艙門開度分別為30°、60°、90°和120°的四種工況進(jìn)行準(zhǔn)定常流動(dòng)模擬,分析其流動(dòng)時(shí)均特性和流場動(dòng)態(tài)特性。從時(shí)均流場(中截面時(shí)均流場和底部壓力分布)、流場的動(dòng)態(tài)特性(總聲壓級、脈動(dòng)壓強(qiáng)功率譜)和瞬時(shí)場湍流結(jié)構(gòu)三個(gè)方面分析艙門對空腔流場產(chǎn)生的影響。本節(jié)采用的數(shù)值方法、計(jì)算方式與干凈空腔模擬方法一致。

    5.1 時(shí)均場特性分析

    將空腔底部的靜壓測點(diǎn)時(shí)均壓力與干凈空腔的結(jié)果進(jìn)行對比分析,圖12 給出了不同狀態(tài)下底部壓力系數(shù)的分布。從腔內(nèi)底部壓力分布來看,艙門打開30°時(shí),底部壓力系數(shù)的分布與干凈空腔的差別較大,前部壓強(qiáng)有所升高,后部壓強(qiáng)明顯降低。在艙門打開60°、90°、120°時(shí),腔內(nèi)前半部分底部壓力系數(shù)基本與干凈空腔的計(jì)算結(jié)果一致,后半部分壓力系數(shù)隨著艙門打開角度的增加逐漸趨于干凈空腔。開啟角度為90°和120°時(shí),時(shí)均壓強(qiáng)系數(shù)分布與干凈空腔基本一致。

    圖12 四種艙門開度(30°、60°、90°、120°)對空腔底部壓力系數(shù)的影響Fig. 12 Impacts of the door-opening angle on mean pressure coefficients at the cavity floor

    圖13 給出了干凈空腔和不同艙門開度的中截面時(shí)均流場壓力系數(shù)云圖和流線分布情況。從圖中可以看出在艙門打開角度為30°時(shí),空腔前緣附近的平板上,壓力系數(shù)會升高,這是由于艙門的阻礙使得流動(dòng)減速造成的。這種壓強(qiáng)升高作用也使得空腔前部的負(fù)壓區(qū)域的壓強(qiáng)系數(shù)相對干凈空腔要偏高。同時(shí),對于30°開啟角度,外流與空腔附近的流動(dòng)受到艙門的阻擋,只能通過中截面附近的開口相互作用。特別是阻礙外流對剪切層的能量供給,使得后壁附近的壓強(qiáng)相對干凈空腔明顯降低,導(dǎo)致30°時(shí),空腔底面后部壓強(qiáng)系數(shù)水平整體降低。此外,后壁附近的高壓也受到艙門阻礙,在縱向方向只能在中截面開口處泄壓,后部流線有更加明顯地向上彎曲。

    圖13 干凈空腔和四種艙門開度下中截面平均流場壓力系數(shù)與流線分布Fig. 13 Impacts of door-opening angles on the streamlines and mean pressure contours in the middle section of cavity C201

    艙門開度為60°、90°和120°時(shí),空腔前緣平板同一位置的壓力系數(shù)逐漸降低,且后壁高壓區(qū)和后壁拐角處低壓區(qū)的范圍較艙門開度為30°時(shí)逐漸變大,并與干凈空腔的結(jié)果趨于一致。這與“艙門打開角度越大,對流動(dòng)的限制明顯越小”的物理原因是一致的。

    另外,相比于干凈空腔,隨著艙門打開角的不斷增大,主渦位置沿流向不斷地后移并靠近干凈空腔主渦位置。艙門打開120°時(shí),湍流結(jié)構(gòu)基本上與干凈空腔的一致,但是在前壁附近的流場結(jié)構(gòu)與干凈空腔有所區(qū)別。由此可知,隨著空腔打開角度的增加,腔內(nèi)時(shí)均流場的形態(tài)、渦的結(jié)構(gòu)與干凈空腔的結(jié)果逐漸趨于一致。

    5.2 總聲壓級

    圖14 給出了四種艙門開度下空腔底部總聲壓級的分布。從圖中可以看出,艙門開度30°時(shí),底部各個(gè)位置的總聲壓級較干凈空腔均明顯下降,但變化趨勢與干凈空腔的結(jié)果基本一致。艙門打開90°時(shí),總聲壓級與120°的變化趨勢一致,前者的值略高1~2 dB。艙門開度為60°時(shí),總聲壓級腔內(nèi)前半部分的分布特點(diǎn)與90°、120°的結(jié)果基本一致;而在腔內(nèi)后半部分,其總聲壓級明顯有所改變。

    本文認(rèn)為導(dǎo)致噪聲升高的原因在于艙門阻礙了腔內(nèi)噪聲向外傳播。從圖14 可以看出,艙門分別開啟60°、90°、120°時(shí),三種條件下的聲壓級分布幾乎一致,因此可知噪聲升高與艙門開啟角度無關(guān),只與艙門的存在有關(guān)。艙門的存在使得空腔內(nèi)部的脈動(dòng)能量無法以聲波的形式沿展向傳播,只能沿著兩艙門開口方向傳播,第二階模態(tài)的能量無法得到有效釋放,這樣就導(dǎo)致腔內(nèi)噪聲強(qiáng)度升高。本文中艙門開度為90°時(shí),總聲壓級的計(jì)算結(jié)果與空腔標(biāo)模M219 在馬赫數(shù)0.85 下的試驗(yàn)結(jié)果[1]分布特點(diǎn)相似,呈現(xiàn)出“W”形態(tài)的分布。

    圖14 干凈彈艙與四種艙門開度下的彈艙底部總聲壓級分布Fig. 14 Impacts of door-opening angles on OASPL at the cavity floor

    5.3 脈動(dòng)壓強(qiáng)功率譜

    空腔中的噪聲主要由寬頻白噪聲和窄頻噪聲(Rossiter 模態(tài))兩部分組成[13]。白噪聲是由流動(dòng)結(jié)構(gòu)產(chǎn)生,其來源為自由來流的脈動(dòng)、剪切層和湍流等。第二部分也稱為Rossiter 模態(tài),是不同流場結(jié)構(gòu)直接相互作用產(chǎn)生的,其特點(diǎn)為幅值大且集中在某幾個(gè)頻率。圖15 給出了四種不同艙門開啟角度下,腔內(nèi)9 個(gè)脈動(dòng)測點(diǎn)的脈動(dòng)壓強(qiáng)功率譜曲線,柱狀條帶是Heller 公式計(jì)算的峰值頻率。從圖中可以看出:

    圖15 不同艙門開度下腔內(nèi)9 個(gè)脈動(dòng)測點(diǎn)脈動(dòng)壓強(qiáng)功率譜Fig. 15 SPL of fluctuating pressure at nine probes on cavity floors for cases with different door-opening angles

    1)艙門開度為30°時(shí),二階Rossiter 所占優(yōu)勢顯著減弱,并且該模態(tài)頻率偏高,其他高頻(第三、第四)模態(tài)已經(jīng)不再明顯。這是由于艙門打開較小時(shí),空腔流動(dòng)中剪切層與回傳壓縮波的相互作用受到限制,所以流場結(jié)構(gòu)相互作用的強(qiáng)度減弱,致使聲壓幅值降低。

    2)60°、90°和120°艙門開度下,各脈動(dòng)壓力測點(diǎn)的Rossiter 模態(tài)明顯,且模態(tài)頻率與經(jīng)驗(yàn)公式基本一致。這說明此時(shí)流場中的結(jié)構(gòu)與干凈空腔一致。二階Rossiter 模態(tài)占主導(dǎo),幅值明顯比其他模態(tài)高10 dB以上。

    3)圖中虛線、實(shí)線分別表示腔內(nèi)前、后半部分采樣點(diǎn)的結(jié)果。從圖中可以看出,四種不同艙門開啟角度下,后半部分的脈動(dòng)強(qiáng)度明顯高于前半部分。

    為研究艙門對腔內(nèi)脈動(dòng)特性的影響,進(jìn)一步選取后壁(p11)的脈動(dòng)壓強(qiáng)測點(diǎn)上的脈動(dòng)壓強(qiáng)功率譜進(jìn)行比較分析,如圖16 所示。各圖中黑色曲線為干凈空腔試驗(yàn)曲線,褐色曲線為干凈空腔計(jì)算曲線,其他彩色曲線為各角度情況下的計(jì)算曲線。

    圖16 不同艙門開度腔內(nèi)后壁測點(diǎn)p11 脈動(dòng)壓強(qiáng)功率譜Fig. 16 SPL of fluctuating pressure at probe p11 on cavity floors for cases with different door-opening angles

    對于后壁測點(diǎn)p11 來說,當(dāng)開啟角為30°時(shí),功率水平明顯低于干凈空腔,其主峰頻率比干凈空腔的主峰頻率偏高。這可能是由于小角度開啟情況下,空腔構(gòu)型更接近封閉腔體,其展向或者縱向的共振成分加強(qiáng)。而在其他開啟角度情況下,變化基本一致,峰值頻率大致和干凈空腔一致,但是第二模態(tài)峰值都顯著升高,三、四模態(tài)峰值沒有升高甚至幅值降低。這也印證了第二階模態(tài)被強(qiáng)化的觀點(diǎn)。同時(shí)這也說明,當(dāng)打開角度較大時(shí),在各處,艙門對主模態(tài)的影響是一致的,是通過改變流場的主結(jié)構(gòu)來影響噪聲的產(chǎn)生。在前壁和腔內(nèi)底部測點(diǎn)的脈動(dòng)壓強(qiáng)功率譜中,出現(xiàn)了相似的變化和影響,不在此贅述。

    整體來看,出現(xiàn)不同的變化的原因是,在艙門打開小角度時(shí),空腔內(nèi)的法向(y向)和展向(z向)流動(dòng)均受到限制,流場脈動(dòng)會受到抑制,噪聲也會減小。艙門打開角度較大(60°以后)時(shí),艙門的限制主要在展向,對法向的限制已經(jīng)不明顯;但是此時(shí)艙門的存在會增加空腔流動(dòng)的復(fù)雜性,阻斷空腔近壁面噪聲沿展向向外傳播,從而會增強(qiáng)腔內(nèi)的噪聲。

    5.4 湍流結(jié)構(gòu)

    為了進(jìn)一步分析不同艙門開度下空腔流動(dòng)的特點(diǎn),通過Q等值面比較瞬時(shí)場中湍流結(jié)構(gòu)的異同。正的Q值區(qū)域表示渦張量比應(yīng)變率張量大的流動(dòng)區(qū)域,是與旋渦有關(guān)的部分,是常被用作判斷旋渦區(qū)域的標(biāo)準(zhǔn)[34]。圖17 給出了干凈空腔和四種不同艙門開度下,瞬時(shí)流場中的Q= 3 000 等值面,采用來流方向的無量綱速度u進(jìn)行染色。

    圖17 瞬時(shí)流場的Q = 3 000 等值面(用u 進(jìn) 行染色)Fig. 17 Instantaneous iso-surfaces (extracted by Q = 3 000) colored by u for clean cavity and the cavity with stationary doors

    從瞬時(shí)圖中可以看出,由于艙門的存在,空腔內(nèi)的流動(dòng)會受到不同程度的影響,具體表現(xiàn)為:

    1)30°的湍流結(jié)構(gòu)明顯受到艙門的限制,從而使得流場的湍流結(jié)構(gòu)較干凈空腔要小,艙內(nèi)渦結(jié)構(gòu)很稀疏,尤其是空腔前部區(qū)域。

    2)在艙門打開60°時(shí),渦結(jié)構(gòu)變得密集,可以明顯看到艙內(nèi)的流動(dòng)結(jié)構(gòu)沿法向(y向)向外運(yùn)動(dòng),但是空腔前部仍然比較稀疏。

    3)到90°和120°時(shí),渦結(jié)構(gòu)的豐富程度變得相似。在艙門開度為120°時(shí),艙門仍然會限制空腔展向的流動(dòng)。這從側(cè)面反映了艙門的存在會阻礙空腔周圍噪聲的傳播,從而使得腔內(nèi)噪聲一定程度地增大。

    4)相比干凈空腔的瞬時(shí)Q準(zhǔn)則等值面圖,可以看出,有艙門時(shí)空腔上方的剪切層渦結(jié)構(gòu)運(yùn)動(dòng)高度更高,這仍然是剪切層在展向受艙門限制的結(jié)果。而隨著艙門開度的不斷增大,展向的空腔流動(dòng)會受到不同程度的限制。

    5)流動(dòng)經(jīng)過前緣平板后形成的剪切層,在流場結(jié)構(gòu)中會出現(xiàn)大尺度的剪切層渦結(jié)構(gòu)。艙門對腔內(nèi)前半部分湍流結(jié)構(gòu)的相互作用較弱;對后腔內(nèi)半部分不同流動(dòng)結(jié)構(gòu)之間相互作用的影響較大,從而影響空腔流動(dòng)的脈動(dòng)和噪聲特性。

    6 運(yùn)動(dòng)艙門空腔模擬

    為進(jìn)一步模擬真實(shí)艙門打開過程對空腔流動(dòng)和聲場的影響,本節(jié)采用前文的空腔和艙門,在考慮艙門運(yùn)動(dòng)的情況下,對空腔的湍流流動(dòng)進(jìn)行模擬。對考慮運(yùn)動(dòng)的非穩(wěn)態(tài)非定常流動(dòng)問題,采用經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)方法提取非穩(wěn)定特征,定量分析其噪聲特性的變化。

    6.1 艙門運(yùn)動(dòng)設(shè)計(jì)

    艙門運(yùn)動(dòng)規(guī)律的選取由條件相似準(zhǔn)則確定,而斯特勞哈爾數(shù)(St)是非定常流動(dòng)中須考慮的相似準(zhǔn)則,物理上代表非定常運(yùn)動(dòng)與定常運(yùn)動(dòng)的慣性力之比,其定義為:

    式中,LD為艙門長度, ΔT為艙門每次開啟的總時(shí)間,V∞表示來流速度。戰(zhàn)斗機(jī)武器艙艙門開啟過程中,需1 ~2 s 打開100°[35]。為使得艙門運(yùn)動(dòng)規(guī)律對流動(dòng)的影響與真實(shí)飛行情況一致,需要保證St數(shù)相同。因此,在來流條件相同的情況下,當(dāng)模型縮比為1∶10 時(shí),在縮比模型模擬過程中,艙門需要0.1~0.2s左右完成一次動(dòng)作。本文選取0.1 s,艙門運(yùn)動(dòng)的角速度為900 rad/s;從艙門完全關(guān)閉到打開120°,模擬的時(shí)間步長取 Δt=2×10?5。

    6.2 非穩(wěn)態(tài)非定常分析方法

    首先對艙門關(guān)閉狀態(tài)進(jìn)行非定常計(jì)算,保證流動(dòng)充分發(fā)展;然后艙門開始運(yùn)動(dòng),進(jìn)行完整的非穩(wěn)態(tài)非定常流動(dòng)模擬。圖18 給出了艙門從0°打開到100°過程中,空腔內(nèi)動(dòng)態(tài)測點(diǎn)的壓力系數(shù)變化。從圖中可以看出,整個(gè)流動(dòng)大致分為三個(gè)階段:初期的小幅穩(wěn)定階段(≤30°)、過渡階段(30°~60°)、充分發(fā)展階段(≥60°)。同時(shí)空腔前半?yún)^(qū)域各動(dòng)態(tài)采集點(diǎn)(p3~p7)的壓力系數(shù)呈對稱擺動(dòng)趨勢,而空腔后半?yún)^(qū)域動(dòng)態(tài)采集點(diǎn)的壓力系數(shù)呈增長趨勢。以p11 測點(diǎn)為例,圖19展示了其壓力值的變化情況。

    圖18 艙門打開過程中各動(dòng)態(tài)采集點(diǎn)的壓力系數(shù)變化Fig. 18 Variations of pressure coefficients at different probes during the opening of cavity doors

    圖19 艙門開啟過程中彈艙后壁測點(diǎn)p11 處壓力變化與EMD 模態(tài)曲線Fig. 19 Variation of pressure and its EMD modes at probe p11 during the opening of cavity doors

    由于艙門運(yùn)動(dòng)使得空腔內(nèi)流動(dòng)呈現(xiàn)非穩(wěn)態(tài)的變化,即壓力脈動(dòng)(p′)需要使用局部平均值而不是簡單的代數(shù)平均。為了準(zhǔn)確描述艙門運(yùn)動(dòng)過程的非穩(wěn)態(tài)變化,本文提出采用希爾伯特-黃轉(zhuǎn)換(Hilbert-Huang Transform) 中 的 經(jīng) 驗(yàn) 模 態(tài) 分 解( Empirical Mode Decomposition,EMD)進(jìn)行分析。該方法的核心在于尋找當(dāng)前數(shù)據(jù)中的所有局部極大值以及局部極小值,然后分別采用三次樣條構(gòu)建極大值包絡(luò)線和極小值包絡(luò)線,兩者的平均值即為EMD 模態(tài);再對扣除EMD 模態(tài)的剩余數(shù)據(jù)進(jìn)行分析,直到剩余數(shù)據(jù)為單調(diào)函數(shù)為止。

    為證明該方法的有效性,采用艙門打開90°的準(zhǔn)定常流動(dòng)進(jìn)行分析,并用代數(shù)平均的計(jì)算結(jié)果進(jìn)行對比驗(yàn)證。以測壓點(diǎn)p11 為例,圖20 給出了不同EMD 模態(tài)的分布曲線。分析結(jié)果表明,低階EMD 模態(tài)分析包含較多高頻信息,與原始數(shù)據(jù)分布具有一致性,無法直接用于捕捉非穩(wěn)態(tài)的信息;高階的EMD 模態(tài)多數(shù)與低頻流場信息有關(guān),可用于描述艙門運(yùn)動(dòng)引起的非穩(wěn)態(tài)流場的變化。從圖20 中,EMD5 和EMD6模態(tài)的分布差異可以觀察到此特征。而更高階的EMD 模態(tài)(如EMD9)基本與代數(shù)平均值一致。從而說明EMD 方法具有與代數(shù)平均計(jì)算方式相同的功能。

    圖20 艙門打開90°時(shí)p11 點(diǎn)壓力變化與不同EMD 模態(tài)曲線Fig. 20 Time histories of pressure and its EMD modes at probe p11 in the cavity with vertical doors

    圖21 中給出了采用不同EMD 模態(tài)作為時(shí)均值時(shí),p11 點(diǎn)總聲壓級的變化。結(jié)果表明,采用低階EMD 模態(tài)所得的OASPL 遠(yuǎn)低于采用代數(shù)平均值方式計(jì)算的值。而提高采用的EMD 模態(tài),OASPL 會逐漸收斂于代數(shù)平均計(jì)算的值(相差在0.5%以內(nèi))。這說明采用EMD 模態(tài)進(jìn)行OASPL 計(jì)算是有效可行的。盡管EMD 模態(tài)的選取無統(tǒng)一標(biāo)準(zhǔn),但可以采用OASPL 的收斂性判斷選取的EMD 模態(tài)。

    圖21 艙門打開90°時(shí)考慮不同EMD 模態(tài)的p11 點(diǎn)OASPL 計(jì)算結(jié)果Fig. 21 OASPL computed using different EMD modes at probe p11 in the cavity with vertical doors

    圖22 給出了EMD 模態(tài)對p11 測點(diǎn)功率譜的影響,主要特點(diǎn)為低頻區(qū)域有差別,高頻區(qū)域無明顯差異。由于EMD 本身就是對高頻濾波的結(jié)果,不同的EMD 曲線代表不同的低頻模態(tài),所以會對低階頻率產(chǎn)生影響。計(jì)算結(jié)果中的典型Rossiter 頻率與采用平均值計(jì)算的方式一致,可見采用EMD 模態(tài)對壓力功率譜的計(jì)算方式是可行的。

    6.3 運(yùn)動(dòng)艙門的噪聲分析

    采用經(jīng)驗(yàn)?zāi)B(tài)分解方法,對亞聲速空腔艙門開啟過程中的非穩(wěn)態(tài)動(dòng)態(tài)特性進(jìn)行分析。圖23 給出了空腔后壁p11 測點(diǎn)的OASPL 隨著選取不同EMD 模態(tài)所產(chǎn)生的變化。當(dāng)取到EMD5 模態(tài)時(shí),計(jì)算得到的OASPL 開始趨于穩(wěn)定,收斂值為艙門運(yùn)動(dòng)過程中的OASPL。值得注意的是,該值低于靜態(tài)艙門的計(jì)算值。這是因?yàn)?,在艙門的開啟過程中,當(dāng)艙門打開角度小時(shí),流動(dòng)的脈動(dòng)較小,艙門打開后期的脈動(dòng)變大(見圖19),空腔內(nèi)流動(dòng)與聲場逐漸產(chǎn)生相互作用并增強(qiáng),進(jìn)而產(chǎn)生了流動(dòng)的振蕩現(xiàn)象。圖24 給出了采用艙門開啟到90°附近的數(shù)據(jù)進(jìn)行分段動(dòng)態(tài)分析的OASPL 計(jì)算結(jié)果。此時(shí)的結(jié)果與靜態(tài)艙門計(jì)算值相當(dāng),說明對靜態(tài)艙門打開90°的準(zhǔn)定常計(jì)算跟真實(shí)運(yùn)動(dòng)過程中的狀態(tài)吻合,這對工程問題的分析具有指導(dǎo)意義。

    圖23 艙門開啟過程中不同EMD 模態(tài)的p11 點(diǎn)OASPL 計(jì)算結(jié)果Fig. 23 OASPL computed using different EMD modes at probe p11 during the opening of cavity doors

    圖24 艙門開啟到90°時(shí)不同EMD 模態(tài)的p11 點(diǎn)的OASPL 分析結(jié)果Fig. 24 OASPL computed using different EMD modes at probe p11 when cavity doors are opened to 90°

    圖25 給出了采用不同EMD 模態(tài)時(shí),對應(yīng)的脈動(dòng)壓力功率譜分析結(jié)果。艙門運(yùn)動(dòng)過程中,二階Rossiter 模態(tài)仍為主導(dǎo)模態(tài),幅值明顯高于其他模態(tài)10 dB 以上,但低于靜態(tài)艙門打開90°的分析結(jié)果,此處與總聲壓級降低的原因一致。一階Rossiter 模態(tài)強(qiáng)度降低,而三、四階Rossiter 模態(tài)略有升高。

    圖25 艙門開啟過程中不同EMD 模態(tài)的p11 點(diǎn)的功率譜比較Fig. 25 SPL computed using different EMD modes at probe p11 during the opening of cavity doors

    采用EMD 方法對空腔內(nèi)全部脈動(dòng)測點(diǎn)進(jìn)行OASPL 分析,如圖26 所示。計(jì)算結(jié)果表明,艙門開啟過程中,OASPL 分布狀態(tài)仍呈現(xiàn)“W”形狀,但其值比靜態(tài)大開度艙門整體降低5 dB 左右。這是由于艙門開啟過程中,存在由小幅脈動(dòng)到強(qiáng)脈動(dòng)的發(fā)展過程。小角度的艙門開啟過程中,流場的脈動(dòng)較小,還沒有完全發(fā)展;強(qiáng)脈動(dòng)的完全發(fā)展需要一定的時(shí)間。在艙門的整個(gè)運(yùn)動(dòng)過程中流動(dòng)呈現(xiàn)過渡發(fā)展的特點(diǎn),因此出現(xiàn)OASPL 的降低。

    圖26 艙門開啟過程中空腔底部總聲壓級分布Fig. 26 OASPL at cavity floor during the opening of cavity doors

    6.4 流動(dòng)分析

    圖27 給出了艙門開啟過程中,艙門打開到不同角度時(shí),物面壓力系數(shù)的分布情況。從圖中可以看出,壓力分布具有一定對稱性,兩側(cè)艙門承受的壓力相當(dāng)。在艙門打開小角度時(shí),艙門限制了空腔內(nèi)部與外部流場的摻混與交換,壓力變化較小。隨著艙門的逐漸打開,空腔與外流場發(fā)生混合,空腔內(nèi)壓力變化增大。艙門打開角度較大(60°以后)時(shí),空腔內(nèi)流動(dòng)出現(xiàn)大幅度振蕩現(xiàn)象,流場與聲場出現(xiàn)強(qiáng)相互作用,從而造成流場脈動(dòng)增強(qiáng),噪聲特性變得顯著。與前文OASPL 分析和功率譜分析結(jié)果一致。

    圖27 艙門開啟過程中物面壓力系數(shù)變化Fig. 27 Instantaneous wall pressure coefficients during the opening of cavity doors

    7 結(jié) 論

    針對考慮艙門的多尺度、強(qiáng)湍流、高雷諾數(shù)、強(qiáng)脈動(dòng)的空腔流動(dòng)問題,采用本文發(fā)展的嵌套重疊網(wǎng)格方法對帶運(yùn)動(dòng)艙門的空腔生成網(wǎng)格,用IDDES 湍流模擬方法對亞聲速(Ma= 0.6)空腔流動(dòng)進(jìn)行模擬,分析了靜、動(dòng)態(tài)艙門對空腔復(fù)雜流動(dòng)和噪聲特性的影響。具體結(jié)論如下:

    1)發(fā)展的動(dòng)態(tài)嵌套重疊網(wǎng)格方法,可自動(dòng)高效生成大幅度運(yùn)動(dòng)物體的計(jì)算網(wǎng)格,用于運(yùn)動(dòng)艙門的空腔復(fù)雜流場模擬。

    2)采用四階對稱重構(gòu)和三階MUSCL 混合的格式,在不明顯增加計(jì)算量的前提下,明顯降低了數(shù)值耗散,提升了計(jì)算的準(zhǔn)確度,可以捕捉空腔中的湍流結(jié)構(gòu)。對標(biāo)模M219 的計(jì)算結(jié)果與文獻(xiàn)給的結(jié)果吻合。對空腔標(biāo)模C201 腔內(nèi)噪聲的預(yù)測和對脈動(dòng)壓力功率譜的分析結(jié)果與試驗(yàn)數(shù)據(jù)一致。

    3)亞聲速空腔標(biāo)模C201 腔內(nèi)底部壓強(qiáng)分布表現(xiàn)出三段分布:前段平緩負(fù)壓區(qū);中段壓力增長過渡區(qū);后段高度增加區(qū),形成后壁的高壓區(qū)域。腔內(nèi)主導(dǎo)峰值頻率為第二峰值。

    4)與干凈空腔相比,當(dāng)艙門為靜態(tài)時(shí),艙門小開度(30°)時(shí),時(shí)均流場中截面上腔內(nèi)的流態(tài)變得更加復(fù)雜,并且腔內(nèi)底部壓力系數(shù)的值均會明顯減小。艙門開度較大(60°以后)時(shí),時(shí)均流場中截面的腔內(nèi)渦結(jié)構(gòu)與干凈空腔一致,底部壓力分布趨于干凈空腔的計(jì)算結(jié)果;但是艙門的存在,會使得底部后半部分的壓力較干凈空腔出現(xiàn)一定程度上的降低。艙門開度為90°時(shí),腔內(nèi)總聲壓級較干凈空腔有所升高,最大升高值近9 dB,呈現(xiàn)出“W”形態(tài)的分布;二階Rossiter模態(tài)占主導(dǎo)。艙門與空腔后半部分的不同流動(dòng)結(jié)構(gòu)之間相互作用,從而對空腔流動(dòng)的脈動(dòng)和噪聲特性產(chǎn)生影響。

    5)針對艙門運(yùn)動(dòng)的非穩(wěn)態(tài)非定常特點(diǎn),本文采用經(jīng)驗(yàn)?zāi)B(tài)分解方法對噪聲特性進(jìn)行分析,計(jì)算結(jié)果表明該方法有效可行。艙門開啟過程中,整體的OASPL有所下降,主導(dǎo)頻率依然為二階Rossiter 模態(tài)。對于小開度的艙門,靜態(tài)艙門模擬與運(yùn)動(dòng)艙門結(jié)果相差較大,因而此類情況需要考慮運(yùn)動(dòng)艙門的影響。在艙門開度較大(90°)時(shí),兩者的分析結(jié)論一致,說明可采用準(zhǔn)定常來分析此狀態(tài)的流場與噪聲特性,指導(dǎo)內(nèi)埋彈艙的設(shè)計(jì)。

    猜你喜歡
    艙門空腔脈動(dòng)
    新學(xué)期,如何“脈動(dòng)回來”?
    家教世界(2023年25期)2023-10-09 02:11:56
    RBI在超期服役脈動(dòng)真空滅菌器定檢中的應(yīng)用
    基于邊光滑有限元法的二維復(fù)合彈性空腔聲振特性分析
    飛機(jī)艙門失效乘客減載計(jì)算方法的探討
    運(yùn)輸機(jī)尾艙門收放液壓控制系統(tǒng)的改進(jìn)設(shè)計(jì)
    基于虛擬鉸鏈打開機(jī)構(gòu)的艙門提升機(jī)構(gòu)研究
    民用飛機(jī)復(fù)合材料艙門優(yōu)化設(shè)計(jì)
    地球脈動(dòng)(第一季)
    空腔參數(shù)對重力壩穩(wěn)定的影響分析
    前置污水去油池
    亚洲国产欧美网| 午夜福利在线免费观看网站| 亚洲熟女精品中文字幕| 久久影院123| 超碰成人久久| 搡老岳熟女国产| 黑丝袜美女国产一区| 日本午夜av视频| 成人影院久久| 另类精品久久| 久久天堂一区二区三区四区| 日韩大码丰满熟妇| 制服人妻中文乱码| 丝袜脚勾引网站| 少妇精品久久久久久久| 男女之事视频高清在线观看 | 伦理电影免费视频| 在线 av 中文字幕| tube8黄色片| 性高湖久久久久久久久免费观看| 18禁黄网站禁片午夜丰满| 校园人妻丝袜中文字幕| 亚洲精品成人av观看孕妇| 无遮挡黄片免费观看| 成人亚洲欧美一区二区av| 欧美日韩福利视频一区二区| 91精品伊人久久大香线蕉| 悠悠久久av| 亚洲第一青青草原| 国产精品一国产av| 国产免费视频播放在线视频| av天堂在线播放| 亚洲欧美精品综合一区二区三区| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲七黄色美女视频| 国产97色在线日韩免费| 亚洲国产最新在线播放| 在线观看人妻少妇| 波野结衣二区三区在线| 精品久久蜜臀av无| 国产精品99久久99久久久不卡| 日韩电影二区| av网站免费在线观看视频| 大码成人一级视频| 狂野欧美激情性bbbbbb| 亚洲精品一二三| 纵有疾风起免费观看全集完整版| 国产高清videossex| 亚洲午夜精品一区,二区,三区| 国产精品久久久av美女十八| 中文字幕人妻丝袜制服| 肉色欧美久久久久久久蜜桃| 老司机在亚洲福利影院| 在线 av 中文字幕| 一边摸一边抽搐一进一出视频| 久久精品久久久久久久性| 18禁国产床啪视频网站| 男人操女人黄网站| 国产免费一区二区三区四区乱码| 国产精品麻豆人妻色哟哟久久| 一级毛片黄色毛片免费观看视频| 成年人黄色毛片网站| 丝袜美腿诱惑在线| 亚洲熟女毛片儿| 在线 av 中文字幕| 婷婷色av中文字幕| 国产成人一区二区在线| 日本色播在线视频| 亚洲精品久久午夜乱码| av在线app专区| 大陆偷拍与自拍| 午夜老司机福利片| 日韩免费高清中文字幕av| 91成人精品电影| 久久中文字幕一级| 国产精品欧美亚洲77777| 老司机在亚洲福利影院| 99国产精品免费福利视频| 一区在线观看完整版| 久久青草综合色| 国产无遮挡羞羞视频在线观看| 男女床上黄色一级片免费看| 国产国语露脸激情在线看| 丝袜美足系列| 精品视频人人做人人爽| 两性夫妻黄色片| 人人妻人人添人人爽欧美一区卜| 悠悠久久av| 女人久久www免费人成看片| 18禁黄网站禁片午夜丰满| 国产免费一区二区三区四区乱码| 国产又爽黄色视频| 国产成人影院久久av| 国产一区二区 视频在线| 精品国产超薄肉色丝袜足j| 日本黄色日本黄色录像| 久久久久久亚洲精品国产蜜桃av| a 毛片基地| 少妇被粗大的猛进出69影院| 青草久久国产| 精品卡一卡二卡四卡免费| 极品少妇高潮喷水抽搐| 色精品久久人妻99蜜桃| 亚洲 欧美一区二区三区| 久久久久久人人人人人| 精品国产一区二区三区久久久樱花| 成人国语在线视频| 亚洲精品第二区| 尾随美女入室| 亚洲伊人久久精品综合| 亚洲自偷自拍图片 自拍| 精品国产一区二区三区四区第35| 国产精品成人在线| 国产99久久九九免费精品| 久久精品国产a三级三级三级| 色网站视频免费| 18禁裸乳无遮挡动漫免费视频| 欧美人与性动交α欧美软件| 深夜精品福利| 菩萨蛮人人尽说江南好唐韦庄| 日韩免费高清中文字幕av| 亚洲国产看品久久| 成人手机av| av在线播放精品| av天堂在线播放| 久久99一区二区三区| 国产男女内射视频| 纯流量卡能插随身wifi吗| 亚洲av电影在线进入| 91麻豆av在线| 老汉色∧v一级毛片| 视频区欧美日本亚洲| 桃花免费在线播放| 国产日韩欧美视频二区| 人体艺术视频欧美日本| 久久ye,这里只有精品| 黄色 视频免费看| 久久女婷五月综合色啪小说| 丰满人妻熟妇乱又伦精品不卡| 一区二区三区激情视频| 成人亚洲精品一区在线观看| 黑人猛操日本美女一级片| 美女视频免费永久观看网站| 亚洲av日韩在线播放| 国产高清videossex| 欧美成人精品欧美一级黄| 人体艺术视频欧美日本| 久久久久精品国产欧美久久久 | 亚洲欧美清纯卡通| 亚洲国产精品999| 操出白浆在线播放| 老司机在亚洲福利影院| 在线精品无人区一区二区三| 国产精品亚洲av一区麻豆| 国产精品免费视频内射| 黑丝袜美女国产一区| 一边亲一边摸免费视频| 久久人人97超碰香蕉20202| 免费女性裸体啪啪无遮挡网站| 国产精品三级大全| 好男人电影高清在线观看| 国产伦人伦偷精品视频| 黄色怎么调成土黄色| 嫩草影视91久久| 久久人妻熟女aⅴ| 搡老乐熟女国产| 999久久久国产精品视频| 人体艺术视频欧美日本| 久久久久久久久免费视频了| 午夜福利视频在线观看免费| 免费高清在线观看视频在线观看| 黄网站色视频无遮挡免费观看| 在线天堂中文资源库| 日韩大片免费观看网站| 啦啦啦啦在线视频资源| 男女免费视频国产| 日日夜夜操网爽| 亚洲成人免费av在线播放| 老鸭窝网址在线观看| 国产精品一区二区免费欧美 | 国产欧美日韩精品亚洲av| 国产一区二区三区av在线| 午夜日韩欧美国产| 久久精品aⅴ一区二区三区四区| 欧美人与性动交α欧美精品济南到| 天天躁日日躁夜夜躁夜夜| avwww免费| 日本午夜av视频| 日韩免费高清中文字幕av| 两个人免费观看高清视频| 欧美日韩一级在线毛片| 久热爱精品视频在线9| 亚洲精品美女久久av网站| 中文字幕色久视频| 中文字幕制服av| 在线看a的网站| 两个人看的免费小视频| 亚洲国产看品久久| 黄色片一级片一级黄色片| 一区二区三区乱码不卡18| av天堂在线播放| 国产午夜精品一二区理论片| 老司机靠b影院| 在线观看人妻少妇| 男女下面插进去视频免费观看| 高清av免费在线| 亚洲熟女毛片儿| 大香蕉久久网| 亚洲欧美成人综合另类久久久| 一区二区三区四区激情视频| 狠狠精品人妻久久久久久综合| 亚洲久久久国产精品| 精品福利永久在线观看| 久久99热这里只频精品6学生| 免费少妇av软件| 亚洲天堂av无毛| 久久精品亚洲熟妇少妇任你| www.999成人在线观看| 桃花免费在线播放| 在线观看www视频免费| 中文字幕另类日韩欧美亚洲嫩草| 国精品久久久久久国模美| 一区福利在线观看| 成人国语在线视频| 亚洲精品国产区一区二| 亚洲精品日韩在线中文字幕| 又黄又粗又硬又大视频| 看免费成人av毛片| 久久午夜综合久久蜜桃| 美女中出高潮动态图| 99精品久久久久人妻精品| 女性被躁到高潮视频| 成年人免费黄色播放视频| 日韩av免费高清视频| 伊人久久大香线蕉亚洲五| 日韩一本色道免费dvd| 亚洲成人手机| 人人妻人人爽人人添夜夜欢视频| 国产成人啪精品午夜网站| av线在线观看网站| 亚洲天堂av无毛| 老司机在亚洲福利影院| 国产伦人伦偷精品视频| 亚洲国产欧美网| 亚洲,欧美精品.| 爱豆传媒免费全集在线观看| 丝瓜视频免费看黄片| 亚洲精品日韩在线中文字幕| 下体分泌物呈黄色| 中文字幕人妻丝袜制服| 午夜日韩欧美国产| 黄色一级大片看看| 高潮久久久久久久久久久不卡| 国产xxxxx性猛交| 欧美日韩一级在线毛片| 亚洲精品国产区一区二| 亚洲精品久久久久久婷婷小说| 欧美 亚洲 国产 日韩一| 久久久久精品人妻al黑| 亚洲av电影在线观看一区二区三区| 亚洲成人免费av在线播放| 69精品国产乱码久久久| 黄网站色视频无遮挡免费观看| 国产免费福利视频在线观看| 咕卡用的链子| 一区二区日韩欧美中文字幕| 黑丝袜美女国产一区| 搡老岳熟女国产| 啦啦啦 在线观看视频| 国产极品粉嫩免费观看在线| 啦啦啦中文免费视频观看日本| 婷婷色av中文字幕| 男女之事视频高清在线观看 | 热re99久久精品国产66热6| 你懂的网址亚洲精品在线观看| 一级毛片女人18水好多 | 丁香六月欧美| 欧美激情 高清一区二区三区| av网站在线播放免费| av天堂在线播放| 免费观看a级毛片全部| 777久久人妻少妇嫩草av网站| 丝袜美腿诱惑在线| 99九九在线精品视频| 亚洲,欧美精品.| 一边亲一边摸免费视频| 日韩 欧美 亚洲 中文字幕| 国产成人91sexporn| 欧美成人精品欧美一级黄| 久久精品成人免费网站| 老司机在亚洲福利影院| 天天操日日干夜夜撸| netflix在线观看网站| 91精品三级在线观看| 久久久精品94久久精品| 日本欧美视频一区| 日韩大片免费观看网站| 91成人精品电影| 黄色片一级片一级黄色片| 妹子高潮喷水视频| 免费人妻精品一区二区三区视频| 日韩一区二区三区影片| 日韩大码丰满熟妇| 一边摸一边做爽爽视频免费| 婷婷色麻豆天堂久久| 国产片特级美女逼逼视频| 叶爱在线成人免费视频播放| 久久精品国产综合久久久| 国产欧美日韩精品亚洲av| 制服人妻中文乱码| 国产91精品成人一区二区三区 | 两人在一起打扑克的视频| 悠悠久久av| 大香蕉久久成人网| 国产精品国产三级国产专区5o| 精品福利永久在线观看| 丝袜脚勾引网站| 国产成人精品久久久久久| 大片免费播放器 马上看| 国产一区二区三区综合在线观看| 搡老乐熟女国产| 午夜免费观看性视频| 大码成人一级视频| 国产真人三级小视频在线观看| 国产97色在线日韩免费| 一个人免费看片子| 一本色道久久久久久精品综合| 亚洲av在线观看美女高潮| 观看av在线不卡| 免费看av在线观看网站| 两性夫妻黄色片| 亚洲,一卡二卡三卡| 国产又色又爽无遮挡免| av又黄又爽大尺度在线免费看| 涩涩av久久男人的天堂| 久久久久久亚洲精品国产蜜桃av| 亚洲精品一区蜜桃| 老司机午夜十八禁免费视频| videos熟女内射| 国产麻豆69| 午夜视频精品福利| 人人妻人人澡人人看| 午夜视频精品福利| 视频区图区小说| 国产欧美日韩一区二区三区在线| 一区在线观看完整版| 国产精品免费视频内射| 亚洲精品久久久久久婷婷小说| 狂野欧美激情性bbbbbb| 国产精品国产三级专区第一集| 久久国产精品影院| 欧美日韩成人在线一区二区| 超碰97精品在线观看| 欧美精品一区二区免费开放| 欧美成人精品欧美一级黄| 波多野结衣一区麻豆| 婷婷色av中文字幕| 操出白浆在线播放| 精品亚洲乱码少妇综合久久| 在线av久久热| 国产日韩欧美在线精品| 亚洲国产中文字幕在线视频| 人妻 亚洲 视频| 一边亲一边摸免费视频| 国产午夜精品一二区理论片| 青春草亚洲视频在线观看| 1024香蕉在线观看| 亚洲欧美日韩高清在线视频 | avwww免费| 国产男女内射视频| 国产精品国产av在线观看| av一本久久久久| 免费看不卡的av| 国产视频首页在线观看| 欧美黑人欧美精品刺激| 啦啦啦中文免费视频观看日本| 久久久国产欧美日韩av| 两性夫妻黄色片| 国产精品久久久人人做人人爽| 久久国产精品大桥未久av| 少妇的丰满在线观看| 国产黄色视频一区二区在线观看| 大香蕉久久网| 中文精品一卡2卡3卡4更新| 一边摸一边做爽爽视频免费| 中文字幕人妻丝袜制服| 亚洲国产精品成人久久小说| 亚洲欧美一区二区三区黑人| 色播在线永久视频| 久久亚洲精品不卡| 成人国产av品久久久| 久久这里只有精品19| 久久国产亚洲av麻豆专区| 精品一品国产午夜福利视频| 国产一区二区在线观看av| 一级黄片播放器| 看十八女毛片水多多多| 久久久精品国产亚洲av高清涩受| 高清欧美精品videossex| www.999成人在线观看| 亚洲国产中文字幕在线视频| 美女视频免费永久观看网站| 日韩大码丰满熟妇| 婷婷色综合www| 午夜福利,免费看| 啦啦啦在线观看免费高清www| 亚洲av在线观看美女高潮| 观看av在线不卡| 999精品在线视频| 国产精品一国产av| 日韩一区二区三区影片| 亚洲激情五月婷婷啪啪| 欧美国产精品va在线观看不卡| 在线天堂中文资源库| 久久久久久久精品精品| 午夜福利乱码中文字幕| 亚洲国产欧美在线一区| 国产亚洲精品第一综合不卡| 欧美日韩视频精品一区| 成年人午夜在线观看视频| 免费观看人在逋| 天天躁日日躁夜夜躁夜夜| 视频区欧美日本亚洲| 高清不卡的av网站| 看免费成人av毛片| 日日夜夜操网爽| 黄片小视频在线播放| 欧美中文综合在线视频| 精品人妻1区二区| 国产国语露脸激情在线看| 999久久久国产精品视频| 操出白浆在线播放| 欧美日韩黄片免| 日韩av在线免费看完整版不卡| 久久精品亚洲熟妇少妇任你| 亚洲欧美一区二区三区黑人| 日本欧美国产在线视频| 色网站视频免费| 咕卡用的链子| 国产欧美日韩一区二区三 | 日韩伦理黄色片| 黄色 视频免费看| 午夜av观看不卡| 亚洲精品一卡2卡三卡4卡5卡 | 人妻 亚洲 视频| 国产男人的电影天堂91| 国产在线观看jvid| av福利片在线| 免费在线观看视频国产中文字幕亚洲 | 人人妻,人人澡人人爽秒播 | 一本—道久久a久久精品蜜桃钙片| 超碰97精品在线观看| 精品福利永久在线观看| av电影中文网址| 美女高潮到喷水免费观看| 中文字幕最新亚洲高清| 一区福利在线观看| 国产一区亚洲一区在线观看| 日本五十路高清| 日本一区二区免费在线视频| 国产日韩欧美视频二区| 亚洲国产精品一区二区三区在线| 无限看片的www在线观看| 女警被强在线播放| 在线观看人妻少妇| av一本久久久久| 精品人妻熟女毛片av久久网站| 亚洲,一卡二卡三卡| 久久中文字幕一级| 侵犯人妻中文字幕一二三四区| 亚洲天堂av无毛| 新久久久久国产一级毛片| 国产精品免费大片| 久久精品国产亚洲av涩爱| 亚洲欧美一区二区三区久久| 日韩一本色道免费dvd| 丝袜在线中文字幕| 爱豆传媒免费全集在线观看| 18禁国产床啪视频网站| 丝袜美足系列| 亚洲国产欧美一区二区综合| 黄色 视频免费看| 桃花免费在线播放| 热99久久久久精品小说推荐| 精品高清国产在线一区| 精品福利永久在线观看| 悠悠久久av| 国产精品亚洲av一区麻豆| 国产午夜精品一二区理论片| 91精品国产国语对白视频| 国产免费福利视频在线观看| 亚洲精品在线美女| 久久久久久亚洲精品国产蜜桃av| 免费看不卡的av| 91精品国产国语对白视频| 亚洲av日韩在线播放| 日韩中文字幕欧美一区二区 | 成年动漫av网址| 91老司机精品| 你懂的网址亚洲精品在线观看| 亚洲精品美女久久av网站| 色播在线永久视频| 伊人亚洲综合成人网| 国产一级毛片在线| 亚洲成国产人片在线观看| 久久久精品免费免费高清| 亚洲精品久久午夜乱码| 黄频高清免费视频| 色婷婷久久久亚洲欧美| 99国产精品一区二区蜜桃av | 亚洲七黄色美女视频| 久久99热这里只频精品6学生| 欧美老熟妇乱子伦牲交| 久久午夜综合久久蜜桃| 午夜福利一区二区在线看| 欧美中文综合在线视频| 午夜福利乱码中文字幕| 国产成人啪精品午夜网站| av欧美777| 国产精品三级大全| 两个人免费观看高清视频| 欧美国产精品一级二级三级| 亚洲视频免费观看视频| 亚洲情色 制服丝袜| 深夜精品福利| 脱女人内裤的视频| 五月开心婷婷网| 国产免费一区二区三区四区乱码| 成在线人永久免费视频| 又黄又粗又硬又大视频| 亚洲欧美精品自产自拍| a级毛片黄视频| 天天躁夜夜躁狠狠久久av| 亚洲欧美精品综合一区二区三区| 1024视频免费在线观看| 热re99久久精品国产66热6| 一本—道久久a久久精品蜜桃钙片| 国产不卡av网站在线观看| 亚洲欧美一区二区三区久久| 欧美日韩成人在线一区二区| 丝袜喷水一区| 我的亚洲天堂| 高清黄色对白视频在线免费看| 免费在线观看影片大全网站 | 日韩制服丝袜自拍偷拍| 久久久精品94久久精品| 色94色欧美一区二区| 日韩 亚洲 欧美在线| 黄色片一级片一级黄色片| 久久久精品免费免费高清| 大片免费播放器 马上看| 亚洲精品日韩在线中文字幕| videos熟女内射| 人人澡人人妻人| 亚洲七黄色美女视频| 国产高清视频在线播放一区 | 欧美日韩av久久| bbb黄色大片| 另类亚洲欧美激情| 免费一级毛片在线播放高清视频 | 成在线人永久免费视频| 国产亚洲一区二区精品| 亚洲精品乱久久久久久| 日本黄色日本黄色录像| 精品福利观看| 国产成人91sexporn| 国产伦理片在线播放av一区| 热99国产精品久久久久久7| 麻豆国产av国片精品| 精品国产一区二区三区四区第35| 国产成人免费无遮挡视频| 1024视频免费在线观看| 国产野战对白在线观看| 日本av手机在线免费观看| 亚洲欧美一区二区三区黑人| 中文字幕制服av| 国产片内射在线| 你懂的网址亚洲精品在线观看| 亚洲精品中文字幕在线视频| 国产男女内射视频| 老司机影院毛片| 久久久精品区二区三区| 日韩中文字幕视频在线看片| 岛国毛片在线播放| 久久人人爽人人片av| 亚洲成国产人片在线观看| 成年人免费黄色播放视频| 日韩 亚洲 欧美在线| 一级毛片 在线播放| 一二三四社区在线视频社区8| 老司机在亚洲福利影院| 99热网站在线观看| 久久精品国产亚洲av涩爱| 水蜜桃什么品种好| 桃花免费在线播放| 欧美日本中文国产一区发布| 丝袜脚勾引网站| 两人在一起打扑克的视频| 老汉色av国产亚洲站长工具| 黄色视频在线播放观看不卡| 极品人妻少妇av视频| 在线精品无人区一区二区三| 777久久人妻少妇嫩草av网站| 女人精品久久久久毛片| www.999成人在线观看| 九色亚洲精品在线播放| 国产成人免费观看mmmm| 国产欧美亚洲国产| 侵犯人妻中文字幕一二三四区| 久久久久视频综合| 日韩熟女老妇一区二区性免费视频| 午夜两性在线视频| 一级片'在线观看视频| 久久女婷五月综合色啪小说| 超碰97精品在线观看|