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

    串列雙圓柱繞流問(wèn)題數(shù)值模擬的多尺度分析

    2017-10-13 03:36:04崔雪揚(yáng)曹博超徐弘一
    關(guān)鍵詞:雷諾數(shù)邊界層圓柱

    崔雪揚(yáng),曹博超,徐弘一

    (復(fù)旦大學(xué) 航空航天系,上海 200433)

    串列雙圓柱繞流問(wèn)題數(shù)值模擬的多尺度分析

    崔雪揚(yáng),曹博超,徐弘一

    (復(fù)旦大學(xué) 航空航天系,上海 200433)

    運(yùn)用有限體積法,對(duì)串列放置的雙圓柱二維不可壓縮流動(dòng)進(jìn)行了直接數(shù)值計(jì)算.在分析Strouhal數(shù)及升阻力系數(shù)等積分量的基礎(chǔ)上,本文從流動(dòng)多尺度層面研究了流場(chǎng)分布及渦結(jié)構(gòu).通過(guò)對(duì)速度場(chǎng)的流動(dòng)顯示和頻譜分析,確定了渦脫落的多個(gè)頻率,以及分別受上游圓柱和下游圓柱邊界層擾動(dòng)形成的多尺度渦的相互作用,并且發(fā)現(xiàn)由于多尺度渦的相互作用形成了更小尺度渦的過(guò)程及機(jī)理.最后,進(jìn)一步將這些渦結(jié)構(gòu)與流場(chǎng)模態(tài)對(duì)應(yīng)起來(lái),使得流場(chǎng)結(jié)構(gòu)更加清晰地展現(xiàn)出來(lái).

    串列雙圓柱; 直接數(shù)值模擬; 多尺度; 渦脫落; 脫落頻率

    流體繞流問(wèn)題是工程中的常見問(wèn)題,各種高大建筑物設(shè)計(jì)、飛行器設(shè)計(jì)、海洋石油平臺(tái)及輸油管道建設(shè)均涉及此類流動(dòng).此外,工業(yè)設(shè)備中也普遍存在繞流現(xiàn)象,如各類管殼式換熱器等.因此研究鈍體繞流的特性對(duì)工程實(shí)際和工業(yè)設(shè)備設(shè)計(jì)尤其重要,相關(guān)課題也一直是流體研究的熱點(diǎn).一個(gè)多世紀(jì)以來(lái),圓柱繞流一直是眾多理論分析、實(shí)驗(yàn)研究及數(shù)值模擬的對(duì)象,但迄今對(duì)該流動(dòng)現(xiàn)象物理本質(zhì)的理解仍未完整.

    大多數(shù)實(shí)驗(yàn)研究的雷諾數(shù)高于104量級(jí),而數(shù)值模擬大都在102量級(jí)進(jìn)行.盡管計(jì)算中雷諾數(shù)遠(yuǎn)小于實(shí)際情況,但實(shí)驗(yàn)結(jié)果表明,伴隨渦脫落的大尺寸尾流在一定程度上呈高、低雷諾數(shù)(102量級(jí))相似的特征[1].因此,通過(guò)對(duì)低雷諾數(shù)問(wèn)題的研究亦可基本反映高雷諾數(shù)流動(dòng)的主要現(xiàn)象.雖然實(shí)際情形為3維流動(dòng),但2維問(wèn)題的研究亦可反映該流動(dòng)的一些主要特征.

    在實(shí)驗(yàn)研究方面,Roshoko[2]最早在實(shí)驗(yàn)室中發(fā)現(xiàn)圓柱繞流存在與流場(chǎng)雷諾數(shù)相關(guān)的3個(gè)尾流狀態(tài);Zdravkovich[3]曾對(duì)兩圓柱串列和交錯(cuò)放置的繞流問(wèn)題進(jìn)行過(guò)實(shí)驗(yàn)研究,針對(duì)兩圓柱中心間距小于5倍圓柱直徑的一系列情況,他研究了兩圓柱間的流動(dòng)相互作用,發(fā)現(xiàn)中心間距小于3倍圓柱直徑時(shí),沒(méi)有明顯的渦自上游圓柱脫落.

    隨著電子計(jì)算機(jī)的飛速發(fā)展,數(shù)值計(jì)算方面的研究也層出不窮,Phuocloc和Bouard[4]數(shù)值模擬繞流圓柱初期流動(dòng)的2次渦結(jié)構(gòu),將得到的結(jié)果和實(shí)驗(yàn)比對(duì),非常接近.鄧見等[5]使用分塊耦合方法,對(duì)單圓柱和串列雙圓柱繞流進(jìn)行了數(shù)值模擬,研究分析了雙柱中心間距變化對(duì)上下游圓柱升阻力系數(shù)的影響,計(jì)算結(jié)果與試驗(yàn)結(jié)果非常吻合.廖俊等[6]使用表面渦法研究了高雷諾數(shù)下排列方式不同的雙圓柱繞流,計(jì)算了圓柱在串列、并列以及錯(cuò)置情況下的各種流動(dòng)結(jié)構(gòu)、渦街變化及作用在圓柱上的力,所得的阻力系數(shù)與試驗(yàn)結(jié)果符合很好.

    然而到目前為止,對(duì)圓柱繞流問(wèn)題的研究大多集中在整體平均量層面,如同一雷諾數(shù)下探究阻力系數(shù)、升力系數(shù)和Strouhal數(shù)與間距直徑比的關(guān)系,或者將數(shù)值模擬得到的流場(chǎng)信息,如流函數(shù)、渦量、渦脫落頻率等與實(shí)驗(yàn)結(jié)果進(jìn)行比較.隨著科學(xué)研究的不斷深入,對(duì)多尺度下流動(dòng)規(guī)律的研究變得越來(lái)越重要.目前國(guó)內(nèi)外對(duì)多尺度下低雷諾數(shù)圓柱繞流的研究還相對(duì)初步,有必要對(duì)低雷諾數(shù)下的圓柱繞流開展深入研究.

    本文用有限體積法將非定常Navier-Stokes(N-S)方程進(jìn)行離散,在求解域中采用直接數(shù)值模擬(Direct Numerical Simulation, DNS)方法計(jì)算得到流場(chǎng)信息,應(yīng)用后處理軟件Fieldview進(jìn)行流場(chǎng)顯示,根據(jù)流動(dòng)結(jié)構(gòu)和圓柱后渦脫落尾跡分布,取相關(guān)位置點(diǎn)進(jìn)行流動(dòng)的多尺度分析,以研究渦脫落頻率、邊界層擾流等流動(dòng)基本特征,并結(jié)合模態(tài)和渦量圖深入解析渦脫落機(jī)理.

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

    1.1控制方程及其差分離散格式

    計(jì)算機(jī)的快速發(fā)展,使得DNS方法可在相對(duì)簡(jiǎn)單的流動(dòng)幾何構(gòu)型中得以應(yīng)用.但此方法要求網(wǎng)格點(diǎn)的數(shù)量多,計(jì)算量巨大.由于本文計(jì)算研究的是相對(duì)較低的雷諾數(shù)(Re=800),使得DNS方法得以應(yīng)用并得到較為精確的結(jié)果.

    研究中的物理模型為2維不可壓縮非定常流動(dòng)的N-S方程:

    式中:u,v是x,y方向的速度;μ為動(dòng)力粘度;ρ為流體密度;t為計(jì)算總時(shí)長(zhǎng);p為作用在研究對(duì)象上的壓力.

    為獲取具有普遍應(yīng)用價(jià)值的流動(dòng)數(shù)據(jù),需要對(duì)方程進(jìn)行無(wú)量綱化處理,所取特征物理量及相關(guān)無(wú)量綱物理量如下:

    無(wú)量綱化后的2維不可壓縮非定常流動(dòng)的N-S方程可改寫為:

    在對(duì)動(dòng)量方程的離散中,對(duì)流項(xiàng)和擴(kuò)撒項(xiàng)采用空間2階中心格式.在時(shí)間方向,采用2階精度的時(shí)間推進(jìn),對(duì)流項(xiàng)為2階Adams-Bashforth顯示格式,擴(kuò)散項(xiàng)則采用2階Adams-Moulton隱式格式.則離散方程可寫為:

    空間離散網(wǎng)格采用交錯(cuò)網(wǎng)格布局,如圖1(見第386頁(yè))所示,即速度分量(u,v)和壓力p各占據(jù)一套網(wǎng)格,這樣就可以避免壓力速度定義在同一點(diǎn)產(chǎn)生的壓力棋盤格式,也可以避免連續(xù)方程通量計(jì)算中的差值守恒性和Poisson方程邊界條件問(wèn)題.動(dòng)量與壓力的解耦求解采用分步方法[7].

    圖1 u方向的交錯(cuò)網(wǎng)格及u分量控制體Fig.1 Staggered grid in u direction and u component control body

    由上可知,得到u方向的動(dòng)量方程:

    其中φ與壓力項(xiàng)p有關(guān),即

    1.2計(jì)算域

    圖2 計(jì)算域展示圖Fig.2 Computational domain sketch map

    在本次模擬計(jì)算中,經(jīng)無(wú)量綱化處理將圓柱直徑設(shè)為1,選取長(zhǎng)l=16,寬d=4的長(zhǎng)方形計(jì)算域,串列雙圓柱的間距比為2,上游圓柱圓心位置與進(jìn)口邊界的距離為4,如圖2所示.

    在以往的計(jì)算中,研究者大多將計(jì)算區(qū)域網(wǎng)格采用分塊的結(jié)構(gòu)化網(wǎng)格,圓柱周圍計(jì)算區(qū)域采用O型網(wǎng)格,由于渦的產(chǎn)生、發(fā)展和脫落的整個(gè)過(guò)程均在圓柱表面進(jìn)行,圓柱表面附近區(qū)域的計(jì)算網(wǎng)格劃分和網(wǎng)格質(zhì)量對(duì)圓柱表面的計(jì)算結(jié)果會(huì)產(chǎn)生較大的影響,因此采用O型網(wǎng)格的優(yōu)勢(shì)是能夠更好地模擬圓柱體表面上渦的產(chǎn)生和發(fā)展過(guò)程.但是在圓柱表面區(qū)域生成正交性較好的高質(zhì)量網(wǎng)格卻比較困難.

    由于模擬中的求解域和圓柱截面的網(wǎng)格均為比較規(guī)整的形狀(Cartesian正交網(wǎng)格),且網(wǎng)格點(diǎn)數(shù)較高(1026×514),使得圓柱表面計(jì)算精度比較高.整體網(wǎng)格在空間為非均勻分布,在圓柱周圍進(jìn)行了加密,如圖3(a)所示的網(wǎng)格整體效果,圖3(b)為將網(wǎng)格放大后的效果.可以看到圓柱附近網(wǎng)格的不均勻性.由于計(jì)算的雷諾數(shù)較低,使得圓柱附近的邊界層不需要太密的網(wǎng)格分辨率,因此使得計(jì)算結(jié)果具有物理可靠性.

    圖3 計(jì)算網(wǎng)格和圓柱邊界層附近計(jì)算網(wǎng)格的放大圖Fig.3 Computational grid and large version of computational grid near the cylinder boundary layer

    1.3邊界條件

    入口邊界條件: 入流為均勻流,u=u0=1,pi+1,j-pi,j=0;

    出口邊界條件: 出口充分發(fā)展ui+1,j-ui,j=0,vi+1,j-vi,j=0;

    圓柱表面無(wú)滑移邊界,給定速度和壓力,u=0,v=0;

    外壁面無(wú)滑移邊界,給定速度,(u,v)=0.

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

    本論文主要研究Re=800時(shí),來(lái)流為層流時(shí)串列雙圓柱的繞流情況,分析展示了流動(dòng)分離、漩渦生成和脫落、漩渦之間以及漩渦-壁面的互相干擾等流體力學(xué)現(xiàn)象,得到非定常流動(dòng)失穩(wěn)后的一系列固有頻率與本征正交分解(Proper Orthogonal Decomposition, POD)模態(tài)分析等物理量的對(duì)應(yīng)關(guān)系.計(jì)算考慮來(lái)流為均勻流動(dòng),總時(shí)間推進(jìn)長(zhǎng)度為300,時(shí)間步長(zhǎng)為Δt=0.0006,推進(jìn)時(shí)間步為500000步.圖4展示了充分發(fā)展時(shí)刻的流場(chǎng)瞬間結(jié)構(gòu),以作為對(duì)串列雙圓柱流場(chǎng)的感性認(rèn)知.

    圖4 均勻來(lái)流的串列雙圓柱繞流的兩個(gè)瞬時(shí)刻流場(chǎng)結(jié)構(gòu)Fig.4 Two moments of instantaneous flow structure about double cylinders in tandem as incoming flow is uniform

    圖5 串列雙圓柱求解域布點(diǎn)情況Fig.5 Domain distribution of double cylinders in tandem

    本文選取的圓柱間距比為2.由于求解域?yàn)樯舷聦?duì)稱,因此只需在對(duì)稱軸以上布點(diǎn)(如圖5),從多尺度層面去探究圓柱繞流脫落事件與頻率的關(guān)系.從左到右依次記為1,2,3,…,12列,共12×12個(gè)研究點(diǎn).

    2.1來(lái)流為均勻流時(shí)的流場(chǎng)頻譜分析

    2.1.1 上下游圓柱之間流場(chǎng)信息分析

    由圖4的流動(dòng)顯示圖可知,因?yàn)樯嫌螆A柱邊界層擾動(dòng)影響,會(huì)在兩圓柱之間產(chǎn)生旋渦.文獻(xiàn)[8]顯示,當(dāng)Re=60~200時(shí),雙柱串列繞流存在一個(gè)臨界間距比(l/d=3.6),小于臨界間距時(shí),上游圓柱不存在渦脫落.在本算例中根據(jù)流動(dòng)顯示圖及后續(xù)的證明,得到當(dāng)Re=800,兩圓柱間距比為2時(shí),上游圓柱也不存在渦脫落.在兩圓柱的共同作用下,該區(qū)域形成有規(guī)律的渦結(jié)構(gòu).對(duì)兩圓柱中間的3列點(diǎn)進(jìn)行分析得到速度場(chǎng)信息和速度場(chǎng)的頻譜分析.選擇其中3個(gè)典型位置: 第4列第1個(gè)點(diǎn)(0.68,0)、第5列第2個(gè)點(diǎn)(0.99,0.14)、第6列第3個(gè)點(diǎn)(1.29,0.29),得到如圖6(見第388頁(yè))所示的u和v方向的速度頻譜分析圖.x軸為頻率,記為f(=t-1);y軸為振幅,u方向的振幅記為Au,v方向的振幅記為Av.

    圖6 點(diǎn)(0.68,0),(0.99,0.14),(1.29,0.29)u和v方向的頻譜分析Fig.6 (0.68,0),(0.99,0.14),(1.29,0.29) spectrum analysis in u and v direction

    以上的頻譜分析結(jié)果可解釋圖4中關(guān)于兩圓柱之間的流動(dòng)顯示觀察.從頻譜分析圖中可明顯看出0.22為主導(dǎo)頻率,0.44為次主導(dǎo)頻率.再結(jié)合u,v方向的頻譜分析,發(fā)現(xiàn)0.22,0.44的頻率會(huì)在不同的位置起不同的主導(dǎo)作用.上游圓柱的邊界層擾動(dòng)是上下游圓柱之間的渦形成的主要原因,但是渦在此位置的擺動(dòng)則是上下游圓柱共同作用的結(jié)果.

    從圖4中可以看出,上下游圓柱之間的渦同時(shí)在轉(zhuǎn)動(dòng)和上下擺動(dòng),因此生成的渦不是簡(jiǎn)單一個(gè)大尺度渦,而是在震蕩的過(guò)程中也產(chǎn)生了小尺度的旋渦.大渦頻率為0.22,小渦頻率為0.44.圖4呈現(xiàn)出的結(jié)果與頻譜分析得出的結(jié)果基本相同.另外,從頻譜分析圖上也可以看到還有一個(gè)比較明顯的頻率0.65,能量不多但也不足以忽視,說(shuō)明除此兩個(gè)渦之外仍有一些細(xì)碎的小渦形成,頻率在0.65附近.

    2.1.2 下游圓柱尾流流場(chǎng)信息分析

    文獻(xiàn)[9]顯示,雷諾數(shù)達(dá)到150和200時(shí),圓柱尾流區(qū)出現(xiàn)明顯的渦脫落,在圓柱后的尾流區(qū)形成了穩(wěn)定的渦街結(jié)構(gòu),渦街排列規(guī)則.本文的雷諾數(shù)為800,圖4的流動(dòng)顯示圖明顯展現(xiàn)了下游圓柱后面有渦的脫落且排列規(guī)則.由流體力學(xué)知識(shí)可知,流體繞圓柱流動(dòng)時(shí),過(guò)流斷面收縮,流速沿程增加,壓強(qiáng)沿程減小,而由于粘性力的存在,柱體周圍形成附面層分離,隨雷諾數(shù)增加,在柱體后半部分形成不同的繞流尾跡現(xiàn)象,可推測(cè)下游渦會(huì)同時(shí)受到上下游圓柱邊界層擾動(dòng)的影響.

    對(duì)下游圓柱后面4列點(diǎn)做出u-t,v-t,p-t和對(duì)應(yīng)的頻譜分析,如圖7所示.結(jié)合流動(dòng)顯示圖4,近下游圓柱尾流區(qū)的渦已經(jīng)形成但是還未脫落.從圖7中可以看到,0.22為最明顯的主導(dǎo)頻率,0.44和0.65為次要頻率.下游圓柱后面的渦發(fā)生了脫落,說(shuō)明0.22為脫落頻率.0.44和0.65分別是受上游圓柱和下游圓柱邊界層擾動(dòng)形成的套在大渦里的小渦.理由如下: 1) 上下游圓柱之間的主頻出現(xiàn)了0.44,但并沒(méi)有出現(xiàn)0.65;2) 下游圓柱上方的位置主頻出現(xiàn)0.44而沒(méi)有0.65;3) 下游圓柱主頻同時(shí)包括有0.44和0.65.兩圓柱之間位置的渦只受上游圓柱的影響,而下游圓柱后面的渦同時(shí)受到上下游圓柱邊界層擾動(dòng)的影響.

    遠(yuǎn)離下游圓柱的尾渦區(qū),渦已經(jīng)脫落,渦經(jīng)過(guò)的位置都有3個(gè)主要的頻率,分別為0.22,0.44,0.65.同樣0.22為渦的脫落頻率,由前面的分析可知,脫落渦中有兩個(gè)頻率為0.44和0.65的小渦,分別是受上游圓柱和下游圓柱的影響.

    由于渦不僅產(chǎn)生了脫落,而且受到上下游圓柱的影響也產(chǎn)生了小渦結(jié)構(gòu).同時(shí)在圖7中發(fā)現(xiàn)小渦頻率除0.44和0.65之外,還有0.91比較突出.說(shuō)明在這些事件的共同影響下使得流動(dòng)現(xiàn)象更加復(fù)雜,出現(xiàn)了其他頻率的小渦結(jié)構(gòu).選擇其中3個(gè)典型位置: (2.72,0.73),(3.92,0.14),(6.17,0.14),得到如圖7所示的u和v方向的速度頻譜分析圖.

    圖7 點(diǎn)(2.72,0.73),(3.92,0.14),(6.17,0.14)u和v方向的頻譜分析Fig.7 (2.72,0.73),(3.92,0.14),(6.17,0.14) spectrum analysis in u and v direction

    圖7展示了下游圓柱尾流區(qū)特定位置渦的特點(diǎn),在研究中也發(fā)現(xiàn),同一列12個(gè)點(diǎn)的位置渦的能量呈現(xiàn)一定的變化規(guī)律,如圖8展示.以x=3.92列為例,比較同一列上的位置,發(fā)現(xiàn)u方向上主頻率的能量從標(biāo)記1處向上開始逐漸增加,到3處到達(dá)最大之后衰減;v方向上主頻率的能量從標(biāo)記1處達(dá)到最大,向上開始逐漸衰減.再以x=6.17列為例,u方向上主頻率的能量從標(biāo)記1處向上開始逐漸增加,到5處到達(dá)最大之后衰減;v方向上主頻率的能量從標(biāo)記1處達(dá)到最大,向上開始逐漸衰減.說(shuō)明渦脫落之后的軌跡并不是一條水平線,且其水平方向的能量最大值隨渦的移動(dòng)而變動(dòng),豎直方向的能量隨位置的升高而降低.可推測(cè)水平方向能量的最大值在渦的中心位置.在x=6.17列的折線圖中看到最低點(diǎn)之后有回升,結(jié)合流動(dòng)顯示圖,說(shuō)明求解域的邊界層在此處已經(jīng)變寬,x=6.17列標(biāo)記12已經(jīng)觸及求解域的邊界層.在這篇論文中不對(duì)邊界層做深入的研究.

    圖8 x=3.92和x=6.17列主頻能量的變化Fig.8 Line x=3.92 and x=6.17 energy change of dominant frequency

    2.1.3 結(jié)果對(duì)比

    在文獻(xiàn)[1]中,劉松等在Re=200,間距比為2,串列雙圓柱繞流的情況下,數(shù)值模擬得到Strouhal數(shù)為0.19,同等條件下文獻(xiàn)[10]的實(shí)驗(yàn)結(jié)果中Strouhal數(shù)為0.13,而在本文中Re=800,間距比為2,由St=fD/u0(式中St為斯特勞哈爾數(shù),f為頻率,D為圓柱直徑,u0為初速度)[10]得到Strouhal數(shù)為0.22,在同一數(shù)量級(jí),這說(shuō)明本文的模擬計(jì)算結(jié)果與其他類似研究的可比性及可靠性.

    2.2均勻來(lái)流的流場(chǎng)模態(tài)分析

    圖9 1~20階模態(tài)能量衰減圖Fig.9 1—20 mode energy attenuation map

    上文在多尺度層面解釋了串列雙圓柱在Re=800,來(lái)流為均勻流時(shí)的流動(dòng)結(jié)構(gòu)以及影響原因.由于湍流多尺度的復(fù)雜性,為了研究湍流的物理機(jī)制,研究人員常用模態(tài)分解的辦法來(lái)提取對(duì)流場(chǎng)發(fā)展有重要影響的相干結(jié)構(gòu).POD法是按照能量的大小對(duì)各模態(tài)進(jìn)行排列,提取流場(chǎng)中較大能量的結(jié)構(gòu).POD目前已用于研究多種流動(dòng)問(wèn)題,例如低雷諾數(shù)下圓柱繞流的動(dòng)態(tài)特性.在本篇論文中,將多尺度微觀層面流場(chǎng)結(jié)構(gòu)分析和POD法模態(tài)分析結(jié)合,來(lái)解釋串列雙圓柱繞流的流動(dòng)現(xiàn)象.

    將計(jì)算得到的前20階模態(tài)能量作線形圖,如圖9所示,可以看到前2階能量大幅衰減.

    根據(jù)能量占比計(jì)算得到,前8階能量已經(jīng)占到總能量的94%,所以取前8階模態(tài)為POD分析的對(duì)象,如圖10所示.

    圖10 1~8階模態(tài)渦量云圖Fig.10 1—8 vortex cloud chart

    圖10是串列雙圓柱繞流前8階模態(tài)的渦量云圖,根據(jù)圖10,1,2,4,6,7階模態(tài)在垂直于流動(dòng)方向上為正對(duì)稱結(jié)構(gòu),3,5,8階模態(tài)在垂直于流動(dòng)方向上為反對(duì)稱結(jié)構(gòu).和圖4的流場(chǎng)結(jié)構(gòu)顯示圖結(jié)合可以從整體上對(duì)該流場(chǎng)有直觀的感受,可以看到流場(chǎng)渦尺度的多樣化和渦結(jié)構(gòu)的復(fù)雜.下面將通過(guò)頻譜分析將模態(tài)、渦結(jié)構(gòu)以及流場(chǎng)的影響因素聯(lián)系起來(lái).圖11為1~8階模態(tài)的速度頻譜圖.

    圖11 1~8階模態(tài)u方向速度頻譜圖Fig.11 1—8 mode in u direction velocity frequency spectrum

    1,2階模態(tài)中只有主頻率0.22的渦,且第1階模態(tài)包含88.5%的能量,前2階模態(tài)包含90.8%的能量,說(shuō)明渦脫落發(fā)生在第1階模態(tài),2階模態(tài)相對(duì)于1階模態(tài)多一些擾動(dòng).3階模態(tài)中的主頻為0.44和0.91,說(shuō)明上游圓柱邊界層擾動(dòng)的影響主要發(fā)生在3階模態(tài)中.4階模態(tài)包含4個(gè)主頻,說(shuō)明4階模態(tài)涵蓋了下游圓柱尾流渦結(jié)構(gòu)的所有信息,也體現(xiàn)了上下游渦共同作用的結(jié)果.除4階模態(tài)之外,7階模態(tài)和8階模態(tài)也都包含了4個(gè)主頻,但明顯7階模態(tài)中頻率0.91的影響比較小,說(shuō)明渦脫落時(shí)結(jié)構(gòu)中主要包含3個(gè)渦,而在8階模態(tài)中頻率0.65的能量相對(duì)頻率0.44對(duì)應(yīng)的能量更大些,說(shuō)明在此階模態(tài)中,下游圓柱邊界層影響效果更大.5階模態(tài)體現(xiàn)了上下游圓柱中間位置渦的結(jié)構(gòu),主要受上游圓柱邊界層擾動(dòng)的影響,頻率0.22能量占比并非最大.6階模態(tài)中體現(xiàn)了下游圓柱邊界層擾動(dòng)對(duì)渦結(jié)構(gòu)的影響,主要頻率為0.22,0.65.這8階模態(tài)中由不同原因造成的渦結(jié)構(gòu)的不同,其流場(chǎng)信息結(jié)構(gòu)也可由圖10的云圖定性地看出.表1為1~8階模態(tài)的能量占比信息及發(fā)生渦的主要頻率.

    表1 1~8階模態(tài)的能量占比信息及發(fā)生渦的主要頻率

    結(jié)合表1以及各階模態(tài)系數(shù)頻譜分析圖11,由于2~8階模態(tài)的能量占比非常小,可知一階模態(tài)的能量即為渦脫落的能量,即主頻0.22的能量占88.5%.上下游圓柱擾動(dòng)影響造成的次主要渦結(jié)構(gòu)總能量占比不到5.5%,其中受下游圓柱層流擾動(dòng)影響形成的渦能量占比最大,即頻率0.65的能量占比排第二,頻率為0.44和0.91的能量占比相近.

    3 結(jié)論與展望

    本文從多尺度層面研究了串列雙圓柱的繞流問(wèn)題,分析了Re=800,來(lái)流為均勻流時(shí)的復(fù)雜流場(chǎng)結(jié)構(gòu),并采用POD方法將能量絕對(duì)占優(yōu)的前8階模態(tài)進(jìn)行剝離,從而探求渦流場(chǎng)結(jié)構(gòu)及其成因以及與相關(guān)模態(tài)之間的關(guān)系.總體而言,本文做了以下的工作并得出相應(yīng)結(jié)論:

    1) 間距比為2時(shí),上下游圓柱之間的渦不會(huì)脫落,且通過(guò)頻譜分析此時(shí)除了主頻為0.22的大渦之外,內(nèi)部還套有主頻為0.44的小渦,且由上游圓柱邊界層擾動(dòng)產(chǎn)生.

    2)Re=800來(lái)流為均勻流時(shí)流場(chǎng)結(jié)構(gòu)已經(jīng)呈現(xiàn)多尺度的復(fù)雜性.在下游圓柱的后面,出現(xiàn)了渦的脫落且脫落頻率為0.22;除此之外,還出現(xiàn)了多尺度的小渦,受上游圓柱邊界層擾動(dòng)產(chǎn)生的頻率為0.44的小渦和受下游圓柱邊界層擾動(dòng)產(chǎn)生的頻率為0.65的小渦,以及相互之間作用產(chǎn)生的更微尺度的頻率為0.91的渦.

    3) 將模態(tài)與流場(chǎng)結(jié)構(gòu)對(duì)應(yīng)起來(lái),發(fā)現(xiàn)由于不同原因產(chǎn)生的渦有些只發(fā)生在某一個(gè)模態(tài)中,有些則在多個(gè)模態(tài)中發(fā)生.前8階模態(tài)的能量占據(jù)總能量的94%.

    本文只是對(duì)Re=800,均勻來(lái)流的情況加以分析,對(duì)于此雷諾數(shù)下來(lái)流為湍流的以及高雷諾數(shù)比如Re=1600,來(lái)流分別為均勻流和充分發(fā)展湍流的情況,后續(xù)仍需要大量的工作對(duì)其進(jìn)行分析以及對(duì)比,才能更深入地從微尺度層面將串列圓柱的繞流問(wèn)題研究清楚.

    [1] 劉松,符松.串列雙圓柱繞流問(wèn)題的數(shù)值模擬 [J].計(jì)算力學(xué)學(xué)報(bào),2000,17(3): 260-266.

    [2] ROSHOKO A. On the development of turbulent wakes from vortex streets [R]. Kitty Hawk, North Carolina, USA: National Advisory Committee for Aeronautics, 1954.

    [3] ZDRAVKOVICH M. Review of flow interference between two circular cylinders in various arrangements [J].ASMEJournalofFluidsEngineering, 1977,99(4): 618-633.

    [4] PHUOCLOC T, BOUARD R.Numerical solution of the early stage of unsteady viscous flow around a circular cylinder: A comparison with experimental visualization and measurements [J].JournalofFluidMechanics,1985,160: 93-117.

    [5] 鄧見,黃鈺期,任安祿.分塊法研究圓柱繞流升阻力 [J].力學(xué)與實(shí)踐,2004,26(1): 24-26.

    [6] 廖俊,景思睿.高雷諾數(shù)下雙圓柱繞流的數(shù)值模擬 [J].水動(dòng)力學(xué)研究與進(jìn)展,2001,16(1): 101-110.

    [7] KIM J, MOIN P. Application of a fractional-step method to incompressible Navier-Stokes equations [J].JournalofComputationalPhysics, 1985,59(2): 308-323.

    [8] CARMO B S, MENEGHINI J R. Numerical investigation of the flow around two circular cylinders in tandem [J].JournalofFluidsandStructures, 2006,22(6): 979-988.

    [9] 洪文鵬,周云龍,劉巍.圓柱繞流漩渦脫落的數(shù)值模擬 [C]∥朱德祥.第七屆全國(guó)水動(dòng)力學(xué)學(xué)術(shù)會(huì)議暨第十九屆全國(guó)水動(dòng)力學(xué)研討會(huì)文集(上冊(cè)).北京: 海洋出版社,2005: 6.

    [10] SLAOUTI A, STANSBY P K. Flow around two circular cylinders by the random-vortex method [J].JournalofFluidsandStructures, 1992,6(6): 641-633.

    Abstract: Two-dimensional flow around double cylinders in tandem was calculated, using the finite volume method. Based on the comparisons of the Strouhal number and the lift and drag coefficients with the relevant researches, the flow field and its vortex structure were further investigated. Flow visualizations were conducted and the spectrum and POD analyses were performed to study the vortex shedding and their respective interactions.The formation of the smaller vortex scales were found due to the vortex interactions. The vortex structure and the modes of flow-field exhibited their distinctive structures in terms of their symmetrically and unsymmetrically characteristic distributions in space.Keywords: double cylinders in tandem; direct numerical simulation; multi-scale flow; vortex shedding; shedding frequency

    Multi-scaleAnalysisofNumericalSimulationaboutFlowaroundDoubleCylindersinTandem

    CUI Xueyang, CAO Bochao, XU Hongyi

    (DepartmentofAeronauticsandAstronautics,FudanUniversity,Shanghai200433,China)

    F830.9

    A

    0427-7104(2017)03-0384-09

    2016-11-24

    國(guó)家自然科學(xué)基金重大研究計(jì)劃培育項(xiàng)目(91434112);上海千人計(jì)劃啟動(dòng)項(xiàng)目(EZH2126503);上海教育委員會(huì)、中航商業(yè)發(fā)動(dòng)機(jī)、復(fù)旦大學(xué)數(shù)值仿真聯(lián)合創(chuàng)新中心項(xiàng)目(AR908.D1RW.002)

    崔雪揚(yáng)(1991—),女,碩士研究生;徐弘一(1963—),男,研究員,通信聯(lián)系人,E-mail: hongyi_xu@fudan.edu.cn.

    猜你喜歡
    雷諾數(shù)邊界層圓柱
    工程學(xué)和圓柱
    圓柱的體積計(jì)算
    基于HIFiRE-2超燃發(fā)動(dòng)機(jī)內(nèi)流道的激波邊界層干擾分析
    基于Transition SST模型的高雷諾數(shù)圓柱繞流數(shù)值研究
    削法不同 體積有異
    失穩(wěn)初期的低雷諾數(shù)圓柱繞流POD-Galerkin 建模方法研究
    基于轉(zhuǎn)捩模型的低雷諾數(shù)翼型優(yōu)化設(shè)計(jì)研究
    民機(jī)高速風(fēng)洞試驗(yàn)的阻力雷諾數(shù)效應(yīng)修正
    一類具有邊界層性質(zhì)的二次奇攝動(dòng)邊值問(wèn)題
    非特征邊界的MHD方程的邊界層
    久久人人爽av亚洲精品天堂 | 国产v大片淫在线免费观看| 最近的中文字幕免费完整| 老熟女久久久| 联通29元200g的流量卡| 亚洲成色77777| 亚洲精品自拍成人| 亚洲精品国产色婷婷电影| 国产精品久久久久久精品电影小说 | 在线观看免费视频网站a站| 高清视频免费观看一区二区| 成人综合一区亚洲| 热re99久久精品国产66热6| 成年av动漫网址| 国产v大片淫在线免费观看| 97在线人人人人妻| 男人添女人高潮全过程视频| 久久亚洲国产成人精品v| 亚洲人成网站高清观看| 国产精品爽爽va在线观看网站| 国产精品国产三级国产av玫瑰| 久久6这里有精品| 最近2019中文字幕mv第一页| 国产一区有黄有色的免费视频| 色5月婷婷丁香| 肉色欧美久久久久久久蜜桃| 亚洲电影在线观看av| 国产爱豆传媒在线观看| 亚洲欧美精品自产自拍| 日产精品乱码卡一卡2卡三| 亚洲av国产av综合av卡| 成人免费观看视频高清| 少妇的逼水好多| 亚洲最大成人中文| 22中文网久久字幕| 久久久久精品久久久久真实原创| 一级毛片 在线播放| 多毛熟女@视频| 日本爱情动作片www.在线观看| 久久人人爽av亚洲精品天堂 | 日韩亚洲欧美综合| 国产精品一区二区在线不卡| 高清午夜精品一区二区三区| 看免费成人av毛片| 人人妻人人看人人澡| 国产又色又爽无遮挡免| 欧美人与善性xxx| 亚洲无线观看免费| 看免费成人av毛片| 熟女电影av网| 寂寞人妻少妇视频99o| 精品一区二区免费观看| 在线观看三级黄色| 久久综合国产亚洲精品| 看免费成人av毛片| 女人久久www免费人成看片| 欧美xxⅹ黑人| 亚洲av国产av综合av卡| 最近最新中文字幕免费大全7| 一区二区三区乱码不卡18| 精品人妻一区二区三区麻豆| 国产免费一级a男人的天堂| 国产一区亚洲一区在线观看| 国产美女午夜福利| 赤兔流量卡办理| 亚洲一区二区三区欧美精品| 嫩草影院新地址| 特大巨黑吊av在线直播| 亚洲色图av天堂| 草草在线视频免费看| 久久韩国三级中文字幕| 成人美女网站在线观看视频| 2022亚洲国产成人精品| 免费观看性生交大片5| 中文字幕亚洲精品专区| 亚洲电影在线观看av| 大码成人一级视频| 午夜激情福利司机影院| 伊人久久精品亚洲午夜| 深夜a级毛片| 国产一区二区三区av在线| 99久久人妻综合| 精品人妻一区二区三区麻豆| 伊人久久国产一区二区| 日韩 亚洲 欧美在线| 韩国av在线不卡| 国内精品宾馆在线| 纯流量卡能插随身wifi吗| 国产成人精品一,二区| 日韩,欧美,国产一区二区三区| 午夜福利网站1000一区二区三区| 精品国产乱码久久久久久小说| 国产精品免费大片| 亚洲婷婷狠狠爱综合网| 久久精品夜色国产| av在线app专区| 黑丝袜美女国产一区| 各种免费的搞黄视频| 国产片特级美女逼逼视频| 久久这里有精品视频免费| 国产黄色免费在线视频| 国产高清国产精品国产三级 | 国产爱豆传媒在线观看| 18禁在线播放成人免费| 老熟女久久久| 制服丝袜香蕉在线| 噜噜噜噜噜久久久久久91| 国产精品久久久久久精品古装| 黄色一级大片看看| 日韩三级伦理在线观看| 国产黄色视频一区二区在线观看| 免费看不卡的av| 纵有疾风起免费观看全集完整版| 国产成人午夜福利电影在线观看| 天堂俺去俺来也www色官网| 久久久久久久大尺度免费视频| 久久久久久久大尺度免费视频| 久久久久人妻精品一区果冻| 七月丁香在线播放| 亚洲av.av天堂| 黄色怎么调成土黄色| 久久国内精品自在自线图片| 男人和女人高潮做爰伦理| 熟妇人妻不卡中文字幕| 欧美人与善性xxx| 黄片无遮挡物在线观看| 久久久国产一区二区| 大香蕉97超碰在线| 欧美极品一区二区三区四区| 97超碰精品成人国产| 国产av码专区亚洲av| 欧美精品一区二区免费开放| 日日撸夜夜添| av卡一久久| 亚洲av综合色区一区| 成人午夜精彩视频在线观看| 亚洲精品亚洲一区二区| 欧美国产精品一级二级三级 | 舔av片在线| 超碰av人人做人人爽久久| 亚洲精品第二区| 一级av片app| 精品亚洲成国产av| 国产伦在线观看视频一区| 精品国产露脸久久av麻豆| 一区二区三区精品91| 高清日韩中文字幕在线| 日韩成人av中文字幕在线观看| 狂野欧美激情性bbbbbb| 国产高清三级在线| av.在线天堂| 一本—道久久a久久精品蜜桃钙片| 一本久久精品| 国产高清三级在线| 在线天堂最新版资源| 我的女老师完整版在线观看| 午夜福利影视在线免费观看| 晚上一个人看的免费电影| 免费看光身美女| 深夜a级毛片| 国产av国产精品国产| 99热全是精品| 精品亚洲成国产av| 不卡视频在线观看欧美| 亚洲成人av在线免费| 在线观看免费日韩欧美大片 | 国产深夜福利视频在线观看| 国产色婷婷99| 人妻制服诱惑在线中文字幕| 交换朋友夫妻互换小说| 网址你懂的国产日韩在线| 成人国产麻豆网| 99久久综合免费| 伊人久久国产一区二区| 色网站视频免费| 少妇人妻 视频| 亚洲国产最新在线播放| 一级爰片在线观看| 寂寞人妻少妇视频99o| 精品久久久久久久久av| av不卡在线播放| h视频一区二区三区| 国产爱豆传媒在线观看| 色婷婷av一区二区三区视频| 国产成人免费观看mmmm| 99九九线精品视频在线观看视频| 99久久精品国产国产毛片| 国产精品久久久久久精品古装| 日韩一区二区视频免费看| 夜夜看夜夜爽夜夜摸| 91aial.com中文字幕在线观看| 亚洲天堂av无毛| 精品久久久久久久久av| 一级毛片我不卡| 婷婷色av中文字幕| 免费人成在线观看视频色| 女性生殖器流出的白浆| 久久久午夜欧美精品| 最近最新中文字幕免费大全7| 国产成人精品婷婷| 激情五月婷婷亚洲| 中文乱码字字幕精品一区二区三区| 少妇精品久久久久久久| 国产成人一区二区在线| 日本与韩国留学比较| 黄色怎么调成土黄色| 精品一区二区三区视频在线| 成年av动漫网址| 亚洲性久久影院| 日本-黄色视频高清免费观看| 亚洲美女黄色视频免费看| 国产男人的电影天堂91| 天天躁日日操中文字幕| 人妻少妇偷人精品九色| av黄色大香蕉| 久久久久国产网址| 午夜福利在线在线| 国产男人的电影天堂91| videossex国产| 亚洲精品456在线播放app| 欧美变态另类bdsm刘玥| 两个人的视频大全免费| 日日撸夜夜添| 天美传媒精品一区二区| 亚洲欧美中文字幕日韩二区| 纯流量卡能插随身wifi吗| 日韩成人伦理影院| 国产乱人视频| 777米奇影视久久| 一本久久精品| 欧美成人午夜免费资源| 国产探花极品一区二区| 亚洲精品乱码久久久久久按摩| 草草在线视频免费看| 国产精品爽爽va在线观看网站| 六月丁香七月| 国产精品久久久久久精品古装| 日日啪夜夜爽| 看免费成人av毛片| 自拍偷自拍亚洲精品老妇| 我的女老师完整版在线观看| 青春草亚洲视频在线观看| 黄色视频在线播放观看不卡| 久久精品国产亚洲网站| 亚洲人成网站在线观看播放| 日韩成人伦理影院| 日韩av不卡免费在线播放| 精品国产三级普通话版| 国产视频首页在线观看| 女性被躁到高潮视频| 亚洲人成网站高清观看| 国产精品久久久久久久电影| 精品国产三级普通话版| 亚洲精品,欧美精品| 丰满人妻一区二区三区视频av| 中国国产av一级| 99热全是精品| 午夜视频国产福利| 午夜福利网站1000一区二区三区| 国产一区二区在线观看日韩| 国产高清有码在线观看视频| 成人二区视频| 插逼视频在线观看| 国产精品一及| 五月天丁香电影| 简卡轻食公司| 在线观看一区二区三区激情| 久久久久久久久久人人人人人人| 亚洲av成人精品一二三区| 亚洲va在线va天堂va国产| 91精品一卡2卡3卡4卡| av播播在线观看一区| 国产在视频线精品| 国产av一区二区精品久久 | 国国产精品蜜臀av免费| 99热这里只有精品一区| 国产免费一级a男人的天堂| 欧美少妇被猛烈插入视频| 乱系列少妇在线播放| 国产成人午夜福利电影在线观看| 亚洲精品日韩av片在线观看| 亚洲成人手机| 免费人成在线观看视频色| 免费人妻精品一区二区三区视频| 好男人视频免费观看在线| 亚洲av成人精品一区久久| 亚洲欧美精品自产自拍| 亚洲精品国产成人久久av| 免费看av在线观看网站| 久久毛片免费看一区二区三区| 免费高清在线观看视频在线观看| 国产av国产精品国产| 乱码一卡2卡4卡精品| 日韩人妻高清精品专区| 男女国产视频网站| 亚洲色图av天堂| 国产欧美另类精品又又久久亚洲欧美| 亚洲欧美一区二区三区国产| 国产伦在线观看视频一区| 国产精品av视频在线免费观看| 国产精品嫩草影院av在线观看| 国产又色又爽无遮挡免| 最新中文字幕久久久久| 色婷婷久久久亚洲欧美| 久久亚洲国产成人精品v| 色婷婷av一区二区三区视频| 下体分泌物呈黄色| 人人妻人人爽人人添夜夜欢视频 | 深爱激情五月婷婷| 午夜免费男女啪啪视频观看| 亚洲内射少妇av| 国产日韩欧美在线精品| 国产一区亚洲一区在线观看| 在线观看一区二区三区激情| 欧美精品国产亚洲| 国产免费又黄又爽又色| 国产男人的电影天堂91| 国产淫语在线视频| 少妇的逼好多水| 国产色爽女视频免费观看| 日本wwww免费看| 亚洲av不卡在线观看| 国产精品福利在线免费观看| 亚洲成人一二三区av| 亚洲真实伦在线观看| 国产 精品1| 1000部很黄的大片| 成年人午夜在线观看视频| 夜夜看夜夜爽夜夜摸| 欧美日韩国产mv在线观看视频 | 久久久久精品性色| 欧美精品一区二区大全| 狂野欧美激情性bbbbbb| 久久久午夜欧美精品| 国产黄色视频一区二区在线观看| 亚洲精品乱码久久久久久按摩| 下体分泌物呈黄色| 91午夜精品亚洲一区二区三区| 国产片特级美女逼逼视频| 国产免费福利视频在线观看| 亚洲国产精品成人久久小说| 国产乱人视频| 国产精品国产三级国产专区5o| 91午夜精品亚洲一区二区三区| 国产人妻一区二区三区在| 纵有疾风起免费观看全集完整版| 亚洲欧美清纯卡通| 欧美日本视频| 久久久久久久久久久丰满| 黄色配什么色好看| 爱豆传媒免费全集在线观看| 亚洲国产精品成人久久小说| 一个人看视频在线观看www免费| 国产成人免费观看mmmm| 秋霞在线观看毛片| 五月开心婷婷网| 丰满乱子伦码专区| 久久精品国产鲁丝片午夜精品| 午夜日本视频在线| 卡戴珊不雅视频在线播放| 亚洲三级黄色毛片| 国产一区二区在线观看日韩| 国产精品嫩草影院av在线观看| 不卡视频在线观看欧美| av又黄又爽大尺度在线免费看| a 毛片基地| 夜夜看夜夜爽夜夜摸| 一本久久精品| 1000部很黄的大片| 丝袜脚勾引网站| 久久热精品热| 欧美xxxx性猛交bbbb| 亚洲av欧美aⅴ国产| 国产精品一区www在线观看| 久久热精品热| 国产乱来视频区| 午夜福利在线在线| 99热全是精品| 日本午夜av视频| 亚洲精品国产av蜜桃| 只有这里有精品99| 嘟嘟电影网在线观看| 欧美日韩视频精品一区| 韩国高清视频一区二区三区| 欧美bdsm另类| 久久久久久久国产电影| 久久久成人免费电影| 亚洲av在线观看美女高潮| 久久久久久伊人网av| av在线老鸭窝| 日日撸夜夜添| tube8黄色片| 七月丁香在线播放| 日韩电影二区| 男女国产视频网站| 精品少妇黑人巨大在线播放| 国产高清有码在线观看视频| 国产精品99久久久久久久久| 亚洲色图综合在线观看| 欧美日韩综合久久久久久| 国模一区二区三区四区视频| 亚洲电影在线观看av| 春色校园在线视频观看| 99精国产麻豆久久婷婷| 在线 av 中文字幕| 亚洲精品亚洲一区二区| 欧美3d第一页| 久久久久精品久久久久真实原创| 高清在线视频一区二区三区| 国产欧美另类精品又又久久亚洲欧美| 一个人免费看片子| 丰满乱子伦码专区| 久久亚洲国产成人精品v| 欧美性感艳星| 网址你懂的国产日韩在线| 美女福利国产在线 | 高清日韩中文字幕在线| 亚洲av成人精品一区久久| 日本猛色少妇xxxxx猛交久久| 成人午夜精彩视频在线观看| 免费大片18禁| 人人妻人人爽人人添夜夜欢视频 | 王馨瑶露胸无遮挡在线观看| 免费av不卡在线播放| 卡戴珊不雅视频在线播放| 99视频精品全部免费 在线| 亚洲va在线va天堂va国产| 亚洲欧美精品自产自拍| 人妻少妇偷人精品九色| 婷婷色麻豆天堂久久| 秋霞在线观看毛片| 五月伊人婷婷丁香| 直男gayav资源| 成人国产av品久久久| 中文字幕精品免费在线观看视频 | 国产v大片淫在线免费观看| 欧美bdsm另类| 观看免费一级毛片| 97在线视频观看| 中文资源天堂在线| 欧美zozozo另类| 网址你懂的国产日韩在线| 99久久综合免费| 日本免费在线观看一区| 亚洲精品日本国产第一区| 亚洲av欧美aⅴ国产| 久久精品国产亚洲网站| 看非洲黑人一级黄片| 蜜桃久久精品国产亚洲av| 色婷婷av一区二区三区视频| 国产日韩欧美在线精品| 国国产精品蜜臀av免费| 97超碰精品成人国产| 日韩成人av中文字幕在线观看| 国产精品.久久久| 插阴视频在线观看视频| 91久久精品电影网| a级毛色黄片| 看非洲黑人一级黄片| 一个人免费看片子| 激情五月婷婷亚洲| 久久鲁丝午夜福利片| 一区二区三区免费毛片| 国产免费福利视频在线观看| 免费在线观看成人毛片| 99热国产这里只有精品6| 丝瓜视频免费看黄片| av免费观看日本| 国产中年淑女户外野战色| .国产精品久久| 亚洲av电影在线观看一区二区三区| 在线精品无人区一区二区三 | 亚洲无线观看免费| 丰满人妻一区二区三区视频av| 中文字幕制服av| 亚洲真实伦在线观看| 欧美精品亚洲一区二区| 黑丝袜美女国产一区| 午夜免费鲁丝| 国产精品不卡视频一区二区| 爱豆传媒免费全集在线观看| 蜜桃久久精品国产亚洲av| 免费高清在线观看视频在线观看| 久久久久网色| 日本-黄色视频高清免费观看| 亚洲国产精品成人久久小说| 亚洲成人中文字幕在线播放| 国产色爽女视频免费观看| 欧美国产精品一级二级三级 | 精品人妻视频免费看| 亚洲国产成人一精品久久久| 黑人高潮一二区| 久久毛片免费看一区二区三区| 亚洲av电影在线观看一区二区三区| 大香蕉久久网| 久久久色成人| av在线观看视频网站免费| 免费黄网站久久成人精品| 爱豆传媒免费全集在线观看| 久久亚洲国产成人精品v| 好男人视频免费观看在线| 中文字幕av成人在线电影| 午夜精品国产一区二区电影| 在线观看三级黄色| 不卡视频在线观看欧美| 女性被躁到高潮视频| 亚洲熟女精品中文字幕| 欧美一级a爱片免费观看看| 美女xxoo啪啪120秒动态图| 国产精品一区二区在线观看99| 精品人妻熟女av久视频| 最新中文字幕久久久久| 国产色爽女视频免费观看| 大码成人一级视频| 联通29元200g的流量卡| 亚洲内射少妇av| 亚洲丝袜综合中文字幕| 夫妻午夜视频| 尤物成人国产欧美一区二区三区| 少妇人妻久久综合中文| 老熟女久久久| 日日撸夜夜添| 新久久久久国产一级毛片| 中国三级夫妇交换| 99热这里只有是精品50| 亚洲国产欧美在线一区| 精品一区二区三卡| 久久人人爽人人爽人人片va| 国产成人精品一,二区| 精品久久久久久久久av| 熟女av电影| 久久久久久久精品精品| 久久女婷五月综合色啪小说| 国模一区二区三区四区视频| 99视频精品全部免费 在线| 国模一区二区三区四区视频| 国产男女超爽视频在线观看| 美女cb高潮喷水在线观看| 一级毛片黄色毛片免费观看视频| 亚洲va在线va天堂va国产| 国产爱豆传媒在线观看| 久久精品国产亚洲网站| 街头女战士在线观看网站| 男人狂女人下面高潮的视频| 日本色播在线视频| 国产精品蜜桃在线观看| 视频区图区小说| 在线免费观看不下载黄p国产| 国产乱人视频| 午夜激情福利司机影院| 在线观看av片永久免费下载| 下体分泌物呈黄色| 国产永久视频网站| 久久久久性生活片| 久久热精品热| 三级国产精品片| 日本色播在线视频| 亚洲不卡免费看| 一级毛片 在线播放| 伊人久久国产一区二区| 免费观看a级毛片全部| 久久久久久久国产电影| 七月丁香在线播放| 联通29元200g的流量卡| 搡老乐熟女国产| 久久鲁丝午夜福利片| 丰满少妇做爰视频| 两个人的视频大全免费| 韩国av在线不卡| 男女下面进入的视频免费午夜| 亚洲精品乱码久久久v下载方式| 中国三级夫妇交换| 久热这里只有精品99| 一级二级三级毛片免费看| 国产精品一区二区三区四区免费观看| 精品久久久久久久末码| 亚洲欧美一区二区三区黑人 | 亚洲精品国产成人久久av| 精品国产露脸久久av麻豆| 亚洲精品456在线播放app| 国产国拍精品亚洲av在线观看| 色综合色国产| 三级国产精品片| 国产亚洲91精品色在线| 亚洲经典国产精华液单| 日本黄色片子视频| 久久人人爽人人爽人人片va| 黄色日韩在线| 美女高潮的动态| 日韩av免费高清视频| 一个人看的www免费观看视频| 国产精品一二三区在线看| 啦啦啦在线观看免费高清www| 国产黄片美女视频| 欧美丝袜亚洲另类| 国产亚洲91精品色在线| 国产精品一区二区在线观看99| 久久精品久久久久久久性| 亚洲av男天堂| 成年女人在线观看亚洲视频| 青青草视频在线视频观看| 国产av精品麻豆| 嫩草影院入口| 欧美高清性xxxxhd video| 国产精品一及| 久久婷婷青草| 一级二级三级毛片免费看| 成人高潮视频无遮挡免费网站| 好男人视频免费观看在线| 1000部很黄的大片| 国产69精品久久久久777片| 亚洲伊人久久精品综合| 欧美97在线视频| 国产综合精华液| 国产av码专区亚洲av|