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

    基于噴流擬序結(jié)構(gòu)預(yù)測的SGS 模型比較研究1)

    2021-11-09 06:26:10劉琪麟賴煥新
    力學(xué)學(xué)報(bào) 2021年7期
    關(guān)鍵詞:模態(tài)模型

    劉琪麟 賴煥新

    (華東理工大學(xué)機(jī)械與動(dòng)力工程學(xué)院,上海 200237)

    引言

    擬序結(jié)構(gòu)普遍存在于湍流中[1-3],并影響湍流的發(fā)生、輸運(yùn)與耗散等過程,對(duì)湍流擬序結(jié)構(gòu)的研究、對(duì)于理解湍流機(jī)理、分析流動(dòng)相互作用有重要的意義.從大渦模擬來看,正確預(yù)測湍流擬序結(jié)構(gòu)的關(guān)鍵就在于算法和亞網(wǎng)格尺度(sub-grid-scale,SGS)模型的準(zhǔn)確性.

    在算法方面,傳統(tǒng)的大渦模擬方法求解網(wǎng)格濾波的控制方程,并采用SGS 模型?;瘉喚W(wǎng)格尺度流動(dòng)與網(wǎng)格尺度流動(dòng)之間的相互作用.雖然近年來部分學(xué)者利用數(shù)值離散的特點(diǎn),也發(fā)展出不需要添加模型的隱式大渦模擬方法[4-5],但為了針對(duì)性地討論應(yīng)用中亞網(wǎng)格模型的選擇問題,本文仍采用高精度有限差分結(jié)合亞網(wǎng)格模型的大渦模擬經(jīng)典算法,并討論SGS 模型的性能.對(duì)大渦模擬算法更全面的介紹可以參考Sagaut 等[6]的綜述.

    目前文獻(xiàn)中廣泛出現(xiàn)的亞網(wǎng)格尺度模型是經(jīng)典的常系數(shù)Smagorinsky 模型(Smagorinsky model,SM),已有眾多文獻(xiàn)報(bào)道該模型過度耗散了脈動(dòng)動(dòng)能[7-9].另一類是基于Germano 恒等式和Lilly 最小二乘法來動(dòng)態(tài)調(diào)整模型系數(shù)的動(dòng)態(tài)Smagorinsky 模型(dynamic Smagorinsky model,DSM),這類模型的局限性在于理論上要求在流場均勻的方向上計(jì)算平均的模型系數(shù),而這一平均過程缺少物理上合理的解釋,且造成分布式并行計(jì)算中不同處理器之間數(shù)據(jù)的頻繁調(diào)用與交換,降低計(jì)算效率.從計(jì)算的角度,較為理想的模型是局部模型,計(jì)算效率高且能較好地克服SM 對(duì)湍流黏性估計(jì)過高的缺點(diǎn).Piomelli和Liu[9]提出了局部動(dòng)態(tài)模型(localized dynamic Smagorinsky model,LDSM),去除了DSM 對(duì)流動(dòng)中需存在各向同性方向的要求,而是利用流動(dòng)的歷史信息來避免平均過程,計(jì)算效率得到提高.Lenormand等[10]提出的選擇多尺度模型(selective mixed-scale model,SMSM),在邊界層流動(dòng)[10]、翼型繞流[11]和空腔流動(dòng)[12]的預(yù)測中具有良好應(yīng)用.這一模型借助渦矢量方向的變化來判斷流動(dòng)的狀態(tài)是否發(fā)生了轉(zhuǎn)捩,同時(shí)反映湍流的間歇性.Kobayashi 等[13-14]則利用流場中黏性耗散與擬序結(jié)構(gòu)的關(guān)聯(lián)性,提出了擬序結(jié)構(gòu)模型(coherent-structure Smagorinsky model,CSM),并應(yīng)用于各向同性湍流、湍流通道流以及橫向射流的計(jì)算中,其模型簡單且具有類似于動(dòng)態(tài)模型的性能.Koyabashi[13]同時(shí)改進(jìn)動(dòng)能模型,提出擬序結(jié)構(gòu)動(dòng)能模型(coherent-structure kinetic-energy model,CKM).這些局部模型都有望克服傳統(tǒng)模型的不足,但對(duì)于預(yù)測可壓縮噴流中的擬序結(jié)構(gòu),目前尚未看到有關(guān)它們預(yù)測性能的分析與評(píng)估的文獻(xiàn).

    典型的噴流擬序結(jié)構(gòu)有流向渦、高低速流體間的條帶[15]、剪切層邊沿的卷吸流動(dòng)結(jié)構(gòu)[16]等,表現(xiàn)為長時(shí)間保持一定形狀的大尺度質(zhì)量團(tuán),且在時(shí)間和空間上具有強(qiáng)相關(guān)性[17].Lumley[18]最早提出采用本征正交分解(proper orthogonal decomposition,POD)方法提取擬序結(jié)構(gòu),且在自由剪切流、分離流、邊界層轉(zhuǎn)捩以及湍流通道流中得到廣泛的應(yīng)用.POD 方法的目標(biāo)是針對(duì)物理場尋找一組空間正交的基函數(shù),并按照基函數(shù)對(duì)物理量方差的貢獻(xiàn)排序,例如對(duì)于脈動(dòng)速度場的研究就是基函數(shù)對(duì)湍動(dòng)能的貢獻(xiàn)[19-20].本征正交分解的特點(diǎn)使得只需要少量的POD 模態(tài)就能夠提取出流場中能量占主導(dǎo)地位的流動(dòng)模式.此外,POD 分解的基函數(shù)直接由湍流場確定,不同于傅里葉變換或者小波變換需要人為指定基函數(shù),因此表征的擬序結(jié)構(gòu)更為客觀[17].在POD 算法方面,雖然快照POD 算法能夠顯著減少計(jì)算量,但是有部分高階模態(tài)沒有參與計(jì)算.基于矩陣奇異值分解(singular value decomposition,SVD)的經(jīng)典POD 算法包含完整的模態(tài)信息,比快照POD算法對(duì)舍入誤差的魯棒性高[19-20].此外數(shù)學(xué)軟件Matlab 中提供了SVD 的庫函數(shù),也便于編制后處理程序.

    本文對(duì)亞聲速可壓縮平行噴流進(jìn)行大渦模擬.在驗(yàn)證計(jì)算結(jié)果的基礎(chǔ)上,研究SM,CKM,SMSM,LDSM 和CSM 模型的特性,分析各個(gè)模型預(yù)測的瞬時(shí)渦結(jié)構(gòu).同時(shí)采用基于SVD 的POD 分解方法提取噴流擬序結(jié)構(gòu)并分析各模型預(yù)測結(jié)果的差異,為可壓縮噴流的大渦模擬提供參考.

    1 數(shù)值計(jì)算與方法

    1.1 亞網(wǎng)格尺度模型

    控制方程組為Favre 質(zhì)量加權(quán)濾波的三維非穩(wěn)態(tài)可壓縮N-S 方程組[21].方程中的亞網(wǎng)格尺度未知量由亞網(wǎng)格模型給出,包括亞網(wǎng)格尺度應(yīng)力τij

    以及亞網(wǎng)格尺度熱通量

    其中ρ,T,νsgs以及ui(i=1,2,3) 分別為密度、溫度、模化的亞網(wǎng)格黏性系數(shù)以及速度矢量的3 個(gè)笛卡爾分量.γ 是量熱完全氣體的比熱比,空氣取γ=1.4.M為馬赫數(shù),湍流普朗特?cái)?shù)常數(shù)Prt=1[22],而,頂標(biāo)“-”表示網(wǎng)格濾波,“~”表示Favre 質(zhì)量加權(quán)的網(wǎng)格濾波,所有計(jì)算變量均為網(wǎng)格濾波后的變量,因此除非特殊說明以下省去頂部符號(hào).

    Smagorinsky 模型(SM)中的亞網(wǎng)格尺度應(yīng)力為

    δij=0.對(duì)于噴流,取模型系數(shù)C1=0.17[23].τkk的?;鶕?jù)文獻(xiàn)[7]的建議得出.其中=Li j?δi jLkk/3,頂標(biāo)“^”表示測試濾波運(yùn)算,Ct為上一時(shí)刻的模型系數(shù),初始值可以根據(jù)SM 模型取Ct=C1,當(dāng)前時(shí)刻的模型系數(shù)C2是未知量,且式(5)對(duì)C2是超定的.C2的估計(jì)值由式(5)兩邊取αij的縮并得出[9],如式(6)所示.LDSM 模型的亞網(wǎng)格尺度應(yīng)力與亞網(wǎng)格黏性系數(shù)通過將式(4)中的替換為得出

    選擇多尺度模型(slective mixed-scale model,SMSM)是根據(jù)局部渦矢量方向的脈動(dòng)來調(diào)整亞網(wǎng)格黏性系數(shù)

    其中C3=0.06,α 取0.5[24].為了反映湍流的多尺度特性,引入小尺度脈動(dòng)的動(dòng)能=(1/2),其中脈動(dòng)速度=?.進(jìn)一步,對(duì)網(wǎng)格濾波尺度的瞬時(shí)渦量進(jìn)行測試濾波,得出.和兩個(gè)向量之間的夾角 θ 反映了渦運(yùn)動(dòng)的脈動(dòng)強(qiáng)弱.對(duì)于均勻各向同性湍流,在夾角 θ 的概率密度函數(shù)的峰值位置處有θ0=20°[24],據(jù)此,可判斷局部流動(dòng)的狀態(tài).模型引入了選擇函數(shù)fθ0,當(dāng) θ> θ0時(shí),取fθ0=1,表示渦方向在兩個(gè)濾波尺度上變化較大,是無規(guī)則的完全湍流狀態(tài);否則,當(dāng) θ<θ0時(shí),取fθ0=tan(θ/2)/tan(θ0/2),意味著湍流逐步減弱,亞網(wǎng)格黏性也應(yīng)相應(yīng)減弱,直到在層流區(qū)僅保留分子黏性.因此,SMSM 模型的亞網(wǎng)格黏性系數(shù)為

    擬序結(jié)構(gòu)模型是利用擬序結(jié)構(gòu)與湍流中耗散的相關(guān)性來調(diào)整亞網(wǎng)格黏性系數(shù)vsgs=CΔ2|S(u)|,模型系數(shù)C=C4|FCS|3/2中常數(shù)C4=0.05[13],FCS=Q/E為擬序結(jié)構(gòu)函數(shù),其中

    為速度導(dǎo)數(shù)Dij=?ui/?xj的第二不變量,它主要由Di j中反映當(dāng)?shù)匦D(zhuǎn)的反對(duì)稱部分Aij和反映當(dāng)?shù)刈冃蔚膶?duì)稱部分Bi j組成,而為Di j的范數(shù).其中

    基于直接數(shù)值模擬(direct numerical simulation,DNS),流場的相關(guān)性分析表明,Q的絕對(duì)值與擬序小尺度渦附近能量耗散的強(qiáng)弱存在關(guān)聯(lián)[13],因此|FCS|反映了流場中黏性與耗散的局部變化,而|FCS|的冪次“3/2”是根據(jù)不可壓縮壁面流動(dòng)中C∝y3,Q∝y2以及E∝ 常數(shù)(y是壁面的法向)而確定的.因此,CSM 模型的亞網(wǎng)格黏性系數(shù)為

    擬序結(jié)構(gòu)動(dòng)能模型是在擬序結(jié)構(gòu)模型的基礎(chǔ)上,在亞網(wǎng)格黏性系數(shù)中顯式地包含亞網(wǎng)格尺度動(dòng)能ksgs,并做量綱調(diào)整后得到

    至此,可以通過式(2)、式(4)來計(jì)算亞網(wǎng)格尺度熱通量和亞網(wǎng)格尺度應(yīng)力.雖然文獻(xiàn)[12]指出亞網(wǎng)格尺度應(yīng)力中各向同性的部分 τkk相對(duì)于熱力學(xué)壓力可以忽略,但根據(jù)Vreman[7]的建議,在基于SM 的LDSM 以及CSM 模型中仍然?;?τkk.

    1.2 計(jì)算對(duì)象與數(shù)值離散

    本文采用Hu 等[22]的亞聲速平行噴流為考核算例.噴流出口雷諾數(shù)ReJ=DJρ∞UJ/μ∞=2000,馬赫數(shù)MJ=UJ/c∞=0.9,其中DJ=1、UJ=1 分別為歸一化的噴口寬度和噴口速度,參考值ρ∞,μ∞和c∞分別為來流密度、動(dòng)力黏度以及環(huán)境聲速.計(jì)算域與網(wǎng)格劃分如圖1 所示,流向(x)、橫向(y)與展向(z)的尺寸分別為Lx×Ly×Lz=24×15×3,采用笛卡爾網(wǎng)格進(jìn)行離散,對(duì)應(yīng)方向的結(jié)點(diǎn)數(shù)分別為 181×181×16.為了提高近噴口區(qū)域和噴流唇線附近的網(wǎng)格分辨率,采用雙曲正切函數(shù)對(duì)x和y方向的網(wǎng)格進(jìn)行拉伸,而噴流在z方向采用周期邊界條件,劃分均勻的網(wǎng)格.基于Celik 等[25]提出的大渦模擬分辨率指標(biāo),由于課題組前期工作已證明該網(wǎng)格能夠分辨湍動(dòng)噴流中75%以上的湍動(dòng)能,考慮到噴流的雷諾數(shù)較低,該計(jì)算網(wǎng)格可以滿足大渦模擬的網(wǎng)格質(zhì)量要求[26].

    圖1 計(jì)算域與網(wǎng)格(3 個(gè)方向都每隔3 條網(wǎng)格線畫一條線)Fig.1 Computional domain and the Grid(One in every three lines is shown in each direction)

    采用有限差分方法開展數(shù)值模擬,本文的計(jì)算源程序以Sandham 課題組[27-28]的(shock-boundary layer interaction,SBLI)代碼為基礎(chǔ),計(jì)算精度和可靠性經(jīng)過了充分的考核驗(yàn)證[22,29].在此基礎(chǔ)上,本文進(jìn)一步采用OpenMP 并行化,使該代碼可用于內(nèi)存共享式的小型計(jì)算工作站.

    空間離散的1 階與2 階導(dǎo)數(shù)均采用5 點(diǎn)4 階中心差分格式計(jì)算,邊界點(diǎn)采用Carpenter 提出的單側(cè)4 階穩(wěn)定差分格式[30],從而保持空間離散的整體高精度.為了抑制差分計(jì)算引起的非物理振蕩,算法上采用了熵分裂方法,通過將對(duì)流項(xiàng)分裂為守恒形式與非守恒形式的加權(quán)求和來實(shí)現(xiàn)[31].同時(shí),在x和y方向的計(jì)算邊界上,采用無反射邊界條件來抑制非物理振蕩波在計(jì)算域邊界上的反射,算法上通過保留走出計(jì)算域的特征波而將走入計(jì)算域的特征波幅值置零來實(shí)現(xiàn)[32].在z方向上采用周期邊界條件來模擬平行噴流在展向的無限延伸.

    時(shí)間離散采用低存儲(chǔ)顯式3 階龍格庫塔算法,這一算法完成3 個(gè)子時(shí)間步的推進(jìn)僅需要兩個(gè)存儲(chǔ)流場的數(shù)組,對(duì)內(nèi)存資源的消耗低[21].為便于模型特性的研究比較,時(shí)間步長取定值Δt=0.003 075,且滿足數(shù)值穩(wěn)定性條件CFL<1.在初始條件方面,給定噴口平面上流向速度u、溫度T以及壓強(qiáng)p的形線,并用噴口平面上的物理量來初始化全場.噴口的速度分布采用雙曲正切函數(shù)給定,并施加伴流速度U0=0.1,可得

    在數(shù)值結(jié)果方面,雖然DNS 能夠準(zhǔn)確給出流場結(jié)構(gòu)的信息,但是DNS 對(duì)計(jì)算資源的消耗量高,而且流場后處理的數(shù)據(jù)量大.以文獻(xiàn)[22-33]中提供的DNS 數(shù)據(jù)為基礎(chǔ),課題組前期工作已開展了計(jì)算代碼的考核[34-35]、LES 的網(wǎng)格分辨率分析[26]以及局部模型的研究[21,36],算法的可靠性與模型都已得到驗(yàn)證.本文在已有工作的基礎(chǔ)上,進(jìn)一步分析局部模型揭示的噴流擬序結(jié)構(gòu).

    1.3 POD 方法

    POD 方法采用一組正交基函數(shù) ψ(x) 的線性組合來表示物理場 φ(x,t),組合系數(shù)為c(t),即

    所得的基函數(shù)也稱為POD 模態(tài).POD 分解的目的是使這種表示的誤差在平方平均的意義下最小,并且在給定誤差下最小化表示所需要的模態(tài)個(gè)數(shù),從而得到最佳的基函數(shù)[20].這組最佳基函數(shù)能夠最優(yōu)地表征原物理量的方差,對(duì)于脈動(dòng)速度場,這種最優(yōu)是模態(tài)反映的湍動(dòng)能隨模態(tài)階數(shù)增大而迅速的減小[19],脈動(dòng)速度的第一階POD 模態(tài)對(duì)湍動(dòng)能的貢獻(xiàn)最大.詳細(xì)的理論介紹可參考文獻(xiàn)[20].

    在計(jì)算中,可對(duì)流場的速度分量進(jìn)行時(shí)間采樣,得到關(guān)于標(biāo)量 φ 的N個(gè)流場樣本,對(duì)于網(wǎng)格點(diǎn)數(shù)為M=Mx×My的二維樣本,可將每個(gè)樣本都寫成M×1的列向量,由此得到M×N的樣本矩陣

    在離散求解時(shí),尋找最佳基函數(shù)的過程在數(shù)學(xué)上等價(jià)于求樣本矩陣F的SVD 問題[20].采用Matlab 中低存儲(chǔ)的“svd()”庫函數(shù)計(jì)算F的奇異值分解

    其中,左矩陣 Ψ 的各列是正交的,Ψ=(ψ1,ψ2,···,ψN),它的第i列 ψi即為物理場的第i階POD 模態(tài),右矩陣C的各列也是正交的,C=(c1,c2,···,cN),它的第i列ci即為模態(tài) ψi對(duì)應(yīng)的系數(shù),反映POD 模態(tài)與時(shí)間有關(guān)的特征[19],而S為奇異值矩陣,S=dig(s1,s2,···,sN),S的對(duì)角元為從大到小排列的奇異值si(i=1,2,···,N),其他矩陣元為零,λi=代表了與第i階POD 模態(tài)對(duì)應(yīng)的能量,且對(duì)各階模態(tài)可以定義能量貢獻(xiàn)率

    以及累計(jì)能量貢獻(xiàn)率

    另外,由于各階POD 模態(tài) ψi是相互正交的,所以λ衰減越快,當(dāng)前模態(tài)所含有的能量相比下一階模態(tài)的差就越大,當(dāng)前模態(tài)所表征的物理場的相關(guān)性就越強(qiáng).

    2 SGS 模型的特性

    2.1 平均流場

    在分析湍流統(tǒng)計(jì)特征時(shí),為排除初場對(duì)統(tǒng)計(jì)結(jié)果的影響,取無量綱時(shí)間t=48~ 216 的流場數(shù)據(jù),對(duì)應(yīng)流體以噴口速度流過2~ 9 倍的計(jì)算域長度,以計(jì)算時(shí)間步長Δt為采樣間隔,得到54 640 個(gè)流場樣本.圖2 以LDSM 模型為例,顯示了在唇線上無量綱位置x=2,6,8,10,20 處的流向脈動(dòng)速度的功率譜,譜線具有一段固定的斜率且能觀察到慣性區(qū),反映了湍流的能量級(jí)串現(xiàn)象,也說明流動(dòng)已經(jīng)達(dá)到了完全湍動(dòng)狀態(tài).

    圖2 LDSM 模型預(yù)測的脈動(dòng)速度無量綱功率譜Fig.2 Dimensionless power spectra of velocity signals predicted by the LDSM

    噴流主流流體與環(huán)境流體存在速度差,形成剪切層.圖3 所示為噴流剪切層的發(fā)展.考慮噴流中心線速度Uc與伴流速度U0,圖3(a)采用速度半值寬δ0.5,即流向速度u= 0.5(Uc+U0) 時(shí)的橫向位置y,顯示出剪切層的發(fā)展呈現(xiàn)兩個(gè)線性階段,即初期的緩慢增長段與下游的快速增長段,這對(duì)應(yīng)噴流勢流核末端x≈6 處的流動(dòng)轉(zhuǎn)捩現(xiàn)象.對(duì)比5 種模型的結(jié)果,SM 模型預(yù)測的剪切層橫向發(fā)展程度較其他模型大,表明該模型過度耗散的特點(diǎn).圖3(b)以LDSM 模型為例,在相似坐標(biāo)中繪出了噴流不同流向位置的速度形線,相似坐標(biāo)的橫向位置采用速度半值寬歸一化,平均噴流速度差〈u〉?Uc采用主、伴流的速度差ΔUc=Uc?U0歸一化,〈·〉 表示時(shí)間平均.圖3(b)中,速度形線在入口呈平頂狀分布,而在下游則快速發(fā)展為自相似的速度形線,可以看到x=4 時(shí),自相似形線還沒有建立,但當(dāng)x≥8 時(shí),速度形線與實(shí)驗(yàn)噴流[37-39]的自相似數(shù)據(jù)點(diǎn)吻合良好,也驗(yàn)證了計(jì)算結(jié)果.

    圖3 噴流剪切層的發(fā)展Fig.3 Developments of the jet shear layer

    噴流流體與環(huán)境流體的強(qiáng)烈摻混伴隨著能量的耗散.對(duì)比5 種SGS 模型的亞網(wǎng)格尺度平均黏性耗散率,定義為.圖4 為在展向中心平面(z=1.5)上的分布.由于噴流入口施加的平均來流條件,流體噴出后在近噴口區(qū)基本為層流狀態(tài),隨后在勢流核末端(x≈6)附近發(fā)生流動(dòng)轉(zhuǎn)捩,并在下游發(fā)展為充分湍流.圖2 中,近噴口的脈動(dòng)水平比勢流核及其下游的區(qū)域低1~ 2 個(gè)數(shù)量級(jí),進(jìn)一步表明來流的層流特征,此處與湍流脈動(dòng)相關(guān)的亞網(wǎng)格尺度黏性耗散應(yīng)趨近于零.但在圖4(a)中,SM 模型預(yù)測的在近噴口區(qū)內(nèi)呈兩條平行的條帶,表明SM 模型對(duì)層流擾動(dòng)敏感,不能正確反映湍流黏性耗散率的分布區(qū)域.圖4(b)~ 圖4(e)分別為CKM,SMSM,LDSM 以及CSM 的預(yù)測結(jié)果,這4 種局部模型雖然能夠正確地反映近噴口區(qū)的分布特征,在轉(zhuǎn)捩區(qū)(x=6~9)內(nèi)呈現(xiàn)一對(duì)峰值,但是以CSM 模型為參考,CKM,LDSM,SMSM 以及SM 模型預(yù)測的峰值耗散不同,分別為CSM 模型預(yù)測值的3.0,1.8,1.3 和0.9 倍.其中CKM 預(yù)測的峰值接近其他模型的1~ 2 倍,稍后分析的分布對(duì)轉(zhuǎn)捩區(qū)內(nèi)渦結(jié)構(gòu)的影響.

    圖4 對(duì)比在中心平面(z=1.5)上的分布Fig.4 Comparisons of the distributi on of ?in the center plane of the jet (z=1.5)

    2.2 瞬時(shí)流場的渦結(jié)構(gòu)

    分析無量綱時(shí)間t=207 的瞬時(shí)流場,此時(shí)噴流以噴口速度流過8 倍以上的計(jì)算域長度,已經(jīng)充分排除初場的影響.參考Hu 等[22]的部分DNS 結(jié)果,以噴流中心線速度〈u〉 衰減為0.99 時(shí)的流向位置得出平均的勢流核長度為x=6.43.圖5 為Q=0.1 的等值面,圖5 中紅色切片位于Hu 等DNS 結(jié)果[22]的勢流核末端.在預(yù)測流動(dòng)轉(zhuǎn)捩方面,圖5(c)~ 圖5(e)中,SMSM,LDSM 以及CSM 模型都能預(yù)測大尺度展向渦在平均勢流核末端附近的破碎,以及下游區(qū)域內(nèi)湍流多尺度渦結(jié)構(gòu)的生成.然而對(duì)于SM 模型,圖5(a)在切片位置下游仍然能夠觀察到流向的大尺度渦結(jié)構(gòu)(如圖5 中箭頭所示),結(jié)合圖4(a)可知,湍流中小尺度渦的發(fā)展可能受到近噴口過度的亞網(wǎng)格尺度耗散的抑制,這與文獻(xiàn)[7-9]對(duì)SM 模型過度耗散特點(diǎn)的報(bào)道一致.圖5(b)中,CKM 模型的預(yù)測結(jié)果與SM 類似,同樣可能由于該模型中局部亞網(wǎng)格尺度耗散較強(qiáng).圖6 為展向渦 ωz的等值面,同樣采用紅色切片標(biāo)記DNS 結(jié)果[22]的勢流核末端位置.可以看到,ωz的等值面在勢流段內(nèi)呈平行的渦片,隨著流動(dòng)向下游的發(fā)展,ωz等值面在勢流核末端附近快速卷起,隨后渦結(jié)構(gòu)被拉伸并在下游發(fā)展出復(fù)雜的三維結(jié)構(gòu).圖6(a)中,SM 模型預(yù)測的渦卷起位置位于切片之后,渦片偏長,對(duì)應(yīng)圖4(a)中近噴口區(qū)的不合理耗散,SM 模型抑制了流場中展向渦的發(fā)展.圖6(b),圖6(d),圖6(e) 中,CKM,LDSM 以及CSM 模型在勢流核末端位置附近都能夠觀察到展向渦的卷起(如圖6 中箭頭標(biāo)記)以及下游的渦破碎,但LDSM 與CSM 模型在切片后預(yù)測的小尺度渦結(jié)構(gòu)比CKM 模型豐富,再次表明CKM 模型的局部強(qiáng)耗散可對(duì)小尺度渦結(jié)構(gòu)產(chǎn)生抑制.圖6(c) 中,SMSM 預(yù)測的渦卷起不顯著,但切片下游的預(yù)測結(jié)果與LDSM,CSM類似.相比SM 與CKM 模型,SMSM,LDSM 以及CSM 模型更能反映勢流核下游展向渦的多尺度特性.

    圖5 對(duì)比 Q=0.1 的等值面,采用流向速度著色,紅色切片位置標(biāo)記了文獻(xiàn)[22]中DNS 結(jié)果的勢流核末端:x=6.43Fig.5 Comparisons of the iso-surface of Q=0.1,colored by streamwize velocity.The red slice marks the end of the potential core x=6.43 of the DNS results in Ref.[22]

    圖6 對(duì)比 ωz 的等值面,紅色切片位置標(biāo)記了文獻(xiàn)[22]中DNS 結(jié)果的勢流核末端:x=6.43Fig.6 Comparisons of the iso-surface of ωz,the red slice marks the end of the potential core x=6.43 of the DNS results in Ref.[22]

    2.3 速度場的POD 模態(tài)分析

    對(duì)5 種SGS 模型的預(yù)測結(jié)果進(jìn)行POD 分解.平行噴流主要表現(xiàn)為x-y平面內(nèi)的二維流動(dòng)特征,因此在湍流統(tǒng)計(jì)的流場數(shù)據(jù)中,截取三維流場的中心平面(z=1.5)為樣本,以10Δt(=0.030 75) 為采樣間隔,得到4685 個(gè)樣本.針對(duì)5 種SGS 模型,提取脈動(dòng)速度分量u′(=u?〈u〉),v′(=v?〈v〉) 和w′(=w?〈w〉)的POD 模態(tài).脈動(dòng)速度u′,v′和w′的第一階模態(tài)分別對(duì)原始脈動(dòng)速度場的脈動(dòng)強(qiáng)度:〈u′u′〉,〈v′v′〉 和〈w′w′〉的貢獻(xiàn)最大,并在這一意義下脈動(dòng)速度的第一階模態(tài)與原始脈動(dòng)速度場最為接近.圖7 為流向脈動(dòng)速度u′模態(tài)的累計(jì)能量貢獻(xiàn)率 η,可以看到,累計(jì)能量貢獻(xiàn)率曲線在前200 階內(nèi)就達(dá)到了90%以上,這200 階模態(tài)僅為總模態(tài)數(shù)的4%,因此曲線 η 反映了POD 的快速收斂特性.但5 種SGS 模型的 η 具有不同的增長速度,SM 模型最快,而SMSM 模型最慢.脈動(dòng)速度v′與w′的 η 曲線具有與圖7 相似的快速收斂特征,因此不再繪出.圖8 為脈動(dòng)速度分量在前20 階模態(tài)的能量貢獻(xiàn)率衰減曲線 ξ.可以看到,5 種SGS 模型得出的 ξ 曲線都在前5 階內(nèi)迅速衰減,說明低階模態(tài)占主導(dǎo)且模態(tài)表征的流場具有強(qiáng)的相關(guān)性.對(duì)比圖8 中第1 階模態(tài)的能量貢獻(xiàn)率,圖8(a)中u′模態(tài)的衰減曲線顯示SM 模型的第1 階模態(tài)貢獻(xiàn)率最高,SMSM 模型最低;圖8(b)中橫向速度脈動(dòng)v′模態(tài)的衰減曲線顯示CKM 模型的貢獻(xiàn)率最高,而SMSM 模型依然排最低;圖8(c)展向速度脈動(dòng)w′模態(tài)的衰減曲線中SMSM 模型的能量貢獻(xiàn)率排序有所提高.圖8 表明,第1 階模態(tài)的能量貢獻(xiàn)率已經(jīng)體現(xiàn)出了不同模型之間的差異,接下來針對(duì)第1 階模態(tài)展開分析.

    圖7 u′ 模態(tài)的累計(jì)能量貢獻(xiàn)率的收斂曲線Fig.7 Convergency curves of the cumulative energy contribution for the modes of u′

    噴流主要表現(xiàn)為流向的動(dòng)量傳遞.對(duì)比5 種SGS 模型預(yù)測的流向速度脈動(dòng)u′的第一階POD 模態(tài),圖9 顯示了模態(tài)幅值在中心平面(z=1.5)上的分布,采用接近零的模態(tài)幅值等值線u′=±0.003 來顯示模態(tài)的輪廓.圖9(a)中,SM 得到兩條沿流向的帶狀模態(tài),且在噴流下游x=12~24 的范圍內(nèi)一直保持高的模態(tài)幅值,對(duì)應(yīng)圖5(a)中大尺度流向渦結(jié)構(gòu)對(duì)流場的攪動(dòng).圖9(b)中,CKM 模型預(yù)測的模態(tài)也在下游保持沿流向的高幅值帶狀模態(tài),在轉(zhuǎn)捩區(qū)(x=6~9)內(nèi)的模態(tài)強(qiáng)度呈現(xiàn)一對(duì)低谷區(qū),位置恰好匹配圖4(b)中的峰值耗散區(qū)位置.這表明流向速度脈動(dòng)的POD 模態(tài)對(duì)亞網(wǎng)格尺度的黏性耗散率敏感,局部的強(qiáng)耗散導(dǎo)致脈動(dòng)強(qiáng)度〈u′u′〉 的降低,從而改變了POD 模態(tài)的形狀.圖9(d)、圖9(e)分別為LDSM與CSM 模型預(yù)測的模態(tài),兩者都沿流向兩側(cè)發(fā)散開來,并且在噴流下游呈現(xiàn)分岔或破碎的形態(tài),模態(tài)幅值也逐漸減小,反映下游的多尺度結(jié)構(gòu)以及脈動(dòng)能量的逐漸消散.圖9(c)為SMSM 模型預(yù)測的模態(tài),模態(tài)輪廓在流向上呈斷續(xù)的條狀,局部的模態(tài)幅值低,這可能與圖8(a)中SMSM 模型第一階POD 模態(tài)的能量貢獻(xiàn)率低有關(guān).而在下游,SMSM 的模態(tài)呈現(xiàn)與CSM,LDSM 模型預(yù)測結(jié)果相似的分岔形態(tài).因此,與SM 和CKM 模型相比,CSM,LDSM 以及SMSM 模型都合理地反映了u′模態(tài)的形態(tài).

    圖8 脈動(dòng)速度分量的POD 模態(tài)的能量貢獻(xiàn)率Fig.8 Energy contribution rate of POD modes for the fluctuate velocity components

    圖9 流向脈動(dòng)速度 u′ 的第一階POD 模態(tài)云圖,實(shí)、虛等值線分別表示 u′=+0.003和u′=?0.003Fig.9 Contours of the first-order POD modes of u′,with solid contour lines for u′=+0.003 and dashed contour lines foru′=?0.003

    圖10 為橫向脈動(dòng)速度v′的第一階POD 模態(tài),等值線v′=±0.003 顯示出模態(tài)的輪廓.同時(shí),圖10中采用面內(nèi)矢量 (u′,v′) 的POD 模態(tài)來顯示噴流中的流動(dòng)模式.觀察圖10,v′的POD 模態(tài)呈正負(fù)相間的肋狀,并沿流向排列.在肋狀模態(tài)上下端部的附近,矢量箭頭顯示出環(huán)形的流動(dòng)模式,在噴流上部的紅圓標(biāo)出了逆時(shí)針環(huán)流,下部的藍(lán)圓標(biāo)出了順時(shí)針環(huán)流.環(huán)流模式穿過了主流與環(huán)境流體,在圓環(huán)的下游側(cè),高速的主流流體被送入兩側(cè)的低速環(huán)境流體中,而在環(huán)流的上游側(cè),低速的環(huán)境流體被帶入主流流體中.同時(shí),v′模態(tài)輪廓的橫向尺寸也沿著流向逐漸增大.可見,(u′,v′) 的POD 模態(tài)反映了噴流的流動(dòng)卷吸現(xiàn)象[16].對(duì)比圖10(a)~ 圖10(e),SM,SMSM,LDSM 以及CSM模型預(yù)測的模態(tài)都反映了噴流的卷吸現(xiàn)象,其中SMSM 模型預(yù)測的旋渦尺度小,這與圖8(a) 和圖8(b) 中的低模態(tài)貢獻(xiàn)率對(duì)應(yīng).而圖10(b)中,CKM模型預(yù)測的模態(tài)基本沒有反映該區(qū)域內(nèi)的流動(dòng)卷吸,結(jié)合圖4(b),脈動(dòng)速度分量u′和v′可能受到CKM模型在x=6 附近局部的強(qiáng)耗散抑制,導(dǎo)致未能預(yù)測出明顯的流動(dòng)卷吸.

    圖10 橫向脈動(dòng)速度v'的第一階POD 模態(tài)云圖與面內(nèi)矢量 (u′,v′) 的第一階POD 模態(tài)矢量圖(實(shí)、虛等值線分別表示v′=+0.003和 v′=?0.003,紅圈標(biāo)示逆時(shí)針環(huán)流、藍(lán)圈標(biāo)示順時(shí)針環(huán)流)Fig.10 Contours of the first-order POD modes of v' and the in plane vector (u′,v′)(Solid contour lines for v′=+0.003 and dashed contour lines for v′=?0.003.Red circle marks counterclockwise circular flow and blue circle indicates clockwise circular flow)

    圖11 為展向脈動(dòng)速度w′的POD 模態(tài)云圖,等值線w′=±0.003 顯示出模態(tài)的輪廓.圖11(b)~圖11(e)分別為CKM,SMSM,LDSM 以及CSM 模型預(yù)測的結(jié)果,w′的POD 模態(tài)在轉(zhuǎn)捩區(qū)(x=3~9)內(nèi)呈現(xiàn)規(guī)則的脊?fàn)钆帕?紅色表示脈動(dòng)速度w′指向平面外,藍(lán)色表示脈動(dòng)速度w′指向平面內(nèi),強(qiáng)弱相間的分布反映了w′模態(tài)對(duì)渦結(jié)構(gòu)的展向拉伸,如圖11(f)所示.圖11 中的脊?fàn)钅B(tài)顯示渦拉伸主要發(fā)生在轉(zhuǎn)捩區(qū),而在下游的模態(tài)輪廓呈碎塊狀,強(qiáng)度也逐漸減弱,對(duì)應(yīng)渦結(jié)構(gòu)逐漸發(fā)展為三維結(jié)構(gòu),展向的拉伸現(xiàn)象不再顯著.圖11(a)中,SM 模型預(yù)測的模態(tài)沒有反映轉(zhuǎn)捩區(qū)內(nèi)的脊?fàn)钆帕?在下游呈大團(tuán)塊,位置與圖4(a) 中SM 模型噴口附近的過度耗散區(qū)域相對(duì)應(yīng).圖11(b)中,CKM 模型預(yù)測的模態(tài)雖然呈現(xiàn)了脊?fàn)钆帕?但是模態(tài)輪廓在下游x=21 處仍然呈大尺度塊狀,與SM 模型的預(yù)測結(jié)果相似.相比而言,SMSM,LDSM 以及CSM 模型的預(yù)測結(jié)果則能夠合理地反映展向拉伸的流動(dòng)模式.

    圖11 展向脈動(dòng)速度 w′ 的第一階POD 模態(tài)云圖Fig.11 Contours of the first-order POD modes of w′,with solid contour lines for w′=+0.003and dashed contour lines forw′=?0.003

    圖11 展向脈動(dòng)速度 w′ 的第一階POD 模態(tài)云圖 (續(xù))Fig.11 Contours of the first-order POD modes of w′,with solid contour lines for w′=+0.003and dashed contour lines for w′=?0.003 (continued)

    在計(jì)算效率方面,表1 對(duì)比了5 種SGS 模型的CPU 時(shí)間,以SM 模型為基準(zhǔn),LDSM,SMSM,CKM,CSM 的計(jì)算時(shí)間分別為SM 模型的:1.64,1.28,1.22 和1.00 倍,其中LDSM,SMSM 以及CKM 模型需要進(jìn)行濾波運(yùn)算,因此需要更多的CPU 時(shí)間,CSM 模型的計(jì)算時(shí)間最接近SM 模型.

    表1 每時(shí)間步內(nèi)每網(wǎng)格點(diǎn)所占用的CPU 時(shí)間Table 1 CPU time per grid point per time step

    3 結(jié)論

    本文開展了亞聲速可壓縮平行噴流的大渦模擬.在分析平均流動(dòng)與耗散的基礎(chǔ)上,針對(duì)瞬時(shí)渦流以及脈動(dòng)速度POD 模態(tài)表征的擬序結(jié)構(gòu),對(duì)SM,CKM,SMSM,LDSM 和CSM,5 種SGS 模型的性能進(jìn)行研究,得出以下主要結(jié)論:

    (1) 初步揭示了u′模態(tài)從流向的帶狀到下游分岔或破碎形態(tài)反映的多尺度特征、v′模態(tài)沿流向的肋狀排列與尺寸增長、矢量 (u′,v′) 模態(tài)的環(huán)流模式表征的流動(dòng)卷吸、以及w′模態(tài)表征的流場結(jié)構(gòu)在展向上受拉伸的模式.

    (2) SMSM,LDSM 以及CSM 均較好地反映出湍流的多尺度特性,清晰地分辨了小尺度結(jié)構(gòu)的發(fā)展過程.而SM 與CKM 模型則未能有效地預(yù)測湍流中的小尺度渦結(jié)構(gòu)及部分流場特征模態(tài).

    (3) 脈動(dòng)速度場的POD 模態(tài)對(duì)SGS 模型的亞網(wǎng)格尺度耗散敏感,表現(xiàn)為CKM 預(yù)測的峰值耗散區(qū)對(duì)應(yīng)u′模態(tài)的低谷區(qū)因而環(huán)流模式不顯著,以及SM 未能反映轉(zhuǎn)捩區(qū)內(nèi)w′模態(tài)的脊?fàn)罾炷J?SMSM,LDSM 和CSM 克服了以上模型的不足之處,其中CSM 模型同時(shí)兼有計(jì)算效率較高的優(yōu)勢.

    猜你喜歡
    模態(tài)模型
    一半模型
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    3D打印中的模型分割與打包
    車輛CAE分析中自由模態(tài)和約束模態(tài)的應(yīng)用與對(duì)比
    國內(nèi)多模態(tài)教學(xué)研究回顧與展望
    FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
    高速顫振模型設(shè)計(jì)中顫振主要模態(tài)的判斷
    基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
    由單個(gè)模態(tài)構(gòu)造對(duì)稱簡支梁的抗彎剛度
    国产伦一二天堂av在线观看| 亚洲成人久久性| 久久精品久久久久久噜噜老黄 | 精品久久久久久久末码| 国产精华一区二区三区| 超碰av人人做人人爽久久| 舔av片在线| 午夜精品久久久久久毛片777| 欧美成人a在线观看| 欧美黑人巨大hd| 禁无遮挡网站| 亚洲最大成人手机在线| 真人一进一出gif抽搐免费| 国产人妻一区二区三区在| 好男人在线观看高清免费视频| 日韩有码中文字幕| 久久久色成人| 国产又黄又爽又无遮挡在线| 九九在线视频观看精品| 亚洲18禁久久av| 国产成人影院久久av| 精品久久久久久成人av| 免费av不卡在线播放| 亚洲精品日韩av片在线观看| 久久6这里有精品| 欧美黄色淫秽网站| 伊人久久精品亚洲午夜| 五月伊人婷婷丁香| 国产精品久久久久久人妻精品电影| 麻豆av噜噜一区二区三区| 国产在线男女| 麻豆久久精品国产亚洲av| 精品乱码久久久久久99久播| 国产精品嫩草影院av在线观看 | 脱女人内裤的视频| 欧美一区二区精品小视频在线| 一二三四社区在线视频社区8| 最近最新免费中文字幕在线| 日本在线视频免费播放| 白带黄色成豆腐渣| 国内揄拍国产精品人妻在线| а√天堂www在线а√下载| 97超级碰碰碰精品色视频在线观看| 免费观看精品视频网站| 日韩免费av在线播放| 九色成人免费人妻av| 高清毛片免费观看视频网站| 欧美高清性xxxxhd video| 日韩国内少妇激情av| 最好的美女福利视频网| 日韩欧美免费精品| 久久草成人影院| 91狼人影院| 国产男靠女视频免费网站| 欧美日韩中文字幕国产精品一区二区三区| 欧美色欧美亚洲另类二区| 亚洲av中文字字幕乱码综合| 人妻夜夜爽99麻豆av| 亚洲,欧美精品.| 夜夜爽天天搞| 亚洲不卡免费看| 成人国产综合亚洲| 在线免费观看的www视频| 久久亚洲真实| 亚洲天堂国产精品一区在线| 又黄又爽又免费观看的视频| 一进一出好大好爽视频| 亚洲成av人片免费观看| 国产在线男女| 精华霜和精华液先用哪个| 性色av乱码一区二区三区2| 日本黄色片子视频| h日本视频在线播放| 久久亚洲精品不卡| 亚洲最大成人手机在线| 在线播放无遮挡| 亚洲自拍偷在线| 少妇裸体淫交视频免费看高清| 美女高潮的动态| 国产精品亚洲美女久久久| 亚洲自拍偷在线| 69人妻影院| 人妻久久中文字幕网| 少妇高潮的动态图| 久久久久久久久久黄片| 非洲黑人性xxxx精品又粗又长| 18+在线观看网站| 99国产极品粉嫩在线观看| 黄色视频,在线免费观看| 欧美3d第一页| 日韩精品中文字幕看吧| 久久久久久久久久黄片| aaaaa片日本免费| 91狼人影院| 国产不卡一卡二| 免费搜索国产男女视频| 国内精品久久久久精免费| 悠悠久久av| 精品久久国产蜜桃| 女人被狂操c到高潮| 欧美激情久久久久久爽电影| 在线免费观看的www视频| 麻豆国产av国片精品| 国产精品野战在线观看| 色噜噜av男人的天堂激情| 色精品久久人妻99蜜桃| 免费人成在线观看视频色| 男人舔奶头视频| 俄罗斯特黄特色一大片| 直男gayav资源| 国语自产精品视频在线第100页| 男女做爰动态图高潮gif福利片| 十八禁人妻一区二区| 欧美一级a爱片免费观看看| 成人av一区二区三区在线看| 亚洲七黄色美女视频| 内地一区二区视频在线| 欧美一区二区精品小视频在线| 国产一区二区三区在线臀色熟女| 色av中文字幕| 亚洲成人久久性| 国产精品自产拍在线观看55亚洲| 一卡2卡三卡四卡精品乱码亚洲| .国产精品久久| 精品乱码久久久久久99久播| xxxwww97欧美| 激情在线观看视频在线高清| 十八禁网站免费在线| 狂野欧美白嫩少妇大欣赏| 91九色精品人成在线观看| 五月伊人婷婷丁香| 欧美乱妇无乱码| 中文字幕av在线有码专区| 久久国产乱子免费精品| 中国美女看黄片| 99热精品在线国产| 久久香蕉精品热| 少妇被粗大猛烈的视频| 亚洲18禁久久av| 高清毛片免费观看视频网站| 国产色爽女视频免费观看| 午夜两性在线视频| 俄罗斯特黄特色一大片| 久久人妻av系列| 在线观看av片永久免费下载| 亚洲黑人精品在线| 丁香六月欧美| 久久亚洲精品不卡| 国产人妻一区二区三区在| 成人永久免费在线观看视频| 免费电影在线观看免费观看| 亚洲人与动物交配视频| 熟妇人妻久久中文字幕3abv| 亚洲熟妇中文字幕五十中出| 如何舔出高潮| 国产爱豆传媒在线观看| 免费看美女性在线毛片视频| 高清在线国产一区| 国产精品久久久久久亚洲av鲁大| 69av精品久久久久久| 亚洲在线自拍视频| 婷婷精品国产亚洲av在线| 日日干狠狠操夜夜爽| 国产在视频线在精品| 12—13女人毛片做爰片一| 99国产综合亚洲精品| 亚洲avbb在线观看| 欧美+亚洲+日韩+国产| 亚洲七黄色美女视频| 真实男女啪啪啪动态图| 日韩人妻高清精品专区| 亚洲精品影视一区二区三区av| 淫妇啪啪啪对白视频| av天堂在线播放| 一个人观看的视频www高清免费观看| 日本一本二区三区精品| 午夜福利欧美成人| 国产欧美日韩一区二区精品| 国产精品乱码一区二三区的特点| АⅤ资源中文在线天堂| 午夜福利成人在线免费观看| 我的老师免费观看完整版| 俄罗斯特黄特色一大片| 啦啦啦观看免费观看视频高清| 国产三级在线视频| 毛片女人毛片| 黄色一级大片看看| 免费高清视频大片| 亚洲精品乱码久久久v下载方式| 12—13女人毛片做爰片一| 男女之事视频高清在线观看| 可以在线观看毛片的网站| eeuss影院久久| 国产精品嫩草影院av在线观看 | 婷婷精品国产亚洲av| 亚洲av成人av| 校园春色视频在线观看| 国产精品影院久久| 毛片女人毛片| .国产精品久久| 欧美性猛交黑人性爽| 亚洲国产精品久久男人天堂| 午夜免费成人在线视频| 99在线视频只有这里精品首页| 日韩国内少妇激情av| 国产高清视频在线观看网站| 免费看日本二区| 精品一区二区三区视频在线观看免费| 亚洲欧美日韩卡通动漫| 在线观看舔阴道视频| 黄色女人牲交| 亚洲成a人片在线一区二区| 变态另类成人亚洲欧美熟女| or卡值多少钱| 成人一区二区视频在线观看| 国产欧美日韩一区二区三| 97超级碰碰碰精品色视频在线观看| 国产老妇女一区| 俄罗斯特黄特色一大片| 国产精品野战在线观看| 99热只有精品国产| 亚洲不卡免费看| 亚洲av成人不卡在线观看播放网| 久久天躁狠狠躁夜夜2o2o| 看片在线看免费视频| 国产伦精品一区二区三区视频9| 国产成人aa在线观看| 精品福利观看| 一级毛片久久久久久久久女| 夜夜躁狠狠躁天天躁| 亚洲国产精品sss在线观看| 俺也久久电影网| 欧美日本视频| 少妇的逼好多水| 中文字幕免费在线视频6| 日韩免费av在线播放| 最近最新免费中文字幕在线| 欧美精品啪啪一区二区三区| 亚洲av日韩精品久久久久久密| 日韩欧美在线二视频| 亚洲 国产 在线| 性插视频无遮挡在线免费观看| 小蜜桃在线观看免费完整版高清| www.色视频.com| 国产亚洲欧美98| av在线蜜桃| 亚洲欧美日韩卡通动漫| 国产精品亚洲美女久久久| 久久精品影院6| 国产久久久一区二区三区| 97碰自拍视频| 国产亚洲av嫩草精品影院| 两性午夜刺激爽爽歪歪视频在线观看| 成年女人永久免费观看视频| 午夜精品一区二区三区免费看| 男女下面进入的视频免费午夜| 午夜激情福利司机影院| 悠悠久久av| 午夜福利在线观看吧| 中文资源天堂在线| 90打野战视频偷拍视频| 免费人成在线观看视频色| 久久精品国产亚洲av香蕉五月| 国内精品久久久久久久电影| 在线观看av片永久免费下载| 亚洲五月婷婷丁香| 黄色日韩在线| 一区二区三区四区激情视频 | 亚洲黑人精品在线| 国产伦一二天堂av在线观看| 国产一区二区在线av高清观看| 精品久久久久久久久久久久久| 日韩欧美精品v在线| 99精品久久久久人妻精品| 黄色日韩在线| 一区二区三区四区激情视频 | 18美女黄网站色大片免费观看| 9191精品国产免费久久| 亚洲av.av天堂| 蜜桃亚洲精品一区二区三区| 亚洲成人精品中文字幕电影| 夜夜夜夜夜久久久久| 国产一级毛片七仙女欲春2| 欧美zozozo另类| 91九色精品人成在线观看| 变态另类丝袜制服| 啦啦啦观看免费观看视频高清| 一区福利在线观看| 欧美3d第一页| 久久人人精品亚洲av| 成人av在线播放网站| 久久香蕉精品热| 亚洲黑人精品在线| 亚洲av五月六月丁香网| 国产黄色小视频在线观看| 国产精品野战在线观看| 国产高清激情床上av| 禁无遮挡网站| 天堂av国产一区二区熟女人妻| 精品人妻偷拍中文字幕| 制服丝袜大香蕉在线| 日本成人三级电影网站| 久久人人爽人人爽人人片va | 亚洲精品一区av在线观看| 精品欧美国产一区二区三| 又紧又爽又黄一区二区| 成人无遮挡网站| 伦理电影大哥的女人| 午夜亚洲福利在线播放| 老司机深夜福利视频在线观看| 国产一区二区三区在线臀色熟女| 日韩欧美在线乱码| 嫩草影院精品99| 成人欧美大片| 久久人人爽人人爽人人片va | 男女做爰动态图高潮gif福利片| 一级av片app| 最好的美女福利视频网| 久久久久性生活片| av中文乱码字幕在线| 18禁黄网站禁片午夜丰满| 亚洲精品在线美女| 亚洲精品影视一区二区三区av| 国内精品一区二区在线观看| 欧美性猛交╳xxx乱大交人| 午夜视频国产福利| 久久久久久久午夜电影| 小蜜桃在线观看免费完整版高清| 综合色av麻豆| 在线观看舔阴道视频| 亚洲av日韩精品久久久久久密| 99精品在免费线老司机午夜| 国产精品一区二区三区四区久久| 日韩大尺度精品在线看网址| 国产成人福利小说| 可以在线观看毛片的网站| 欧美精品国产亚洲| 看免费av毛片| 99在线人妻在线中文字幕| 欧美一级a爱片免费观看看| 黄片小视频在线播放| 中文在线观看免费www的网站| 国产精品av视频在线免费观看| 欧美3d第一页| 日韩成人在线观看一区二区三区| 女人被狂操c到高潮| 自拍偷自拍亚洲精品老妇| 1024手机看黄色片| 日韩免费av在线播放| 成人性生交大片免费视频hd| 精品人妻偷拍中文字幕| 很黄的视频免费| 国产精品久久久久久亚洲av鲁大| 国产视频内射| 亚洲自拍偷在线| 亚洲一区高清亚洲精品| 亚洲欧美日韩无卡精品| 精品日产1卡2卡| 亚洲精品日韩av片在线观看| 亚洲avbb在线观看| 丰满的人妻完整版| 观看免费一级毛片| АⅤ资源中文在线天堂| 窝窝影院91人妻| 国产一区二区在线av高清观看| 免费在线观看日本一区| 国产一区二区激情短视频| 精品国内亚洲2022精品成人| 精品久久国产蜜桃| netflix在线观看网站| 亚洲无线在线观看| netflix在线观看网站| 黄色视频,在线免费观看| 国产熟女xx| 在线观看美女被高潮喷水网站 | a级毛片免费高清观看在线播放| 脱女人内裤的视频| 露出奶头的视频| 国产av麻豆久久久久久久| 国产黄a三级三级三级人| 精品一区二区三区av网在线观看| 久久人人爽人人爽人人片va | 午夜福利欧美成人| 亚洲精品乱码久久久v下载方式| 男人舔女人下体高潮全视频| 国产成人啪精品午夜网站| 男女下面进入的视频免费午夜| 成人亚洲精品av一区二区| 神马国产精品三级电影在线观看| 久久久久九九精品影院| 国产成人a区在线观看| 精品久久久久久成人av| 1000部很黄的大片| 国内毛片毛片毛片毛片毛片| 亚洲人成网站高清观看| 国产极品精品免费视频能看的| 成人午夜高清在线视频| 天堂网av新在线| 大型黄色视频在线免费观看| 午夜亚洲福利在线播放| 国产成人a区在线观看| 成人三级黄色视频| 精品熟女少妇八av免费久了| av福利片在线观看| 极品教师在线免费播放| 男女那种视频在线观看| 一级a爱片免费观看的视频| 国产真实乱freesex| 国产黄片美女视频| 噜噜噜噜噜久久久久久91| 国产综合懂色| 国产一区二区在线观看日韩| 午夜激情欧美在线| 动漫黄色视频在线观看| 日本免费a在线| 免费在线观看日本一区| 男女做爰动态图高潮gif福利片| 久久99热6这里只有精品| 三级国产精品欧美在线观看| 男女视频在线观看网站免费| 亚洲国产日韩欧美精品在线观看| 日韩欧美三级三区| 如何舔出高潮| 久久久久亚洲av毛片大全| 一本久久中文字幕| 成人国产综合亚洲| 一个人免费在线观看电影| 日日摸夜夜添夜夜添小说| 国产蜜桃级精品一区二区三区| 日本黄色片子视频| 麻豆久久精品国产亚洲av| 欧美xxxx黑人xx丫x性爽| 成人永久免费在线观看视频| 久久久久久久午夜电影| av在线老鸭窝| 男女做爰动态图高潮gif福利片| 欧美黄色淫秽网站| 日韩人妻高清精品专区| 给我免费播放毛片高清在线观看| 99国产精品一区二区三区| 在线免费观看的www视频| 成年女人毛片免费观看观看9| 欧美不卡视频在线免费观看| 欧美黑人欧美精品刺激| 免费看日本二区| 国产三级中文精品| 成人无遮挡网站| 757午夜福利合集在线观看| 九九热线精品视视频播放| 国产大屁股一区二区在线视频| 搡老岳熟女国产| 很黄的视频免费| 三级国产精品欧美在线观看| 十八禁人妻一区二区| 国产伦在线观看视频一区| 久久精品综合一区二区三区| 国产探花在线观看一区二区| 我的女老师完整版在线观看| 精品人妻一区二区三区麻豆 | 18+在线观看网站| 免费观看人在逋| av福利片在线观看| 国产真实乱freesex| 亚洲乱码一区二区免费版| 日韩欧美在线二视频| 日本精品一区二区三区蜜桃| 午夜日韩欧美国产| 国产精品久久久久久人妻精品电影| 宅男免费午夜| 日韩国内少妇激情av| 免费搜索国产男女视频| 色哟哟哟哟哟哟| 无人区码免费观看不卡| 3wmmmm亚洲av在线观看| 欧洲精品卡2卡3卡4卡5卡区| 最新中文字幕久久久久| 精品午夜福利视频在线观看一区| 怎么达到女性高潮| 又粗又爽又猛毛片免费看| 亚洲国产精品成人综合色| xxxwww97欧美| 日本五十路高清| 久久久国产成人免费| 在线天堂最新版资源| 欧美成人一区二区免费高清观看| 色5月婷婷丁香| 国产高清视频在线观看网站| 国产aⅴ精品一区二区三区波| 精华霜和精华液先用哪个| 桃色一区二区三区在线观看| 我的老师免费观看完整版| 成年人黄色毛片网站| 黄色一级大片看看| 成人特级av手机在线观看| 精品国内亚洲2022精品成人| 免费电影在线观看免费观看| 夜夜夜夜夜久久久久| ponron亚洲| 成人性生交大片免费视频hd| 动漫黄色视频在线观看| 69av精品久久久久久| 黄色配什么色好看| 免费大片18禁| 成人精品一区二区免费| 99在线人妻在线中文字幕| 99久久成人亚洲精品观看| 亚洲精品色激情综合| 国产熟女xx| 18禁裸乳无遮挡免费网站照片| 国产精品久久久久久亚洲av鲁大| 亚洲第一电影网av| 久久久久久久午夜电影| 国产精品1区2区在线观看.| 日韩 亚洲 欧美在线| 精华霜和精华液先用哪个| 亚洲在线自拍视频| 精品免费久久久久久久清纯| 精品一区二区三区人妻视频| 婷婷丁香在线五月| 九色成人免费人妻av| 久久久久国产精品人妻aⅴ院| 午夜福利视频1000在线观看| 日韩欧美精品v在线| 国产又黄又爽又无遮挡在线| 免费大片18禁| 亚洲精品一卡2卡三卡4卡5卡| 一进一出抽搐动态| 一区福利在线观看| 少妇的逼水好多| 日本 欧美在线| 又粗又爽又猛毛片免费看| 日本在线视频免费播放| 国产精品爽爽va在线观看网站| 亚洲精品一卡2卡三卡4卡5卡| 热99re8久久精品国产| 国产免费av片在线观看野外av| 免费一级毛片在线播放高清视频| 久久久久久久精品吃奶| 日韩人妻高清精品专区| 亚洲不卡免费看| 长腿黑丝高跟| 黄色一级大片看看| 一本一本综合久久| 少妇丰满av| 日日夜夜操网爽| av视频在线观看入口| ponron亚洲| 日本在线视频免费播放| 我要看日韩黄色一级片| 夜夜看夜夜爽夜夜摸| 免费看日本二区| 国产国拍精品亚洲av在线观看| 伊人久久精品亚洲午夜| 精品熟女少妇八av免费久了| 黄色女人牲交| 亚洲一区二区三区不卡视频| 男人的好看免费观看在线视频| 国产真实伦视频高清在线观看 | 9191精品国产免费久久| 男人的好看免费观看在线视频| 无遮挡黄片免费观看| 国产色婷婷99| 日韩av在线大香蕉| 成人国产一区最新在线观看| 国产在视频线在精品| 少妇裸体淫交视频免费看高清| 精品一区二区三区视频在线观看免费| 91麻豆av在线| 国产一区二区亚洲精品在线观看| 伦理电影大哥的女人| 非洲黑人性xxxx精品又粗又长| 宅男免费午夜| 国产精品乱码一区二三区的特点| 亚洲成人久久爱视频| 亚洲激情在线av| 欧美日韩国产亚洲二区| 久久性视频一级片| 在线观看一区二区三区| 欧美高清成人免费视频www| 亚洲欧美清纯卡通| 少妇人妻一区二区三区视频| 成人欧美大片| 日韩高清综合在线| 国产午夜精品久久久久久一区二区三区 | 欧美不卡视频在线免费观看| 狂野欧美白嫩少妇大欣赏| 宅男免费午夜| 久久久久国产精品人妻aⅴ院| 51国产日韩欧美| 精品一区二区免费观看| 国产精品不卡视频一区二区 | 日韩精品青青久久久久久| 久久人妻av系列| 欧美zozozo另类| 国产精品久久久久久亚洲av鲁大| 亚洲国产精品久久男人天堂| 国产精品自产拍在线观看55亚洲| 熟女电影av网| 免费看日本二区| 亚洲无线在线观看| 在线天堂最新版资源| 午夜福利成人在线免费观看| 嫩草影院精品99| 1000部很黄的大片| 精品人妻一区二区三区麻豆 | 国产一区二区在线观看日韩| 午夜日韩欧美国产| 日韩亚洲欧美综合| 久久久成人免费电影| 色5月婷婷丁香| 精品人妻熟女av久视频| 黄片小视频在线播放| 成年版毛片免费区| 麻豆av噜噜一区二区三区| 国产三级黄色录像|