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

    譜單元法及其在多圓柱繞流分析中的應(yīng)用

    2014-04-30 02:29:12韓兆龍陳亞楠桂曉瀾李俊龍
    空氣動力學(xué)學(xué)報 2014年1期
    關(guān)鍵詞:斯托旋渦尾流

    韓兆龍,周 岱,陳亞楠,桂曉瀾,李俊龍

    (上海交通大學(xué) 船舶海洋與建筑工程學(xué)院,上海 200240)

    譜單元法及其在多圓柱繞流分析中的應(yīng)用

    韓兆龍,周 岱,陳亞楠,桂曉瀾,李俊龍

    (上海交通大學(xué) 船舶海洋與建筑工程學(xué)院,上海 200240)

    系統(tǒng)闡明譜單元方法,基于譜單元方法對低雷諾數(shù)Re=200時不同間距下的順排兩圓柱和Re=150正方形排列的四圓柱繞流及其阻力系數(shù)、升力系數(shù)等進(jìn)行數(shù)值模擬。研究比較不同間距比L/D(兩圓柱圓心距離與圓柱直徑之比)對兩圓柱和四圓柱繞流的影響,計算分析了渦量圖分布、平均阻力系數(shù)和斯托羅哈數(shù)隨間距比的變化。研究表明,間距比對順排兩圓柱和正方形四圓柱繞流影響顯著;順排兩圓柱繞流存在臨界間距比,在Re=200時臨界間距比約為3.6。正方形排列四圓柱存在三種流態(tài)。當(dāng)流場從一種流動形態(tài)變成另外一種流動形態(tài)時,力學(xué)參數(shù)發(fā)生顯著變化,在某些間距比區(qū)間內(nèi)出現(xiàn)驟升或驟降現(xiàn)象。

    譜單元法;順排兩圓柱;正方形四圓柱;間距比;流動形態(tài)

    0 引 言

    自從1977年Gottlieb和Orszag[1]系統(tǒng)性地從數(shù)學(xué)方面對譜方法進(jìn)行了理論闡述以來,譜方法被廣泛地應(yīng)用于更多的領(lǐng)域。1984年,Gottlieb和Hussaini開始將譜方法向計算流體動力學(xué)推廣[2]。20世紀(jì)80年代初期,Patera結(jié)合譜方法的精度和有限元的思想提出了所謂的譜單元法[3]。譜單元法是譜方法與有限元方法結(jié)合起來,將求解區(qū)域剖分成單元,采用等參元,在每個單元內(nèi)用N階多項式展開來近似變量,對方程采用Galerkin方法,用Gauss點求積分,得到離散化的方程組。當(dāng)單元數(shù)目固定時,計算精度隨多項式階數(shù)N的增加而增加。因此,譜單元法是求解偏微分方程的一種高階加權(quán)余量法,它繼承了有限元方法對復(fù)雜幾何區(qū)域的適應(yīng)性和譜方法的高精度與收斂特性等優(yōu)點。

    圓柱繞流的理論研究和工程應(yīng)用意義重大。順排兩圓柱或多圓柱是常見的排列形式,例如發(fā)電廠冷卻塔群等。圓柱尾流中的旋渦脫落對圓柱繞流的周期性升力、阻力產(chǎn)生影響,引起壓力脈動、切應(yīng)力脈動和升阻力系數(shù)的脈動。相對于單圓柱的繞流,兩圓柱或多圓柱繞流及其流致受力更為復(fù)雜?,F(xiàn)今,有限體積法、Boltzmann等多種數(shù)值模擬方法可用于分析圓柱繞流,獲取柱體表面上的作用力和關(guān)鍵參數(shù)[4-6]。

    本文較系統(tǒng)闡述譜單元方法及其表達(dá),運(yùn)用譜單元方法數(shù)值模擬低雷諾數(shù)Re=200下的順排兩圓柱、Re=150正方形排列四圓柱的繞流及其互擾效應(yīng),分析不同間距比對流場特性的影響,揭示平均升阻力系數(shù)和斯托羅哈數(shù)(旋渦脫落頻率)等隨間距比的變化,并解釋其形成機(jī)制。

    本文的數(shù)值計算采用由Monash大學(xué)的Blackburn教授所提供的Semtex程序,對二維計算,該程序采用Gauss-Lobatto-Legendre(GLL)多項式作為譜單元基函數(shù),時間和空間的數(shù)值精度可由使用人員自由控制。該方法已成功應(yīng)用于多項科學(xué)研究[7-10]。

    1 基本理論和計算公式

    二維粘性不可壓縮非定常流動的控制方程為Navier-Stokes方程:

    采用基于混合剛性穩(wěn)定格式的高階時間分裂算法[7],可將Navier-Stokes方程的求解轉(zhuǎn)化為如下三部分:

    流體壓力滿足泊松方程:

    對流體域的空間離散,采用譜單元法,且選取Legendre多項式[7]的極值點作為插值點。

    Φ(Ω)(Ω上任意連續(xù)函數(shù)構(gòu)成的空間)中的任何連續(xù)實函數(shù)u(x)可表示為:

    插值多項式

    (1)在[-1,1]以外為0;(2)在[-1,1]內(nèi),φi(xj)=δij。其中LN(x)為Legendre多項式。

    若將式(4)、(5a)、(5b)、(6a)寫成變分形式,則有:

    對計算區(qū)域Ω采取等參變換,且取φ=φj,Navier-Stokes方程的離散格式變?yōu)椋?/p>

    其中

    選取Legendre多項式,并利用Gauss積分,各系數(shù)計算公式為

    其中

    其中,[φb]表示相應(yīng)變量(壓力或速度)在計算區(qū)域邊界點上的未知值,[φi]表示相應(yīng)變量是在計算區(qū)域的內(nèi)點上的未知值。

    先由下式解出邊界上的未知物理量值:

    再由下式解出計算區(qū)域內(nèi)點的未知物理量值:

    若求解的方程組維數(shù)較高,可采用對單元(而不是對計算區(qū)域)求解的策略,即先采用類似于式(19)的方程求得各個單元邊界上的未知量值,而后再對各個單元分別求解其內(nèi)點上的未知物理量值。

    2 單圓柱繞流模擬與數(shù)值驗證

    為驗證計算的準(zhǔn)確性和可靠性,茲選用圖1示的單圓柱繞流問題進(jìn)行數(shù)值方法驗證??紤]直徑為D的圓柱受未經(jīng)擾動的均勻來流作用,基于圓柱直徑和來流流速的雷諾數(shù)取Re=200。選定計算域為50D×40D,圓柱位于坐標(biāo)系原點(0,0)。入口邊界和出口邊界分別位于圓柱中心上游20D和下游30D處,流域頂部和底部離圓柱中心20D。相應(yīng)的邊界條件為:進(jìn)口處自由來流速度為特征速度,即u=U∞=1,v=0.0;上下邊界條件與進(jìn)口邊界條件相同;出口邊界處壓力均為0.0,即p=0.0;圓柱表面處為不可滑移邊界條件,即u=0.0,v=0.0。計算域網(wǎng)格劃分采用四邊形非結(jié)構(gòu)化譜單元網(wǎng)格,網(wǎng)格單元數(shù)為354,如圖1(a)所示。在靠近圓柱壁面的地方進(jìn)行幾層細(xì)化網(wǎng)格加密,離圓柱壁面最近的一層網(wǎng)格厚度為0.1D,如圖1(b)示。同時,在圓柱尾流區(qū)域亦進(jìn)行加密處理。

    圖1 單圓柱繞流的譜單元網(wǎng)格示意Fig.1 Mesh system of spectral element for cross-flow around a single cylinder

    研究圓柱繞流流場常用的無量綱化系數(shù)包括阻力系數(shù)CD、升力系數(shù)CL和斯托羅哈數(shù)St,定義為:

    其中,F(xiàn)D為阻力,與來流方向一致,主要由流體圓柱前后壓力差和繞圓柱柱表面摩擦阻力造成;FL為升力,與來流方向垂直,主要由渦交替從圓柱上下表面脫落產(chǎn)生上下表面壓力脈動造成;fs為旋渦脫落頻率,D為圓柱直徑。

    為分析插值函數(shù)的階數(shù)對計算結(jié)果的影響,對單圓柱繞流的三種計算工況進(jìn)行基于三種不同階數(shù)的插值函數(shù)的數(shù)值模擬。計算工況1中,譜函數(shù)插值采用N=5階Gauss-Lobatto-Legendre多項式形函數(shù);計算工況2中,N=7;工況3中,N=9。計算時間步長Δi=0.005。從表1中看出,采取高階譜插值的譜單元法的計算精度和效率較高,即使對較粗糙網(wǎng)格劃分,當(dāng)N=5、7、9時,亦可得到滿意的結(jié)果,且與既有文獻(xiàn)中的實驗結(jié)果和數(shù)值模擬結(jié)果吻合。為取得計算精度和計算效率的平衡,本文取N=7。這樣,計算的空間精度為7階,時間上的精度為2階。

    需要說明的是,盡管在本單圓柱擾流的算例中Re=200有弱三維效應(yīng),但是本研究采用二維計算,仍能得到與其他文獻(xiàn)接近的結(jié)果。

    表1 平均阻力系數(shù)和斯托羅哈數(shù)(雷諾數(shù)Re=200)Table 1 Average drag coefficients and Strouhal numbers(Re=200)

    3 順排兩圓柱繞流數(shù)值模擬

    將兩圓柱順排排列放置在流體中,通過改變圓柱軸線之間的距離L,分析流場的變化。取Re=200,計算時間步長Δi=0.005,采用N=7階插值函數(shù),數(shù)值模擬間距比L/D=1.2、1.5、2.0、3.6、4.0、6.0、8.0和10.0情形下的流場。

    3.1 計算域和邊界條件

    流場計算域和邊界條件與上述單個圓柱繞流數(shù)值模擬情形相同,圖2為計算域和邊界條件示意圖。圖3(a)為兩圓柱圓心間距為2.5D時的網(wǎng)格劃分示意圖,共劃分432個網(wǎng)格,在靠近圓柱壁面的幾層進(jìn)行細(xì)化網(wǎng)格加密處理,離圓柱壁面最近的一層網(wǎng)格厚度為0.1D,如圖3(b)所示。同時,在圓柱尾流區(qū)域進(jìn)行加密處理。

    圖2 雙圓柱順排排列的計算域和邊界條件示意圖Fig.2 Schematic diagram of the computational domain and boundary conditions

    圖3 兩圓柱繞流的譜單元網(wǎng)格示意圖Fig.3 Mesh system of spectral element for cross-flow around two cylinders

    3.2 不同間距比下的流態(tài)

    圖4為不同間距比下的順排兩圓柱繞流渦量圖,從圖4中看出,當(dāng)兩圓柱串列時存在一個臨界間距,小于臨界間距時,上游圓柱后不存在旋渦脫落;大于臨界間距時,上游圓柱后發(fā)生周期性旋渦脫落。據(jù)圖4可定性判定出此臨界間距比大約為L/D≈3.6;即當(dāng)L/D<3.6時,上游圓柱后沒有產(chǎn)生渦脫落,上游圓柱繞流的分離分界層附著在下游圓柱上,而下游圓柱的后側(cè)產(chǎn)生渦脫落,尾流模式類似于單圓柱繞流尾流模式,當(dāng)L/D=3.6時,上游圓柱后開始形成渦脫落,當(dāng)L/D=4時,上游圓柱的后側(cè)已形成Kármán渦街,且旋渦脫落后撞擊下游圓柱前表面,干擾下圓柱的旋渦脫落模式,隨著間距比的繼續(xù)增大,上游圓柱后的脫落渦繼續(xù)影響下游圓柱的尾流。Re=200時,圓柱繞流發(fā)生三維效應(yīng);可以看到上游圓柱后的渦脫落穩(wěn)定,但下游圓柱后的尾流渦脫落形態(tài)變得復(fù)雜,因此二維數(shù)值模擬不能精確確定臨界間距比周圍的互擾區(qū)域。

    圖4 雙圓柱不同間距比下的渦量云圖Fig.4 Vorticity contours with different spacing ratios for the two cylinders

    3.3 圓柱表面受力特性隨間距比的變化

    圖5為兩圓柱繞流平均阻力系數(shù)隨間距比的變化。從圖5可以看到,兩圓柱的平均阻力系數(shù)均小于單個圓柱的平均阻力系數(shù),圓柱1(上游圓柱)的平均阻力系數(shù)要大于圓柱2(下游圓柱)的平均阻力系數(shù)。對于小間距比下圓柱1,先觀察到其平均阻力系數(shù)隨著間距比的增大而減小,間距比L/D=3.6的時取得最小值;然后隨著間距比增大,平均阻尼增大,在L/D=3.6和L/D=4之間劇增后,其值緩慢增大,逐漸接近于單個圓柱的平均阻力系數(shù)。對于圓柱2,間距比L/D<3.6時,其值隨著間距比增大而增大,在L/D=3.6和L/D=4之間劇增,從負(fù)值變到相當(dāng)大的正值,之后又緩慢減小,趨于穩(wěn)定。從中易見臨界間距比的存在(3.6<L/D<4.0)。當(dāng)間距比小于臨界值時,上游圓柱的分離剪切層附著在下游圓柱上,渦脫落現(xiàn)象只在下游圓柱出現(xiàn),此時下游圓柱受到較小的阻力,且為負(fù)值。當(dāng)間距比大于臨界值時,上下游圓柱均產(chǎn)生渦脫落,且升力振幅均驟然增大(圖4g),平均阻力系數(shù)突然變大,下游圓柱變化更劇烈。之后,隨著間距比繼續(xù)增大,上下游圓柱的渦脫落均呈現(xiàn)出穩(wěn)定性,變化緩慢。

    圖5 雙圓柱平均阻力系數(shù)隨間距比的變化Fig.5 Variation of mean drag coefficients with different spacing ratios for the two cylinders

    3.4 斯托羅哈數(shù)隨間距比的變化

    本文所有變量均經(jīng)過無量綱化處理,因此斯托羅哈數(shù)即等于渦脫落頻率。渦脫落頻率可通過對升力曲線作FFT分析求得。

    圖6為斯托羅哈數(shù)隨間距比的變化??傮w上,兩圓柱的斯托羅哈數(shù)基本保持相等,均隨著間距比的增大先減小再增大,并趨向于單個圓柱繞流時的斯托羅哈數(shù)。但間距比在某些區(qū)間內(nèi)對斯托羅哈數(shù)的影響顯著,即間距比從3.6增大到4.0時,斯托羅哈數(shù)劇烈增大。當(dāng)間距比小于臨界值時,上游圓柱后無渦脫落,僅下游圓柱后產(chǎn)生渦脫落;當(dāng)間距比大于臨界值時,上游圓柱和下游圓柱均存在周期性的渦脫落,從而使得斯托羅哈數(shù)在臨界間距之間發(fā)生顯著跳躍。同樣,在雙圓柱擾流的算例中Re=200有弱三維效應(yīng),二維計算,仍能得到與其他文獻(xiàn)較好擬合的結(jié)果。

    圖6 斯托羅哈數(shù)隨間距比的變化Fig.6 Variation of the Strouhal numbers,St,with different spacing ratios,L/D for the two tandem cylinders

    4 正方形順排排列四圓柱繞流數(shù)值模擬

    把四圓柱繞正方形順排排列放置在流體中,通過改變圓柱軸線之間的距離L,分析流場的變化。取Re=150,計算時間步長Δi=0.005,采用N=7階插值函數(shù),數(shù)值模擬L=1.4D、2.5D、3.5D和4.0D時的流場。在Re=150這樣的雷諾數(shù)下,流體無三維效應(yīng)。同樣,計算的空間精度為7階,時間上的精度為2階。

    4.1 計算域和邊界條件

    圖7 四圓柱排列的計算域和邊界條件示意圖Fig.7 Schematic diagram of the computational domain and boundary conditions for four cylinders

    流場計算域和邊界條件與上述單圓柱繞流數(shù)值模擬情形相同,圖7為計算域和邊界條件示意圖。圖8(a)為兩圓柱圓心間距為2.5D時的網(wǎng)格劃分示意圖,共劃分778個網(wǎng)格,在靠近圓柱壁面的地方進(jìn)行幾層網(wǎng)格加密,離圓柱壁面最近的一層網(wǎng)格厚度為0.1D,如圖8(b)所示。同時,在圓柱尾流區(qū)域進(jìn)行加密處理。

    圖8 四圓柱繞流的譜單元網(wǎng)格示意圖Fig.8 Mesh system of spectral element for cross-flow around four cylinders

    4.2 不同間距比下的流態(tài)

    圖9為不同間距比下的四圓柱繞流的渦量云圖。觀察發(fā)現(xiàn)三種不同的流動形態(tài):間距比較?。↙/D=1.4)時,圓柱互擾以臨近效應(yīng)為主,上游兩個圓柱后面靠內(nèi)側(cè)的自由剪切層附著在下游兩個圓柱表面,但是靠外側(cè)的兩個自由剪切層并未附著在下游兩個圓柱表面,而是覆蓋下游兩個圓柱,并從下游兩個圓柱后面流出,內(nèi)外兩側(cè)的自由剪切層未出現(xiàn)擺動,圓柱尾流呈現(xiàn)為單鈍體繞流流動形態(tài);在中等間距比(L/D=2.5)條件下,結(jié)構(gòu)互擾以臨近效應(yīng)和剪切層干擾為主,上游兩個圓柱后面靠內(nèi)側(cè)的自由剪切層再附著在下游兩個圓柱表面,但是靠外側(cè)的兩個自由剪切層并未再附著在下游兩個圓柱表面,而是在下游兩個圓柱表面交替出現(xiàn)擺動,此時下游圓柱尾部形成旋渦脫落流態(tài),且干擾上游圓柱尾部旋渦發(fā)展并抵制其脫落;圓柱間距比較大(L/D=3.5/4.0)時,結(jié)構(gòu)互擾以尾流效應(yīng)為主,此時所有圓柱尾流均充分發(fā)展,上游圓柱的自由剪切層形成成熟的渦脫落,然后撞擊到下游圓柱,且剪切層不斷向下游泄送Kármán旋渦,圓柱互擾表現(xiàn)為相鄰剪切層之間的相互干擾。由上可見,流動形態(tài)明顯影響圓柱表面的壓力分布,尤其是當(dāng)從一種流動形態(tài)轉(zhuǎn)變到另一種流動形態(tài),壓力分布變化顯著。

    圖9 不同間距比下的渦量云圖Fig.9 Vorticity contours with different spacing ratios for the four cylinders

    4.3 圓柱表面受力特性隨間距比的變化

    圖10為四圓柱繞流情形下的柱體平均阻力系數(shù)隨間距比的變化,同時與相關(guān)文獻(xiàn)[11-12]在不同雷諾數(shù)下的結(jié)果對比。從圖10中看出,圓柱1和2的平均阻力系數(shù)CD1和CD2基本上相等,圓柱3和4CD3和CD4也呈現(xiàn)相同的特征,這是對稱性造成的結(jié)果,但是在間距比較小時,則有細(xì)微差異。另外,比較不同雷諾數(shù)的結(jié)果,發(fā)現(xiàn)雷諾數(shù)對阻力系數(shù)的影響較大,但是數(shù)值在同一數(shù)量級范圍內(nèi)。

    圖10 四圓柱平均阻力系數(shù)隨間距比的變化Fig.10 Variation of mean drag coefficients with different spacing ratios for the four cylinders

    圖11 四圓柱平均升力系數(shù)隨間距比的變化Fig.11 Variation of mean lift coefficients with different spacing ratios for the four cylinders

    圖11為四圓柱繞流情形下的柱體平均升力系數(shù)隨間距比的變化,同時與相關(guān)文獻(xiàn)[14-15]在不同雷諾數(shù)下的結(jié)果對比。從圖中可見,圓柱1和2平均升力系數(shù)CL1和CL2相反,圓柱3和4平均升力系數(shù)CL3和CL4也相反,這反映出兩對圓柱間互相排斥。隨著間距比的增大,4個圓柱的平均升力系數(shù)的絕對值均減小,且趨于0。上游兩個圓柱的平均升力系數(shù)絕對值大于下游兩個圓柱。在層流范圍內(nèi)不同雷諾數(shù)的情況下,柱體的平均升力系數(shù)接近。但是在湍流情況下,間距越小,平均升力系數(shù)會更大。

    4.4 斯托羅哈數(shù)隨間距比的變化

    圖12為四圓柱繞流情形下的斯托羅哈數(shù)St1-St4隨間距比的變化。因本文變量均經(jīng)無量綱化處理,所以斯托羅哈數(shù)即等于旋渦脫落頻率。旋渦脫落頻率可通過對升力曲線作傅立葉變換(FFT)分析求得。從圖9中看出,四個圓柱旋渦脫落頻率近乎相同。從L/D=1.4到L/D=2.5是從一種尾流形態(tài)轉(zhuǎn)變成另一種完全不同尾流形態(tài),旋渦脫落頻率有所減小,從L/D=2.5到L/D=3.5尾流形態(tài)再發(fā)生轉(zhuǎn)變,旋渦脫落頻率增大。L/D=1.4時,尾流形態(tài)類似于單鈍體繞流時的尾流,尾流只有一個占主要地位的旋渦脫落;L/D=2.5時,上下側(cè)兩對圓柱尾流中分別形成旋渦脫落,類似于兩對反對稱的串列圓柱,但是每一對串列圓柱后面的尾流又不同于兩串列圓柱時的尾流,尾流中有一個占主要地位的旋渦,還有一個小渦伴隨,而并非兩排相同大小的旋渦;之后隨著間距比繼續(xù)增大,旋渦脫落頻率隨之增大,并逐漸趨近于單圓柱繞流時的旋渦脫落頻率,上下兩側(cè)尾流形態(tài)也類似于兩對串列圓柱尾流形態(tài)。

    圖12 四圓柱斯托羅哈數(shù)隨間距比的變化Fig.12 Variation of the Strouhal numbers with different spacing ratios for the four cylinders

    比較雷諾數(shù)的影響,發(fā)現(xiàn)從層流到湍流范圍內(nèi),旋渦脫落頻率在0.05~0.2之間,在同一數(shù)量級。

    5 結(jié) 論

    本文系統(tǒng)闡明了譜單元方法的基本原理,基于譜單元方法對低雷諾數(shù)和不同間距比下的順排雙圓柱繞流和正方形排列的四圓柱擾流及其阻力系數(shù)、升力系數(shù)等關(guān)鍵參數(shù)進(jìn)行數(shù)值模擬和比較分析。

    對順排雙圓柱擾流Re=200,間距比選擇包括1.2、1.5、2.0、3.6、4.0、6.0、8.0和10.0等八種。分析顯示,間距比對兩圓柱繞流流場特性影響明顯,當(dāng)流場從一種流動形態(tài)變成另外一種流動形態(tài)時,關(guān)鍵力學(xué)參數(shù)變化顯著,在某些間距比區(qū)間內(nèi)出現(xiàn)驟升或驟降的現(xiàn)象。順排雙圓柱繞流存在臨界間距比,在Re=200時臨界間距比約為3.6;當(dāng)間距比較小時(L/D≤3.6),上游圓柱后不發(fā)生渦脫落,而僅在下游圓柱后發(fā)生渦脫落;間距比增大時(L/D≥4.0),兩圓柱之間產(chǎn)生渦脫落,且上游圓柱后的渦脫落對下游圓柱產(chǎn)生較大的影響。在臨界間距比區(qū)間(3.6≤L/D≤4.0),平均阻力系數(shù)、平均阻力系數(shù)均方根、升力系數(shù)均方根和斯托羅哈數(shù)變化顯著。

    對于正方形排列四圓柱繞流問題在Re=150,發(fā)現(xiàn)了三種不同的流動形態(tài):間距比較?。↙/D=1.4)時,圓柱互擾以臨近效應(yīng)為主,圓柱尾流呈現(xiàn)為單鈍體繞流流動形態(tài);中等間距比(L/D=2.5)時,結(jié)構(gòu)互擾以臨近效應(yīng)和剪切層干擾為主,此時下游圓柱尾部形成旋渦脫落流態(tài),且干擾上游圓柱尾部旋渦發(fā)展并抑制其脫落;圓柱間距較大(L/D=3.5/4.0)時,結(jié)構(gòu)互擾以尾流效應(yīng)為主,此時所有圓柱尾流均充分發(fā)展,且剪切層不斷向下游泄送Kármán旋渦,圓柱互擾表現(xiàn)為相鄰剪切層之間的相互干擾。圓柱1和圓柱2的平均阻力系數(shù)在L/D=2.5處發(fā)生跳躍,四個圓柱的斯托羅哈數(shù)均在L/D=2.5處有一個彎折,之后隨著間距比的增大而增大。

    [1]GOTTLIEB D,ORGSZAG S A.Numerical analysis of spectral method:theory and application[A].CB Ms-NFS Monograph No 26.Society for Industry and Applied Mathematics[C].Philadelphia,1977.

    [2]GOTTLIEB D,HUSSAINI M Y,ORGSZAG S A.Theory and application of spectral method[A].In:Voigt RG.Gottlieb D.Hussaini M Y.Spectral Method for Partial Differential Equation[C].SIAM Philadelphia,1984:1-54.

    [3]PATERA A T.A spectral element method for fluid dynamics:laminar flow in a channel expansion[J].Journal of computattonal Phystcs,1984,154:468-488.

    [4]TAN Lingyan.Numerical simulation of flow past a main cylinder with an affiliated cylinder[J].Journal of Jtltn Untverstty(Sctence Edttton),2012,50(1):69-72.(in Chinese)

    譚玲燕.?dāng)?shù)值模擬放置附屬圓柱的主圓柱繞流[J].吉林大學(xué)學(xué)報(理學(xué)版),2012,50(1):69-72.

    [5]WANG Tong,CAO Shuyang,ZHOU Qiang.Investigation of aerodynamic lift of a circular cylinder in linear shear flows at subcritical Reynolds number[J].Chtnese Quarterly of Mechantcs,2011,32(4):473-479.(in Chinese)

    王通,曹曙陽,周強(qiáng).亞臨界雷諾數(shù)線性剪切流場中圓柱的氣動升力特性研究[J].力學(xué)季刊,2011,32(4):473-479.

    [6]GONG Shuai,GUO Zhaoli.Lattice Boltzmann simulation of flow over a transversly oscillation circular cylinder[J].Chtnese Journal of Theorettcal and Applted Mechantcs,2011,43(5):809-818.(in Chinese)

    龔帥,郭照立.橫向振蕩圓柱繞流的格子Boltzmann方法模擬[J].力學(xué)學(xué)報,2011,43(5):809-818.

    [7]BLACKBURN H M,SHERWIN S J.Formulation of a Galerkin spectral element-Fourier method for three-dimensional incompressible flows in cylindrical geometries[J].Journal of Computattonal Phystcs,2004,197:759-778.

    [8]BLACKBURN H M,LOPEZ J M.Modulated waves in a periodically driven annular cavity[J].Journal of Flutd Mechantcs,2011,667:336-357.

    [9]BLACKBURN H M,SHEARD G J.On quasi-periodic and subharmonic Floquet wake instabilities[J].Phystcs of Flutds,2010,22:031701-1-4.

    [10]BLACKBURN H M,BARKLEY D,SHERWIN S J.Convective instability and transient growth in flow over a backward-facing step[J].Journal of Flutd Mechantcs,2008,603:271-304.

    [11]FRANKE R,RODI W,SCHONUNG B.Numerical calculation of laminar vortex shedding flow past cylinders[J].Journal of Wtnd Engtneertng and Industrtal Aerodynamtcs,1990,35:237-257.

    [12]MENEGHINI J R,SALTARA F,SIQUEIRA C L R,et al.Numerical simulation of flow interference between two circular cylinders in tandem and side-by-side arrangements[J].Journal of Flutds and Structures,2001,15:327-350.

    [13]FARRANT T,TAN M,PRICE W G.A cell boundary element method applied to laminar vortex-shedding from arrays of cylinders in various arrangements[J].Journal of Flutds and Structures,2000,14:375-402.

    [14]LAM K,LI J Y,SO R M C.Force coefficient and strouhal numbers of four cylinders in cross flow[J].Journal of Flutds and Structures,2003,18:305-324.

    [15]LAM K,GONG W Q,SO R M C.Numerical simulation of cross-flow around four cylinders in an in-line square configuration[J].Journal of Flutds and Structures,2008,24:34-57.

    Numerical simulation of cross-flow around multiple circular cylinders by spectral element method

    HAN Zhaolong,ZHOU Dai,CHEN Yanan,GUI Xiaolan,LI Junlong
    (School of Naval architecture,Ocean and Civil Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

    The details of the spectral element method are presented.Flow past two tandem circular cylinders at the Reynolds numberRe=200 and four square-arranged circular cylinders with various spacing ratios atRe=150 are numerically investigated.The effects of spacing ratiosL/D(ratio of center-to-center distance to diameter of cylinder)for the two cases above are main concerned.The relationship between the spacing ratio and the flow patterns,mean drag and lift coefficients as well as the Strouhal numbers are studied.The results reveal that the spacing ratio has a great impact on the flow characteristics of both the two-cylinder andthe four-cylinder cases.There is a critical spacing ratio of around 3.6 atRe=200 for the case of two tandem cylinders.Three different flow patterns are observed in the case of four-cylinders.Among different flow patterns,a transformation from one pattern to another will change the drag/lift forces and vortex shedding frequencies.The key mechanical parameters jump suddenly in some critical spacing ratio intervals.

    spectral element method;two tandem cylinders;four square-arranged cylinders;spacing ratio;flow patterns

    O353.4

    A doi:10.7638/kqdlxxb-2012.0070

    0258-1825(2014)01-0021-10

    2012-08-07;

    2012-11-22

    國家自然科學(xué)基金(51078230,11172174,51278297)、上海市優(yōu)秀學(xué)術(shù)帶頭人計劃(13XD1402100)項目資助

    韓兆龍(1983-),男,工學(xué)博士,專業(yè):結(jié)構(gòu)工程,研究方向:風(fēng)工程與流固耦合.E-mail:han.a(chǎn)rkey@gmail.com

    周岱,男(1963-),工學(xué)博士,教授,研究方向:大跨空間結(jié)構(gòu)與流固耦合.E-mail:zhoudai@sjtu.edu.cn

    韓兆龍,周 岱,陳亞楠,等.譜單元法及其在多圓柱繞流分析中的應(yīng)用[J].空氣動力學(xué)學(xué)報,2014,32(1):21-29.

    10.7638/kqdlxxb-2012.0070.HAN Z L,ZHOU D,CHEN Y N,et al.Numerical simulation of cross-flow around multiple circular cylinders by spectral element method[J].ACTA Aerodynamica Sinica,2014,32(1):21-29.

    猜你喜歡
    斯托旋渦尾流
    小心,旋渦來啦
    大班科學(xué)活動:神秘的旋渦
    旋渦笑臉
    山間湖
    飛機(jī)尾流的散射特性與探測技術(shù)綜述
    錐形流量計尾流流場分析
    化繁為簡,費斯托工具2步法拋光工藝
    雷貝斯托:技術(shù)革新讓OE與售后兩市場相得益彰
    水面艦船風(fēng)尾流效應(yīng)減弱的模擬研究
    芬蘭坎尼斯托學(xué)校
    神州·校長(2013年7期)2013-04-29 19:10:31
    久久久精品区二区三区| 日本wwww免费看| 免费久久久久久久精品成人欧美视频| 国产1区2区3区精品| 精品少妇一区二区三区视频日本电影| 精品一区二区三卡| 国产有黄有色有爽视频| 国产精品亚洲av一区麻豆| 亚洲欧美日韩高清在线视频 | 18禁裸乳无遮挡动漫免费视频| 脱女人内裤的视频| 美女视频免费永久观看网站| 久久av网站| cao死你这个sao货| 自拍欧美九色日韩亚洲蝌蚪91| 99久久综合免费| 国产成人欧美| 日本黄色日本黄色录像| 中文字幕人妻熟女乱码| 久久99热这里只频精品6学生| 9色porny在线观看| 国产精品免费视频内射| 精品亚洲成a人片在线观看| 一区二区日韩欧美中文字幕| 19禁男女啪啪无遮挡网站| 亚洲一码二码三码区别大吗| 男女边摸边吃奶| 国产一区亚洲一区在线观看| a级毛片黄视频| videos熟女内射| 国产精品九九99| 青草久久国产| 亚洲欧洲精品一区二区精品久久久| 五月天丁香电影| av天堂在线播放| 免费看不卡的av| 纵有疾风起免费观看全集完整版| 国产成人欧美在线观看 | 中文字幕人妻丝袜制服| 中文字幕色久视频| 国产片特级美女逼逼视频| 丝袜喷水一区| 伊人亚洲综合成人网| 自拍欧美九色日韩亚洲蝌蚪91| 欧美人与性动交α欧美软件| 亚洲精品一区蜜桃| 精品国产乱码久久久久久男人| 亚洲欧美成人综合另类久久久| 久久影院123| 水蜜桃什么品种好| a级片在线免费高清观看视频| 美女大奶头黄色视频| 久久久国产一区二区| 欧美国产精品va在线观看不卡| 涩涩av久久男人的天堂| 亚洲七黄色美女视频| 一本—道久久a久久精品蜜桃钙片| 脱女人内裤的视频| 亚洲少妇的诱惑av| 国产成人精品久久久久久| 黄色片一级片一级黄色片| 亚洲 欧美一区二区三区| 免费看av在线观看网站| 90打野战视频偷拍视频| 操美女的视频在线观看| 在现免费观看毛片| 巨乳人妻的诱惑在线观看| 久久中文字幕一级| 少妇的丰满在线观看| 一本综合久久免费| 搡老岳熟女国产| 免费观看av网站的网址| √禁漫天堂资源中文www| 国产片特级美女逼逼视频| 青春草视频在线免费观看| 国产成人精品久久久久久| 色综合欧美亚洲国产小说| 国产精品久久久久成人av| 丁香六月天网| 亚洲黑人精品在线| 久久热在线av| 伊人久久大香线蕉亚洲五| www.av在线官网国产| 99香蕉大伊视频| 看免费av毛片| 亚洲国产成人一精品久久久| 香蕉丝袜av| 在线av久久热| 水蜜桃什么品种好| 成年人免费黄色播放视频| 久9热在线精品视频| xxxhd国产人妻xxx| av在线app专区| 精品国产一区二区三区久久久樱花| 亚洲av男天堂| 9热在线视频观看99| 在线观看免费视频网站a站| av网站免费在线观看视频| 一级毛片电影观看| 国产精品熟女久久久久浪| 亚洲av在线观看美女高潮| 亚洲人成77777在线视频| 日韩一卡2卡3卡4卡2021年| 热re99久久国产66热| 国产在线一区二区三区精| 大陆偷拍与自拍| 七月丁香在线播放| 国产精品一区二区精品视频观看| 午夜福利,免费看| 成人手机av| 午夜免费鲁丝| 男人舔女人的私密视频| 精品免费久久久久久久清纯 | 午夜福利乱码中文字幕| 色精品久久人妻99蜜桃| 久久影院123| 在线观看www视频免费| 欧美人与善性xxx| 日本午夜av视频| 热re99久久国产66热| 十八禁高潮呻吟视频| 香蕉丝袜av| 青青草视频在线视频观看| 中文乱码字字幕精品一区二区三区| 国产淫语在线视频| 亚洲精品久久久久久婷婷小说| 久久久久精品国产欧美久久久 | 晚上一个人看的免费电影| 丝袜喷水一区| 19禁男女啪啪无遮挡网站| 国产黄色免费在线视频| 亚洲天堂av无毛| 精品一区在线观看国产| 精品国产国语对白av| 亚洲国产av新网站| 久久人妻福利社区极品人妻图片 | 黄色怎么调成土黄色| 9热在线视频观看99| 丝瓜视频免费看黄片| 欧美黄色片欧美黄色片| 嫁个100分男人电影在线观看 | 另类亚洲欧美激情| 黄频高清免费视频| 成人国产av品久久久| 欧美变态另类bdsm刘玥| 日韩av不卡免费在线播放| 久久久久视频综合| 在线 av 中文字幕| 在线观看人妻少妇| 国产精品一区二区在线观看99| 久久这里只有精品19| 午夜免费鲁丝| 999久久久国产精品视频| 成年人免费黄色播放视频| 一区在线观看完整版| 激情五月婷婷亚洲| 下体分泌物呈黄色| 国产野战对白在线观看| 国产精品久久久人人做人人爽| 美女脱内裤让男人舔精品视频| 中文字幕人妻熟女乱码| 国产高清国产精品国产三级| 秋霞在线观看毛片| 青春草视频在线免费观看| 大型av网站在线播放| 各种免费的搞黄视频| 99国产精品99久久久久| 麻豆乱淫一区二区| 久久国产精品男人的天堂亚洲| 少妇裸体淫交视频免费看高清 | 多毛熟女@视频| 精品卡一卡二卡四卡免费| 晚上一个人看的免费电影| 蜜桃在线观看..| 丝袜美腿诱惑在线| 男女边摸边吃奶| 91麻豆av在线| 亚洲中文日韩欧美视频| 悠悠久久av| 久久久欧美国产精品| 欧美xxⅹ黑人| av又黄又爽大尺度在线免费看| 亚洲欧美清纯卡通| 18禁观看日本| 亚洲久久久国产精品| 亚洲精品日韩在线中文字幕| 午夜激情久久久久久久| 女性被躁到高潮视频| 国产日韩欧美亚洲二区| 看十八女毛片水多多多| 久9热在线精品视频| 丁香六月欧美| 亚洲精品一卡2卡三卡4卡5卡 | 欧美黄色片欧美黄色片| 汤姆久久久久久久影院中文字幕| 国产精品二区激情视频| 成在线人永久免费视频| 国产精品一国产av| av一本久久久久| 中文乱码字字幕精品一区二区三区| 欧美av亚洲av综合av国产av| 国产精品99久久99久久久不卡| 国产熟女午夜一区二区三区| 又粗又硬又长又爽又黄的视频| 夜夜骑夜夜射夜夜干| 国产野战对白在线观看| 伊人久久大香线蕉亚洲五| 99精国产麻豆久久婷婷| 日韩制服骚丝袜av| 2021少妇久久久久久久久久久| 久久亚洲精品不卡| 美女脱内裤让男人舔精品视频| 午夜老司机福利片| 久久 成人 亚洲| 亚洲图色成人| 黄色一级大片看看| a级毛片在线看网站| 黑人猛操日本美女一级片| av国产久精品久网站免费入址| 热re99久久精品国产66热6| 亚洲欧美色中文字幕在线| 免费看不卡的av| 国产精品偷伦视频观看了| 9色porny在线观看| 亚洲av综合色区一区| 国产精品二区激情视频| 亚洲人成电影免费在线| 丁香六月欧美| 日韩制服丝袜自拍偷拍| 成人免费观看视频高清| 久久精品久久久久久噜噜老黄| 精品少妇内射三级| 日韩制服丝袜自拍偷拍| 精品国产一区二区三区久久久樱花| 亚洲免费av在线视频| av国产久精品久网站免费入址| av在线老鸭窝| 亚洲中文av在线| 免费久久久久久久精品成人欧美视频| 男人舔女人的私密视频| 欧美日韩视频精品一区| 99热全是精品| 美女大奶头黄色视频| 精品卡一卡二卡四卡免费| 欧美老熟妇乱子伦牲交| av网站在线播放免费| 午夜精品国产一区二区电影| 欧美日韩成人在线一区二区| 亚洲欧美精品综合一区二区三区| 大陆偷拍与自拍| 黄色a级毛片大全视频| 久久狼人影院| 日本vs欧美在线观看视频| 国产亚洲av高清不卡| 你懂的网址亚洲精品在线观看| 久久人人爽av亚洲精品天堂| 在线观看免费视频网站a站| 国产一区二区在线观看av| 韩国高清视频一区二区三区| 色精品久久人妻99蜜桃| 亚洲国产精品成人久久小说| 国产有黄有色有爽视频| 大片免费播放器 马上看| 赤兔流量卡办理| 亚洲精品乱久久久久久| 国产日韩欧美亚洲二区| 国产成人一区二区三区免费视频网站 | svipshipincom国产片| 久久精品久久久久久久性| 国产女主播在线喷水免费视频网站| 91精品国产国语对白视频| 欧美黑人精品巨大| av又黄又爽大尺度在线免费看| 自线自在国产av| 亚洲五月婷婷丁香| 久久人人爽人人片av| 久久精品国产亚洲av高清一级| 亚洲,一卡二卡三卡| 91麻豆精品激情在线观看国产 | av不卡在线播放| 91精品三级在线观看| 国产精品免费大片| 亚洲国产毛片av蜜桃av| 国产一区二区三区av在线| 国产精品99久久99久久久不卡| 久久人妻福利社区极品人妻图片 | 男女免费视频国产| 欧美乱码精品一区二区三区| 国产亚洲欧美在线一区二区| 日本av手机在线免费观看| 午夜免费观看性视频| 性少妇av在线| www.自偷自拍.com| 人人妻人人澡人人看| 另类精品久久| 视频在线观看一区二区三区| 亚洲欧美中文字幕日韩二区| 欧美xxⅹ黑人| 不卡av一区二区三区| 成人亚洲精品一区在线观看| 99国产精品99久久久久| 成人免费观看视频高清| 美女福利国产在线| 日韩熟女老妇一区二区性免费视频| 中国国产av一级| 丝袜喷水一区| 少妇的丰满在线观看| 精品亚洲乱码少妇综合久久| 国产国语露脸激情在线看| 国产精品久久久久久精品古装| 色精品久久人妻99蜜桃| 在线观看66精品国产| 51午夜福利影视在线观看| 精品久久久久久久人妻蜜臀av| 88av欧美| 欧美一级a爱片免费观看看 | 欧美成狂野欧美在线观看| 亚洲欧美精品综合一区二区三区| a级毛片在线看网站| 欧美日韩瑟瑟在线播放| 18禁美女被吸乳视频| 欧美乱码精品一区二区三区| 99国产综合亚洲精品| 亚洲av中文字字幕乱码综合 | 一本大道久久a久久精品| a级毛片a级免费在线| 人妻丰满熟妇av一区二区三区| 久久中文看片网| 亚洲男人的天堂狠狠| 亚洲第一av免费看| 欧美日韩亚洲国产一区二区在线观看| 国产不卡一卡二| 成年免费大片在线观看| 亚洲真实伦在线观看| 热99re8久久精品国产| 老熟妇仑乱视频hdxx| tocl精华| 日韩欧美 国产精品| 一本一本综合久久| 亚洲国产高清在线一区二区三 | 又黄又粗又硬又大视频| 日韩欧美免费精品| 欧美国产精品va在线观看不卡| 91av网站免费观看| 一本大道久久a久久精品| 俄罗斯特黄特色一大片| 搞女人的毛片| 满18在线观看网站| 岛国在线观看网站| av超薄肉色丝袜交足视频| 色婷婷久久久亚洲欧美| 免费搜索国产男女视频| 一级黄色大片毛片| 国产精品 国内视频| 欧美日韩瑟瑟在线播放| 级片在线观看| 美女午夜性视频免费| aaaaa片日本免费| 男人舔女人的私密视频| 精品无人区乱码1区二区| www.自偷自拍.com| 国产亚洲欧美98| 18禁观看日本| 欧美色欧美亚洲另类二区| 婷婷精品国产亚洲av在线| 国产日本99.免费观看| 人人妻,人人澡人人爽秒播| 日本一区二区免费在线视频| 成人国产一区最新在线观看| 久久香蕉激情| 丝袜在线中文字幕| 99久久国产精品久久久| 最好的美女福利视频网| 欧美一级毛片孕妇| 男女下面进入的视频免费午夜 | 色在线成人网| 少妇的丰满在线观看| 国产在线观看jvid| 亚洲国产精品久久男人天堂| 日韩精品免费视频一区二区三区| 法律面前人人平等表现在哪些方面| 中文字幕人成人乱码亚洲影| 久久伊人香网站| 在线永久观看黄色视频| 精品久久久久久久久久免费视频| 国产一区二区在线av高清观看| 精品国内亚洲2022精品成人| 亚洲中文av在线| 色综合婷婷激情| 成在线人永久免费视频| 国产精品乱码一区二三区的特点| 国产精品久久久久久人妻精品电影| 女性被躁到高潮视频| 国产精品久久久久久人妻精品电影| 欧美不卡视频在线免费观看 | 亚洲国产精品久久男人天堂| 在线免费观看的www视频| 欧美一级a爱片免费观看看 | 国产日本99.免费观看| 2021天堂中文幕一二区在线观 | 欧美激情 高清一区二区三区| 国产精品一区二区免费欧美| 天堂动漫精品| 成人一区二区视频在线观看| 久久久久久大精品| 91九色精品人成在线观看| 精品久久久久久久毛片微露脸| av电影中文网址| 91av网站免费观看| 国产精品一区二区精品视频观看| 国产亚洲av嫩草精品影院| 国产成人欧美| 很黄的视频免费| 欧美av亚洲av综合av国产av| 欧美日韩亚洲国产一区二区在线观看| 久久精品国产清高在天天线| 黄色a级毛片大全视频| 久久香蕉国产精品| 黑人操中国人逼视频| 国产亚洲精品久久久久5区| 好看av亚洲va欧美ⅴa在| 亚洲七黄色美女视频| 婷婷精品国产亚洲av| 亚洲精品久久国产高清桃花| 成人av一区二区三区在线看| 桃色一区二区三区在线观看| av欧美777| 在线播放国产精品三级| 欧美午夜高清在线| 日日夜夜操网爽| 亚洲无线在线观看| 中文字幕精品免费在线观看视频| 中亚洲国语对白在线视频| 亚洲 欧美 日韩 在线 免费| 三级毛片av免费| 波多野结衣av一区二区av| 国产精品久久电影中文字幕| 亚洲成av片中文字幕在线观看| 亚洲av中文字字幕乱码综合 | 亚洲人成网站在线播放欧美日韩| 久99久视频精品免费| 成人国产一区最新在线观看| 亚洲第一欧美日韩一区二区三区| 国产高清激情床上av| xxx96com| 亚洲人成77777在线视频| 欧美性长视频在线观看| 免费看十八禁软件| 日韩免费av在线播放| aaaaa片日本免费| 99re在线观看精品视频| 中文字幕最新亚洲高清| 女性生殖器流出的白浆| 岛国在线观看网站| 两性夫妻黄色片| 好男人电影高清在线观看| 在线av久久热| 91麻豆精品激情在线观看国产| 一级毛片女人18水好多| 嫩草影院精品99| 不卡av一区二区三区| 天天一区二区日本电影三级| 国产伦人伦偷精品视频| 亚洲欧美精品综合一区二区三区| 国产真实乱freesex| 午夜影院日韩av| 精品久久久久久成人av| 88av欧美| 亚洲一区二区三区不卡视频| 一边摸一边抽搐一进一小说| 99久久久亚洲精品蜜臀av| 亚洲免费av在线视频| 一区福利在线观看| 脱女人内裤的视频| 日本免费a在线| 侵犯人妻中文字幕一二三四区| 成年人黄色毛片网站| 人人妻,人人澡人人爽秒播| 欧美一级毛片孕妇| 日本精品一区二区三区蜜桃| 美国免费a级毛片| 色精品久久人妻99蜜桃| 日韩免费av在线播放| 久久久久国产一级毛片高清牌| 久久婷婷人人爽人人干人人爱| 怎么达到女性高潮| 国产伦人伦偷精品视频| 精品国产亚洲在线| 国产亚洲欧美98| 国产又色又爽无遮挡免费看| 观看免费一级毛片| 视频区欧美日本亚洲| 中文字幕精品免费在线观看视频| 黄片大片在线免费观看| 此物有八面人人有两片| 久久伊人香网站| 国产一区二区在线av高清观看| 亚洲自偷自拍图片 自拍| 1024视频免费在线观看| 国内少妇人妻偷人精品xxx网站 | aaaaa片日本免费| 久热这里只有精品99| 夜夜躁狠狠躁天天躁| 国内毛片毛片毛片毛片毛片| 色老头精品视频在线观看| 真人一进一出gif抽搐免费| 日本免费一区二区三区高清不卡| 亚洲av电影不卡..在线观看| 亚洲成av片中文字幕在线观看| 好看av亚洲va欧美ⅴa在| 亚洲精品久久国产高清桃花| 久久久久亚洲av毛片大全| 可以免费在线观看a视频的电影网站| 国产又色又爽无遮挡免费看| 精品国内亚洲2022精品成人| 成人午夜高清在线视频 | 国产国语露脸激情在线看| 午夜福利欧美成人| 夜夜躁狠狠躁天天躁| 欧美日韩乱码在线| 999久久久精品免费观看国产| 最新在线观看一区二区三区| 久久中文看片网| 搡老熟女国产l中国老女人| 这个男人来自地球电影免费观看| 婷婷精品国产亚洲av在线| 国产av一区在线观看免费| 亚洲无线在线观看| 在线观看66精品国产| 久久久精品国产亚洲av高清涩受| 久久亚洲精品不卡| 一区二区三区国产精品乱码| 亚洲五月天丁香| 好男人电影高清在线观看| svipshipincom国产片| 成人亚洲精品一区在线观看| 侵犯人妻中文字幕一二三四区| 很黄的视频免费| 久久天堂一区二区三区四区| 一级毛片高清免费大全| 精品久久久久久久末码| 成人特级黄色片久久久久久久| www.www免费av| 亚洲精品中文字幕在线视频| 中文字幕人妻熟女乱码| 一区二区三区激情视频| 最近最新中文字幕大全免费视频| 亚洲一区二区三区不卡视频| 欧美成人性av电影在线观看| 欧美日韩精品网址| 极品教师在线免费播放| 一级毛片精品| 99国产综合亚洲精品| 成人国产综合亚洲| 日韩大码丰满熟妇| 1024视频免费在线观看| 在线观看免费视频日本深夜| 好男人在线观看高清免费视频 | 成人手机av| 欧美+亚洲+日韩+国产| 国产单亲对白刺激| 国产av一区二区精品久久| 亚洲天堂国产精品一区在线| 久久99热这里只有精品18| 伦理电影免费视频| 天堂动漫精品| 国产激情欧美一区二区| 国产一区二区激情短视频| netflix在线观看网站| 少妇裸体淫交视频免费看高清 | 人妻丰满熟妇av一区二区三区| 女性生殖器流出的白浆| 好看av亚洲va欧美ⅴa在| 制服丝袜大香蕉在线| 亚洲精品一卡2卡三卡4卡5卡| 禁无遮挡网站| 一进一出抽搐gif免费好疼| 久久久精品国产亚洲av高清涩受| 成人三级黄色视频| 高清在线国产一区| 免费高清视频大片| 国产亚洲精品av在线| 大型av网站在线播放| 亚洲中文字幕日韩| 国产成人av教育| 亚洲成人精品中文字幕电影| 欧美中文综合在线视频| 久久精品影院6| 手机成人av网站| 精品一区二区三区四区五区乱码| 国产一卡二卡三卡精品| 欧美性长视频在线观看| 女警被强在线播放| 国产免费av片在线观看野外av| 国产伦在线观看视频一区| 亚洲一卡2卡3卡4卡5卡精品中文| 99re在线观看精品视频| 色综合婷婷激情| 亚洲成av人片免费观看| av在线播放免费不卡| 99国产精品99久久久久| 大香蕉久久成人网| 国产精品香港三级国产av潘金莲| 精品卡一卡二卡四卡免费| 麻豆av在线久日| 国产成人系列免费观看| 国产激情欧美一区二区| 99riav亚洲国产免费| 美女午夜性视频免费| 亚洲国产欧美网| 亚洲色图 男人天堂 中文字幕| 一区二区日韩欧美中文字幕| 国产伦一二天堂av在线观看| 好男人在线观看高清免费视频 | 国产激情偷乱视频一区二区|