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

    SAS與PANS模型在圓柱繞流中的應(yīng)用比較

    2019-08-29 09:14:40躍管小榮徐
    關(guān)鍵詞:實(shí)驗(yàn)模型

    劉 躍管小榮徐 誠

    (南京理工大學(xué) 機(jī)械工程學(xué)院,南京 210094)

    0 引 言

    廣泛存在于航空航天、風(fēng)工程等空氣動(dòng)力學(xué)領(lǐng)域的鈍體繞流問題通常為無序、非線性、多尺度的湍流分離運(yùn)動(dòng),在使用數(shù)值模擬方法(CFD)對(duì)其進(jìn)行分析計(jì)算時(shí),湍流的求解成為影響CFD計(jì)算精度的瓶頸問題[1]。在諸多的湍流模型中,兼顧計(jì)算精度與資源消耗的RANS/LES混合方法受到人們的青睞[2],混合模型在近壁區(qū)及遠(yuǎn)場(chǎng)使用RANS模型,而在分離區(qū)及尾跡區(qū)切換為更合理的類LES濾波模型,極大地提高了計(jì)算精度,且資源消耗遠(yuǎn)小于LES僅略大于RANS。近年來,不同學(xué)者發(fā)展了種類繁多的混合方法,其中包括DES(Detached Eddy Simulation)方法[3]、SAS(Scale Adaptive Simulation)方 法[4-7]、PANS(Partially Averaged Navier-Stokes)方法[8]、PITM(Partially Integrated Transport Model)方法、LNS(Limited Numerical Scales)方法等。其中DES方法是最早出現(xiàn)的混合方法,經(jīng)過一系列改進(jìn)(DDES、IDDES)已發(fā)展較為成熟[9-13];SAS模型為能夠解析寬頻非定常湍流脈動(dòng)而對(duì)計(jì)算網(wǎng)格依賴較小的新一代非定常湍流預(yù)測(cè)方法(URANS),通過加入?yún)⒄债?dāng)?shù)亓鲃?dòng)的馮卡門(von Karman)長(zhǎng)度尺度L vk,合理地釋放更多流場(chǎng)信息[14-15];PANS模型同樣是最近發(fā)展的一種RANS-DNS橋接模型,模型構(gòu)造相對(duì)簡(jiǎn)單,通過加入?;膭?dòng)能比例參數(shù)f k調(diào)整湍流的?;縼碚{(diào)整湍流求解尺度[16];DES模型的提出者Spalart在2000年時(shí)曾預(yù)測(cè),2045年左右計(jì)算機(jī)的計(jì)算能力才能滿足利用LES解決90%的湍流問題,因此,在現(xiàn)階段,驗(yàn)證與發(fā)展RANS/LES方法仍是非常必要的[17-18]。

    鈍體繞流中的圓柱繞流計(jì)算模型簡(jiǎn)單,同時(shí)又包含如湍流壁面、剪切層失穩(wěn)、流動(dòng)分離、周期性渦脫落等豐富的流動(dòng)現(xiàn)象,是考察湍流模型計(jì)算性能的標(biāo)準(zhǔn)算例。Re=3900的圓柱繞流是典型的亞臨界雷諾數(shù)流動(dòng),流動(dòng)結(jié)構(gòu)復(fù)雜且耗費(fèi)資源較少,人們對(duì)該雷諾數(shù)的圓柱繞流進(jìn)行了大量的實(shí)驗(yàn)[19-22]及數(shù)值研究[23-29],有可靠豐富的數(shù)據(jù)可供參考。本文選用SAS模型及兩種變f k構(gòu)造函數(shù)的PANS模型計(jì)算了Re=3900的圓柱繞流流場(chǎng),并重點(diǎn)從流場(chǎng)特征捕捉、氣動(dòng)力特性計(jì)算、渦黏性控制等方面對(duì)SAS模型及PANS模型進(jìn)行了比較,通過對(duì)比標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格后的計(jì)算結(jié)果,比較了兩類混合方法對(duì)網(wǎng)格的適應(yīng)性,并探討了兩類混合模型的尺度函數(shù)對(duì)流場(chǎng)的調(diào)控機(jī)制。

    1 計(jì)算模型及數(shù)值方法

    1.1 控制方程

    控制方程為三維N-S方程,在一般計(jì)算坐標(biāo)系(τ,ξ,η,ζ)下的 守恒形 式為:

    式中:為守恒變量為無粘通量為黏性通量,各物理量的具體表述參見文獻(xiàn)[1]。

    1.2 湍流模型

    1.2.1 SST-k-ω模型

    SST模型由標(biāo)準(zhǔn)k-ω模型與k-ε模型混合修改而成,該模型兼顧了k-ω模型的近壁性能及k-ε模型的遠(yuǎn)場(chǎng)準(zhǔn)確性,SST模型中湍動(dòng)能k和比耗散率ω的輸運(yùn)方程可表示如下:

    式中模型參數(shù)σk、α、β、σω、σω2均由對(duì)應(yīng)的k-ω模型參數(shù)φ1和k-ε模型參數(shù)φ2通過一定函數(shù)計(jì)算得到,其對(duì)應(yīng)的函數(shù)關(guān)系為:

    SST模型中湍動(dòng)黏度μt可表示為:

    式中S為應(yīng)變率張量S ij的模,SST模型中F1、F2函數(shù)構(gòu)造方式及各模型參數(shù)值參考文獻(xiàn)[30]。

    1.2.2 SST-SAS模型

    尺度自適應(yīng)SAS(Scale-Adaptive Simulation)模型是一種較新的創(chuàng)新型湍流模型,它通過在湍流模型中加入刻畫當(dāng)?shù)亓鲃?dòng)拓?fù)涞膙on Karman長(zhǎng)度尺度L vk來作為湍流模型的第二長(zhǎng)度尺度[6]。基于SST模型的SAS模型構(gòu)造方法為在式(3)SST模型的ω方程中加入發(fā)揮自適應(yīng)作用的源項(xiàng)QSAS,其定義式為:

    由式(6)可以看出:在SAS模型中出現(xiàn)了由一階及二階速度導(dǎo)數(shù)決定的長(zhǎng)度尺度項(xiàng)L vk,該長(zhǎng)度尺度基于當(dāng)?shù)亓鲃?dòng)而與網(wǎng)格尺度無關(guān),在湍流邊界層內(nèi),L vk能?;瘧T性子區(qū)所有湍流脈動(dòng),同時(shí),在非穩(wěn)態(tài)區(qū)域L vk能根據(jù)當(dāng)?shù)鼐W(wǎng)格分辨的湍流渦動(dòng)態(tài)調(diào)整RANS長(zhǎng)度尺度。在SAS模型中,湍動(dòng)黏度μt仍按照式(5)計(jì)算。

    1.2.3 SST-PANS模型

    本文中所用PANS模型同樣是基于SST模型修改而成[8],其原理是在SST模型中引入兩個(gè)模型控制參數(shù)f k和fω:

    式中:f k為模化湍動(dòng)能與總的湍動(dòng)能之比,fω是?;群纳⒙逝c總的比耗散率的比值。如果假定f k和fω是不隨時(shí)間和空間變化的常數(shù),則k和ω的輸運(yùn)方程(2)、(3)變?yōu)?

    將式(8)、(9)與式(2)、(3)對(duì)比可以發(fā)現(xiàn)在PANS模型中主要是模型參數(shù)發(fā)生了變化,PANS模型中模型參數(shù)σku、σωu、σω2u、β表示為:

    通過以上變化,將SST模型轉(zhuǎn)變?yōu)镾ST-PANS模型,由k、ε、ω的關(guān)系可得式(10)中fω與f k和fε的關(guān)系,在SST-PANS模型中,fε一般取為1,這樣PANS中的參數(shù)可以表示為f k的函數(shù)。通過調(diào)節(jié)f k的值(1-0),改變模化的湍動(dòng)能與比耗散率的比值,從而實(shí)現(xiàn)RANS-DNS的轉(zhuǎn)化。這時(shí)湍動(dòng)黏度μtu變?yōu)?

    在早期PANS模型的應(yīng)用中,f k通常被設(shè)置為空間常數(shù),通過計(jì)算發(fā)現(xiàn)對(duì)于整個(gè)流場(chǎng)很難確定合適統(tǒng)一的f k值。對(duì)于復(fù)雜的分離流動(dòng)問題,計(jì)算時(shí)分離區(qū)與尾跡區(qū)需要較小的f k值(f k<1)釋放湍流信息,而近壁區(qū)與遠(yuǎn)場(chǎng)區(qū)采用RANS模型(f k=1)更合理,因此近年來人們逐漸嘗試用適當(dāng)?shù)暮瘮?shù)來計(jì)算f k分布,本文選取其中兩種相對(duì)簡(jiǎn)單、計(jì)算相對(duì)穩(wěn)定的變f k函數(shù)來進(jìn)行對(duì)比。

    (a)類DES型f k函數(shù)[31-32]。該函數(shù)參考DES模型中模型尺度的構(gòu)造方法,將f k表示為網(wǎng)格尺度Δ與湍流尺度L u的函數(shù):

    式中模型常數(shù)CPANS取0.3。用該函數(shù)進(jìn)行計(jì)算時(shí),在分離區(qū)與尾跡區(qū),網(wǎng)格尺度要小于湍流尺度,因此f k<1;在壁面及遠(yuǎn)場(chǎng)區(qū)域,網(wǎng)格尺度大于湍流尺度,因此f k=1。文中將使用類DES型f k分布函數(shù)的PANS模型稱為PANS(des)模型。

    (b)Tanh型f k函數(shù)[33]。在該函數(shù)中,f k可表示為:

    對(duì)比式(12)與式(13)可知,Tanh型f k函數(shù)由湍流尺度、網(wǎng)格尺度及tanh函數(shù)修改而成,tanh函數(shù)的取值范圍為(-1,1),這樣根據(jù)流動(dòng)與網(wǎng)格關(guān)系得到不同位置的f k取值范圍為(0,1)。文中將使用Tanh型f k函數(shù)的PANS模型稱為PANS(tanh)模型。

    2 計(jì)算模型及數(shù)值方法

    2.1 計(jì)算域、網(wǎng)格及邊界設(shè)置

    文中圓柱繞流計(jì)算雷諾數(shù)為Re=3900(基于圓柱直徑D與遠(yuǎn)場(chǎng)來流速度U∞),圖1為計(jì)算域及邊界條件示意圖。為減小上游來流對(duì)計(jì)算的影響并使尾跡區(qū)流動(dòng)充分發(fā)展,上游來流區(qū)域(x方向)及橫向區(qū)域(y方向)延展至20D,下游計(jì)算區(qū)域設(shè)置為30D。根據(jù)近尾跡區(qū)流向結(jié)構(gòu)的展向波長(zhǎng)設(shè)置計(jì)算域展向長(zhǎng)度(z方向)為πD[2]。圓柱體設(shè)置為無滑移壁面,上游入口流動(dòng)各物理量由無擾動(dòng)自由來流條件給定,下游出口設(shè)置為無反射邊界,計(jì)算域的上下界面設(shè)置為自由滑移壁面,計(jì)算域的前后邊界設(shè)置為周期邊界。文中坐標(biāo)系設(shè)置如圖1右下角所示,軸向從上游到下游方向?yàn)閤正方向,坐標(biāo)系符合右手螺旋法則,橫向?yàn)閥方向,展向?yàn)閦方向,方位角從上游圓柱前駐點(diǎn)為0°開始,順時(shí)針旋轉(zhuǎn)為壁面方位角θ的正方向,下游背風(fēng)點(diǎn)為180°。

    圖1 計(jì)算域及邊界條件Fig.1 Computational domain and boundary conditions

    圖2為圓柱繞流的局部網(wǎng)格示意,計(jì)算網(wǎng)格為分塊結(jié)構(gòu)化網(wǎng)格。為了保證計(jì)算精度,在壁面附近及尾跡分離區(qū)域?qū)W(wǎng)格進(jìn)行適當(dāng)加密,壁面第一層網(wǎng)格高度為0.002D,保證壁面y+<1。為了探討不同模型的尺度函數(shù)對(duì)網(wǎng)格的敏感性,文中設(shè)置了兩套計(jì)算網(wǎng)格,其中標(biāo)準(zhǔn)網(wǎng)格保證近壁增長(zhǎng)率為1.06,分離區(qū)網(wǎng)格固定為0.1D;加密網(wǎng)格為保證近壁增長(zhǎng)率為1.04,分離區(qū)網(wǎng)格固定為0.05D。參考文獻(xiàn)[24,25]中的網(wǎng)格設(shè)置見表1。

    圖2 圓柱繞流局部網(wǎng)格Fig.2 Computational grid for the cylinder

    表1 計(jì)算網(wǎng)格Table 1 Computational grids

    2.2 方程求解設(shè)置

    文中對(duì)所求方程組空間離散采用有限體積法(FVM),對(duì)流項(xiàng)中速度項(xiàng)的離散采用三階限制(TVD)中心差分格式,湍流項(xiàng)離散采用二階限制中心差分格式,梯度項(xiàng)及黏性項(xiàng)的離散采用四階中心差分格式,壓力與速度耦合選用PISO算法。為了保證計(jì)算的收斂性及穩(wěn)定性,時(shí)間項(xiàng)離散采用二階隱式后退歐拉(Backward)格式,并采用自適應(yīng)時(shí)間步長(zhǎng)方法[24],時(shí)間步長(zhǎng)設(shè)置為庫朗數(shù)(Courant number)小于0.5,并且保證在一個(gè)渦脫落周期內(nèi)包含近500個(gè)時(shí)間步。

    3 計(jì)算結(jié)果分析

    3.1 數(shù)據(jù)采集及流場(chǎng)宏觀結(jié)構(gòu)

    文中分別計(jì)算了SST、SAS、PANS(des)、PANS(tanh)模型的圓柱繞流流場(chǎng)。周期性渦脫落是該計(jì)算雷諾數(shù)下的顯著特征。為了獲得合理的統(tǒng)計(jì)數(shù)據(jù),根據(jù)Franke等[34]的研究結(jié)果,圓柱繞流的統(tǒng)計(jì)周期應(yīng)大于40個(gè)脫落周期。因此參照文獻(xiàn)[34,25],本文中總的計(jì)算時(shí)間約為60個(gè)渦脫落周期,在對(duì)非定常結(jié)果的統(tǒng)計(jì)時(shí)去掉前10個(gè)脫落周期,將后50個(gè)周期內(nèi)的統(tǒng)計(jì)數(shù)據(jù)作為有效數(shù)據(jù)進(jìn)行模型的對(duì)比。計(jì)算結(jié)果進(jìn)行時(shí)間與空間平均。圖3中給出了以標(biāo)準(zhǔn)網(wǎng)格計(jì)算結(jié)果為例的圓柱繞流升力系數(shù)(C l)、阻力系數(shù)(C d)的時(shí)程(U∞/D無量綱化)變化曲線??梢钥闯鯮e=3900時(shí)的升力幅值已不再穩(wěn)定,但仍呈現(xiàn)出周期性變化趨勢(shì),同時(shí)阻力系數(shù)表現(xiàn)出強(qiáng)烈的脈動(dòng)現(xiàn)象,并在一平均值附近上下擺動(dòng)。

    圖3 升力、阻力系數(shù)隨時(shí)間變化及升力頻譜Fig.3 Time history of lift C l,drag C d coef ficients

    圖4中給出了由升力系數(shù)時(shí)程曲線經(jīng)過傅里葉變換(FFT)獲得的渦脫落頻率功率譜,由此可以計(jì)算渦脫落頻率f及斯特勞哈數(shù)S t(S t=f D/U∞)。各模型阻力系數(shù)(C d)及斯特勞哈數(shù)S t見表2。

    圖4 升力頻譜Fig.4 Spectra of lift

    圖5為四種模型基于標(biāo)準(zhǔn)網(wǎng)格計(jì)算得到的圓柱繞流流場(chǎng)二維瞬時(shí)(t=300)渦量分布及用軸向速度染色的三維Q分布,對(duì)比可知:由于RANS模型湍流黏性較大,抹去了很多重要的流場(chǎng)細(xì)節(jié),計(jì)算結(jié)果只捕捉到大尺度的脫落渦結(jié)構(gòu),Q圖也只先顯示出近似二維的流動(dòng)結(jié)構(gòu),這說明RANS模型基本上沒有求解小尺度流動(dòng)的能力;SAS模型及PANS模型均捕捉到豐富的流動(dòng)結(jié)構(gòu),二維渦量及三維Q值均顯示出復(fù)雜的周期性渦脫落現(xiàn)象。分析三者的Q圖可清晰地發(fā)現(xiàn)附著于壁面的二維層流剪切層及分離區(qū)流向與展向相嵌套的三維無序渦體結(jié)構(gòu)。SAS模型與PANS模型均允許形成湍流譜,表現(xiàn)出預(yù)期的類似LES模型的求解小尺度運(yùn)動(dòng)的能力。同時(shí),仔細(xì)對(duì)比SAS模型與PANS模型的結(jié)果可以看出,PANS模型能捕捉到比SAS模型更豐富的湍流結(jié)構(gòu)。

    圖5 不同模型計(jì)算得到的流場(chǎng)瞬時(shí)渦量與Q值圖(t=300)Fig.5 Instantaneous vorticity and Q of different models(t=300)

    由實(shí)驗(yàn)及以往計(jì)算結(jié)果[23]可知,圓柱繞流壁面剪切層的發(fā)展對(duì)外界影響因素非常敏感,微小外部干擾很容易造成邊界層的提前分離,使得剪切層變短,進(jìn)而影響回流區(qū)及整個(gè)流場(chǎng)的測(cè)量或計(jì)算結(jié)果。圖6為用渦量等值線表示的標(biāo)準(zhǔn)網(wǎng)格條件下四種模型的瞬時(shí)(t=300)剪切層,通過對(duì)比可知剪切層長(zhǎng)度關(guān)系為:SAS>PANS(des)>PANS(tanh)>SST。

    圖6 不同模型瞬時(shí)分離剪切層示意圖(t=300)Fig.6 Instantaneous separating shear layers of different models(t=300)

    3.2 回流區(qū)速度分布

    圖7為以SAS模型為例的y=0平面軸向速度時(shí)均圖,可以清晰地看出在圓柱近尾跡區(qū)域形成一個(gè)流速與來流相反的區(qū)域,即圓柱繞流的回流區(qū)。圖8為圓柱下游尾跡區(qū)中心線上時(shí)均軸向速度分布。由圖可知,圓柱背風(fēng)點(diǎn)處速度為0,速度為逐漸增大的負(fù)值,到達(dá)回流區(qū)中心時(shí)速度負(fù)值最大。背風(fēng)點(diǎn)與空間流動(dòng)的再附點(diǎn)即流向速度恢復(fù)到0的點(diǎn)之間的距離即為回流區(qū)長(zhǎng)度(L r/D)。四種模型均模擬出了圓柱后的速度分區(qū),但彼此間差別較大,特別是對(duì)回流區(qū)的模擬。圖8(b)圖是對(duì)回流區(qū)中心線時(shí)均軸向速度的放大圖。首先需要說明的是用空心圓圈表示的實(shí)驗(yàn)數(shù)據(jù)是Lourenco&Shih(1993)用PIV測(cè)量獲得[20],用實(shí)心圓表示的實(shí)驗(yàn)數(shù)據(jù)是Parnaudeau等(2008)同樣用PIV測(cè)量獲得[21],可以看到兩者實(shí)驗(yàn)測(cè)得的回流區(qū)長(zhǎng)度差別較大。最近幾年的文獻(xiàn)也出現(xiàn)了與Parnaudeau實(shí)驗(yàn)接近或與Lourenco實(shí)驗(yàn)接近的兩種數(shù)值計(jì)算結(jié)果,有意思的是Kravchenko等(2000)在用B-Spline伽遼金方法對(duì)Re=3900進(jìn)行的LES數(shù)值計(jì)算中[23],最先發(fā)現(xiàn)了數(shù)值結(jié)果與Lourenco實(shí)驗(yàn)數(shù)據(jù)在回流區(qū)與近尾跡區(qū)速度型的差別。Kravchenko計(jì)算的回流區(qū)長(zhǎng)度大于Lourenco實(shí)驗(yàn)值,其當(dāng)時(shí)推測(cè)可能有外加擾動(dòng)導(dǎo)致Lourenco實(shí)驗(yàn)中圓柱壁面剪切層過早分離,造成剪切層變短,進(jìn)而影響回流區(qū)長(zhǎng)度。Parnaudeau實(shí)驗(yàn)數(shù)據(jù)在回流區(qū)及速度型方面均驗(yàn)證了Kravchenko的數(shù)值計(jì)算結(jié)果。由于實(shí)驗(yàn)時(shí)間較晚,實(shí)驗(yàn)方法更先進(jìn)、嚴(yán)謹(jǐn),本文中我們認(rèn)為Parnaudeau的實(shí)驗(yàn)值是更可信的。分析圖8(b)可知,SST模型計(jì)算的圓柱回流區(qū)長(zhǎng)度誤差較大,SAS與PANS計(jì)算的回流區(qū)長(zhǎng)度關(guān)系為SAS>PANS(des)>PANS(tanh)。三種模型計(jì)算的回流區(qū)區(qū)域位于兩個(gè)實(shí)驗(yàn)值之間,與Parnaudeau實(shí)驗(yàn)值仍存在一定誤差,其中SAS模型與Parnaudeau實(shí)驗(yàn)最接近。以上結(jié)論與圖6中得出的剪切層長(zhǎng)度結(jié)論也是一致的。

    圖7 y=0平面軸向速度時(shí)均分布Fig.7 Mean streamwise velocity along the y=0 surface

    圖8 圓柱尾跡中心線上時(shí)均軸向速度分布Fig.8 Mean streamwise velocity along the wake center line

    圖9~圖11為圓柱繞流下游不同截面處的時(shí)均軸向速度分布(圖例見圖8)??梢钥吹捷S向速度流型從近尾跡區(qū)到遠(yuǎn)尾跡區(qū)由U型向V型轉(zhuǎn)變,結(jié)合前部分析可知,剪切層過早分離同樣會(huì)導(dǎo)致速度型的過早轉(zhuǎn)變,圖中兩種實(shí)驗(yàn)值的對(duì)比也說明了這一點(diǎn)。圖9近壁x/D=0.58時(shí),四種模型的速度型均為U型,彼此差別不大,而到回流區(qū)x/D=1.06時(shí),四種模型出現(xiàn)差別,SST模型為完全的V型,SAS模型與兩種PANS模型同樣沒有計(jì)算出更加正確的U型,而是介于U型與V型之間,其中SAS模型與Parnaudeau實(shí)驗(yàn)得到的U型更接近,同時(shí)PANS(des)要稍好于PANS(tanh);在近尾跡回流區(qū)后部x/D=1.54、2.02位置處,除了SST沒有計(jì)算出正確的速度型,SAS模型與兩種PANS模型的計(jì)算結(jié)果只是速度峰值與實(shí)驗(yàn)值存在一定誤差,整體均與實(shí)驗(yàn)值吻合較好;在遠(yuǎn)尾跡區(qū)x/D=4、7、10位置處,四種模型計(jì)算的速度型均較好,速度峰值與實(shí)驗(yàn)存在較小誤差。以上分析說明近尾跡區(qū)域特別是回流區(qū)流動(dòng)較為復(fù)雜、敏感,計(jì)算控制更困難,同時(shí)遠(yuǎn)尾跡區(qū)受近壁區(qū)影響較小,使得遠(yuǎn)尾跡區(qū)計(jì)算誤差小于近尾跡區(qū)。

    圖9 不同截面的平均軸向速度分布(x/D=0.58,1.06)Fig.9 Mean streamwise velocity at different surfaces(x/D=0.58,1.06)

    圖10 不同截面的平均軸向速度分布(x/D=1.54,2.02)Fig.10 Mean streamwise velocity at different surfaces(x/D=1.54,2.02)

    圖11 不同截面的平均軸向速度分布(x/D=4,7,10)Fig.11 Mean streamwise velocity at different surfaces(x/D=4,7,10)

    圖12為圓柱繞流下游近尾跡回流區(qū)附近不同截面處的時(shí)均橫向速度分布(圖例見圖8)。由x/D=1.06處的時(shí)均橫向速度與實(shí)驗(yàn)值可以看到,SAS與兩種PANS模型得到的速度型均與Parnaudeau實(shí)驗(yàn)值較吻合,與SST模型比較顯示出較高的計(jì)算精度。x/D=1.54、2.02處的對(duì)比結(jié)果可以看出四種模型在速度型與峰值上均存在一定誤差,與軸向速度計(jì)算結(jié)果類似。

    圖12 不同截面的平均橫向速度分布(x/D=1.06,1.54,2.02)Fig.12 Mean cross-flow velocity at different surface(x/D=1.06,1.54,2.02)

    3.3 加密網(wǎng)格計(jì)算結(jié)果

    根據(jù)混合模型的構(gòu)造原理可知,其湍流尺度求解能力是基于對(duì)渦黏性的控制作用。圖13為三種混合模型計(jì)算的圓柱繞流時(shí)均湍動(dòng)黏度μt/μ(湍動(dòng)黏度與來流黏度之比)分布,上部為標(biāo)準(zhǔn)網(wǎng)格計(jì)算結(jié)果,下部為加密網(wǎng)格計(jì)算結(jié)果。僅從數(shù)值來看,SAS模型與兩種PANS模型的湍動(dòng)黏度分布關(guān)系為SAS>PANS(des)>PANS(tanh)。PANS模型所求湍動(dòng)黏度要小于SAS模型,這也解釋了圖4中PANS模型對(duì)流場(chǎng)的捕捉要好于SAS模型。但需要注意的是,較低的湍動(dòng)黏度不一定獲得較準(zhǔn)確的數(shù)值計(jì)算結(jié)果。由文中氣動(dòng)力的比較可知,SAS模型要好于PANS模型。通過比較網(wǎng)格加密后的湍流黏度可知,兩種PANS模型在近壁網(wǎng)格密集區(qū)域出現(xiàn)了湍流黏度較低的區(qū)域,這是由于PANS模型的f k構(gòu)造函數(shù)仍是由當(dāng)?shù)赝牧鞒叨群途W(wǎng)格尺度來控制,在近壁網(wǎng)格較密的條件下,PANS模型顯現(xiàn)出類似DES模型對(duì)網(wǎng)格敏感的缺點(diǎn),這時(shí)近壁RANS區(qū)域被破壞,而網(wǎng)格尺度還沒有小到能夠進(jìn)行大渦模擬的程度。

    圖13 時(shí)均湍流黏度比(μt/μ)分布(標(biāo)準(zhǔn)網(wǎng)格:上,加密網(wǎng)格:下)Fig.13 Time-averaged distributions ofμt/μ(baseline grid:top,fine grid:bottom)

    圖14為三種混合模型在不同網(wǎng)格條件下計(jì)算的局部(x/D=1.06)時(shí)均軸向雷諾應(yīng)力分布(圖例見圖15)??梢钥吹絊AS模型及兩種PANS模型計(jì)算的雷諾應(yīng)力與實(shí)驗(yàn)值均有較大差距,這可能是由于渦黏性模型的內(nèi)在缺陷造成的。我們?cè)诖岁P(guān)注不同的模型在不同網(wǎng)格條件下的計(jì)算差異??梢钥吹絊AS模型在標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格條件下均高估了雷諾應(yīng)力值,但兩種網(wǎng)格下的計(jì)算結(jié)果差別不大;兩種PANS模型對(duì)不同的網(wǎng)格表現(xiàn)出相似的規(guī)律,在加密網(wǎng)格時(shí)雷諾應(yīng)力均出現(xiàn)了明顯的降低,結(jié)合圖13中的黏度分布可知,由于近壁區(qū)域網(wǎng)格的加密,湍流黏度降低,導(dǎo)致?;睦字Z應(yīng)力減少,同時(shí)有沒有求解足夠的雷諾應(yīng)力,導(dǎo)致了雷諾應(yīng)力不足現(xiàn)象的發(fā)生,這也驗(yàn)證了PANS構(gòu)造函數(shù)對(duì)網(wǎng)格較敏感的缺陷。同時(shí)兩種PANS模型函數(shù)比較可以看出,PANS(des)模型計(jì)算的雷諾應(yīng)力及氣動(dòng)力分布要稍優(yōu)于PANS(tanh)模型,這也說明了PANS(tanh)模型構(gòu)造函數(shù)對(duì)網(wǎng)格的敏感程度要大于PANS(des)模型的構(gòu)造函數(shù)。

    圖14 局部時(shí)均軸向雷諾應(yīng)力分布Fig.14 Time-averaged distributions of streamwise Reynolds stress

    圖15為三種模型在不同網(wǎng)格條件下的時(shí)均軸向速度分布對(duì)比。結(jié)合前文分析可知,SAS模型網(wǎng)格獨(dú)立性較好,標(biāo)準(zhǔn)網(wǎng)格與加密網(wǎng)格下的計(jì)算結(jié)果沒有體現(xiàn)出太大差別,且由于對(duì)剪切層的計(jì)算更準(zhǔn)確,計(jì)算得到的回流區(qū)大小更接近實(shí)驗(yàn)值。兩種PANS模型由于對(duì)網(wǎng)格較敏感,在加密網(wǎng)格時(shí)延緩了剪切層的轉(zhuǎn)捩,導(dǎo)致回流區(qū)變長(zhǎng)。這在PANS(tanh)模型的計(jì)算結(jié)果中表現(xiàn)更明顯。

    圖15 時(shí)均軸向速度分布Fig.15 Time-averaged distributions of streamwise velocity

    3.4 積分統(tǒng)計(jì)量比較

    表2為四種模型計(jì)算的圓柱繞流平均積分統(tǒng)計(jì)量與相關(guān)實(shí)驗(yàn)值及數(shù)值結(jié)果對(duì)比??梢钥吹皆谂c剪切層運(yùn)動(dòng)密切相關(guān)的參數(shù)中,普通網(wǎng)格下SAS與兩種PANS模型均取得了與實(shí)驗(yàn)值較吻合的回流區(qū)長(zhǎng)度(L r/D)計(jì)算結(jié)果,與構(gòu)造相對(duì)完善的LES-THE、SA-IDDES、V2-f-DES模型的計(jì)算結(jié)果也是相近的,其中SAS模型的計(jì)算結(jié)果與實(shí)驗(yàn)值最接近;加密網(wǎng)格時(shí),兩種PANS模型均過大地預(yù)測(cè)了回流區(qū)長(zhǎng)度。在阻力系數(shù)(C d)及分離角(θSep)的計(jì)算中,在標(biāo)準(zhǔn)網(wǎng)格下三種混合模型計(jì)算誤差均較小,C d滿足實(shí)驗(yàn)值的0.98±0.05區(qū)間,在加密網(wǎng)格時(shí),SAS模型的計(jì)算結(jié)果變化不大,兩種PANS模型獲得了較小的分離角及阻力系數(shù)。在其余的比較項(xiàng)中,可以看到四種模型對(duì)圓柱尾跡旋渦脫落頻率的計(jì)算及壁面壓力的計(jì)算均較準(zhǔn)確,所計(jì)算的斯特勞哈數(shù)(S t)與背風(fēng)點(diǎn)壓力系數(shù)(-C pb)差別不大。對(duì)于SST模型,由于其構(gòu)造?;怂械耐牧髁?導(dǎo)致重要的流場(chǎng)細(xì)節(jié)被抹平,其計(jì)算的阻力系數(shù)、背風(fēng)點(diǎn)壓力系數(shù)、回流區(qū)長(zhǎng)度、分離角均存在較大誤差。同時(shí)可以看到兩種網(wǎng)格條件下PANS(des)計(jì)算結(jié)果稍好于PANS(tanh)。比較文獻(xiàn)[25]中兩種f k為常數(shù)的PANS模型計(jì)算結(jié)果可以看出,f k=0.1時(shí),流場(chǎng)f k值已處于較低水平,但計(jì)算的回流區(qū)長(zhǎng)度要小于本文變f k計(jì)算結(jié)果,而f k=0.3時(shí)的計(jì)算結(jié)果誤差已較大,這說明了變f k在流場(chǎng)計(jì)算中更合理,同時(shí)也說明了流場(chǎng)對(duì)f k值是比較敏感的。

    表2 圓柱繞流平均積分量對(duì)比Table 2 Summary of the mean integrated quantities for the cylinder flow

    3.5 混合模型調(diào)控機(jī)制比較

    由前述分析可知,SAS及PANS等渦黏性混合模型均能通過降低流場(chǎng)湍動(dòng)黏度來求解更多湍流尺度,但SAS及PANS對(duì)求解尺度的調(diào)控機(jī)制是完全不同的。由前文模型構(gòu)造我們可知,SAS模型中引入了von Karman長(zhǎng)度尺度L vk,通過當(dāng)?shù)赝牧鞒叨萀與L vk的比例關(guān)系來控制ω輸運(yùn)方程中源項(xiàng)QSAS的生成量,因此SAS模型在分離區(qū)是通過L/L vk的調(diào)控作用來增大ω的值,從而降低湍動(dòng)黏度的值;PANS模型是通過加入一個(gè)f k參數(shù)來調(diào)控流場(chǎng)可解/不可解湍流量的比例,f k值越小,可解尺度就釋放的越多,文中PANS模型使用了兩種變f k函數(shù),其核心是當(dāng)?shù)赝牧鞒叨扰c網(wǎng)格尺度的比例關(guān)系,兩者都是通過減小ω方程中耗散項(xiàng)的量,從而使求解的ω值增大,進(jìn)而使湍動(dòng)黏度減小。

    圖16為SAS模型計(jì)算的瞬時(shí)及時(shí)均L/L vk分布??梢钥吹皆谖槽E分離區(qū)L/L vk的值較大,而在遠(yuǎn)場(chǎng)區(qū)L/L vk基本為0。圖17和圖18是兩種變f k的PANS模型計(jì)算得到的瞬時(shí)及時(shí)均f k分布,可以看到在分離區(qū)f k的值處在較低水平,而在遠(yuǎn)場(chǎng)范圍f k=1湍流信息均被模化。通過比較兩類模型中對(duì)渦黏度起關(guān)鍵調(diào)控作用的L/L vk與f k分布可以看出,立足于當(dāng)?shù)亓鲃?dòng)、對(duì)網(wǎng)格無依賴的SAS模型局部分辨率更高,可以反映出流場(chǎng)結(jié)構(gòu)的無序與豐富狀態(tài)。PANS模型中網(wǎng)格尺度的關(guān)聯(lián)性使得其局部分辨率降低,且計(jì)算結(jié)果對(duì)網(wǎng)格及計(jì)算初值設(shè)置非常敏感。從以上來看,SAS模型對(duì)流場(chǎng)渦黏性的調(diào)控機(jī)制是有優(yōu)勢(shì)的。

    圖16 SAS模型計(jì)算的瞬時(shí)及時(shí)均L/L vk分布Fig.16 Time-averaged and instantaneous distributions of L/L vk for SAS

    圖17 PANS(des)模型計(jì)算得到的瞬時(shí)及時(shí)均f k分布Fig.17 Time-averaged and instantaneous distributions of f k for PANS(des)

    圖18 PANS(tanh)模型計(jì)算得到的瞬時(shí)及時(shí)均f k分布Fig.18 Time-averaged and instantaneous distributions of f k for PANS(tanh)

    對(duì)比圖17和圖18,顯示了PANS模型的兩種函數(shù)計(jì)算的f k分布差別不大,在近尾跡區(qū)f k均處在較低水平,這有效地調(diào)整了流場(chǎng)的求解尺度。同時(shí)可以看到PANS(des)模型相比PANS(tanh)模型在近尾跡區(qū)的分辨率更高,其在回流區(qū)后部及遠(yuǎn)尾跡區(qū)的f k值要稍小于后者。分析兩種模型的構(gòu)造可知,PANS(des)模型f k函數(shù)構(gòu)造相對(duì)簡(jiǎn)單,但其對(duì)流場(chǎng)的調(diào)控作用更加直接和靈敏,而PANS(tanh)模型中f k函數(shù)中借助tanh函數(shù),tanh函數(shù)整體從-1~1變化,這使得函數(shù)的整體靈敏度下降,從而影響f k的分布。同時(shí)圖13中可以看到PANS(des)的湍動(dòng)黏度稍高于PANS(tanh),其獲得的氣動(dòng)力分布也稍好于后者。這也說明了兩種變f k函數(shù)比較,類DES的f k調(diào)控函數(shù)更好一些,計(jì)算的湍流耗散也更合理。

    4 結(jié) 論

    本文以Re=3900的圓柱繞流作為算例比較了SAS、PANS(des)、PANS(tanh)湍流模型在大分離流動(dòng)中的計(jì)算能力,主要從流場(chǎng)特征捕捉、氣動(dòng)力特性計(jì)算、流場(chǎng)渦黏性控制等方面對(duì)所用湍流模型進(jìn)行了比較,并通過不同的網(wǎng)格計(jì)算條件對(duì)比了不同模型調(diào)控函數(shù)的網(wǎng)格敏感性。計(jì)算結(jié)果表明:三種混合模型均具有求解三維小尺度渦運(yùn)動(dòng)的能力,并能較完整地展現(xiàn)出亞臨界雷諾數(shù)下圓柱繞流層流剪切層形成、發(fā)展、失穩(wěn)、分離,在下游近尾跡區(qū)形成回流區(qū)直至遠(yuǎn)尾跡區(qū)形成剪切分離區(qū)的過程。同時(shí),通過對(duì)比計(jì)算結(jié)果得到以下結(jié)論:

    1)在對(duì)湍流結(jié)構(gòu)的顯示方面,PANS模型相比SAS模型能捕捉到更豐富的湍流結(jié)構(gòu),兩類模型對(duì)流場(chǎng)氣動(dòng)力的計(jì)算均存在一定誤差,回流區(qū)計(jì)算誤差要大于遠(yuǎn)尾跡區(qū),三種模型氣動(dòng)力計(jì)算與實(shí)驗(yàn)吻合關(guān)系為SAS>PANS(des)>PANS(tanh)。

    2)SAS模型L vk尺度立足當(dāng)?shù)亓鲃?dòng),局部流場(chǎng)分辨度更高,調(diào)控作用更強(qiáng);PANS模型對(duì)網(wǎng)格設(shè)置較敏感。兩類模型相比,SAS模型的流場(chǎng)信息調(diào)控模式更優(yōu)秀,獲得的流動(dòng)渦黏性更合理。

    3)兩種PANS模型中,類DES可變f k分布函數(shù)構(gòu)造相對(duì)簡(jiǎn)單,所得f k分布更準(zhǔn)確;借助tanh函數(shù)的f k分布函數(shù)調(diào)節(jié)靈敏度較低,計(jì)算得到的尾跡區(qū)f k值偏低,對(duì)流場(chǎng)調(diào)控作用稍差。

    猜你喜歡
    實(shí)驗(yàn)模型
    一半模型
    記一次有趣的實(shí)驗(yàn)
    微型實(shí)驗(yàn)里看“燃燒”
    重要模型『一線三等角』
    重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
    做個(gè)怪怪長(zhǎng)實(shí)驗(yàn)
    3D打印中的模型分割與打包
    NO與NO2相互轉(zhuǎn)化實(shí)驗(yàn)的改進(jìn)
    實(shí)踐十號(hào)上的19項(xiàng)實(shí)驗(yàn)
    太空探索(2016年5期)2016-07-12 15:17:55
    FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
    亚洲欧美一区二区三区国产| 亚洲av二区三区四区| 亚洲自拍偷在线| 亚洲精品色激情综合| 精品久久久久久久末码| 嘟嘟电影网在线观看| 美女被艹到高潮喷水动态| 熟女电影av网| 草草在线视频免费看| 狂野欧美激情性xxxx在线观看| 国产老妇伦熟女老妇高清| 99久久人妻综合| 久久久欧美国产精品| 欧美三级亚洲精品| 国产 亚洲一区二区三区 | 91久久精品国产一区二区三区| 久久久久九九精品影院| 国产久久久一区二区三区| 高清视频免费观看一区二区 | 一级二级三级毛片免费看| 国产亚洲av嫩草精品影院| 亚洲av不卡在线观看| 欧美变态另类bdsm刘玥| 男人舔奶头视频| 免费不卡的大黄色大毛片视频在线观看 | 欧美日韩视频高清一区二区三区二| 亚洲成人久久爱视频| 观看免费一级毛片| 肉色欧美久久久久久久蜜桃 | 七月丁香在线播放| 国产av国产精品国产| 亚洲图色成人| 国精品久久久久久国模美| 国产乱来视频区| 可以在线观看毛片的网站| 国内精品美女久久久久久| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 精品午夜福利在线看| 少妇猛男粗大的猛烈进出视频 | 黄色配什么色好看| 国产精品av视频在线免费观看| 日日干狠狠操夜夜爽| 亚洲精品中文字幕在线视频 | 国产精品日韩av在线免费观看| 国产精品嫩草影院av在线观看| 国产 一区 欧美 日韩| 国产精品1区2区在线观看.| 又大又黄又爽视频免费| 欧美xxxx黑人xx丫x性爽| 高清毛片免费看| 白带黄色成豆腐渣| 免费电影在线观看免费观看| 国产片特级美女逼逼视频| 夜夜爽夜夜爽视频| 亚洲欧美精品自产自拍| 久久久久久国产a免费观看| 国产人妻一区二区三区在| 日韩视频在线欧美| 亚洲精品自拍成人| 国国产精品蜜臀av免费| 51国产日韩欧美| 精品久久久久久久人妻蜜臀av| 成人毛片a级毛片在线播放| 精品久久久久久久久av| 午夜激情福利司机影院| 黄色一级大片看看| 丰满人妻一区二区三区视频av| 一区二区三区免费毛片| 国产免费又黄又爽又色| 日韩精品青青久久久久久| 国产精品日韩av在线免费观看| 精品一区二区三区视频在线| 午夜免费男女啪啪视频观看| 久久99蜜桃精品久久| 偷拍熟女少妇极品色| 国产亚洲一区二区精品| 伦理电影大哥的女人| 久久精品国产自在天天线| 69av精品久久久久久| 蜜桃久久精品国产亚洲av| 国产精品1区2区在线观看.| 99热网站在线观看| 亚洲国产色片| 国产成人一区二区在线| 国产精品av视频在线免费观看| 国模一区二区三区四区视频| 免费人成在线观看视频色| 国产高潮美女av| 国产国拍精品亚洲av在线观看| 日韩亚洲欧美综合| 久久这里有精品视频免费| 中文字幕av成人在线电影| 国产精品三级大全| 天堂俺去俺来也www色官网 | 寂寞人妻少妇视频99o| 亚洲精品色激情综合| 男人爽女人下面视频在线观看| 亚洲av中文字字幕乱码综合| 久久国产乱子免费精品| 一区二区三区免费毛片| 99久久精品一区二区三区| freevideosex欧美| 深爱激情五月婷婷| 亚洲图色成人| 22中文网久久字幕| 一夜夜www| 中文字幕免费在线视频6| 不卡视频在线观看欧美| 国产一区二区在线观看日韩| 五月天丁香电影| 国产午夜精品一二区理论片| 欧美日韩精品成人综合77777| 亚洲精品影视一区二区三区av| 国产91av在线免费观看| 最近中文字幕2019免费版| 在线免费观看的www视频| 在线 av 中文字幕| 亚洲精品一二三| 久久国产乱子免费精品| 久久久久九九精品影院| 久久久欧美国产精品| 日韩欧美精品v在线| 成人午夜高清在线视频| 午夜精品一区二区三区免费看| 天天躁夜夜躁狠狠久久av| 狂野欧美激情性xxxx在线观看| 女人被狂操c到高潮| 伦理电影大哥的女人| 亚洲一级一片aⅴ在线观看| 中文乱码字字幕精品一区二区三区 | 热99在线观看视频| 精品一区二区三区视频在线| 99re6热这里在线精品视频| 最新中文字幕久久久久| 国内精品宾馆在线| 中国美白少妇内射xxxbb| 久久精品国产自在天天线| 亚洲综合精品二区| 国产伦理片在线播放av一区| 国产 亚洲一区二区三区 | 亚洲四区av| 久久久久久久久中文| 亚洲真实伦在线观看| 久久综合国产亚洲精品| 久久久久久久亚洲中文字幕| 亚洲av中文av极速乱| 天美传媒精品一区二区| 尤物成人国产欧美一区二区三区| 欧美精品一区二区大全| 两个人视频免费观看高清| 伊人久久国产一区二区| 国产极品天堂在线| 一区二区三区乱码不卡18| 成人一区二区视频在线观看| 精品国产三级普通话版| 中文乱码字字幕精品一区二区三区 | 亚洲一级一片aⅴ在线观看| 日日摸夜夜添夜夜添av毛片| 天堂影院成人在线观看| 国产成人精品久久久久久| 男女边吃奶边做爰视频| 久热久热在线精品观看| 男的添女的下面高潮视频| 欧美性感艳星| 久久久久久久亚洲中文字幕| 五月伊人婷婷丁香| 亚洲成人久久爱视频| 日韩,欧美,国产一区二区三区| 亚洲av.av天堂| 国产美女午夜福利| 成人亚洲精品一区在线观看 | 国产伦一二天堂av在线观看| av免费在线看不卡| 成人鲁丝片一二三区免费| 三级男女做爰猛烈吃奶摸视频| 九色成人免费人妻av| 天堂俺去俺来也www色官网 | 最后的刺客免费高清国语| 天堂av国产一区二区熟女人妻| 国产午夜精品久久久久久一区二区三区| 91久久精品国产一区二区成人| 国产色爽女视频免费观看| 亚洲精品久久久久久婷婷小说| 一二三四中文在线观看免费高清| 成人美女网站在线观看视频| 久久99精品国语久久久| 国产精品熟女久久久久浪| 国产亚洲精品av在线| www.色视频.com| 亚洲精品第二区| 女人久久www免费人成看片| 午夜激情久久久久久久| 亚洲久久久久久中文字幕| 国产高潮美女av| 国产黄片视频在线免费观看| 亚洲人成网站在线观看播放| 我的女老师完整版在线观看| 国产v大片淫在线免费观看| 亚洲电影在线观看av| 听说在线观看完整版免费高清| 亚洲国产精品国产精品| 日本色播在线视频| 最近2019中文字幕mv第一页| 一级毛片黄色毛片免费观看视频| 美女xxoo啪啪120秒动态图| 亚洲国产精品国产精品| 亚洲精品日韩在线中文字幕| 2021少妇久久久久久久久久久| 欧美日本视频| 亚洲av福利一区| 又爽又黄无遮挡网站| 国产精品久久久久久久电影| av在线播放精品| 国产免费又黄又爽又色| 亚洲精品成人久久久久久| 一本一本综合久久| 成人亚洲欧美一区二区av| 国国产精品蜜臀av免费| 国产成人a区在线观看| 欧美日韩在线观看h| 亚洲综合精品二区| 三级经典国产精品| 搡老乐熟女国产| 日韩欧美国产在线观看| eeuss影院久久| 在线免费十八禁| 黄色配什么色好看| 青青草视频在线视频观看| 亚洲在线自拍视频| 国内揄拍国产精品人妻在线| 国产69精品久久久久777片| 亚洲国产高清在线一区二区三| 国产精品精品国产色婷婷| 中文字幕人妻熟人妻熟丝袜美| 国产精品国产三级专区第一集| 能在线免费看毛片的网站| 97超碰精品成人国产| 高清日韩中文字幕在线| 国产乱来视频区| 日本免费在线观看一区| 国产成人a∨麻豆精品| 日韩欧美一区视频在线观看 | 老女人水多毛片| 亚洲乱码一区二区免费版| 日本一二三区视频观看| 波野结衣二区三区在线| 日日干狠狠操夜夜爽| 久久久久久久亚洲中文字幕| 色综合亚洲欧美另类图片| 麻豆精品久久久久久蜜桃| 全区人妻精品视频| 一个人观看的视频www高清免费观看| 国产激情偷乱视频一区二区| 日韩三级伦理在线观看| 午夜久久久久精精品| 亚洲精品,欧美精品| 99久国产av精品国产电影| 一个人看视频在线观看www免费| 777米奇影视久久| 欧美极品一区二区三区四区| 在线观看av片永久免费下载| 国产老妇伦熟女老妇高清| 亚洲av二区三区四区| 免费看美女性在线毛片视频| 午夜视频国产福利| 在线免费观看的www视频| 久久久久久久大尺度免费视频| 成人毛片a级毛片在线播放| 日韩欧美精品免费久久| 免费无遮挡裸体视频| 精品一区在线观看国产| 免费人成在线观看视频色| 国产精品久久久久久精品电影| 国产精品伦人一区二区| av播播在线观看一区| 日韩亚洲欧美综合| 亚洲第一区二区三区不卡| 91久久精品国产一区二区成人| 久久亚洲国产成人精品v| xxx大片免费视频| 欧美高清性xxxxhd video| 久久99精品国语久久久| 国产人妻一区二区三区在| 美女脱内裤让男人舔精品视频| 久久99热这里只频精品6学生| 天堂av国产一区二区熟女人妻| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 日韩欧美国产在线观看| 精品久久久久久电影网| 久热久热在线精品观看| 91av网一区二区| 一级黄片播放器| 国产成人精品久久久久久| 亚洲高清免费不卡视频| 成人亚洲精品一区在线观看 | 五月伊人婷婷丁香| 欧美丝袜亚洲另类| 久久久久久九九精品二区国产| 亚洲综合精品二区| 熟女电影av网| 最近手机中文字幕大全| 美女脱内裤让男人舔精品视频| 国精品久久久久久国模美| 91精品一卡2卡3卡4卡| 国产亚洲91精品色在线| 日韩欧美国产在线观看| 神马国产精品三级电影在线观看| 亚洲国产欧美人成| 男女边摸边吃奶| 婷婷色综合www| 国内精品一区二区在线观看| 国产综合精华液| 中文字幕av在线有码专区| 日韩三级伦理在线观看| 国产成人a区在线观看| 草草在线视频免费看| 99久国产av精品| 国产精品国产三级专区第一集| 午夜爱爱视频在线播放| 国产高潮美女av| 2021少妇久久久久久久久久久| .国产精品久久| 女人十人毛片免费观看3o分钟| 国产伦在线观看视频一区| 色综合站精品国产| 小蜜桃在线观看免费完整版高清| 国产亚洲91精品色在线| 免费人成在线观看视频色| 日韩三级伦理在线观看| 午夜精品一区二区三区免费看| 高清日韩中文字幕在线| 真实男女啪啪啪动态图| 色吧在线观看| 丝袜喷水一区| 97超视频在线观看视频| 国产精品麻豆人妻色哟哟久久 | 午夜福利成人在线免费观看| 亚洲精品日韩av片在线观看| 国产v大片淫在线免费观看| 日本与韩国留学比较| 少妇人妻精品综合一区二区| 亚洲精品视频女| 国产精品久久久久久精品电影小说 | 色尼玛亚洲综合影院| 男女下面进入的视频免费午夜| 婷婷色av中文字幕| 国产 一区精品| 亚洲一区高清亚洲精品| 日韩欧美精品v在线| 韩国高清视频一区二区三区| 亚洲自拍偷在线| 色视频www国产| 婷婷色麻豆天堂久久| 国产黄色视频一区二区在线观看| 男女边摸边吃奶| 最近2019中文字幕mv第一页| 神马国产精品三级电影在线观看| 国产亚洲一区二区精品| 亚洲国产精品sss在线观看| 人妻系列 视频| 国产精品国产三级专区第一集| 在线免费观看不下载黄p国产| 天天躁夜夜躁狠狠久久av| 成人av在线播放网站| 1000部很黄的大片| 51国产日韩欧美| 日韩av不卡免费在线播放| 中文乱码字字幕精品一区二区三区 | 成人午夜高清在线视频| 国产精品.久久久| 欧美激情在线99| 色综合色国产| 好男人视频免费观看在线| 亚洲av二区三区四区| av卡一久久| 成年版毛片免费区| 午夜福利在线观看吧| 国产老妇女一区| 成人无遮挡网站| 插逼视频在线观看| 美女cb高潮喷水在线观看| 九九在线视频观看精品| 国语对白做爰xxxⅹ性视频网站| 精品国产露脸久久av麻豆 | 大香蕉久久网| 美女内射精品一级片tv| 日本-黄色视频高清免费观看| 禁无遮挡网站| 久久久久久九九精品二区国产| 亚洲成人久久爱视频| 久久精品人妻少妇| 国产69精品久久久久777片| 亚洲自拍偷在线| 国产一级毛片在线| 日本午夜av视频| 精品99又大又爽又粗少妇毛片| 亚洲美女视频黄频| 国产亚洲av片在线观看秒播厂 | 中文字幕制服av| 美女主播在线视频| av黄色大香蕉| 欧美精品国产亚洲| 亚洲av免费在线观看| 国产精品麻豆人妻色哟哟久久 | av在线亚洲专区| 一个人免费在线观看电影| 成人高潮视频无遮挡免费网站| 两个人视频免费观看高清| 亚洲色图av天堂| 欧美日韩国产mv在线观看视频 | 欧美不卡视频在线免费观看| 一级毛片我不卡| 亚洲欧美精品自产自拍| 51国产日韩欧美| av网站免费在线观看视频 | 你懂的网址亚洲精品在线观看| 婷婷色综合www| 1000部很黄的大片| 三级国产精品片| 亚洲成人精品中文字幕电影| 九九爱精品视频在线观看| 菩萨蛮人人尽说江南好唐韦庄| 国产乱人偷精品视频| 免费不卡的大黄色大毛片视频在线观看 | 黄片wwwwww| 内地一区二区视频在线| 在线免费观看的www视频| 啦啦啦韩国在线观看视频| 成人特级av手机在线观看| 91狼人影院| 狂野欧美激情性xxxx在线观看| 亚洲,欧美,日韩| 国精品久久久久久国模美| 国产男女超爽视频在线观看| 精品欧美国产一区二区三| 男女边摸边吃奶| 干丝袜人妻中文字幕| 看非洲黑人一级黄片| 直男gayav资源| 国产欧美另类精品又又久久亚洲欧美| 国产精品久久视频播放| 最近中文字幕2019免费版| 日韩亚洲欧美综合| 汤姆久久久久久久影院中文字幕 | 国内精品美女久久久久久| 最近2019中文字幕mv第一页| av免费观看日本| 亚洲精品一区蜜桃| 一级片'在线观看视频| 日韩在线高清观看一区二区三区| 69av精品久久久久久| 国产有黄有色有爽视频| 在线 av 中文字幕| 日日啪夜夜撸| 狠狠精品人妻久久久久久综合| 菩萨蛮人人尽说江南好唐韦庄| 亚洲在线观看片| a级毛片免费高清观看在线播放| 亚洲无线观看免费| 六月丁香七月| 亚洲电影在线观看av| 成人av在线播放网站| av播播在线观看一区| 欧美成人午夜免费资源| 国产亚洲91精品色在线| 国产伦一二天堂av在线观看| 久久午夜福利片| 永久免费av网站大全| 国产精品国产三级国产专区5o| 亚洲一级一片aⅴ在线观看| 一边亲一边摸免费视频| 午夜精品一区二区三区免费看| 国产精品爽爽va在线观看网站| 国产成人freesex在线| 亚洲婷婷狠狠爱综合网| ponron亚洲| 久久草成人影院| 美女国产视频在线观看| 可以在线观看毛片的网站| 亚洲av一区综合| 国产高清不卡午夜福利| 日韩 亚洲 欧美在线| 亚洲欧美日韩无卡精品| 国产 一区 欧美 日韩| 精品一区在线观看国产| 久久久久久久久久久免费av| 在线 av 中文字幕| 深夜a级毛片| 亚洲av在线观看美女高潮| 久久久久久久久久黄片| 国产精品人妻久久久久久| 亚洲欧美日韩卡通动漫| 永久网站在线| 麻豆av噜噜一区二区三区| 亚洲av成人精品一二三区| 国产免费又黄又爽又色| 99九九线精品视频在线观看视频| av福利片在线观看| 97精品久久久久久久久久精品| 欧美bdsm另类| 成人一区二区视频在线观看| 久久99热这里只频精品6学生| 国产精品久久久久久久电影| 国产精品一二三区在线看| 久久久成人免费电影| 日韩欧美一区视频在线观看 | 国产午夜精品论理片| 亚洲精品影视一区二区三区av| 日韩欧美 国产精品| 嘟嘟电影网在线观看| 爱豆传媒免费全集在线观看| 国产成人一区二区在线| 亚洲欧美清纯卡通| 亚洲欧美精品专区久久| 中文字幕人妻熟人妻熟丝袜美| av一本久久久久| 网址你懂的国产日韩在线| 精品一区二区三区视频在线| 美女主播在线视频| 国产淫语在线视频| 欧美日韩视频高清一区二区三区二| 99视频精品全部免费 在线| 亚洲精品第二区| 久久99热这里只有精品18| 亚洲丝袜综合中文字幕| 91久久精品国产一区二区成人| 永久免费av网站大全| 晚上一个人看的免费电影| 国产亚洲91精品色在线| 国内精品一区二区在线观看| 麻豆成人午夜福利视频| 最新中文字幕久久久久| 天堂av国产一区二区熟女人妻| 亚洲国产精品成人综合色| av在线蜜桃| 特大巨黑吊av在线直播| 国产黄a三级三级三级人| 免费看光身美女| 亚洲精品日本国产第一区| 婷婷色综合大香蕉| 成人欧美大片| 天天躁日日操中文字幕| 亚洲av国产av综合av卡| 蜜桃久久精品国产亚洲av| 国产精品日韩av在线免费观看| 三级经典国产精品| 亚洲精品成人av观看孕妇| 日本wwww免费看| 亚洲av中文字字幕乱码综合| 国产精品久久视频播放| 人人妻人人澡欧美一区二区| 最近最新中文字幕大全电影3| 免费看美女性在线毛片视频| 亚洲欧美精品自产自拍| 特大巨黑吊av在线直播| 夫妻午夜视频| 免费av不卡在线播放| 99热6这里只有精品| 日日撸夜夜添| 欧美激情国产日韩精品一区| 好男人在线观看高清免费视频| 少妇被粗大猛烈的视频| 80岁老熟妇乱子伦牲交| 大片免费播放器 马上看| 2022亚洲国产成人精品| av在线老鸭窝| 国产高清有码在线观看视频| 成人毛片a级毛片在线播放| 亚洲美女搞黄在线观看| 国产视频首页在线观看| 亚洲国产高清在线一区二区三| 国产一级毛片在线| 国产精品三级大全| 成人亚洲欧美一区二区av| 久久精品久久久久久久性| 黑人高潮一二区| 久久鲁丝午夜福利片| 亚洲成人久久爱视频| 免费av毛片视频| 国产精品99久久久久久久久| 国产淫片久久久久久久久| 草草在线视频免费看| 中国美白少妇内射xxxbb| 亚洲av免费高清在线观看| 男女下面进入的视频免费午夜| 人人妻人人看人人澡| 男插女下体视频免费在线播放| 久久国产乱子免费精品| 欧美一级a爱片免费观看看| 在线观看免费高清a一片| 亚洲激情五月婷婷啪啪| 真实男女啪啪啪动态图| 国产av国产精品国产| 亚洲av中文av极速乱| 日本一本二区三区精品| 免费观看性生交大片5| 久久久久久九九精品二区国产| 亚洲国产精品国产精品| 一个人免费在线观看电影| 国产永久视频网站| 国产免费一级a男人的天堂| 51国产日韩欧美| 久久精品久久久久久噜噜老黄| 久久精品熟女亚洲av麻豆精品 | 国产一区二区三区综合在线观看 | 国产亚洲av嫩草精品影院| 少妇丰满av| 亚洲精品一二三| 好男人在线观看高清免费视频| 草草在线视频免费看| 特级一级黄色大片| 亚洲成人久久爱视频| 在线免费观看的www视频| 欧美不卡视频在线免费观看| 国产精品一二三区在线看|