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

    基于RANS/LES混合方法的箱梁氣動(dòng)導(dǎo)納函數(shù)研究

    2023-10-10 07:15:54張偉峰張志田
    振動(dòng)與沖擊 2023年18期
    關(guān)鍵詞:湍流脈動(dòng)壁面

    張偉峰, 張志田

    (1. 華北水利水電大學(xué) 土木與交通學(xué)院,鄭州 450045;2. 海南大學(xué) 土木建筑工程學(xué)院,海口 570228)

    大跨度橋梁的主梁斷面形式多樣,主要有整體式鋼箱梁、分離式鋼箱梁、鋼桁梁、混凝土箱梁、混凝土雙主梁等。由于斷面形式的不同,氣流流經(jīng)斷面時(shí)占據(jù)主導(dǎo)地位的流動(dòng)現(xiàn)象也不同 。例如對(duì)于扁平鋼箱梁,斷面周圍的流動(dòng)大部分為附著流,僅在轉(zhuǎn)角和尾部有較弱的分離流。相反,對(duì)于分離式雙箱梁、混凝土雙主梁等會(huì)在斷面周圍和尾部產(chǎn)生顯著的渦脫。CFD(computational fluid dynamics)模擬橋梁氣動(dòng)力的關(guān)鍵是對(duì)這些流動(dòng)現(xiàn)象進(jìn)行準(zhǔn)確的模擬。

    RANS(Reynolds-averaged navier-stokes)由于計(jì)算花費(fèi)少,在當(dāng)前的工程中得到廣泛應(yīng)用。但是RANS由于對(duì)湍流黏度過大的?;?極大的抑制了剪切層失穩(wěn)和小尺度運(yùn)動(dòng)的生成與演化,僅能預(yù)測(cè)到大尺度和低頻運(yùn)動(dòng),這對(duì)于時(shí)均氣動(dòng)力的求解是可行的,卻難以應(yīng)用于非定常氣動(dòng)力的求解,如抖振力。LES(large eddy simulation)是一種真正的非定常模擬方法,它對(duì)大尺度的脈動(dòng)進(jìn)行直接求解,而只對(duì)具有普適性的小尺度脈動(dòng)進(jìn)行模化。在近壁面由于湍流尺度Lt正比于到壁面的距離y,因此越靠近壁面,積分尺度越小,需要的空間、時(shí)間分辨率就越高,現(xiàn)階段的計(jì)算能力很難滿足高雷諾數(shù)壁面流的需求。Spalart[1]指出,對(duì)于一個(gè)雷諾數(shù)1×106左右的簡(jiǎn)單邊界層流動(dòng),網(wǎng)格數(shù)目和時(shí)間步數(shù)也分別達(dá)到了驚人的1×1011和1×107,Spalart估計(jì)在2045年有望實(shí)現(xiàn)LES在工程問題中的廣泛應(yīng)用。需要注意的是,雖然LES已經(jīng)被大量應(yīng)用于橋梁風(fēng)工程中,但已有研究的網(wǎng)格精度并沒有達(dá)到LES的標(biāo)準(zhǔn)而形成VLES(very large eddy simulation)模擬,從而不能準(zhǔn)確模擬含能區(qū)的湍流結(jié)構(gòu)與湍流能量,導(dǎo)致模擬結(jié)果與網(wǎng)格尺寸相關(guān)[2-4]。目前,嚴(yán)格的LES主要用于自由剪切層和低雷諾數(shù)的壁面流。在可見的未來,同時(shí)滿足網(wǎng)格分辨率和時(shí)間分辨率的LES在橋梁風(fēng)工程中的應(yīng)用仍很難普及[5-6]。

    RANS/LES混合方法屬于尺度解析模擬(scale-resolving simulation,SRS),在近壁面的平衡湍流區(qū)域中采用RANS,而在遠(yuǎn)離壁面的非平衡區(qū)域采用LES。這樣既降低了LES在近壁面處的計(jì)算代價(jià),又獲得了相對(duì)準(zhǔn)確且豐富的非定常信息,有效實(shí)現(xiàn)了計(jì)算效率和計(jì)算精度的統(tǒng)一。事實(shí)上,RANS/LES混合方法正是基于在當(dāng)前的計(jì)算條件下,滿足用適中的計(jì)算量提供工程應(yīng)用所需要的非定常信息而提出的[7]。三維的RANS/LES混合方法在航空航天、葉輪機(jī)械、汽車等領(lǐng)域已取得了較為廣泛的應(yīng)用,但在橋梁風(fēng)工程中卻應(yīng)用較少。Bai等[8]采用當(dāng)三維DES(detached eddy simulation)模擬了箱梁斷面的靜三分力系數(shù)并識(shí)別得到了顫振導(dǎo)數(shù),研究表明當(dāng)三維DES用于橋梁氣動(dòng)力分析時(shí)具有明顯的優(yōu)勢(shì)和較高的精度;魏志剛[9]采用IDDES(improved delayed DES)對(duì)箱梁斷面的靜力三分力系數(shù)、表面壓強(qiáng)分布、流場(chǎng)信息進(jìn)行了詳細(xì)的研究;Zhang等[10]采用DDES(delayed DES)研究了展向布置的吸氣孔對(duì)橋梁渦脫的影響。

    氣動(dòng)導(dǎo)納函數(shù)是大跨度橋梁抖振分析的關(guān)鍵參數(shù),用來描述作用在結(jié)構(gòu)物上氣動(dòng)力的非定常特性。RANS/LES混合方法由于具備對(duì)非定常流動(dòng)信息的捕捉能力,因此是研究氣動(dòng)導(dǎo)納函數(shù)非常有利的工具。遺憾的是,相比于RANS/LES混合方法在橋梁時(shí)均流場(chǎng)的研究,在非定常氣動(dòng)力的研究方面更加稀少。大量研究表明橋梁氣動(dòng)導(dǎo)納函數(shù)與來流的風(fēng)場(chǎng)特性,如湍流度、湍流積分尺度等有關(guān)[11-12]。實(shí)橋順風(fēng)向的湍流積分尺度通常位于80~300 m,約為橋梁寬度W的3倍~10倍,豎風(fēng)向的湍流積分尺度通常位于15~50 m,約為橋梁寬度W的0.5倍~2.0倍。傳統(tǒng)的基于格柵湍流場(chǎng)的風(fēng)洞試驗(yàn),對(duì)于積分尺度的模擬嚴(yán)重失真,特別是順風(fēng)向積分尺度,多數(shù)試驗(yàn)位于0.1W~0.5W[13-14]??紤]到積分尺度對(duì)氣動(dòng)導(dǎo)納函數(shù)的影響,把風(fēng)洞試驗(yàn)識(shí)別的氣動(dòng)導(dǎo)納用于實(shí)橋的抖振分析,必然會(huì)對(duì)結(jié)果帶來影響。雖然可以借助于大尺寸風(fēng)洞或者主動(dòng)風(fēng)洞實(shí)現(xiàn)大積分尺度湍流的模擬,但這兩者都需要較大的投入。

    滿足自然風(fēng)場(chǎng)特性的湍流數(shù)值模擬,一直是風(fēng)工程領(lǐng)域的熱點(diǎn)。Huang等[15-20]先后分別提出了不同的各向異性湍流場(chǎng)模擬方法,可以實(shí)現(xiàn)任意特性脈動(dòng)風(fēng)的模擬。通過將脈動(dòng)風(fēng)場(chǎng)模擬方法和RANS/LES混合方法相結(jié)合,可以實(shí)現(xiàn)任意積分尺度湍流場(chǎng)下橋梁非定常氣動(dòng)力的模擬。

    本文選取3種典型的RANS/LES混合方法,即:SAS(scale-adaptive simulation)、SBES(stress-blended eddy simulation)、WMLES(wall-modeled LES)。在前人研究的基礎(chǔ)上提出了一種改進(jìn)的脈動(dòng)風(fēng)合成方法,該方法既滿足無散度條件,同時(shí)又可以實(shí)現(xiàn)任意空間相關(guān)性和時(shí)間相關(guān)性的風(fēng)場(chǎng)。在均勻來流下,選取流線型箱梁斷面和分離式雙箱梁斷面,通過表面均勻壓強(qiáng)分布和周圍流動(dòng)特征的比較,探討了3種混合方法的特點(diǎn)和適用性。分別在與實(shí)際風(fēng)環(huán)境相似的大湍流積分尺度來流下,和風(fēng)洞試驗(yàn)格柵湍流場(chǎng)相似的小尺度來流下識(shí)別了箱梁斷面的氣動(dòng)導(dǎo)納函數(shù),研究了大積分尺度湍流對(duì)箱梁氣動(dòng)導(dǎo)納的影響。

    1 RANS/LES混合方法

    自1997年Spalart等提出DES以來,近20年來研究者們提出并發(fā)展了多種RANS/LES混合模型。不幸的是,并不存在一種普適的可以應(yīng)用于任意流動(dòng)的RANS/LES混合模型。針對(duì)于橋梁風(fēng)工程,由于各類橋梁斷面的形式不同,占據(jù)主導(dǎo)地位的流動(dòng)現(xiàn)象也不同,各種混合模型的表現(xiàn)不一。

    1.1 SAS方法

    SAS方法嚴(yán)格意義上來說不屬于RANS/LES混合模型,但考慮到它在定常流動(dòng)區(qū)域表現(xiàn)出RANS的特性,而在脈動(dòng)顯著的區(qū)域表現(xiàn)出LES的特性,因此通常把它歸類到RANS/LES混合模型。SAS模型最早由Menter等[21]提出,并由Menter和Egorov[22]改進(jìn)。它與RANS模型最主要的區(qū)別是引入馮卡門長(zhǎng)度尺度LvK。LvK包含了速度的一階導(dǎo)數(shù)和二階導(dǎo)數(shù),它只與當(dāng)?shù)亓鲃?dòng)有關(guān)而與網(wǎng)格尺度無關(guān)。當(dāng)流場(chǎng)中存在脈動(dòng)時(shí)LvK減小,而LvK的減小會(huì)使得湍流黏度變小,從而使小的脈動(dòng)得以保留和發(fā)展。

    SAS通過在k-ω模型的ω方程中引入源項(xiàng)QSAS

    QSAS=max(T1-T2,0)

    (1)

    (2)

    (3)

    式中:T1,T2為組成源項(xiàng)的附加項(xiàng);ρ為流體的密度;k為湍動(dòng)能;ω為比耗散率;xj為空間坐標(biāo)(j=1,2,3);S為應(yīng)變率張量Sij的模;η2,σφ,C,κ為常系數(shù),分別為η2=3.51,σφ=2/3,C=2,κ=0.41。

    L為湍流長(zhǎng)度尺度

    (4)

    式中,cμ為常數(shù)0.09。LvK為馮卡門長(zhǎng)度尺度

    (5)

    式中,U′,U″分別為速度的一階、二階導(dǎo)數(shù)

    (6)

    式中:Ui為流動(dòng)速度;xk為空間坐標(biāo)(k=1,2,3)。式(1)中:T1項(xiàng)通過LvK來分辨小尺度的渦;T2項(xiàng)是為了保障在定常流時(shí)模型能回歸到RANS。從以上分析可以看出,SAS模型通過LvK的變化,使得非穩(wěn)態(tài)流域的渦發(fā)展,比如鈍體斷面周圍以及尾部流動(dòng)區(qū)域。因此在均勻來流時(shí),SAS可以應(yīng)用于具有顯著流動(dòng)不穩(wěn)定現(xiàn)象的鈍體斷面,但是對(duì)于流動(dòng)不穩(wěn)定現(xiàn)象較弱的流線型斷面,SAS有可能會(huì)退化為RANS。

    1.2 SBES方法

    Spalart早期提出的DES模型,對(duì)RANS和LES計(jì)算域的判定嚴(yán)重依賴于當(dāng)?shù)鼐W(wǎng)格尺度,即對(duì)網(wǎng)格有強(qiáng)烈的依賴性。當(dāng)進(jìn)行網(wǎng)格的加密操作時(shí),有可能在邊界層內(nèi)激活LES,此時(shí)網(wǎng)格尺度過于粗糙,無法準(zhǔn)確模擬求解渦的脈動(dòng),導(dǎo)致雷諾應(yīng)力耗損,從而誘發(fā)邊界層的提前分離,也就是所謂的網(wǎng)格誘導(dǎo)分離(grid induced separation,GIS) 問題。為了克服GIS問題,Spalart等[23]通過引入延遲函數(shù)來保障在網(wǎng)格加密時(shí),邊界層仍然處于RANS,這就是所謂的DDES。

    DES類模型存在的另外一個(gè)問題是,在RANS和LES的交接處附近存在“灰區(qū)”。在這個(gè)過渡區(qū)由于邊界層內(nèi)部湍流黏度的影響,使得LES區(qū)域的小尺度脈動(dòng)得到抑制,從而解析的雷諾應(yīng)力也較小。在“灰區(qū)”以外,湍流黏度的影響減小,小尺度脈動(dòng)得以發(fā)展,從而解析的雷諾應(yīng)力增大。由于“灰區(qū)”的存在,導(dǎo)致了邊界層內(nèi)部和外部的對(duì)數(shù)率不匹配 (log layer mismatch,LLM) 問題。LLM通常會(huì)導(dǎo)致錯(cuò)誤的壁面摩擦因數(shù)。為此,Shur等[24]提出了IDDES方法來消除LLM的影響。

    近年來,Menter[25]提出了SDES(shielded detached eddy simulation)方法和SBES方法,兩種方法都基于一種新的屏蔽函數(shù),該屏蔽函數(shù)可以有效解決RANS/LES混合方法普遍存在的GIS問題和LLM問題。Menter建議使用SBES方法代替現(xiàn)有的DES、DDES、IDDES方法。

    SBES方法是把RANS和LES通過屏蔽函數(shù)結(jié)合起來,以總應(yīng)力的形式表示為

    τij,SBES=fSBESτij,RANS+(1-fSBES)τij,LES

    (7)

    式中:τij,RANS為總應(yīng)力;fSBES為屏蔽函數(shù),在邊界層內(nèi)fSBES=1,在流動(dòng)分離區(qū)和自由剪切流中fSBES=0;τij,RANS為總應(yīng)力中對(duì)應(yīng)RANS的部分;τij,LES為總應(yīng)力中對(duì)應(yīng)LES的部分。

    Menter等[26]將其應(yīng)用于槽道流、圓柱繞流、汽車?yán)@流、擴(kuò)散器等,認(rèn)為SBES方法在附著流、大尺度分離流、中等尺度分離流等流動(dòng)問題中都有著優(yōu)異的表現(xiàn)。總的來說它具有如下的特點(diǎn):現(xiàn)有的RANS類湍流模型和LES類湍流模型可以任意組合到SBES方法中;在網(wǎng)格細(xì)化時(shí)對(duì)邊界層內(nèi)的RANS提供保護(hù),防止GIS的發(fā)生;在RANS和LES交界處,可以實(shí)現(xiàn)快速過渡,減少“灰區(qū)”,消除LLM;既可以應(yīng)用于顯著流動(dòng)分離的鈍體斷面也可以應(yīng)用于附著流和小尺度分離流的流線型斷面;計(jì)算花費(fèi)小等。

    1.3 WMLES方法

    為了避免LES的計(jì)算花費(fèi)隨雷諾數(shù)Re的指數(shù)增長(zhǎng),可以在靠近壁面的對(duì)數(shù)率層使用RANS湍流模型,而在網(wǎng)格分辨率足夠時(shí)的外層使用LES湍流模型,該方法即為WMLES方法。在平行于壁面的方向,相比于完全求解的LES,WMLES方法的網(wǎng)格尺寸可以顯著增大。

    本文采用的WMLES最早由Shur等提出。它將普朗特混合長(zhǎng)度和Smagorinsky模型結(jié)合起來,對(duì)于湍流黏度系數(shù)νt有

    νt=min[(κdw)2,(CSmagΔ)2]·S·
    {1-exp[-(y+/25)3]}

    (8)

    式中:dw為網(wǎng)格中心到壁面的距離;κ,CSmag為常數(shù),κ=0.418 7,CSmag=0.2;y+為無量綱距離; Δ為亞格子濾波尺度,它同時(shí)考慮了網(wǎng)格的各向異性和到壁面的距離。

    Δ=min(max[Cw·dw,Cw·hmax,hwn];hmax)

    (9)

    式中:hmax對(duì)于六面體網(wǎng)格為網(wǎng)格最大邊長(zhǎng);hwn為垂直壁面方向的網(wǎng)格尺度;Cw為常數(shù),Cw=0.15。從構(gòu)造上來說,WMLES在近壁面采用0方程的普朗特混合長(zhǎng)度、遠(yuǎn)離壁面區(qū)域采用Smagorinsky渦黏模式,具有構(gòu)造簡(jiǎn)單、容易收斂的特點(diǎn),是最簡(jiǎn)單的RANS/LES混合方法。同時(shí),WMLES在計(jì)算中的表現(xiàn)也具有這兩種模型的特點(diǎn),如對(duì)逆壓梯度、分離再附流動(dòng)等現(xiàn)象無法準(zhǔn)確預(yù)測(cè),總體表現(xiàn)不如兩方程的RANS/LES混合方法。

    2 CFD數(shù)值模型

    在以往研究的基礎(chǔ)上,首先提出了一種入口湍流合成方法,該方法能夠模擬任意湍流特性的風(fēng)場(chǎng)。

    2.1 入口湍流模擬

    為了模擬橋梁斷面在湍流來流下的氣動(dòng)力,計(jì)算域的入口需要設(shè)置湍流邊界條件。滿足自然風(fēng)場(chǎng)特性的湍流數(shù)值模擬,一直是風(fēng)工程領(lǐng)域的熱點(diǎn)。基于Huang等提出的DSRFG(discrete and synthesizing random flow generation)方法,Castro等引入了一個(gè)時(shí)間尺度,用來調(diào)整生成風(fēng)速時(shí)程的時(shí)間相關(guān)性,并將此方法命名為MDSRFG(modified DSRFG)。為了克服DSRFG生成風(fēng)速時(shí)程的相干性與頻率無關(guān)的缺點(diǎn),Aboshosha等引入了一個(gè)與頻率有關(guān)的積分尺度,但是由于該參數(shù)的引入,使得CDRFG(consistent DRFG)方法不再滿足無散度條件。綜合MDSRFG方法和CDRFG方法,本文同時(shí)引入時(shí)間尺度因子和與頻率有關(guān)的積分尺度,并通過修正使模擬的風(fēng)速滿足無散度條件,具體過程如下。

    在CDRFG方法的基礎(chǔ)上,同時(shí)借鑒MDSRFG方法,在脈動(dòng)風(fēng)速時(shí)程的表達(dá)式中引入時(shí)間尺度,此時(shí)脈動(dòng)風(fēng)速時(shí)程ui(i=1,2,3)可以表示為

    (10)

    其中

    (11)

    (12)

    為了滿足連續(xù)性方程,波數(shù)k與參數(shù)p,q的關(guān)系應(yīng)該滿足

    (13)

    圖1 空流域中的瞬時(shí)漩渦結(jié)構(gòu)Fig.1 Instantaneous vortex structures in the empty domain

    圖2 模擬功率譜與入口功率譜比較Fig.2 Comparison of the simulated power spectrum and inlet power spectrum

    圖3 模擬風(fēng)場(chǎng)的時(shí)間相關(guān)性Fig.3 Time correlation of the simulated wind velocity components

    圖4 模擬風(fēng)場(chǎng)的空間相關(guān)性Fig.4 Spatial correlation of the simulated wind velocity components

    此外,本文還驗(yàn)證了流場(chǎng)的均勻性,結(jié)果表明橋梁斷面所處位置范圍內(nèi)的流場(chǎng)特征具有很好的均勻性,限于篇幅這里不再給出。

    2.2 數(shù)值計(jì)算工況

    選取扁平鋼箱梁作為研究對(duì)象,扁平鋼箱梁因具有流線化的氣動(dòng)外形,被大跨度橋梁廣泛采用。通常,在其斷面周圍除了棱角處有小范圍的流動(dòng)分離現(xiàn)象外,大部分區(qū)域處于附著流狀態(tài),計(jì)算模型尺寸如圖5所示。圖5中:s為以斷面迎風(fēng)端前緣為起始點(diǎn),沿著順時(shí)針方向各個(gè)壓強(qiáng)監(jiān)測(cè)點(diǎn)到起始點(diǎn)的距離與斷面周長(zhǎng)比值的無量綱距離。計(jì)算域與邊界條件如圖6所示。斷面前緣到入口處的距離為8B。斷面尾緣到出口處的距離為16B,上下邊界距離斷面形心的距離為8B,對(duì)應(yīng)的模型阻塞比為0.64%,小于風(fēng)洞試驗(yàn)?zāi)P妥枞炔粦?yīng)超過5%的要求。參照研究者以往研究成果,計(jì)算域展向長(zhǎng)度選為B[27]。計(jì)算域的選取一方面是為了盡量減少邊界條件對(duì)近壁面流動(dòng)和尾流流動(dòng)的影響;另一方面當(dāng)采用人工合成的速度場(chǎng)作為入口時(shí),減少非物理的流場(chǎng)在計(jì)算域中引起的壓強(qiáng)脈動(dòng)對(duì)斷面壓強(qiáng)的影響。邊界條件的設(shè)定為,入口處采用速度進(jìn)口,平均速度U為8.0 m/s,對(duì)應(yīng)的雷諾數(shù)Re=UB/ν=3.7×105,其中:ν為運(yùn)動(dòng)黏度系數(shù),出口處為紐曼邊界條件,前后邊界和上下邊界都采用周期性邊界。

    圖5 橋梁斷面尺寸(mm)Fig.5 Geometry of the bridge section (mm)

    網(wǎng)格采取混合型的劃分策略,在模型周圍采用結(jié)構(gòu)化的六面體網(wǎng)格,在此六面體網(wǎng)格之外采用四面體網(wǎng)格和六面體網(wǎng)格對(duì)斷面周圍和尾流區(qū)進(jìn)行加密,加密區(qū)之外恢復(fù)為結(jié)構(gòu)化的六面體網(wǎng)格。此外,為了減少湍流脈動(dòng)從入口處到模型所在位置處的衰減,需要對(duì)此區(qū)域的網(wǎng)格進(jìn)行加密,網(wǎng)格大小為0.029B。斷面周圍邊界層的網(wǎng)格尺寸如表1所示。表1中:x,y,z分別為順流向、橫向和豎向。表1中給出了3種不同的網(wǎng)格尺寸,用以驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確性,其中網(wǎng)格3為在網(wǎng)格2的基礎(chǔ)上繼續(xù)加密而來的,用于LES計(jì)算。網(wǎng)格1對(duì)應(yīng)的計(jì)算域網(wǎng)格如圖7所示。數(shù)值計(jì)算格式設(shè)置如下:動(dòng)量方程采用有界中心差分格式,壓力采用二階迎風(fēng)格式,時(shí)間離散采用有界二階隱式格式,時(shí)間步長(zhǎng)保障CFL數(shù)(courant number)小于1.0,壓力、速度耦合采用PISO(pressure-implicit with splitting of operators)算法。

    表1 邊界層網(wǎng)格尺寸Tab.1 Grid size in the boundary layer

    圖7 計(jì)算域網(wǎng)格Fig.7 Computational mesh

    3 數(shù)值計(jì)算結(jié)果

    采用SAS、SBES、WMLES計(jì)算了均勻來流下箱梁斷面的氣動(dòng)力,分別從氣動(dòng)力系數(shù)和流場(chǎng)特征分析了這3種方法的特點(diǎn)。然后采用SBES方法,分別在大尺度湍流和小尺度湍流下,識(shí)別了氣動(dòng)導(dǎo)納函數(shù)。

    3.1 壓強(qiáng)模擬結(jié)果

    在網(wǎng)格2,分別采用SAS、SBES、WMLES在均勻來流下得到壓強(qiáng)系數(shù)Cp沿箱梁斷面的時(shí)均值和脈動(dòng)均方根值σcp,如圖8所示。圖8中同時(shí)給出了LES模擬結(jié)果以及Yan等[28]風(fēng)洞試驗(yàn)均勻流場(chǎng)中的結(jié)果作為比較。LES的網(wǎng)格尺寸如表1中網(wǎng)格3所示。需要注意的是,雖然LES所用網(wǎng)格尺寸遠(yuǎn)小于其他兩種網(wǎng)格尺寸,但平行于壁面方向的網(wǎng)格尺寸仍不足以滿足LES的網(wǎng)格尺寸要求。

    圖8 壓強(qiáng)系數(shù)沿箱梁表面分布圖Fig.8 Distribution of pressure coefficients along the surface of box girder

    由圖8可知,對(duì)于平均壓強(qiáng)系數(shù),箱梁頂板處SAS的結(jié)果明顯大于試驗(yàn)值,相反WMLES的值小于試驗(yàn)值。除了頂板以外,其他位置處SAS和WMLES的值都與試驗(yàn)值吻合較好。整體上,SBES和LES的結(jié)果基本一致,都與試驗(yàn)值吻合較好。對(duì)于壓強(qiáng)的均方根值,SAS和WMLES的結(jié)果小于試驗(yàn)值,LES的結(jié)果大于試驗(yàn)值,SBES的結(jié)果與試驗(yàn)值較為接近。數(shù)值模擬結(jié)果與試驗(yàn)值的差別:一方面由于數(shù)值模擬的監(jiān)測(cè)點(diǎn)布置在靠近斷面的拐點(diǎn)處,而這在試驗(yàn)中是不可能的,因此造成拐點(diǎn)處的數(shù)值大于試驗(yàn)值;另一方面由于風(fēng)洞試驗(yàn)的流場(chǎng)并不是完全的均勻來流,湍流的存在抑制了斷面周圍的脈動(dòng)。

    Q準(zhǔn)則顯示的流場(chǎng)瞬時(shí)漩渦結(jié)構(gòu),如圖9所示。用流向速度V進(jìn)行渲染。由圖9可知,在斷面周圍以及尾部區(qū)域SBES方法都展示出了豐富的三維漩渦結(jié)構(gòu),流動(dòng)結(jié)構(gòu)與LES類似。相比于SBES方法,SAS方法和WMLES方法的流動(dòng)呈現(xiàn)出二維特性。由于SAS方法在建立之初就是為了應(yīng)用于大的流動(dòng)分離現(xiàn)象,當(dāng)流動(dòng)中存在顯著的流動(dòng)不穩(wěn)定時(shí),方程中的長(zhǎng)度尺度LvK會(huì)相應(yīng)調(diào)整以降低湍流黏度,從而釋放更多的渦。由圖9(c)可知,在前緣風(fēng)嘴與上下表面的轉(zhuǎn)折處有漩渦結(jié)構(gòu),其他區(qū)域流動(dòng)都呈現(xiàn)出二維特性。這說明SAS方法并不適用于這種流線型的箱梁斷面,在流動(dòng)的大部分區(qū)域SAS都表現(xiàn)為URANS(unsteady RANS)。WMLES方法在剪切層分離處會(huì)產(chǎn)生較大的湍流黏度,此外WMLES沒有屏蔽函數(shù),會(huì)產(chǎn)生較大的“灰區(qū)”,造成RANS?;耐牧黟ざ韧ㄟ^對(duì)流進(jìn)入LES區(qū)域,對(duì)LES區(qū)域的湍流脈動(dòng)帶來較大的抑制。為了促進(jìn)湍渦的發(fā)展,WMLES方法需要配合湍流入口使用。

    圖9 不同RANS/LES方法的流場(chǎng)瞬時(shí)漩渦結(jié)構(gòu)Fig.9 Instantaneous vortex structures with various RANS/LES models

    為了進(jìn)一步驗(yàn)證RANS/LES混合方法對(duì)復(fù)雜流動(dòng)現(xiàn)象的預(yù)測(cè)能力,分別將3種混合方法用于分離式雙箱梁斷面。分離式雙箱梁由于上下幅箱梁的相互干擾,斷面周圍的流場(chǎng)異常復(fù)雜。分離式雙箱梁的斷面尺寸見圖5(b)。分離式雙箱梁斷面的平均壓強(qiáng)和脈動(dòng)壓強(qiáng),如圖10、圖11所示。圖10、圖11同時(shí)給出了Li等[29]的風(fēng)洞試驗(yàn)結(jié)果以及吳步晨[30]的LES模擬結(jié)果,雷諾數(shù)Re=UD/ν=4.0×104(D為模型高度),與本文的3.85×104接近。由圖10、圖11可知,除了上游間隙的平均壓強(qiáng)和下游箱梁下表面的脈動(dòng)壓強(qiáng),SBES方法的結(jié)果與試驗(yàn)結(jié)果和吳步晨的LES結(jié)果吻合較好,顯示出SBES方法對(duì)復(fù)雜流動(dòng)現(xiàn)象較強(qiáng)的預(yù)測(cè)能力。SAS方法對(duì)下游箱梁平均壓強(qiáng)的預(yù)測(cè)與試驗(yàn)結(jié)果有一定偏差,但是對(duì)間隙壓強(qiáng)的預(yù)測(cè)在3種方法中與試驗(yàn)值卻最為接近。同時(shí)從僅有的下表面脈動(dòng)壓強(qiáng)風(fēng)洞試驗(yàn)結(jié)果可以看出,SAS的預(yù)測(cè)結(jié)果與試驗(yàn)值也最為接近。相比較而言,WMLES方法的模擬結(jié)果整體上與試驗(yàn)值相差較大,特別是對(duì)間隙壓強(qiáng)的預(yù)測(cè)。

    圖10 分離式雙箱梁平均壓強(qiáng)系數(shù)分布圖Fig.10 Distribution of mean pressure coefficients for the separated twin-box girder

    圖11 分離式雙箱梁脈動(dòng)壓強(qiáng)均方根分布Fig.11 Distribution of the root mean square of the fluctuating pressure for the separated twin-box girder

    Q準(zhǔn)則顯示的分離式雙箱梁流場(chǎng)瞬時(shí)漩渦結(jié)構(gòu)圖,如圖12所示。由圖12可知,在斷面周圍、間隙中和斷面尾部SBES方法同樣展示出了非常豐富的三維漩渦結(jié)構(gòu)。由于上游箱梁尾部具有顯著的渦脫,SAS方法在間隙處成功捕捉到了三維的流動(dòng)結(jié)構(gòu),隨著漩渦在下游箱梁的附著,流動(dòng)不穩(wěn)定現(xiàn)象減弱,SAS方法的湍流黏性增大,抑制了小尺度渦的發(fā)展,最終在尾部?jī)H有大尺度的漩渦結(jié)構(gòu)。對(duì)于WMLES方法,上游箱梁的流動(dòng)形態(tài)與流線型箱梁類似,都表現(xiàn)為近似的二維結(jié)構(gòu)。下游箱梁的表面有較大的漩渦,同時(shí)尾部的渦脫也呈現(xiàn)三維特性,但由壓強(qiáng)結(jié)果可知,WMLES方法對(duì)間隙和下游箱梁的流動(dòng)形態(tài)預(yù)測(cè)是不準(zhǔn)確的。

    圖12 分離式雙箱梁瞬時(shí)漩渦截圖Fig.12 Instantaneous vortex structures for the separated twin-box girder

    綜上分析可知,SBES方法無論是對(duì)流線型箱梁還是分離式雙箱梁的氣動(dòng)力和三維流動(dòng)結(jié)構(gòu)都具有非常好的預(yù)測(cè)能力,相反WMLES方法由于采用了混合長(zhǎng)度理論,并且在不同區(qū)域之間沒有采用屏蔽函數(shù),所以無法適用于流動(dòng)分離存在的流場(chǎng)。SAS方法對(duì)于具有顯著流動(dòng)分離的流場(chǎng)具有較好的預(yù)測(cè)能力,比如鈍體的雙主梁斷面、橋塔斷面等,對(duì)于寬高比較大附著流占主導(dǎo)地位的結(jié)構(gòu),預(yù)測(cè)能力和URANS類似。

    3.2 氣動(dòng)導(dǎo)納識(shí)別

    鑒于SBES方法對(duì)非定常氣動(dòng)力模擬的優(yōu)秀表現(xiàn),采用SBES方法對(duì)湍流來流下的箱梁斷面氣動(dòng)力進(jìn)行分析。網(wǎng)格1、網(wǎng)格2兩種不同網(wǎng)格精度下箱梁斷面在均勻來流下的平均壓強(qiáng)系數(shù),如圖13所示。由圖13可知,網(wǎng)格1和網(wǎng)格2的結(jié)果基本一致,并且與試驗(yàn)值吻合較好。壓強(qiáng)系數(shù)的均方根值具有相同的特點(diǎn),不再給出??紤]到計(jì)算效率,本文選取網(wǎng)格1作為計(jì)算網(wǎng)格。

    圖13 不同網(wǎng)格精度下壓強(qiáng)系數(shù)沿箱梁表面分布Fig.13 Mean pressure coefficient distribution of the box-girder section with different grid density

    表2 來流湍流特性Tab.2 Turbulence characteristics of the incoming flow

    圖14 湍流來流下壓強(qiáng)系數(shù)沿箱梁表面分布Fig.14 Pressure coefficient distribution of the box-girder section under turbulent flow

    (14)

    圖15 氣動(dòng)導(dǎo)納擬合結(jié)果Fig.15 Fitting of the aerodynamic admittance function

    在兩種風(fēng)場(chǎng)中識(shí)別的升力氣動(dòng)導(dǎo)納,如圖16所示。圖16同時(shí)給出了Sears函數(shù)和Yan等[32]的風(fēng)洞試驗(yàn)結(jié)果。由圖16可知,在K<3內(nèi),三者差別較大。具體表現(xiàn)為,小尺度湍流風(fēng)場(chǎng)2中的結(jié)果與Sears函數(shù)基本一致,而大尺度湍流風(fēng)場(chǎng)1中的結(jié)果明顯大于Sears函數(shù)和風(fēng)場(chǎng)2中的結(jié)果,這表明氣動(dòng)導(dǎo)納隨積分尺度的增大而增大,這與Larose等的研究結(jié)果一致。Yan等的試驗(yàn)結(jié)果與Sears函數(shù)的差別顯著,數(shù)值遠(yuǎn)小于Sears函數(shù),這種現(xiàn)象普遍存在于風(fēng)洞試驗(yàn)中,主要是由于風(fēng)洞試驗(yàn)中的湍流積分尺度偏小。相反,在折算頻率K>3內(nèi),三者結(jié)果與Sears具有較好的一致性。

    圖16 不同流場(chǎng)中的氣動(dòng)導(dǎo)納函數(shù)Fig.16 Aerodynamic admittance function in different turbulent flows

    從本文大積分尺度下氣動(dòng)導(dǎo)納的識(shí)別結(jié)果可以看出,積分尺度的增大會(huì)導(dǎo)致低頻處氣動(dòng)導(dǎo)納函數(shù)的增大。傳統(tǒng)的利用格柵湍流節(jié)段模型測(cè)力風(fēng)洞試驗(yàn)得到的結(jié)果,在數(shù)值上明顯小于大積分尺度湍流下的結(jié)果,這意味著如果將小尺度格柵湍流場(chǎng)中識(shí)別的氣動(dòng)導(dǎo)納用于實(shí)橋的抖振分析,必然會(huì)導(dǎo)致抖振力顯著偏小。

    4 結(jié) 論

    RANS/LES混合方法,通過在邊界層使用RANS在遠(yuǎn)離壁面區(qū)域使用LES,既避免了LES應(yīng)用于橋梁風(fēng)工程時(shí)所需的大量的計(jì)算代價(jià),也能滿足橋梁風(fēng)工程對(duì)于靜氣動(dòng)力和脈動(dòng)氣動(dòng)力的需求。本文將RANS/LES混合方法和入口湍流合成方法相結(jié)合,實(shí)現(xiàn)了大積分尺度湍流場(chǎng)的數(shù)值模擬,并將此用于箱梁斷面氣動(dòng)導(dǎo)納的識(shí)別。本文主要的研究結(jié)論如下:

    (1) 在以往人工合成湍流方法的基礎(chǔ)上,本文通過引進(jìn)時(shí)間尺度因子和與頻率有關(guān)的積分尺度,使得合成的風(fēng)速時(shí)程既滿足無散度條件,又能正確模擬自然風(fēng)的空間相關(guān)性和時(shí)間相關(guān)性。本文采用此方法,成功實(shí)現(xiàn)了大積分尺度湍流的人工模擬。

    (2) SBES方法無論是對(duì)流動(dòng)分離現(xiàn)象較弱的流線型箱梁還是存在氣動(dòng)干擾效應(yīng)的分離式雙箱梁,都能得到豐富的三維流動(dòng)結(jié)構(gòu)和準(zhǔn)確的氣動(dòng)力。相比較而言,SAS方法對(duì)于具有顯著流動(dòng)分離的流場(chǎng)具有較好的預(yù)測(cè)能力,而WMLES方法整體表現(xiàn)都較差。

    (3) 采用SBES方法分別在小尺度湍流和與實(shí)際風(fēng)環(huán)境相似的大尺度湍流中,計(jì)算了一箱梁斷面受到的非定常氣動(dòng)力,并識(shí)別得到了氣動(dòng)導(dǎo)納函數(shù)。整體上,小尺度湍流場(chǎng)中的結(jié)果與Sears函數(shù)吻合較好,而大尺度湍流場(chǎng)中的結(jié)果在低頻段明顯大于Sears函數(shù)和小尺度湍流場(chǎng)中的結(jié)果,尤其與節(jié)段模型風(fēng)洞試驗(yàn)的結(jié)果差距非常顯著。

    (4) 箱梁斷面的升力氣動(dòng)導(dǎo)納隨積分尺度的增大而增大,這表明如果采用小尺度格柵湍流場(chǎng)中識(shí)別的氣動(dòng)導(dǎo)納進(jìn)行橋梁抖振分析,會(huì)導(dǎo)致低頻處的抖振力小于實(shí)際值。

    猜你喜歡
    湍流脈動(dòng)壁面
    新學(xué)期,如何“脈動(dòng)回來”?
    家教世界(2023年25期)2023-10-09 02:11:56
    二維有限長(zhǎng)度柔性壁面上T-S波演化的數(shù)值研究
    RBI在超期服役脈動(dòng)真空滅菌器定檢中的應(yīng)用
    重氣瞬時(shí)泄漏擴(kuò)散的湍流模型驗(yàn)證
    地球脈動(dòng)(第一季)
    壁面溫度對(duì)微型內(nèi)燃機(jī)燃燒特性的影響
    “青春期”湍流中的智慧引渡(三)
    “青春期”湍流中的智慧引渡(二)
    顆?!诿媾鲎步Ec數(shù)據(jù)處理
    考慮裂縫壁面?zhèn)Φ膲毫丫a(chǎn)能計(jì)算模型
    国产又色又爽无遮挡免| 久久久欧美国产精品| 国产精品二区激情视频| 狠狠婷婷综合久久久久久88av| 国产女主播在线喷水免费视频网站| 免费观看人在逋| 最近最新中文字幕大全免费视频 | 母亲3免费完整高清在线观看| 日韩一本色道免费dvd| 免费在线观看黄色视频的| 日本午夜av视频| 99香蕉大伊视频| 啦啦啦在线观看免费高清www| 青春草亚洲视频在线观看| 色播在线永久视频| 老司机靠b影院| 国产免费又黄又爽又色| 黄色视频不卡| 免费高清在线观看视频在线观看| av网站在线播放免费| 丰满乱子伦码专区| 2021少妇久久久久久久久久久| 国产乱来视频区| 精品亚洲成a人片在线观看| 国产成人午夜福利电影在线观看| 亚洲国产精品一区二区三区在线| 亚洲精品视频女| 国产麻豆69| 黑人猛操日本美女一级片| 丰满迷人的少妇在线观看| 亚洲欧美一区二区三区久久| 亚洲av综合色区一区| 99久久人妻综合| 国产深夜福利视频在线观看| 99国产精品免费福利视频| 18在线观看网站| 天天操日日干夜夜撸| 男女下面插进去视频免费观看| 亚洲国产欧美一区二区综合| 国产亚洲最大av| 成年人午夜在线观看视频| 久久精品久久久久久久性| 久久久欧美国产精品| 精品国产一区二区三区久久久樱花| 久久性视频一级片| 又粗又硬又长又爽又黄的视频| 啦啦啦在线免费观看视频4| 精品国产国语对白av| 亚洲欧美日韩另类电影网站| 久久狼人影院| 欧美xxⅹ黑人| 桃花免费在线播放| 精品国产一区二区三区久久久樱花| 亚洲一码二码三码区别大吗| 久久久久国产一级毛片高清牌| 又黄又粗又硬又大视频| 人人妻,人人澡人人爽秒播 | 精品人妻在线不人妻| 亚洲精华国产精华液的使用体验| 51午夜福利影视在线观看| 午夜av观看不卡| 亚洲国产欧美在线一区| 亚洲精品国产av成人精品| 亚洲国产欧美日韩在线播放| 考比视频在线观看| 交换朋友夫妻互换小说| 七月丁香在线播放| av卡一久久| 侵犯人妻中文字幕一二三四区| 日韩制服丝袜自拍偷拍| 99国产精品免费福利视频| 国产女主播在线喷水免费视频网站| 亚洲美女视频黄频| 亚洲欧洲日产国产| 午夜91福利影院| 国产精品久久久久久久久免| 久久天堂一区二区三区四区| 人成视频在线观看免费观看| av片东京热男人的天堂| 99久久99久久久精品蜜桃| 一级黄片播放器| 老司机靠b影院| 久久毛片免费看一区二区三区| av.在线天堂| 在线观看免费午夜福利视频| 亚洲精品久久午夜乱码| 亚洲精品在线美女| 欧美在线黄色| 日韩熟女老妇一区二区性免费视频| 麻豆av在线久日| 久久这里只有精品19| 久久午夜综合久久蜜桃| 日日爽夜夜爽网站| 三上悠亚av全集在线观看| 久久久久久人妻| 1024香蕉在线观看| 日韩av在线免费看完整版不卡| 中文字幕色久视频| 99久久综合免费| 999精品在线视频| 欧美在线黄色| 女性生殖器流出的白浆| 男的添女的下面高潮视频| 色婷婷av一区二区三区视频| 黄色毛片三级朝国网站| 久久人人97超碰香蕉20202| 午夜福利一区二区在线看| 免费观看av网站的网址| 一本—道久久a久久精品蜜桃钙片| 国产免费视频播放在线视频| www.自偷自拍.com| 国产成人精品在线电影| 桃花免费在线播放| 国产免费一区二区三区四区乱码| 男人爽女人下面视频在线观看| 亚洲av综合色区一区| 91成人精品电影| 婷婷成人精品国产| 美女午夜性视频免费| 欧美在线黄色| 男女床上黄色一级片免费看| av视频免费观看在线观看| 国产成人精品在线电影| 日日撸夜夜添| 99久国产av精品国产电影| 老熟女久久久| 精品国产一区二区三区久久久樱花| 久久久久精品久久久久真实原创| 亚洲精品av麻豆狂野| 一级黄片播放器| 中文字幕亚洲精品专区| 日本91视频免费播放| 18禁动态无遮挡网站| 纯流量卡能插随身wifi吗| 无遮挡黄片免费观看| 两个人免费观看高清视频| 国产成人一区二区在线| 另类亚洲欧美激情| 亚洲久久久国产精品| 国产av精品麻豆| 天天躁夜夜躁狠狠久久av| 久久性视频一级片| 国产亚洲最大av| 丝袜脚勾引网站| 丝袜喷水一区| 免费高清在线观看日韩| 欧美日韩福利视频一区二区| 青青草视频在线视频观看| 狂野欧美激情性bbbbbb| av有码第一页| 又粗又硬又长又爽又黄的视频| 亚洲三区欧美一区| 五月开心婷婷网| 少妇的丰满在线观看| 国产 精品1| 十八禁网站网址无遮挡| 狂野欧美激情性xxxx| 亚洲精品成人av观看孕妇| 午夜免费男女啪啪视频观看| 亚洲国产欧美一区二区综合| 热re99久久国产66热| 亚洲国产看品久久| 久久久久精品久久久久真实原创| 超碰成人久久| 1024香蕉在线观看| 十八禁人妻一区二区| 日本爱情动作片www.在线观看| 性色av一级| 久久毛片免费看一区二区三区| 亚洲av欧美aⅴ国产| 国产精品无大码| 亚洲国产av新网站| 黄片无遮挡物在线观看| 国产黄色免费在线视频| av在线观看视频网站免费| 一本大道久久a久久精品| 国产男人的电影天堂91| 精品午夜福利在线看| 成人18禁高潮啪啪吃奶动态图| 久久综合国产亚洲精品| 美女福利国产在线| a 毛片基地| 国产精品久久久久久久久免| 国产精品香港三级国产av潘金莲 | 国产野战对白在线观看| 免费看不卡的av| 国产精品二区激情视频| 久久免费观看电影| 美女扒开内裤让男人捅视频| 在现免费观看毛片| 国产在线一区二区三区精| 看十八女毛片水多多多| 亚洲免费av在线视频| 亚洲婷婷狠狠爱综合网| 国产xxxxx性猛交| 天堂中文最新版在线下载| 国产精品亚洲av一区麻豆 | 国产日韩欧美视频二区| 国产精品二区激情视频| 欧美日韩福利视频一区二区| 精品国产国语对白av| 男女午夜视频在线观看| 成人毛片60女人毛片免费| 国产精品人妻久久久影院| 青青草视频在线视频观看| 国产精品久久久久久精品电影小说| 欧美在线黄色| 美女大奶头黄色视频| 一级片'在线观看视频| 最近2019中文字幕mv第一页| 免费在线观看完整版高清| 亚洲天堂av无毛| 一区二区三区四区激情视频| 丝袜喷水一区| 中文天堂在线官网| 综合色丁香网| 欧美 日韩 精品 国产| 女人久久www免费人成看片| 亚洲欧洲国产日韩| 国产精品熟女久久久久浪| 国产精品香港三级国产av潘金莲 | 午夜福利在线免费观看网站| 国产成人欧美在线观看 | 两个人免费观看高清视频| 美女脱内裤让男人舔精品视频| 久久久久国产一级毛片高清牌| 久热爱精品视频在线9| 国产成人一区二区在线| 少妇被粗大猛烈的视频| 色婷婷久久久亚洲欧美| 亚洲,欧美精品.| 欧美精品一区二区免费开放| 欧美乱码精品一区二区三区| 男女无遮挡免费网站观看| 欧美xxⅹ黑人| 天美传媒精品一区二区| 十八禁高潮呻吟视频| 欧美日韩一区二区视频在线观看视频在线| 国产成人a∨麻豆精品| 婷婷色av中文字幕| 午夜福利免费观看在线| 亚洲情色 制服丝袜| 超碰97精品在线观看| 美女高潮到喷水免费观看| 亚洲国产av新网站| 欧美日韩国产mv在线观看视频| 涩涩av久久男人的天堂| 老鸭窝网址在线观看| 制服人妻中文乱码| 国产日韩欧美亚洲二区| 亚洲av中文av极速乱| 啦啦啦中文免费视频观看日本| 国产爽快片一区二区三区| 午夜福利免费观看在线| 夫妻午夜视频| 精品人妻熟女毛片av久久网站| 免费在线观看黄色视频的| 国产在视频线精品| av视频免费观看在线观看| 亚洲国产欧美一区二区综合| 久久韩国三级中文字幕| 80岁老熟妇乱子伦牲交| 一区二区av电影网| 黄色视频在线播放观看不卡| 国产一区二区三区综合在线观看| 国产一区二区在线观看av| 国产成人a∨麻豆精品| 国产精品香港三级国产av潘金莲 | 大陆偷拍与自拍| 男女边吃奶边做爰视频| 自线自在国产av| 亚洲七黄色美女视频| 久久影院123| 国产欧美日韩一区二区三区在线| 免费观看人在逋| 男女国产视频网站| 国产乱人偷精品视频| 久久久亚洲精品成人影院| 看免费av毛片| 亚洲精品久久成人aⅴ小说| 91国产中文字幕| 久久 成人 亚洲| 欧美av亚洲av综合av国产av | 亚洲综合精品二区| 午夜福利一区二区在线看| 日韩伦理黄色片| 欧美人与性动交α欧美软件| 日韩视频在线欧美| 日本色播在线视频| 亚洲综合色网址| 在现免费观看毛片| 人人妻人人爽人人添夜夜欢视频| 日韩中文字幕视频在线看片| 亚洲熟女毛片儿| 黑丝袜美女国产一区| 午夜久久久在线观看| 亚洲美女搞黄在线观看| 国产极品粉嫩免费观看在线| av网站在线播放免费| 啦啦啦在线免费观看视频4| 18禁国产床啪视频网站| 欧美国产精品一级二级三级| 久热这里只有精品99| 制服诱惑二区| 黄网站色视频无遮挡免费观看| 亚洲欧美一区二区三区黑人| 大片免费播放器 马上看| 国产在线视频一区二区| www.精华液| 精品亚洲成国产av| 99香蕉大伊视频| 亚洲精品在线美女| 国产精品无大码| 99热全是精品| 中文天堂在线官网| 欧美日韩av久久| 国产有黄有色有爽视频| 亚洲国产欧美网| 国产毛片在线视频| 国产日韩一区二区三区精品不卡| a级片在线免费高清观看视频| 大片电影免费在线观看免费| 波多野结衣av一区二区av| 日韩不卡一区二区三区视频在线| 国产欧美日韩一区二区三区在线| 免费在线观看黄色视频的| 亚洲成人一二三区av| 1024视频免费在线观看| 韩国av在线不卡| 伊人亚洲综合成人网| 亚洲四区av| 欧美日韩一区二区视频在线观看视频在线| 一本大道久久a久久精品| 免费高清在线观看日韩| 女人爽到高潮嗷嗷叫在线视频| 国产一区二区激情短视频 | 午夜免费鲁丝| 久久av网站| 男女国产视频网站| 午夜福利一区二区在线看| 大片电影免费在线观看免费| 久久婷婷青草| 多毛熟女@视频| 男人添女人高潮全过程视频| 熟女少妇亚洲综合色aaa.| 日本欧美视频一区| 亚洲av福利一区| 一个人免费看片子| 亚洲国产欧美在线一区| 国产淫语在线视频| 国产成人免费无遮挡视频| 91aial.com中文字幕在线观看| 老司机在亚洲福利影院| 一级,二级,三级黄色视频| 国产 精品1| 精品福利永久在线观看| 97精品久久久久久久久久精品| 纵有疾风起免费观看全集完整版| 老司机在亚洲福利影院| 亚洲精品视频女| 最近的中文字幕免费完整| 久久久亚洲精品成人影院| 水蜜桃什么品种好| 日韩av不卡免费在线播放| 国产精品99久久99久久久不卡 | 飞空精品影院首页| 香蕉国产在线看| 超碰97精品在线观看| 欧美精品一区二区免费开放| 亚洲国产成人一精品久久久| 久久久精品国产亚洲av高清涩受| 精品亚洲成a人片在线观看| 99久久精品国产亚洲精品| 亚洲av成人不卡在线观看播放网 | 国产精品一国产av| 纯流量卡能插随身wifi吗| 最黄视频免费看| 国产精品久久久久久久久免| 亚洲欧洲国产日韩| 精品亚洲成国产av| 啦啦啦在线观看免费高清www| 日本欧美视频一区| 国产精品 欧美亚洲| 日韩一区二区视频免费看| 日本猛色少妇xxxxx猛交久久| 秋霞伦理黄片| 亚洲国产精品一区三区| 精品人妻熟女毛片av久久网站| xxxhd国产人妻xxx| 国产一区二区激情短视频 | 久久久久久久精品精品| 看免费成人av毛片| 精品少妇内射三级| 欧美日韩一区二区视频在线观看视频在线| 成年人午夜在线观看视频| 欧美精品一区二区大全| 久久久久国产精品人妻一区二区| 一区二区三区精品91| 精品国产超薄肉色丝袜足j| 香蕉国产在线看| 欧美成人精品欧美一级黄| 别揉我奶头~嗯~啊~动态视频 | 久久久久久免费高清国产稀缺| 国产黄频视频在线观看| 国产精品蜜桃在线观看| 日韩一卡2卡3卡4卡2021年| 91成人精品电影| av又黄又爽大尺度在线免费看| 男女无遮挡免费网站观看| 多毛熟女@视频| 蜜桃在线观看..| 中文字幕最新亚洲高清| 日韩,欧美,国产一区二区三区| 久久精品久久精品一区二区三区| 久久国产精品大桥未久av| 精品久久久精品久久久| 高清黄色对白视频在线免费看| 久久久久精品性色| 国产精品偷伦视频观看了| 欧美成人精品欧美一级黄| 亚洲av成人精品一二三区| 新久久久久国产一级毛片| av又黄又爽大尺度在线免费看| 久久婷婷青草| 国产精品女同一区二区软件| 九草在线视频观看| 女人爽到高潮嗷嗷叫在线视频| 国产免费福利视频在线观看| 侵犯人妻中文字幕一二三四区| 成人手机av| 黄色 视频免费看| 亚洲欧美成人综合另类久久久| 精品亚洲成a人片在线观看| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲精品国产av蜜桃| 欧美日韩一级在线毛片| 无遮挡黄片免费观看| 亚洲精品国产一区二区精华液| 日本欧美视频一区| 国产精品一区二区在线不卡| 成年人免费黄色播放视频| 亚洲国产欧美网| 天堂中文最新版在线下载| 亚洲欧美精品综合一区二区三区| 国产精品无大码| 丝瓜视频免费看黄片| 亚洲av男天堂| 国产淫语在线视频| 天堂8中文在线网| 五月开心婷婷网| 夜夜骑夜夜射夜夜干| 欧美精品一区二区大全| 婷婷色av中文字幕| 日日撸夜夜添| av视频免费观看在线观看| 综合色丁香网| 色网站视频免费| www日本在线高清视频| 日韩一本色道免费dvd| 久久精品久久久久久久性| 日韩精品有码人妻一区| 国产在线一区二区三区精| 青草久久国产| 麻豆精品久久久久久蜜桃| av福利片在线| 人人妻,人人澡人人爽秒播 | www.熟女人妻精品国产| 亚洲国产日韩一区二区| 丰满少妇做爰视频| 久久天躁狠狠躁夜夜2o2o | 不卡视频在线观看欧美| 亚洲精品国产av成人精品| 一级片免费观看大全| 在线观看www视频免费| 精品国产一区二区久久| h视频一区二区三区| 黑人欧美特级aaaaaa片| 精品国产一区二区三区久久久樱花| 久久久久精品性色| 免费观看a级毛片全部| 91成人精品电影| 男女国产视频网站| 亚洲av日韩精品久久久久久密 | 久久久精品区二区三区| 99久久精品国产亚洲精品| 午夜福利乱码中文字幕| 一区二区三区四区激情视频| 久久精品国产亚洲av高清一级| 一边摸一边做爽爽视频免费| 日本一区二区免费在线视频| 一级爰片在线观看| 久久久久久久久免费视频了| 成人亚洲欧美一区二区av| 国产老妇伦熟女老妇高清| 黄片无遮挡物在线观看| 免费黄频网站在线观看国产| 丰满迷人的少妇在线观看| 亚洲av电影在线进入| a级毛片黄视频| 国产精品一国产av| 国产精品麻豆人妻色哟哟久久| 搡老乐熟女国产| 久久女婷五月综合色啪小说| 午夜激情av网站| 别揉我奶头~嗯~啊~动态视频 | 国产乱人偷精品视频| av网站在线播放免费| 亚洲精品美女久久久久99蜜臀 | 少妇精品久久久久久久| 黑人欧美特级aaaaaa片| 国产熟女欧美一区二区| 建设人人有责人人尽责人人享有的| 日韩精品有码人妻一区| 在线观看免费视频网站a站| av国产精品久久久久影院| 国产av码专区亚洲av| 少妇 在线观看| 久久久久精品久久久久真实原创| 秋霞伦理黄片| 中文字幕高清在线视频| 男女下面插进去视频免费观看| 亚洲成人一二三区av| 亚洲成国产人片在线观看| 免费在线观看完整版高清| 精品久久久精品久久久| 欧美亚洲 丝袜 人妻 在线| 麻豆精品久久久久久蜜桃| av片东京热男人的天堂| 亚洲,欧美精品.| 久久精品国产综合久久久| 色精品久久人妻99蜜桃| 99久久人妻综合| 亚洲欧洲日产国产| 国产免费现黄频在线看| 一二三四在线观看免费中文在| 国产极品粉嫩免费观看在线| 久热这里只有精品99| 91精品伊人久久大香线蕉| 久久免费观看电影| 亚洲精品一区蜜桃| 精品一区在线观看国产| 国产免费现黄频在线看| 97精品久久久久久久久久精品| 在线观看国产h片| 国产高清不卡午夜福利| 黄网站色视频无遮挡免费观看| 最新在线观看一区二区三区 | 国产av码专区亚洲av| 性色av一级| 国产一区亚洲一区在线观看| 18在线观看网站| 少妇 在线观看| 人人妻人人爽人人添夜夜欢视频| 亚洲综合精品二区| 日本av免费视频播放| 久久国产精品男人的天堂亚洲| xxxhd国产人妻xxx| 欧美精品一区二区大全| 黑人巨大精品欧美一区二区蜜桃| 男女国产视频网站| 丝袜美腿诱惑在线| 婷婷成人精品国产| 我要看黄色一级片免费的| 两性夫妻黄色片| 精品一区二区三卡| 美女福利国产在线| 久久免费观看电影| 久久精品国产a三级三级三级| 久久久精品94久久精品| 国产精品av久久久久免费| 亚洲熟女毛片儿| 国产精品人妻久久久影院| 亚洲国产日韩一区二区| 亚洲第一av免费看| 久久久久久久精品精品| 伊人亚洲综合成人网| 人人妻,人人澡人人爽秒播 | 热re99久久国产66热| 汤姆久久久久久久影院中文字幕| 亚洲成人免费av在线播放| 伦理电影免费视频| 免费女性裸体啪啪无遮挡网站| 亚洲精品国产区一区二| e午夜精品久久久久久久| 亚洲精品久久成人aⅴ小说| 国产免费一区二区三区四区乱码| 亚洲精品自拍成人| 伊人久久国产一区二区| 大码成人一级视频| videosex国产| 久久久精品区二区三区| 日韩av免费高清视频| 夫妻性生交免费视频一级片| 建设人人有责人人尽责人人享有的| 在线观看三级黄色| 成人亚洲精品一区在线观看| 韩国av在线不卡| 亚洲精品乱久久久久久| 19禁男女啪啪无遮挡网站| 欧美精品亚洲一区二区| 亚洲伊人久久精品综合| 尾随美女入室| 在现免费观看毛片| 日本wwww免费看| 亚洲,欧美精品.| 一本大道久久a久久精品| 欧美精品一区二区免费开放| 日韩精品免费视频一区二区三区| 久久精品国产综合久久久| 成人黄色视频免费在线看| 日韩欧美一区视频在线观看| 欧美日韩国产mv在线观看视频| 成人手机av| 久久精品熟女亚洲av麻豆精品| 嫩草影院入口|