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

    亞臨界圓柱繞流的DES方法比較

    2017-11-20 01:53:50唐虎常士楠成竹馬蘭
    航空學(xué)報(bào) 2017年3期
    關(guān)鍵詞:尾跡湍流圓柱

    唐虎, 常士楠,*, 成竹, 馬蘭

    1.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院, 北京 100083 2.中國飛機(jī)強(qiáng)度研究所, 西安 710065

    亞臨界圓柱繞流的DES方法比較

    唐虎1, 常士楠1,*, 成竹2, 馬蘭2

    1.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院, 北京 100083 2.中國飛機(jī)強(qiáng)度研究所, 西安 710065

    針對地面結(jié)冰試驗(yàn)設(shè)施噴霧系統(tǒng)尾跡湍流模擬需求,考察了3種脫體渦模擬(DES)方法對三維亞臨界圓柱繞流預(yù)測的準(zhǔn)確性,比較和分析了瞬時(shí)流動特征和流場統(tǒng)計(jì)量。研究發(fā)現(xiàn):1) 從瞬時(shí)流動特征來看,k-ω雷諾平均Navier-Stokes(RANS)控制的區(qū)域?qū)羟袘?yīng)力輸運(yùn)k-ωDES (SSTk-ωDES)求解的準(zhǔn)確性有顯著影響;2) 從流動參數(shù)統(tǒng)計(jì)量的誤差范圍來看,回流區(qū)長度和流向最小速度的預(yù)測質(zhì)量能在一定程度上反應(yīng)圓柱繞流數(shù)值模擬的準(zhǔn)確性;3) 綜合比較,SSTk-ωDES預(yù)測值與試驗(yàn)值和大渦模擬(LES)數(shù)據(jù)最為吻合,具有應(yīng)用于地面結(jié)冰試驗(yàn)設(shè)施噴霧系統(tǒng)尾跡湍流計(jì)算的潛力。

    湍流擴(kuò)散; 尾跡流動; 脫體渦模擬; 渦結(jié)構(gòu); 統(tǒng)計(jì)量

    地面結(jié)冰試驗(yàn)設(shè)施由于具備可控的結(jié)冰云霧模擬能力而在飛機(jī)結(jié)冰研究和冰防護(hù)系統(tǒng)發(fā)展過程中發(fā)揮重要作用,而液態(tài)水含量均勻性作為地面結(jié)冰試驗(yàn)設(shè)施的一個(gè)重要設(shè)計(jì)指標(biāo),是設(shè)計(jì)階段必須要考慮的。此前的研究[1-7]已經(jīng)表明地面結(jié)冰試驗(yàn)設(shè)施噴霧棒尾跡湍流對水滴擴(kuò)散有顯著影響。要正確地預(yù)測尾跡湍流對水滴的擴(kuò)散行為,就必須合適地捕捉噴霧棒尾跡湍流。

    氣流繞噴霧棒流動屬于亞臨界鈍體繞流,涉及邊界層分離、剪切層發(fā)展、非定常渦的產(chǎn)生與脫落、渦與渦之間相互作用等復(fù)雜物理現(xiàn)象。脫體渦模擬(Detached Eddy Simulation,DES)是一種致力于解決三維非定常大分離流動的湍流求解方法[8],最早由Spalart等[9]提出,該方法在網(wǎng)格分辨率適用于大渦模擬(Large Eddy Simulation,LES)的區(qū)域發(fā)揮類似于LES亞格子尺度模型的作用,而在網(wǎng)格分辨率不適用于LES的區(qū)域發(fā)揮雷諾平均Navier-Stokes(Reynolds Averaged Navier-Stokes,RANS)模型的作用[10],這樣可在保證計(jì)算精度可接受的同時(shí)較大幅度地節(jié)省計(jì)算資源。

    DES方法經(jīng)過逐步改進(jìn)[11-19],已成功應(yīng)用于多項(xiàng)實(shí)際工程中[20-32]。目前構(gòu)建的較為常用的DES方法有Spalart-Allmaras DES(S-A DES)[33]方法、Realizablek-εDES(Rk-εDES)[34]方法和剪切應(yīng)力輸運(yùn)k-ωDES(SSTk-ωDES)[35]方法(k為湍動能;ε為耗散率;ω為比耗散率)。在DES方法比較研究方面,Strelets[8]在將S-A DES推廣至SSTk-ωDES時(shí)以ReD=50 000(ReD為圓柱雷諾數(shù))圓柱繞流案例進(jìn)行了模型測試,比較了瞬時(shí)渦量、阻力系數(shù)、升力系數(shù)、雷諾應(yīng)力和圓柱表面壓力系數(shù)分布,但僅有壓力系數(shù)分布與試驗(yàn)數(shù)據(jù)進(jìn)行了比較。Nichols[36]則以ReD=8×106圓柱繞流案例對S-A DES、SSTk-ωDES和SST多尺度模型進(jìn)行了比較,分析了計(jì)算結(jié)果對時(shí)間步長和網(wǎng)格分辨率的敏感性,認(rèn)為采用隱式算法時(shí)需要對每個(gè)渦脫落周期捕捉至少200步,但Nichols[36]僅對阻力系數(shù)CD和斯特羅哈爾數(shù)St與試驗(yàn)數(shù)據(jù)進(jìn)行了比較,沒有對這3種方法預(yù)測的流動統(tǒng)計(jì)量以及瞬時(shí)流動特征與試驗(yàn)進(jìn)行全面比較,因此難以說明哪種方法更適用于亞臨界流動。

    作為噴霧棒尾跡湍流對水滴的擴(kuò)散研究的第1步,本研究擬采用上述3種DES方法對亞臨界圓柱繞流進(jìn)行計(jì)算,并將預(yù)測的瞬時(shí)流動特征和流動統(tǒng)計(jì)量與可用的試驗(yàn)數(shù)據(jù)進(jìn)行比較與分析,確定可準(zhǔn)確描述亞臨界鈍體尾跡湍流特征的DES方法,為下一步不同噴霧棒外流場數(shù)值計(jì)算和水滴擴(kuò)散計(jì)算奠定基礎(chǔ)。

    1 數(shù)值方法

    1.1 控制方程

    (1)

    將這種形式的速度表達(dá)式代入瞬態(tài)連續(xù)性方程和動量方程,并取時(shí)間平均可得到

    (2)

    式中:ρ為密度;t為時(shí)間。

    (3)

    湍流建模的RANS方法要求式(3)中的雷諾應(yīng)力適當(dāng)模化。一般采用Boussinesq假設(shè)將雷諾應(yīng)力與平均速度梯度聯(lián)系起來,即

    (4)

    式中:μt為湍流黏性。

    下面給出本文考察的3種雷諾應(yīng)力?;疍ES方法。

    1) S -A DES

    (5)

    式中:CDES,S -A為S-A DES的模型常數(shù),取值為0.65;Δmax為網(wǎng)格在x、y或z方向的最大值,即

    Δmax=max(Δx,Δy,Δz)

    (6)

    2) Rk-εDES

    在Rk-εDES模型中,耗散項(xiàng)Yk由式(7)修正。

    (7)

    式中:lDES為DES模型湍流長度,其表達(dá)式為

    lDES=min(lRk ε,lLES)

    (8)

    (9)

    lLES=CDES,Rk εΔmax

    (10)

    其中:lRk ε和lLES分別為Rk-εRANS模型和LES模型湍流長度;CDES,Rk ε為Rk-εDES的模型常數(shù),取值為0.61。

    3) SSTk-ωDES

    采用式(11)對SSTk-ωDES模型中的湍動能耗散項(xiàng)Yk進(jìn)行修正。

    Yk=ρβ*kωFDES

    (11)

    式中:β*為修正系數(shù);FDES的表達(dá)式為

    (12)

    其中:CDES,SST為SSTk-ωDES的模型常數(shù),取值為0.61;Lt為湍流長度尺度,在RANS中的定義為

    (13)

    1.2 求解方法

    采用ANSYS 14.5 Fluent軟件平臺開展數(shù)值計(jì)算。由于研究的流動為低速不可壓縮非定常流動,因此選用基于壓力的求解器進(jìn)行瞬態(tài)模擬,并采用壓力隱式分裂算子(Pressure-Implicit with Splitting of Operators,PISO)進(jìn)行壓力-速度耦合求解。對于瞬態(tài)模擬,控制方程必須在空間和時(shí)間上進(jìn)行離散。在空間離散上,連續(xù)性方程采用有限差分方法,并采用算術(shù)平均處理密度插值;動量方程采用有界中心差分格式(Bounded Central Differencing Scheme)進(jìn)行離散,其壓力梯度項(xiàng)采用基于單元的最小平方法進(jìn)行計(jì)算,壓力項(xiàng)采用動量方程系數(shù)處理表面處的壓力差值;湍動能和耗散率采用二階迎風(fēng)格式。在時(shí)間離散上,采用二階隱式格式。

    2 結(jié)果與討論

    2.1 計(jì)算域、網(wǎng)格與邊界條件

    選用具有代表性的ReD=3 900亞臨界圓柱繞流進(jìn)行模型測試,因?yàn)檫@種流動有Lourenco和Shih[37]、Norberg[37](報(bào)道于文獻(xiàn)[37]中)、Ong和Wallace[38]、Parnaudeau等[39]的試驗(yàn)數(shù)據(jù)用于比較,而且這些試驗(yàn)數(shù)據(jù)也被廣泛地用于模型測試。本文主要采用Parnaudeau等[39]在2008年采用熱線風(fēng)速儀(Hot Wire Anemometry, HWA)和粒子圖像測速(Particle Image Velocimetry,PIV)獲取的ReD=3 900圓柱繞流試驗(yàn)數(shù)據(jù)進(jìn)行模型驗(yàn)證,另外還引用了Parnaudeau等[39]在相同條件下的LES數(shù)據(jù)進(jìn)行比較。文獻(xiàn)[39]中圓柱直徑D=12 mm,長徑比L/D=24,空氣密度ρ=1.205 kg/m3,空氣溫度ta=20 ℃,空氣動力黏性系數(shù)μa=1.81×10-5Pa·s,來流空氣速度Uc=4.8 m/s。

    參照文獻(xiàn)[37]中的圓柱繞流計(jì)算域選取經(jīng)驗(yàn),確定計(jì)算域?yàn)椋毫飨?15≤x/D≤15,側(cè)向-15≤y/D≤15,展向-2≤z/D≤2,圓柱體中心位于x/D=0的位置,如圖 1(a)所示。在網(wǎng)格設(shè)計(jì)時(shí),要考慮圓柱邊界層長度尺度、分離剪切層和渦結(jié)構(gòu)。Spalart[40]和Bunge等[41]對DES方法的使用、網(wǎng)格系統(tǒng)和網(wǎng)格尺寸(Δx,Δy,Δz)設(shè)計(jì)給出了建議。此處參考Travin等[42]開展圓柱繞流DES計(jì)算時(shí)的網(wǎng)格系統(tǒng)設(shè)計(jì)方法,將流體計(jì)算域分為3塊,即RANS域、LES域和歐拉域,3個(gè)域的交點(diǎn)A由剪切層長度確定,約為(1.09,1.03),如圖1(b)所示。

    為了準(zhǔn)確地描述邊界層的發(fā)展與脫落,RANS域第1層網(wǎng)格厚度Δy取0.003D,對應(yīng)的y+≈1,y方向網(wǎng)格增長率取1.15,這樣可保證邊界層內(nèi)(以0.99Uc確定)分布有足夠的節(jié)點(diǎn)。這樣得到RANS域在交點(diǎn)A處的Δy約為0.045D,考慮到DES方法要求LES域具有近似的各向同性網(wǎng)格,因此RANS域第1層網(wǎng)格Δx取0.015D,這樣可以保證交點(diǎn)A處的Δx≈Δy。展向網(wǎng)格分辨率要考慮渦結(jié)構(gòu)尺度,根據(jù)Mansy[43]和Williamson[44]等的試驗(yàn)研究,在近圓柱尾跡區(qū)域的展向波長尺度約為

    (14)

    根據(jù)式(14),ReD=3 900時(shí)的展向波長λz/D≈0.4,故第1層網(wǎng)格展向取Δz=0.05D可保證每個(gè)波長內(nèi)分布約8個(gè)節(jié)點(diǎn)。因此整個(gè)網(wǎng)格系統(tǒng)節(jié)點(diǎn)分布:RANS域?yàn)?95×65×85,LES域?yàn)?46×151×85,歐拉域?yàn)?5×51×85,總網(wǎng)格數(shù)約為200萬。

    根據(jù)Parnaudeau等[39]的試驗(yàn)條件,設(shè)置計(jì)算的入口速度Uc=4.8 m/s,入口湍流強(qiáng)度為0.2%,圓柱表面為無滑移靜止壁面,展向方向設(shè)置為對稱邊界。為了保證方程求解滿足CFL(Courant-Fredrichs-Levy)限制條件,取固定時(shí)間步長Δt=0.02D/Uc進(jìn)行計(jì)算,這樣可使每個(gè)渦脫落周期捕捉約240次,滿足文獻(xiàn)[36]提到的時(shí)間步長要求。

    2.2 瞬時(shí)流動特征

    在對DES模擬數(shù)據(jù)進(jìn)行處理時(shí),發(fā)現(xiàn)Rk-εDES給出的預(yù)測數(shù)據(jù)具有顯著的RANS特征,為了表明Rk-εDES結(jié)果的可靠性,此處同時(shí)引入Rk-εRANS結(jié)果進(jìn)行比較。圖2和圖3 分別給出了瞬態(tài)展向渦分量ΩzD/Uc和側(cè)向渦分量ΩyD/Uc在2個(gè)對稱面(z=0和y=0)上的分布,為了與Parnaudeau等[39]的PIV試驗(yàn)數(shù)據(jù)和LES數(shù)據(jù)進(jìn)行比較,此處坐標(biāo)軸與渦量范圍的選取均與Parnaudeau等[39]的一致。從圖2 可看出,在z=0平面,S-A DES和SSTk-ωDES在圓柱兩側(cè)均預(yù)測出兩條細(xì)長的剪切層,且剪切層幾乎一直延伸到x/D=2,這與試驗(yàn)測量以及LES結(jié)果較為吻合,然而Rk-εDES卻預(yù)測出了兩條相對粗短的剪切區(qū)域,沒能預(yù)測出明顯的剪切層特征。在尾跡區(qū)內(nèi),3種DES方法給出了尺度較為接近的渦結(jié)構(gòu),這是因?yàn)镈ES在該區(qū)域內(nèi)轉(zhuǎn)換為LES,但受上游解析能力的影響,DES在尾跡區(qū)內(nèi)的微小尺度渦結(jié)構(gòu)預(yù)測能力稍弱于LES。Rk-εRANS預(yù)測出的邊界層和剪切層與Rk-εDES得到的結(jié)果極為相似,而Rk-εDES有在近壁區(qū)保持為Rk-εRANS的特點(diǎn),這也表明Rk-εDES的預(yù)測結(jié)果是可靠的。另外,Rk-εRANS在尾跡區(qū)幾乎沒能捕捉到展向渦分量,這是由于RANS構(gòu)架下產(chǎn)生過高的湍流黏性,使得流場的脈動信息被劇烈耗散。

    從圖3可看出,在y=0平面,3種DES方法均預(yù)測出了交替出現(xiàn)的湍流渦以及無旋區(qū),因此3種DES方法均能預(yù)測出圓柱下游區(qū)的卡門渦街,而且預(yù)測的渦街排列與試驗(yàn)及LES較為接近,均與圓柱基本平行,表明了流動結(jié)構(gòu)的三維特征。但3種DES方法中,Rk-εDES預(yù)測的側(cè)向渦量ΩyD/Uc較弱。另外,Rk-εRANS幾乎沒能捕捉到側(cè)向渦分量,這說明RANS在預(yù)測流動的三維特征方面表現(xiàn)欠佳。

    文獻(xiàn)[39]沒有給出總渦量ΩD/Uc在z=0平面分布的試驗(yàn)數(shù)據(jù),這里給出3種DES方法在z=0平面的瞬態(tài)渦量(ΩD/Uc取0~10)以供比較,如圖4所示。從圖4可見,S -A DES和SSTk-ωDES均預(yù)測出了不穩(wěn)定的剪切層從圓柱表面分離,并逐步發(fā)展為周期性的卡門渦街這一過程,而Rk-εDES和Rk-εRANS對流動結(jié)構(gòu)的發(fā)展過程預(yù)測能力較弱。Rk-εDES在預(yù)測下游尾跡區(qū)的渦結(jié)構(gòu)能力方面與S-A DES和SSTk-ωDES相當(dāng),這也體現(xiàn)了DES方法在該區(qū)域保持LES特性的解析能力。Rk-εRANS預(yù)測的尾跡區(qū)渦量偏弱,這是由于RANS產(chǎn)生過高的湍流黏性。

    從構(gòu)造的角度來看,DES方法在邊界層區(qū)域保持RANS,在遠(yuǎn)場轉(zhuǎn)化為LES。而SSTk-ωDES中的RANS部分本身就包含了k-ωRANS和k-εRANS兩部分,僅是采用混合函數(shù)進(jìn)行了處理,使得SSTk-ωRANS在邊界層近壁區(qū)保持為k-ωRANS,在邊界層遠(yuǎn)壁區(qū)切換為k-εRANS,以發(fā)揮出更好的性能[45]。Strelets[8]認(rèn)為影響SSTk-ωDES準(zhǔn)確性的主要是k-εRANS控制的區(qū)域。本文給出的Rk-εDES和SSTk-ωDES結(jié)果在邊界層分離和剪切層發(fā)展區(qū)域存在較大差異,然而Rk-εDES預(yù)測的邊界層分離和剪切層發(fā)展行為又與Rk-εRANS較為接近,這說明Rk-εDES在邊界層區(qū)域的確是采用k-εRANS進(jìn)行求解。以上分析表明k-ωRANS控制的區(qū)域?qū)STk-ωDES求解的準(zhǔn)確性也有較大影響。

    從渦量數(shù)據(jù)來看,流向渦分量ΩxD/Uc和側(cè)向渦分量ΩyD/Uc的量值相當(dāng);展向渦分量ΩzD/Uc較大,是渦量的主要分量,這進(jìn)一步體現(xiàn)了流動的三維特征。

    2.3 流動統(tǒng)計(jì)量

    在比較流動統(tǒng)計(jì)量時(shí),首先要考慮流場統(tǒng)計(jì)量的收斂性。很多研究者在對流場統(tǒng)計(jì)量進(jìn)行處理時(shí)所取的平均時(shí)間都不一致,例如Kravchenko和Moin[37]用7個(gè)渦脫落周期計(jì)算流場統(tǒng)計(jì)量,而Dong等[46]使用40~50個(gè)渦脫落周期,Ma等[47]又采用131個(gè)渦脫離周期。Franke和Frank[48]研究了流場統(tǒng)計(jì)量的收斂性,認(rèn)為至少需要40個(gè)渦脫離周期才能獲取收斂的流動統(tǒng)計(jì)量。目前的數(shù)值模擬均在尾跡湍流形態(tài)規(guī)律穩(wěn)定之后開始對瞬態(tài)流場進(jìn)行采樣。所有的統(tǒng)計(jì)量都是基于T=240D/Uc(約50個(gè)渦脫落周期)取平均得到,平均符號由“〈·〉”表示。

    表1中給出了ReD=3 900圓柱繞流統(tǒng)計(jì)量。圓括號內(nèi)的值表示數(shù)值計(jì)算數(shù)據(jù)與試驗(yàn)數(shù)據(jù)的誤差,“+”表示偏大,“-”表示偏小。從表 1可以看出,3種DES方法得到的阻力系數(shù)CD均稍小于試驗(yàn)測量值,而LES方法得到的CD稍大于試驗(yàn)測量值,這說明DES方法預(yù)計(jì)的分離點(diǎn)比試驗(yàn)測量得到的分離點(diǎn)稍微延后,這也由DES方法計(jì)算得到的分離角θsep均大于試驗(yàn)測量值所證實(shí)。3種DES方法預(yù)計(jì)的阻力系數(shù)與試驗(yàn)最為接近的是SSTk-ωDES,誤差僅為1%。-Cpb為背壓系數(shù),Cpb為圓柱基點(diǎn)(Base Point,θ=180°)處的壓力系數(shù)。3種DES方法計(jì)算的-Cpb均比試驗(yàn)測量值小,精度最高的是SSTk-ωDES模型,誤差為14%,精度不及LES方法(誤差為7%),這和Norberg[37]是在ReD=4 020條件下測量的-Cpb有關(guān)。此前有研究者認(rèn)為,St數(shù)(St=fD/Uc,f為渦脫落頻率)對網(wǎng)格分辨率和時(shí)間步長不敏感,筆者在網(wǎng)格分辨率和時(shí)間步長調(diào)試過程中也發(fā)現(xiàn)了這一點(diǎn),但DES模型仍然對St數(shù)有輕微影響,3種DES方法中SSTk-ωDES模型預(yù)計(jì)的St數(shù)與試驗(yàn)最為接近,誤差為1%,精度與LES相當(dāng)。3種DES方法計(jì)算的分離角θsep基本相同,最接近試驗(yàn)值的是Rk-εDES,誤差為5%。回流區(qū)長度Lrec為圓柱基點(diǎn)與中心線上平均流向速度〈u〉符號轉(zhuǎn)變點(diǎn)之間的距離,如圖5所示。由表1可看出,3種DES方法在回流區(qū)長度和流向最小速度umin預(yù)測方面存在較大差異,S-A DES預(yù)測的Lrec比Parnaudeau等[39]得到的試驗(yàn)值大39%,Rk-εDES預(yù)測的Lrec比Parnaudeau等[39]得到的試驗(yàn)值小50%,SSTk-ωDES預(yù)測的Lrec與Parnaudeau等[39]得到的試驗(yàn)值最為接近,僅偏大8%。S-A DES方法給出了較為準(zhǔn)確的流向最小速度umin,約偏低6%。從流動參數(shù)統(tǒng)計(jì)量的誤差范圍可以看出,3種DES方法得到的流動參數(shù)CD、-Cpb、St和θsep的誤差較為接近,僅回流區(qū)長度Lrec和流向最小速度umin的誤差有明顯差異,因此回流區(qū)長度Lrec和流向最小速度umin的預(yù)測質(zhì)量能在一定程度上反應(yīng)圓柱繞流數(shù)值模擬的準(zhǔn)確性。另外,通過比較Rk-εRANS和Rk-εDES預(yù)計(jì)的流動統(tǒng)計(jì)量可以發(fā)現(xiàn),Rk-εDES預(yù)計(jì)的流動統(tǒng)計(jì)量具有顯著的RANS特征,這表明在Rk-εDES中k-εRANS對邊界層區(qū)域表現(xiàn)出了過分抑制。

    表1 ReD=3 900圓柱繞流統(tǒng)計(jì)量Table 1 Statistics of flow around circular cylinder at ReD=3 900

    圖5給出了尾跡中心線上平均流向速度〈u〉分布。雖然3種DES方法預(yù)測的〈u〉變化趨勢幾乎一致,均是隨著尾跡中心線的延伸先逐步減小,達(dá)到最小速度umin后再增大,最后趨于平穩(wěn),但是僅有SSTk-ωDES模型得到的〈u〉在整個(gè)尾跡區(qū)(x/D≥0.5)與Parnaudeau等[39]的試驗(yàn)數(shù)據(jù)最為吻合(和LES數(shù)據(jù)基本一致),而S -A DES模型和Rk-εDES模型預(yù)測的〈u〉僅在遠(yuǎn)尾跡區(qū)(約x/D≥6)與試驗(yàn)數(shù)據(jù)較為接近。這里引入了Lourenco和Shih[37]以及Ong和Wallace[38]的試驗(yàn)數(shù)據(jù)作為比較。Ong和Wallace[38]僅給出了回流區(qū)下游的〈u〉分布,Lourenco和Shih[37]給出了包含回流區(qū)的〈u〉分布。Lourenco和Shih[37]試驗(yàn)測量得到的Lrec比Parnaudeau等[39]測量得到的Lrec小,Kravchenko和Moin[37]認(rèn)為這是由于內(nèi)流干擾導(dǎo)致剪切層過早轉(zhuǎn)捩形成的。

    將3種DES方法獲取的z=0平面3個(gè)不同流向位置處的平均流向速度〈u〉分布(見圖6)和平均側(cè)向速度〈v〉分布(見圖7)與Parnaudeau等[39]以及Lourenco和Shih[37]的近尾跡區(qū)試驗(yàn)數(shù)據(jù)進(jìn)行比較,并引入Parnaudeau等[39]的LES數(shù)據(jù)作為參考。為了便于與試驗(yàn)數(shù)據(jù)進(jìn)行比較,流向位置的選取與試驗(yàn)一致。

    總體看來,〈u〉分布呈現(xiàn)出明顯的對稱性,〈v〉分布呈現(xiàn)出明顯的反對稱性。從模擬與試驗(yàn)數(shù)據(jù)的接近程度來看,對于〈u〉的預(yù)測,S -A DES數(shù)據(jù)偏大,Rk-εDES數(shù)據(jù)〈u〉偏??;對于〈v〉的預(yù)測,S -A DES模擬數(shù)據(jù)又偏小,而Rk-εDES數(shù)據(jù)偏大,只有SSTk-ωDES與Parnaudeau等[39]的試驗(yàn)數(shù)據(jù)最為接近,與LES僅存在極小的差異。值得注意的是,在位于回流區(qū)的x/D=1.06處(見圖 6(a)),S -A DES、SSTk-ωDES和LES都給出了U型輪廓的〈u〉分布,這和Parnaudeau等[39]試驗(yàn)給出的〈u〉分布輪廓幾乎完全一致。然而Rk-εDES卻給出了V型〈u〉分布,這又和Lourenco和Shih[37]試驗(yàn)給出的〈u〉分布輪廓接近。注意到試驗(yàn)和模擬數(shù)據(jù)在回流區(qū)邊緣區(qū)域x/D=1.54和回流區(qū)下游x/D=2.02處都給出了V型〈u〉分布,這說明〈u〉分布可能是從回流區(qū)的U型輪廓發(fā)展為回流下游尾跡區(qū)的V型輪廓。從圖7的〈v〉分布可明顯看出,除Lourenco和Shih[37]的試驗(yàn)數(shù)據(jù)外,所有的模擬和試驗(yàn)數(shù)據(jù)均給出了反對稱的〈v〉分布,這種異常行為也表明Lourenco和Shih[37]的試驗(yàn)測量可能受到某些外部干擾的影響,導(dǎo)致剪切層過早轉(zhuǎn)捩,影響回流區(qū)尺寸和近尾跡區(qū)的速度輪廓形狀,這也從圖 5中Lourenco和Shih[37]試驗(yàn)得到的回流區(qū)長度Lrec得到證實(shí)。

    圖8給出了圓柱表面z=0處壓力系數(shù)分布。由于Parnaudeau等[39]沒有給出圓柱表面壓力系數(shù)分布的試驗(yàn)數(shù)據(jù),這里用Norberg[37]在ReD=4 020 下的試驗(yàn)數(shù)據(jù)進(jìn)行比較。在DES方法中,盡管3種DES方法均給出了相同的Cp變化趨勢,但S-A DES和Rk-εDES(兩條線幾乎完全重合)方法給出的Cp分布與試驗(yàn)數(shù)據(jù)偏離較遠(yuǎn),SSTk-ωDES方法與試驗(yàn)數(shù)據(jù)較為接近,精度能比前兩種DES方法高1倍,但仍不及LES。另外,DES方法中駐點(diǎn)處的壓力系數(shù)略微大于1,這在Travin等[42]的圓柱繞流DES方法中也曾出現(xiàn),Travin等[42]認(rèn)為這種誤差會隨計(jì)算域的增大而減小,然而背壓卻對計(jì)算域的增大不敏感。研究較為關(guān)注的是背壓系數(shù),而不是駐點(diǎn)處的壓力系數(shù),因此在目前的計(jì)算域條件下得到的壓力系數(shù)分布誤差是可以接受的。

    3 結(jié) 論

    1) 從瞬時(shí)流動特征來看,k-ωRANS控制的區(qū)域?qū)STk-ωDES求解的準(zhǔn)確性有顯著影響。

    2) 從流動參數(shù)統(tǒng)計(jì)量的誤差范圍來看,回流區(qū)長度Lrec和流向最小速度umin的預(yù)測質(zhì)量能在一定程度上反應(yīng)圓柱繞流數(shù)值模擬的準(zhǔn)確性。

    3) 綜合比較,SSTk-ωDES預(yù)測值與試驗(yàn)和LES數(shù)據(jù)最為吻合,具有應(yīng)用于地面結(jié)冰試驗(yàn)設(shè)施噴霧系統(tǒng)尾跡湍流計(jì)算的潛力。

    [1] MAREK J, OLSEN W A. Turbulent dispersion of the icing cloud from spray nozzles used in icing tunnels: NASA TM-87316[R]. Washington, D.C.: NASA, 1986.

    [2] BARTLETT C S. Turbine engine icing spray bar design issues[J]. Journal of Engineering for Gas Turbines and Power, 1995, 117(3): 406-412.

    [3] DEANGELIS B C, LOTH E. Simulations of turbulent droplet dispersion in wind-tunnel icing clouds[J]. Journal of Aircraft, 1997, 34(2): 213-219.

    [4] HANCIR P, ANDERSON A, LOTH E. Computations of droplet distributions in the NASA icing research tunnel: AIAA-2000-0101[R]. Reston: AIAA, 2000.

    [5] BHARGAVA C, LOTH E, POTAPCZUK M. Aerodynamic simulations of the NASA Glenn icing research tunnel: AIAA-2003-0566[R]. Reston: AIAA, 2003.

    [6] LEE A, LOTH E. Droplet dispersion in the NASA Glenn icing research tunnel: AIAA-2010-7533[R]. Reston: AIAA, 2010.

    [7] CLARK K, MALINOWSKI M, LOTH E. Air flow and liquid water concentration simulations of the 2012 NASA Glenn icing research tunnel: AIAA-2012-2936[R]. Reston: AIAA, 2012.

    [8] STRELETS M. Detached eddy simulation of massively separated flows: AIAA-2001-0879[R]. Reston: AIAA, 2001.

    [9] SPALART P R, JOU W H, STRELETS M, et al. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[C]//1st AFOSR International Conference on DNS/LES. Dayton, OH: Greyden Press, 1997: 14-18.

    [10] SPALART P R. Detached-eddy simulation[J]. Annual Review of Fluid Mechanics, 2009, 41: 181-202.

    [11] VATSA V N, SINGER B A. Evaluation of a second accurate Navier-Stokes code for detached eddy simulation past a circular cylinder: AIAA-2003-4085[R]. Reston: AIAA, 2003.

    [12] HANSEN R P, FORSYTHE I R. Large and detached eddy simulation of flow over a circular cylinder using unstructured grids: AIAA-2003-0775[R]. Reston: AIAA, 2003.

    [13] ROY C J, DECHANT L J, PAYNE J L, et al. Bluff-body flow simulations using hybrid RANS/LES: AIAA-2003-3889[R]. Reston: AIAA, 2003.

    [14] LO S C, HOFFMANN K A, DIETIKER J F. Numerical investigation of high Reynolds number flow over square and circular cylinders[J]. Journal of Thermophisics and Heat Transfer, 2005, 19(1): 72-80.

    [15] NISHINO T, ROBERTS G T, ZHANG X. Unsteady RANS and detached-eddy simulations of flow around a circular cylinder in ground effect[J]. Journal of Fluid and Structure, 2008, 24(1): 18-33.

    [16] JEE S, SHARIFF K. Detached-eddy simulation based on theυ2-fmodel[J]. International Journal of Heat and Fluid Flow, 2014, 46(2): 84-101.

    [17] 徐晶磊, 高歌, 楊焱. 基于當(dāng)?shù)亓鲃咏Y(jié)構(gòu)的RANS/LES混合模型[J]. 航空學(xué)報(bào), 2014, 35(11): 2992-2999. XU J L, GAO G, YANG Y. A RANS/LES hybrid model based on local flow structure[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(11): 2992-2999 (in Chinese).

    [18] 王翔宇, 李棟. SST-DES在小分離流動數(shù)值模擬中的改進(jìn)[J]. 北京航空航天大學(xué)學(xué)報(bào), 2014, 40(9): 1245-1249. WANG X Y, LI D. Improved SST-DES in numerical simulation of mild separation[J]. Journal of Beijing University of Aeronautics and Astronautics, 2014, 40(9): 1245-1249 (in Chinese).

    [19] 李斌, 吳頌平. 基于湍流尺度的混合RANS/LES模型[J]. 北京航空航天大學(xué)學(xué)報(bào), 2008, 34(7): 755-758. LI B, WU S P. Hybrid RANS/LES model based on turbulent scale[J]. Journal of Beijing University of Aeronautics and Astronautics, 2008, 34(7): 755-758 (in Chinese).

    [20] 李棟, 焦予秦, SHOV I M, 等. Detached Eddy Simulation方法模擬不同類型翼型的失速特性[J]. 航空學(xué)報(bào), 2005, 26(4): 406-410. LI D, JIAO Y Q, SHOV I M, et al. Detached eddy simulation for airfoil stall[J]. Acta Aeronautica et Astronautica Sinica, 2005, 26(4): 406-410 (in Chinese).

    [21] 白俊強(qiáng), 張揚(yáng), 華俊. 一種濾波SST方法在翼型深失速模擬中的應(yīng)用[J]. 航空學(xué)報(bào), 2013, 34(5): 979-987. BAI J Q, ZHANG Y, HUA J. Application of filter-SST method in airfoil deep stall simulation[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(5): 979-987 (in Chinese).

    [22] 劉周, 楊云軍, 周偉江, 等. 基于RANS-LES混合方法的翼型大迎角非定常分離流動研究[J]. 航空學(xué)報(bào), 2014, 35(2): 372-380. LIU Z, YANG Y J, ZHOU W J, et al. Study of unsteady separation flow around airfoil at high angle of attack using hybrid RANS -LES method [J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(2): 372-380 (in Chinese).

    [23] DECK S. Zonal-detached-eddy simulation of the flow around a high-lift configuration[J]. AIAA Journal, 2005, 43(11): 2372-2384.

    [24] FORSYTHE J R, SQUIRES K D, WURTZLER E, et al. Detached-eddy simulation of the F-15E at high alpha[J]. Journal of Aircraft, 2004, 41(2): 193-200.

    [25] BROCK J M, SUBBAREDDY P K, CANDLER G V. Detached-eddy simulation of hypersonic capsule wake flow[J]. AIAA Journal, 2015, 53(1): 70-80.

    [26] 楊小龍, 林鐵平. 汽車外流場DES/RANS模擬研究[J]. 湖南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2011, 38(1): 29-34. YANG X L, LIN T P. DES and RANS of vehicle external flow field[J]. Journal of Hunan University (Natural Science), 2011, 38(1): 29-34 (in Chinese).

    [27] 謝超, 谷正氣, 楊振東, 等. 不同RANS/LES混合模型的汽車氣動噪聲分析[J]. 汽車工程, 2015, 37(4): 440-445. XIE C, GU Z Q, YANG Z D, et al. Analysis on vehicle aerodynamic noise with different hybrid RANS/LES models[J]. Automotive Engineering, 2015, 37(4): 440-445 (in Chinese).

    [28] 劉學(xué)強(qiáng), 伍貽兆. 用DES數(shù)值模擬具有橫向噴流的紊流流場[J]. 航空學(xué)報(bào), 2004, 25(3): 209-213. LIU X Q, WU Y Z. The computation of the lateral jet turbulence flow using DES method[J]. Acta Aeronautica et Astronautica Sinica, 2004, 25(3): 209-213 (in Chinese).

    [29] 周昊, 岑可法, 樊建人. 分離渦方法模擬濃氮?dú)夤躺淞鲀上喾欠€(wěn)態(tài)流動特性研究[J]. 中國電機(jī)工程學(xué)報(bào), 2005, 25(7): 1-6. ZHOU H, CEN K F, FAN J R. The numerical investigation on the transient characteristics of the gas-solid two-phase fuel rich-lean burn flow[J]. Proceedings of the CSEE, 2005, 25(7): 1-6 (in Chinese).

    [30] 孫明波, 梁劍寒, 王振國. 二維凹腔超聲速流動的混合RANS/LES模擬[J]. 推進(jìn)技術(shù), 2006, 27(2): 119-123. SUN M B, LIANG J H, WANG Z G. Hybrid RANS/LES simulation of the supersonic flow over two-dimentional cavities[J]. Journal of Propulsion Technology, 2006, 27(2): 119-123 (in Chinese).

    [31] 汪洪波, 孫明波, 吳海燕, 等. 超聲速燃燒凹腔質(zhì)量交換特性的混合RANS/LES模擬[J]. 航空動力學(xué)報(bào), 2010, 25(1): 41-46. WANG H B, SUN M B, WU H Y, et al. Hybrid RANS/LES simulation of mass exchange characteristics of cavity for supersonic combustion[J]. Journal of Aerospace Power, 2010, 25(1): 41-46 (in Chinese).

    [32] 薛幫猛, 楊永. 基于兩方程湍流模型的DES方法在超聲速圓柱底部流動計(jì)算中的應(yīng)用[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2006, 24(5): 544-547. XUE B M, YANG Y. Technical details in applying DES method to computing supersonic cylinder-base flow[J]. Journal of Northwestern Polytechnical University, 2006, 24(5): 544-547 (in Chinese).

    [33] SHUR M, SPALART P R, STRELETS M, et al. Detached-eddy simulation of an airfoil at high angle of attack[C]//4th International Symposium on Engineering Turbulence Modeling and Experiments. Amsterdam: Elsevier Science Ltd., 1999: 669-678.

    [34] SHIH T H, LIOU W W, SHABBIR A, et al. A newk-εeddy-viscosity model for high Reynolds number turbulent flow[J]. Computers & Fluids, 1995, 24(3): 227-238.

    [35] MENTER F R, KUNTZ M, LANGTRY R. Ten years of industrial experience with the SST turbulence model[J]. Turbulence, Heat and Mass Transfer, 2003, 4(1): 625-632.

    [36] NICHOLS R H. Comparison of hybrid RANS/LES turbulence models on a circular cylinder at high Reynolds number: AIAA-2005-0498[R]. Reston: AIAA, 2005.

    [37] KRAVCHENKO A G, MOIN P. Numerical studies of flow over a circular cylinder atReD=3 900[J]. Physics of Fluids, 2000, 12(2): 403-417.

    [38] ONG L, WALLACE J. The velocity field of the turbulent very near wake of a circular cylinder[J]. Experiments in Fluids, 1996, 20(6): 441-453.

    [39] PARNAUDEAU P, CARLIER J, HEITZ D, et al. Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[J]. Physics of Fluids, 2008, 20(8): 1-14.

    [40] SPALART P R. Yong-person’s guide to detached-eddy simulation grids: NASA/CR-2001-211032[R]. Washington, D.C.: NASA, 2001.

    [41] BUNGE U, MOCKETT C, THIELE F. Guidelines for implementing detached-eddy simulation using different models[J]. Aerospace Science and Technology, 2007, 11(5): 376-385.

    [42] TRAVIN A, SHUR M, STRELETS M, et al. Detached-eddy simulations past a circular cylinder[J]. Flow, Turbulence and Combustion, 1999, 63(1): 293-313.

    [43] MANSY H, YANG P M, WILLIAMS D R. Quantitative measurements of three-dimensional structures in the wake of a circular cylinder[J]. Journal of Fluid Mechanics, 1994, 270: 227-296.

    [44] WILLIAMSON C H K, WU J, SHERIDAN J. Scaling of streamwise vortices in wakes[J]. Physics of Fluids, 1995, 7(10): 2307-2309.

    [45] MENTER F R. Two-equation eddy-viscosity turbulence model for engineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605.

    [46] DONG S, KARNIADAKIS G E, EKMEKCI A, et al. A combined direct numerical simulation-particle image velocimetry study of the turbulent near wake[J]. Journal of Fluid Mechanics, 2006, 569: 185-207.

    [47] MA X, KARAMANOS G S, KARNIADAKIS G E. Dynamics and low-dimensionality of turbulent near wake[J]. Journal of Fluid Mechanics, 2000, 410: 29-65.

    [48] FRANKE J, FRANK W. Large eddy simulation of the flow past a circular cylinder atRe=3 900[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2002, 90(10): 1191-1206.

    (責(zé)任編輯:鮑亞平)

    *Corresponding author. E-mail: sn_chang@buaa.edu.cn

    Comparison of detached eddy simulation schemes on a subcritical flow around circular cylinder

    TANG Hu1, CHANG Shinan1,*, CHENG Zhu2, MA Lan2

    1.SchoolofAeronauticScienceandEngineering,BeihangUniversity,Beijing100083,China2.AircraftStrengthResearchInstituteofChina,Xi’an710065,China

    Considering the requirements of simulating the flow around spray system in ground based icing test facilities, the accuracy of three detached eddy simulation (DES) schemes applied to three-dimensional subcritical flow around circular cylinder was examined through the comparison and analysis of instantaneous flow characteristics and flow statistic parameters. It is found that, from the point of instantaneous flow characterization, thek-ωReynolds averaged Navier-Stokes (RANS) branch has remarkable influence on the accuracy of shear stress transportk-ωDES (SSTk-ωDES). From the point of error range in flow statistics, the length of recirculation region and streamwise minimum velocity are the key parameters which determine the accuracy of numerical simulation on flow around circular cylinder. Compared comprehensively, the instantaneous flow characteristics and flow statistics obtained from SSTk-ωDES have well agreements with the data of experiments and large eddy simulation (LES), which thus has the potential of being applied to computing the wake turbulence of spray system in ground based icing test facilities.

    turbulent dispersion; wake flow; DES; vortex structure; statistics

    2016-04-07; Revised:2016-05-07; Accepted:2016-05-27; Published online:2016-06-07 08:40

    URL:www.cnki.net/kcms/detail/11.1929.V.20160607.0840.004.html

    s:National Natural Science Foundation of China (11372026); Technique Innovation Foundation of Aviation Industry Corporation of China (2013F62302)

    http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

    10.7527/S1000-6893.2016.0174

    2016-04-07; 退修日期:2016-05-07; 錄用日期:2016-05-27; 網(wǎng)絡(luò)出版時(shí)間:2016-06-07 08:40

    www.cnki.net/kcms/detail/11.1929.V.20160607.0840.004.html

    國家自然科學(xué)基金 (11372026); 中國航空工業(yè)集團(tuán)公司技術(shù)創(chuàng)新基金 (2013F62302)

    *通訊作者.E-mail: sn_chang@buaa.edu.cn

    唐虎, 常士楠, 成竹, 等. 亞臨界圓柱繞流的DES方法比較[J]. 航空學(xué)報(bào), 2017, 38(3): 120294. TANG H, CHANG S N, CHENG Z, et al. Comparison of detached eddy simulation schemes on a subcritical flow around circular cylinder[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(3): 120294.

    V216.8

    A

    1000-6893(2017)03-120294-11

    猜你喜歡
    尾跡湍流圓柱
    工程學(xué)和圓柱
    一種基于Radon 變換和尾跡模型的尾跡檢測算法
    圓柱的體積計(jì)算
    重氣瞬時(shí)泄漏擴(kuò)散的湍流模型驗(yàn)證
    基于EEMD-Hilbert譜的渦街流量計(jì)尾跡振蕩特性
    削法不同 體積有異
    “青春期”湍流中的智慧引渡(三)
    “青春期”湍流中的智慧引渡(二)
    弱分層湍流輸運(yùn)特性的統(tǒng)計(jì)分析
    圓柱殼的聲輻射特性分析
    河南科技(2014年14期)2014-02-27 14:11:56
    国产老妇女一区| 黄色欧美视频在线观看| 免费观看无遮挡的男女| 国产成人freesex在线| 欧美成人一区二区免费高清观看| 欧美bdsm另类| 80岁老熟妇乱子伦牲交| 美女被艹到高潮喷水动态| 婷婷色综合www| 精品久久国产蜜桃| 日日撸夜夜添| 国产精品福利在线免费观看| 久久久久久久久久人人人人人人| 免费少妇av软件| 国产v大片淫在线免费观看| 简卡轻食公司| 99九九线精品视频在线观看视频| 亚洲av.av天堂| 欧美精品人与动牲交sv欧美| 久久鲁丝午夜福利片| 亚洲天堂国产精品一区在线| 性色av一级| 男人爽女人下面视频在线观看| 丝袜美腿在线中文| 高清av免费在线| 国产精品无大码| 日韩一本色道免费dvd| 亚洲国产日韩一区二区| 观看免费一级毛片| 婷婷色av中文字幕| 在线 av 中文字幕| 好男人视频免费观看在线| 一级毛片久久久久久久久女| 亚洲精品国产成人久久av| 亚洲成人精品中文字幕电影| 大话2 男鬼变身卡| 成人免费观看视频高清| 欧美激情国产日韩精品一区| 亚洲精华国产精华液的使用体验| 亚洲人成网站高清观看| 狠狠精品人妻久久久久久综合| 亚洲精品国产av成人精品| xxx大片免费视频| av在线蜜桃| 晚上一个人看的免费电影| 亚洲国产精品成人综合色| 美女xxoo啪啪120秒动态图| 亚洲精品乱久久久久久| 国内揄拍国产精品人妻在线| 亚洲自拍偷在线| 国产永久视频网站| 日本免费在线观看一区| 在线精品无人区一区二区三 | 亚洲精品aⅴ在线观看| 久久久精品欧美日韩精品| 99久久精品热视频| 欧美bdsm另类| 亚洲伊人久久精品综合| 亚洲国产av新网站| 亚洲精品成人久久久久久| 久久久午夜欧美精品| 国产爱豆传媒在线观看| 亚洲国产高清在线一区二区三| 99视频精品全部免费 在线| 三级经典国产精品| 亚洲无线观看免费| 国产精品国产av在线观看| 欧美 日韩 精品 国产| 熟女人妻精品中文字幕| 黄色视频在线播放观看不卡| 少妇 在线观看| 国产乱来视频区| 五月天丁香电影| 一级二级三级毛片免费看| 亚洲内射少妇av| 亚洲欧洲日产国产| 啦啦啦在线观看免费高清www| 国产精品久久久久久精品电影| 深夜a级毛片| 亚洲真实伦在线观看| 日本色播在线视频| 亚洲欧美日韩无卡精品| av播播在线观看一区| 精品人妻视频免费看| 日韩欧美精品v在线| 欧美日韩视频精品一区| 婷婷色综合大香蕉| 哪个播放器可以免费观看大片| 日本与韩国留学比较| 亚洲精品视频女| 中文字幕av成人在线电影| 1000部很黄的大片| 午夜福利网站1000一区二区三区| 永久网站在线| 国产有黄有色有爽视频| 亚洲成人av在线免费| 三级经典国产精品| 久久久久久久久久久免费av| 亚洲av欧美aⅴ国产| 人妻夜夜爽99麻豆av| 夜夜看夜夜爽夜夜摸| 老司机影院毛片| 免费av不卡在线播放| 欧美激情久久久久久爽电影| 人妻一区二区av| 中国美白少妇内射xxxbb| 十八禁网站网址无遮挡 | 免费播放大片免费观看视频在线观看| 久久精品久久久久久噜噜老黄| 婷婷色综合大香蕉| 黄片无遮挡物在线观看| 我的女老师完整版在线观看| 男女边吃奶边做爰视频| 欧美国产精品一级二级三级 | 亚洲欧美一区二区三区国产| 日韩一区二区三区影片| 如何舔出高潮| av黄色大香蕉| 午夜精品一区二区三区免费看| 欧美老熟妇乱子伦牲交| 麻豆成人av视频| 成人免费观看视频高清| 国产黄片美女视频| 九九在线视频观看精品| 你懂的网址亚洲精品在线观看| 成人高潮视频无遮挡免费网站| 成人亚洲欧美一区二区av| 女人久久www免费人成看片| 亚洲成人av在线免费| 成年女人看的毛片在线观看| 午夜福利网站1000一区二区三区| 国产黄色免费在线视频| 亚洲欧美日韩无卡精品| 视频区图区小说| 亚洲怡红院男人天堂| 一本色道久久久久久精品综合| 国产一区亚洲一区在线观看| 色吧在线观看| 18禁裸乳无遮挡动漫免费视频 | 熟妇人妻不卡中文字幕| 中文在线观看免费www的网站| 国产精品久久久久久精品古装| 久久精品夜色国产| 男人狂女人下面高潮的视频| 97热精品久久久久久| 亚洲aⅴ乱码一区二区在线播放| 少妇人妻久久综合中文| 国产欧美另类精品又又久久亚洲欧美| 综合色av麻豆| 久久精品综合一区二区三区| 国产精品熟女久久久久浪| 乱系列少妇在线播放| 日本黄色片子视频| 在现免费观看毛片| 中文乱码字字幕精品一区二区三区| 下体分泌物呈黄色| 久热久热在线精品观看| 在线观看人妻少妇| 久久久a久久爽久久v久久| 午夜视频国产福利| 黄色一级大片看看| 午夜福利在线观看免费完整高清在| 亚洲精品自拍成人| 国产黄色免费在线视频| 日韩制服骚丝袜av| 亚洲美女视频黄频| 亚洲精品,欧美精品| 成人黄色视频免费在线看| 精品国产乱码久久久久久小说| 高清日韩中文字幕在线| 国产国拍精品亚洲av在线观看| 一本色道久久久久久精品综合| 国产精品三级大全| 日本爱情动作片www.在线观看| 亚洲欧美日韩卡通动漫| 国语对白做爰xxxⅹ性视频网站| 简卡轻食公司| 超碰97精品在线观看| 18禁动态无遮挡网站| 最近最新中文字幕大全电影3| 久热这里只有精品99| 人妻制服诱惑在线中文字幕| 美女脱内裤让男人舔精品视频| 熟女电影av网| 国产永久视频网站| 搡女人真爽免费视频火全软件| 精品久久久久久电影网| 一级片'在线观看视频| 欧美日韩视频高清一区二区三区二| 毛片一级片免费看久久久久| 国产精品久久久久久av不卡| 最后的刺客免费高清国语| 久久精品夜色国产| 成年av动漫网址| 大又大粗又爽又黄少妇毛片口| 亚洲精华国产精华液的使用体验| 一个人看视频在线观看www免费| 一个人看的www免费观看视频| 亚洲av欧美aⅴ国产| 精华霜和精华液先用哪个| 777米奇影视久久| 偷拍熟女少妇极品色| 久久97久久精品| 一级毛片久久久久久久久女| 97热精品久久久久久| 亚洲欧美日韩另类电影网站 | 日韩中字成人| 亚洲欧美成人综合另类久久久| 一本色道久久久久久精品综合| 在线观看美女被高潮喷水网站| 一本久久精品| 久久人人爽人人片av| 国产日韩欧美亚洲二区| 国产探花在线观看一区二区| 波野结衣二区三区在线| 欧美精品一区二区大全| 又大又黄又爽视频免费| 国产av不卡久久| 免费观看性生交大片5| 美女xxoo啪啪120秒动态图| 能在线免费看毛片的网站| 丝袜脚勾引网站| 国产视频内射| 日韩欧美一区视频在线观看 | 国产女主播在线喷水免费视频网站| 久久99精品国语久久久| 免费播放大片免费观看视频在线观看| 激情 狠狠 欧美| 国产在视频线精品| 极品教师在线视频| 亚洲丝袜综合中文字幕| 黄色一级大片看看| 熟女电影av网| 少妇人妻一区二区三区视频| 国产男人的电影天堂91| videossex国产| 日韩一区二区三区影片| 久久精品国产亚洲网站| 欧美一级a爱片免费观看看| 亚洲欧美日韩东京热| 午夜福利在线在线| 亚洲真实伦在线观看| 国产精品av视频在线免费观看| 毛片女人毛片| 久久99热6这里只有精品| 久久久久久久久久成人| 国产视频首页在线观看| 亚洲综合精品二区| 99热6这里只有精品| videos熟女内射| 三级男女做爰猛烈吃奶摸视频| 97人妻精品一区二区三区麻豆| 久久久精品94久久精品| 亚洲欧美清纯卡通| 直男gayav资源| 欧美人与善性xxx| 亚洲精品久久久久久婷婷小说| 精品一区二区免费观看| 久久综合国产亚洲精品| 搡老乐熟女国产| 日韩亚洲欧美综合| 免费黄色在线免费观看| av在线老鸭窝| 亚洲人成网站在线观看播放| 免费高清在线观看视频在线观看| 自拍欧美九色日韩亚洲蝌蚪91 | 中文资源天堂在线| h日本视频在线播放| 日韩制服骚丝袜av| av卡一久久| 啦啦啦啦在线视频资源| 麻豆成人av视频| 一区二区三区精品91| 国产黄a三级三级三级人| 久久精品久久精品一区二区三区| 久久久久国产精品人妻一区二区| 男女无遮挡免费网站观看| 国产老妇女一区| 菩萨蛮人人尽说江南好唐韦庄| 色视频在线一区二区三区| 亚洲欧美一区二区三区国产| 国模一区二区三区四区视频| 女人被狂操c到高潮| 最近最新中文字幕大全电影3| 岛国毛片在线播放| 精品人妻视频免费看| 麻豆国产97在线/欧美| 久久午夜福利片| 中文字幕免费在线视频6| 黄色欧美视频在线观看| 精品熟女少妇av免费看| 看非洲黑人一级黄片| 18禁在线播放成人免费| 国产乱来视频区| 国产精品人妻久久久久久| 看免费成人av毛片| 亚洲av二区三区四区| 2021天堂中文幕一二区在线观| 美女xxoo啪啪120秒动态图| 最近的中文字幕免费完整| av在线app专区| 国产精品一区二区三区四区免费观看| 汤姆久久久久久久影院中文字幕| 99热网站在线观看| 午夜福利视频精品| 欧美成人精品欧美一级黄| 日韩国内少妇激情av| 亚洲成人中文字幕在线播放| 国产成人a∨麻豆精品| 国产成人精品一,二区| 男女边摸边吃奶| 欧美激情在线99| 亚洲精品乱码久久久v下载方式| 1000部很黄的大片| 久久综合国产亚洲精品| 精品亚洲乱码少妇综合久久| 91久久精品电影网| 一级二级三级毛片免费看| 狠狠精品人妻久久久久久综合| 国产有黄有色有爽视频| 精品一区在线观看国产| 中文字幕亚洲精品专区| 国产老妇女一区| 九色成人免费人妻av| 女人被狂操c到高潮| 在线a可以看的网站| 精品午夜福利在线看| 成年免费大片在线观看| 天天一区二区日本电影三级| 99久久精品热视频| 热99国产精品久久久久久7| 人人妻人人看人人澡| 精品久久久久久久久亚洲| 午夜福利视频精品| av在线app专区| 亚洲精品aⅴ在线观看| 亚洲国产欧美人成| 内地一区二区视频在线| 免费看光身美女| 老司机影院毛片| 日韩大片免费观看网站| 国产精品久久久久久精品电影小说 | 精品亚洲乱码少妇综合久久| 99久久九九国产精品国产免费| 成人国产麻豆网| 日韩成人av中文字幕在线观看| 国产精品伦人一区二区| 免费看av在线观看网站| 欧美97在线视频| 国产毛片a区久久久久| 亚洲成人中文字幕在线播放| 日韩成人av中文字幕在线观看| 美女高潮的动态| 亚洲在线观看片| 欧美激情国产日韩精品一区| 久久99精品国语久久久| 日韩成人av中文字幕在线观看| 国产乱来视频区| 91精品一卡2卡3卡4卡| 国产真实伦视频高清在线观看| 亚洲av.av天堂| 日韩欧美精品免费久久| 久久精品夜色国产| 亚洲人成网站在线观看播放| 久久精品人妻少妇| 亚洲内射少妇av| 永久网站在线| 日韩中字成人| 最近中文字幕高清免费大全6| 国产91av在线免费观看| 亚洲人成网站在线观看播放| 午夜福利在线在线| 18禁在线播放成人免费| 特大巨黑吊av在线直播| 国产精品.久久久| 亚洲天堂av无毛| 91精品伊人久久大香线蕉| 极品少妇高潮喷水抽搐| 1000部很黄的大片| 色哟哟·www| 亚洲精品成人av观看孕妇| 又爽又黄a免费视频| 插逼视频在线观看| 国产精品人妻久久久影院| 久久久久精品性色| 中文资源天堂在线| 亚洲天堂av无毛| 免费观看无遮挡的男女| 国产色爽女视频免费观看| 亚洲av日韩在线播放| 久久久国产一区二区| 国产精品一二三区在线看| 少妇人妻精品综合一区二区| 日韩一区二区视频免费看| 午夜日本视频在线| 一级二级三级毛片免费看| 久久人人爽人人片av| 精品酒店卫生间| 亚洲欧美成人精品一区二区| 亚洲精品亚洲一区二区| 91久久精品国产一区二区三区| 国内精品宾馆在线| 秋霞伦理黄片| 男女啪啪激烈高潮av片| 成人国产av品久久久| 性插视频无遮挡在线免费观看| 青青草视频在线视频观看| 国产高清有码在线观看视频| 建设人人有责人人尽责人人享有的 | 婷婷色av中文字幕| 免费看光身美女| 国产精品一及| 一本一本综合久久| 看黄色毛片网站| 欧美日韩在线观看h| 六月丁香七月| 在线观看一区二区三区激情| 神马国产精品三级电影在线观看| 国产高清国产精品国产三级 | 免费高清在线观看视频在线观看| 91在线精品国自产拍蜜月| 欧美成人精品欧美一级黄| 天天一区二区日本电影三级| 久久久久久伊人网av| 久久久久国产精品人妻一区二区| 搞女人的毛片| 成人国产麻豆网| 日本一二三区视频观看| 美女被艹到高潮喷水动态| 自拍偷自拍亚洲精品老妇| 国精品久久久久久国模美| 精品久久久噜噜| 汤姆久久久久久久影院中文字幕| 欧美zozozo另类| tube8黄色片| 日韩av不卡免费在线播放| 黄色怎么调成土黄色| 乱码一卡2卡4卡精品| 亚洲精品日韩在线中文字幕| 欧美激情在线99| 国产乱人偷精品视频| 大话2 男鬼变身卡| 中国国产av一级| 欧美成人午夜免费资源| av网站免费在线观看视频| av又黄又爽大尺度在线免费看| 日韩一区二区三区影片| 免费大片黄手机在线观看| 国产亚洲91精品色在线| 久久久久国产精品人妻一区二区| 岛国毛片在线播放| 亚洲av免费在线观看| 久久人人爽人人片av| 国产精品一区二区在线观看99| 日韩中字成人| 久久久久久国产a免费观看| 亚洲人成网站在线播| 自拍偷自拍亚洲精品老妇| 99久国产av精品国产电影| 全区人妻精品视频| av卡一久久| 国产成人精品福利久久| 一级二级三级毛片免费看| 赤兔流量卡办理| 啦啦啦在线观看免费高清www| 菩萨蛮人人尽说江南好唐韦庄| 2022亚洲国产成人精品| 春色校园在线视频观看| 国产探花极品一区二区| 一边亲一边摸免费视频| av在线app专区| 草草在线视频免费看| 久久99热6这里只有精品| 成人毛片60女人毛片免费| 精品久久久精品久久久| 看十八女毛片水多多多| 美女xxoo啪啪120秒动态图| av专区在线播放| 亚州av有码| 新久久久久国产一级毛片| 最新中文字幕久久久久| 又大又黄又爽视频免费| 99热这里只有是精品在线观看| 国产爱豆传媒在线观看| 一本一本综合久久| h日本视频在线播放| 国精品久久久久久国模美| 看黄色毛片网站| 欧美bdsm另类| 18禁动态无遮挡网站| 免费观看a级毛片全部| 水蜜桃什么品种好| av在线蜜桃| 高清欧美精品videossex| 成人亚洲精品一区在线观看 | 亚洲一区二区三区欧美精品 | 国产爱豆传媒在线观看| 日本一二三区视频观看| 国产精品久久久久久久电影| 欧美成人一区二区免费高清观看| 网址你懂的国产日韩在线| 丝袜美腿在线中文| 亚洲性久久影院| 综合色丁香网| 伦精品一区二区三区| 舔av片在线| 天堂网av新在线| 大码成人一级视频| 午夜精品国产一区二区电影 | 日本wwww免费看| 久久精品国产鲁丝片午夜精品| eeuss影院久久| 蜜桃亚洲精品一区二区三区| 亚洲熟女精品中文字幕| 婷婷色av中文字幕| 2021天堂中文幕一二区在线观| 久久久久性生活片| 欧美亚洲 丝袜 人妻 在线| 狂野欧美白嫩少妇大欣赏| 亚洲av日韩在线播放| 成年人午夜在线观看视频| 制服丝袜香蕉在线| 国产在视频线精品| 看十八女毛片水多多多| 国产免费视频播放在线视频| 看十八女毛片水多多多| av国产精品久久久久影院| 日韩,欧美,国产一区二区三区| 久久99热这里只频精品6学生| 91狼人影院| 国国产精品蜜臀av免费| 亚洲丝袜综合中文字幕| 精品久久久久久久久av| 欧美97在线视频| 国产男女内射视频| 精品久久久久久久人妻蜜臀av| 国产成人福利小说| 日韩一区二区视频免费看| 国产一区二区亚洲精品在线观看| 美女被艹到高潮喷水动态| 2018国产大陆天天弄谢| 真实男女啪啪啪动态图| 视频区图区小说| 亚洲图色成人| 亚洲av.av天堂| 丝袜喷水一区| 国产欧美另类精品又又久久亚洲欧美| 丝袜喷水一区| 各种免费的搞黄视频| 亚洲人成网站高清观看| av又黄又爽大尺度在线免费看| 男女无遮挡免费网站观看| 国产精品三级大全| 小蜜桃在线观看免费完整版高清| 激情五月婷婷亚洲| 欧美高清成人免费视频www| 搡老乐熟女国产| 国产毛片在线视频| 少妇裸体淫交视频免费看高清| 成人高潮视频无遮挡免费网站| 丝袜美腿在线中文| 老司机影院毛片| 国产亚洲精品久久久com| 亚洲电影在线观看av| 老司机影院成人| 中文字幕av成人在线电影| 99热这里只有是精品在线观看| av在线观看视频网站免费| 日本黄色片子视频| 少妇人妻一区二区三区视频| 欧美成人精品欧美一级黄| 亚洲av成人精品一二三区| 极品少妇高潮喷水抽搐| 一级毛片 在线播放| 亚洲精品乱码久久久v下载方式| 国产成人免费观看mmmm| 国产一区有黄有色的免费视频| 菩萨蛮人人尽说江南好唐韦庄| 国产日韩欧美亚洲二区| 欧美精品国产亚洲| 免费av观看视频| 偷拍熟女少妇极品色| 黄色配什么色好看| 久久午夜福利片| av.在线天堂| 狂野欧美白嫩少妇大欣赏| 亚洲欧美精品专区久久| 国产91av在线免费观看| 免费电影在线观看免费观看| 欧美97在线视频| 亚洲精品国产色婷婷电影| 国产一区二区三区综合在线观看 | 亚洲天堂av无毛| 直男gayav资源| 五月玫瑰六月丁香| 久久久欧美国产精品| 欧美性猛交╳xxx乱大交人| 亚洲精华国产精华液的使用体验| 精品少妇久久久久久888优播| 国产精品一区二区在线观看99| 亚洲精品日本国产第一区| 精品人妻偷拍中文字幕| 亚洲最大成人手机在线| 成人毛片60女人毛片免费| 少妇熟女欧美另类| 精品人妻视频免费看| 亚洲av在线观看美女高潮| 国产成年人精品一区二区| 精品一区二区三卡| 91久久精品国产一区二区三区| 久久精品国产亚洲av涩爱| 免费观看无遮挡的男女| 看免费成人av毛片| 国产黄片视频在线免费观看| 亚洲成色77777|