亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于Lighthill聲類(lèi)比的流激噪聲三維計(jì)算及驗(yàn)證

        2014-12-05 02:19:18張?jiān)侜t劉繼明
        艦船科學(xué)技術(shù) 2014年9期
        關(guān)鍵詞:模型

        張?jiān)侜t,張 濤,劉繼明,趙 威,李 奇

        (1.華中科技大學(xué) 船舶與海洋工程學(xué)院,湖北 武漢430074;2.船舶和海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室,湖北 武漢430074;3.武漢第二船舶設(shè)計(jì)研究所,湖北 武漢430064;4.上海?;⒃畔⒖萍加邢薰荆虾?00235)

        0 引 言

        隨著艦艇管路系統(tǒng)流激噪聲問(wèn)題的日益突出,流激噪聲數(shù)值計(jì)算方法逐漸受到關(guān)注。從1952年開(kāi)創(chuàng)性的Lighthill 方程提出后,近代聲學(xué)開(kāi)始了蓬勃的發(fā)展。近年來(lái),氣動(dòng)聲學(xué)方面已取得廣泛的研究成果,但是水動(dòng)力聲學(xué)方面的研究卻進(jìn)展緩慢。

        對(duì)于低馬赫數(shù)的流激噪聲數(shù)值計(jì)算,國(guó)內(nèi)外都提出了不同的計(jì)算方法。張楠[1-2]等通過(guò)大渦模擬結(jié)合FW-H 聲類(lèi)比方法對(duì)空腔流激噪聲問(wèn)題進(jìn)行了數(shù)值計(jì)算與驗(yàn)證,并利用該方法預(yù)報(bào)了三維空腔的流激噪聲。耿冬寒、劉正先[3]利用大渦模擬-Lighthill 等效聲源法對(duì)二維空腔的水動(dòng)力噪聲進(jìn)行了預(yù)測(cè)。Moon[4]等利用LES/LPCE 混合方法對(duì)臺(tái)階繞流聲學(xué)問(wèn)題進(jìn)行了計(jì)算,其結(jié)果與試驗(yàn)數(shù)據(jù)相符。Ji和Wang[5]以臺(tái)階模型為研究對(duì)象,利用LES和Lighthill 理論進(jìn)行了低馬赫數(shù)下流場(chǎng)和聲學(xué)求解,并與實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。此外,劉敏[6]等基于邊界元方法對(duì)水下流噪聲問(wèn)題進(jìn)行了分析。

        在艦艇流激噪聲研究方面,張?jiān)剩?]等人在二維模型驗(yàn)證的基礎(chǔ)上,對(duì)三維開(kāi)孔潛體流激噪聲進(jìn)行了數(shù)值預(yù)報(bào)。江文成[8-9]等對(duì)潛體流噪聲和流固耦合產(chǎn)生的結(jié)構(gòu)噪聲進(jìn)行了研究。在管路系統(tǒng)流激噪聲研究方面,袁壽其[10]和趙威[11]等人對(duì)離心泵內(nèi)部流激噪聲進(jìn)行了數(shù)值計(jì)算,并進(jìn)行了一定的驗(yàn)證。Liu[12]等采用大渦模擬對(duì)蒸汽截止閥的流激噪聲進(jìn)行了預(yù)報(bào)。郭濤[13]采用大渦模擬結(jié)合Lighthill 聲類(lèi)比對(duì)直管和彎頭的流激噪聲進(jìn)行了分析。張碩[14]等采用FW- H 聲類(lèi)比對(duì)三通、四通管路流激噪聲進(jìn)行了預(yù)報(bào)并提出了降噪方案。

        針對(duì)閥門(mén)及管路流激噪聲問(wèn)題,本文在前期對(duì)泵[11]及截止閥[12]流激噪聲進(jìn)行預(yù)測(cè)的工作基礎(chǔ)上,對(duì)大渦模擬結(jié)合Lighthill 聲類(lèi)比的混合求解方法進(jìn)行了計(jì)算及驗(yàn)證。鑒于水動(dòng)力聲學(xué)試驗(yàn)費(fèi)用昂貴,試驗(yàn)條件要求較高。本文選取的驗(yàn)證對(duì)象為L(zhǎng)afon[15]在2003年于法國(guó)航空技術(shù)研究所中所使用的蒸氣管路開(kāi)口閥門(mén)模型。計(jì)算時(shí)先采用大渦模擬的方法進(jìn)行低馬赫數(shù)下三維類(lèi)閥空腔模型的非定常流場(chǎng)計(jì)算。然后將流場(chǎng)結(jié)果作為噪聲源,利用ACTRAN 基于有限元/無(wú)限元計(jì)算的Lighthill 聲類(lèi)比法進(jìn)行聲場(chǎng)求解。最終將計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行了對(duì)比,驗(yàn)證了該混合方法可行,計(jì)算結(jié)果可靠。

        1 流激噪聲數(shù)值模擬方法

        1.1 大渦模擬

        本文采用大渦模擬對(duì)非定常流場(chǎng)進(jìn)行仿真。大渦模擬將湍流中的渦按照尺度分成大小兩類(lèi)。大渦的形態(tài)和強(qiáng)度因流動(dòng)而異,是高度各向異性的,大部分質(zhì)量、動(dòng)量、能量的輸運(yùn)是由大渦引起的。而小渦主要是通過(guò)大渦之間的非線(xiàn)性相互作用間接產(chǎn)生的,它與平均運(yùn)動(dòng)或流場(chǎng)邊界形狀幾乎沒(méi)有關(guān)系,因而近似是各向同性的。

        通過(guò)將非定常的N-S 方程進(jìn)行濾波,得到大渦模擬的控制方程。濾波過(guò)程有效地過(guò)濾掉了那些尺度小于濾波寬度(或網(wǎng)格尺度)的小渦[2]。

        濾波函數(shù)G(x,x′)取為:

        連續(xù)性方程可以寫(xiě)為:

        則濾波后的NS 方程如下式所示:

        式中:σij為分子粘性引起的應(yīng)力張量;τij為亞格子雷諾應(yīng)力。本文采用LES WALE 模型來(lái)模擬亞格子尺度效應(yīng)。

        1.2 Lighthill 聲類(lèi)比法

        Lighthill 聲類(lèi)比理論由流體力學(xué)基本N-S 方程導(dǎo)出。由于方程的非線(xiàn)性和流動(dòng)與聲場(chǎng)的耦合性使方程不易求解,將聲場(chǎng)分為近場(chǎng)和遠(yuǎn)場(chǎng)。近場(chǎng)為聲源區(qū),遠(yuǎn)場(chǎng)為輻射區(qū),假定輻射區(qū)的流動(dòng)對(duì)聲場(chǎng)沒(méi)有影響。在該假定下通過(guò)連續(xù)方程和動(dòng)量方程簡(jiǎn)化得到Lighthill 聲類(lèi)比方程:

        式中:c0為等熵條件下的聲速值;ρ′=ρ- ρ0,ρ 與ρ′分別為擾動(dòng)與未擾動(dòng)時(shí)的流體密度;為L(zhǎng)ighthill應(yīng)力張量,定義為

        本文中的聲學(xué)計(jì)算采用軟件ACTRAN,其基于有限元和無(wú)限元的方法不僅可以考慮偶極子噪聲,也可以考慮由湍流引起的四極子噪聲。

        1.3 流激噪聲混合模擬方法

        本文流激噪聲混合計(jì)算步驟可總結(jié)如下:

        1)流場(chǎng)計(jì)算:首先基于不可壓流體假設(shè)采用CFX 中的k-ε 模型對(duì)模型進(jìn)行定常求解;然后,以定常計(jì)算得到的流場(chǎng)為初始條件,通過(guò)大渦模擬對(duì)模型非定常流場(chǎng)進(jìn)行計(jì)算;最后,待測(cè)點(diǎn)壓力變化范圍基本不變后,保存并導(dǎo)出流場(chǎng)計(jì)算結(jié)果;

        2)網(wǎng)格插值:將流場(chǎng)計(jì)算結(jié)果,流場(chǎng)計(jì)算網(wǎng)格信息以及聲學(xué)計(jì)算網(wǎng)格信息導(dǎo)入ACTRAN 的iCFD模塊,對(duì)計(jì)算結(jié)果進(jìn)行插值和傅里葉變換。其中流場(chǎng)計(jì)算結(jié)果插值到聲源區(qū)作為聲學(xué)計(jì)算時(shí)的聲源項(xiàng);

        3)聲場(chǎng)求解:在傅里葉變換后的流場(chǎng)計(jì)算結(jié)果基礎(chǔ)上,采用ACTRAN 的Lighthill 聲類(lèi)比理論和有限元/無(wú)限元計(jì)算方法,對(duì)可壓的流體聲學(xué)特性進(jìn)行計(jì)算。

        2 計(jì)算模型及網(wǎng)格

        針對(duì)艦船管路系統(tǒng)流激噪聲問(wèn)題,本文采用研究較多的類(lèi)閥模型的低馬赫數(shù)氣動(dòng)聲學(xué)試驗(yàn)結(jié)果進(jìn)行驗(yàn)證,即Lafon[15]在法國(guó)航空技術(shù)研究所試驗(yàn)的模型,試驗(yàn)風(fēng)速為U0=62.8 m/s,馬赫數(shù)Ma=0.183。

        本文采用的仿真模型幾何尺寸與試驗(yàn)?zāi)P拖嗤?,測(cè)點(diǎn)位置及模型尺寸如圖1所示。其中:H=0.137 m,d=0.05 m,h=0.02 m,h1=0.008 m,L=0.073 m。仿真空腔外流場(chǎng)上游長(zhǎng)0.075 m,下游長(zhǎng)0.175 m,模型寬度為0.02 m。

        圖1 空腔模型示意圖Fig.1 Characteristic dimensions of cavity model

        在采用混合方法求解流激噪聲時(shí),時(shí)間步長(zhǎng)為0.000 08 s,大渦模擬穩(wěn)定后保存5 000 步作為聲場(chǎng)計(jì)算輸入。流場(chǎng)計(jì)算網(wǎng)格和聲學(xué)計(jì)算網(wǎng)格均為ICEM劃分的結(jié)構(gòu)網(wǎng)格,且在近壁面處對(duì)網(wǎng)格進(jìn)行了加密處理,具體網(wǎng)格參數(shù)如表1所示。在流場(chǎng)計(jì)算時(shí)CFL 數(shù)為1 ~5,所有壁面y+值小于100。

        表1 流場(chǎng)及聲學(xué)計(jì)算網(wǎng)格參數(shù)表Tab.1 CFD and acoustic mesh parameters

        流場(chǎng)計(jì)算邊界條件和三維模型整體及局部網(wǎng)格示意圖如圖2所示。模型網(wǎng)格寬度方向劃分為20層。流場(chǎng)計(jì)算邊界條件如下:

        1)入口采用速度入口邊界條件,速度分布采用試驗(yàn)測(cè)量值;

        2)出口采用壓力出口邊界條件;

        3)壁面1、2 為滑移壁面,壁面3、4 為無(wú)滑移壁面。垂直于壁面1 ~4 的兩側(cè)面采用對(duì)稱(chēng)壁面邊界條件。

        圖2 流場(chǎng)網(wǎng)格及邊界條件示意圖Fig.2 CFD mesh with boundary conditions

        聲學(xué)網(wǎng)格及邊界如圖3所示。聲學(xué)網(wǎng)格中模型寬度方向劃分為10 層網(wǎng)格。

        通過(guò)將大渦模擬的湍流流場(chǎng)計(jì)算結(jié)果插值到圖3中聲源區(qū)聲學(xué)網(wǎng)格來(lái)模擬噪聲聲源。為接近實(shí)際試驗(yàn)邊界條件,聲源區(qū)兩端設(shè)置不導(dǎo)入流場(chǎng)結(jié)果的聲傳播區(qū)。同時(shí)聲傳播區(qū)兩邊界設(shè)為無(wú)限透射邊界,通過(guò)無(wú)限元對(duì)無(wú)網(wǎng)格區(qū)域聲學(xué)進(jìn)行計(jì)算。除無(wú)限透射邊界外,其他邊界均為剛性壁面。

        圖3 聲學(xué)網(wǎng)格及邊界條件示意圖Fig.3 Acoustic mesh with boundary conditions

        3 結(jié)果及分析

        3.1 流場(chǎng)結(jié)果分析

        大渦模擬穩(wěn)定后,空腔內(nèi)觀測(cè)點(diǎn)壓力時(shí)域曲線(xiàn)如圖4所示,圖中橫坐標(biāo)為無(wú)因次化處理的流場(chǎng)時(shí)間,縱坐標(biāo)為觀測(cè)點(diǎn)壓力值。由圖可以看到測(cè)點(diǎn)壓力在一定范圍內(nèi)波動(dòng),壓力變化范圍不隨時(shí)間明顯變化,可以認(rèn)為大渦模擬計(jì)算已趨于穩(wěn)定。

        圖4 測(cè)點(diǎn)壓力時(shí)域變化曲線(xiàn)Fig.4 Pressure fluctuations at the monitor point

        流場(chǎng)穩(wěn)定后某時(shí)刻類(lèi)閥空腔截面的渦量云圖如圖5所示,圖中從左往右為外流場(chǎng)來(lái)流方向。流體流過(guò)空腔后,在空腔內(nèi)形成了大量渦,但對(duì)外流場(chǎng)影響不大。此外,空腔對(duì)邊界層分布產(chǎn)生了一定的影響,特別是在圖中區(qū)域1和區(qū)域2 流場(chǎng)變化較大,上游邊界層在區(qū)域2 處破碎并重新形成新的邊界層,導(dǎo)致此處渦量值較大。

        由圖5 知,渦量較大的區(qū)域主要在空腔內(nèi)和空腔口。由渦聲理論可知,低馬赫數(shù)下流場(chǎng)渦量較大的區(qū)域?qū)⒊蔀橹饕肼曉础?/p>

        圖5 截面流場(chǎng)渦量云圖Fig.5 Vorticity contour at the cross section

        圖6 給出了500 Hz,1 000 Hz和1 500 Hz 頻率下的空腔處聲壓級(jí)云圖。如圖所示,與流場(chǎng)特性相似,主要噪聲源在空腔內(nèi)和區(qū)域2,其中區(qū)域2 在各頻段內(nèi)聲壓級(jí)均較高,區(qū)域1和空腔內(nèi)聲壓級(jí)在低頻時(shí)也有較大值。

        圖6 500 Hz,1 000 Hz和1 500 Hz 聲壓級(jí)云圖Fig.6 Sound pressure level contour at 500 Hz,1 000 Hz and 1 500 Hz

        為與文獻(xiàn)[15]中測(cè)點(diǎn)聲壓級(jí)(見(jiàn)圖7)比較,本文在與試驗(yàn)測(cè)點(diǎn)相同位置提取聲壓,其在100 ~3 200 Hz 頻段內(nèi)聲壓級(jí)曲線(xiàn)如圖8所示。

        圖7 文獻(xiàn)[15]聲壓測(cè)量曲線(xiàn)Fig.7 Experimental data in[15]

        圖8 本文大渦模擬和聲類(lèi)比混合法測(cè)點(diǎn)聲壓級(jí)曲線(xiàn)Fig.8 Sound pressure level obtained by the hybrid method in the paper

        對(duì)比圖7和圖8 可知,在測(cè)點(diǎn)處,聲壓級(jí)從100 Hz 的100 dB 左右降低到3 200 Hz 的70 dB 以下,整體上隨著頻率的增加而逐漸降低。其中在1 200 Hz和2 400 Hz 左右出現(xiàn)共振峰值,共振峰值的試驗(yàn)值與仿真計(jì)算值對(duì)比如表2所示。其中誤差百分?jǐn)?shù)為試驗(yàn)值與計(jì)算值差值占試驗(yàn)值的百分比。

        表2 峰值頻率及聲壓級(jí)計(jì)算值與試驗(yàn)值對(duì)比Tab.1 Frequency and SPL results of experiment and simulation at two peaks

        由表2 可知,試驗(yàn)和計(jì)算得到的2 個(gè)峰值頻率的頻率值和聲壓級(jí)值基本一致。文獻(xiàn)[15]中的仿真結(jié)果第1 峰值聲壓級(jí)偏低10 dB 以上,沒(méi)有明顯的第2 峰值。本文計(jì)算得到的峰值頻率比試驗(yàn)值略為偏小,誤差小于3%,而第1 峰值聲壓級(jí)偏大3.5 dB,約為2.8%。對(duì)于復(fù)雜的空腔湍流致聲問(wèn)題,計(jì)算精度已可為工程應(yīng)用接受。

        4 結(jié) 語(yǔ)

        為對(duì)艦船管路系統(tǒng)水動(dòng)力噪聲進(jìn)行計(jì)算,本文在前期工作的基礎(chǔ)上,采用低馬赫數(shù)下的三維類(lèi)閥空腔模型試驗(yàn)數(shù)據(jù),對(duì)大渦模擬和Lighthill 聲類(lèi)比理論的流激噪聲混合計(jì)算方法進(jìn)行了驗(yàn)證。

        首先采用大渦模擬對(duì)三維類(lèi)閥空腔湍流流場(chǎng)進(jìn)行了模擬。然后將流場(chǎng)信息作為聲源插值到聲學(xué)網(wǎng)格并做傅里葉變換,通過(guò)ACTRAN 中的Lighthill 聲類(lèi)比理論對(duì)空腔流激噪聲進(jìn)行了計(jì)算。最后對(duì)三維類(lèi)閥空腔模型的流動(dòng)特性及聲學(xué)特性進(jìn)行了一定的分析,并將空腔內(nèi)測(cè)點(diǎn)聲壓值與對(duì)應(yīng)的低馬赫數(shù)下試驗(yàn)結(jié)果進(jìn)行了對(duì)比。仿真結(jié)果與試驗(yàn)結(jié)果在測(cè)量頻率范圍內(nèi)量級(jí)相近,趨勢(shì)相同,峰值頻率和聲壓級(jí)大小吻合良好,驗(yàn)證了大渦模擬和Lighthill 聲類(lèi)比混合方法計(jì)算流激噪聲可行,結(jié)果可靠,可用于低馬赫數(shù)下氣體流激噪聲及水動(dòng)力噪聲的預(yù)測(cè)。

        [1]張楠,沈泓萃,姚惠之,等.孔穴流激噪聲的計(jì)算與驗(yàn)證研究[J].船舶力學(xué),2008,12(5):799-805.

        ZHANG Nan,SHEN Hong-cui,YAO Hui-zhi,et al.Validation and calculation of flow-induced noise of cavity[J].Journal of Ship Mechanics,2008,12(5):799-805.

        [2]張楠,沈泓萃,朱錫清,等.三維孔腔流激噪聲的大渦模擬與聲學(xué)類(lèi)比預(yù)報(bào)與驗(yàn)證研究[J].船舶力學(xué),2010,14(1-2):181-190.

        ZHANG Nan,SHEN Hong-cui,ZHU Xi-qing,et al.Validation and prediction of flow induced noise of 3-dimensional cavity with large eddy simulation and acoustic analogy[J].Journal of Ship Mechanics,2010,14(1-2):181-190.

        [3]耿冬寒,劉正先.大渦模擬-Lighthill 等效聲源法的空腔水動(dòng)噪聲預(yù)測(cè)[J].哈爾濱工程大學(xué)學(xué)報(bào),2010,31(2):182-187.

        GENG Dong-han,LIU Zheng-xian.Predicting cavity hydrodynamic noise using a hybrid large eddy simulation-Lighthill′s equivalent acoustic source method[J].Journal of Harbin Engineering University,2010,31(2):182-187.

        [4]MOON Y J,SEO J H,BAE Y M,et al.A hybrid prediction method for low-subsonic turbulent flow noise[J].Comput.Fluids,2010,39(7):1125-1135.

        [5]JI M,WANG M.Sound generation by turbulent boundarylayer flow over small steps[J].J.Fluid Mech.,2010,654:161-193.

        [6]劉敏,張寧,周友明,等.表面開(kāi)孔水下回轉(zhuǎn)體流噪聲數(shù)值模擬[J].艦船科學(xué)技術(shù),2011,33(11):34-38.

        LIU Min,ZHANG Ning,ZHOU You-ming,et al.Flowinduced noise simulation of an underwater fenestrated revolution body[J].Ship Science and Technology,2011,33(11):34-38.

        [7]張?jiān)?,傅慧萍,繆?guó)平.基于大渦模擬的開(kāi)孔潛體流噪聲數(shù)值模擬[J].上海交通大學(xué)學(xué)報(bào),2011,45(12):1868-1873.

        ZHANG Yun,F(xiàn)U Hui-ping,MIAO Guo-ping.LES-based numerical simulation of flow noise for submerged body with cavities[J].Journal of Shanghai Jiaotong University,2011,45(12):1868-1873.

        [8]江文成,張懷新,孟堃宇.基于邊界元理論求解水下潛艇流噪聲的研究[J].水動(dòng)力學(xué)研究與進(jìn)展,2013,28(4):453-459.

        JIANG Wen-cheng,ZHANG Huai-xin,MENG Kun-yu.Research on the flow noise of underwater submarine based on boundary element method[J].Chinese Journal of Hydrodynamics,2013,28(4):453-459.

        [9]江文成.潛艇流噪聲與流固耦合作用下流激噪聲的數(shù)值模擬[D].上海:上海交通大學(xué),2013.

        JIANG Wen-cheng.Numerical simulation on flow noise and flow-excited noise of submarine under fluid structure interaction[D].Shanghai:Shanghai Jiaotong University,2013.

        [10]袁壽其,司喬瑞,薛菲,等.離心泵蝸殼內(nèi)部流動(dòng)誘導(dǎo)噪聲的數(shù)值計(jì)算[J].排灌機(jī)械工程學(xué)報(bào),2011,29(2):93-98.

        YUAN Shou-qi,SI Qiao-rui,XUE Fei,et al.Numerical calculation of internal flow-induced noise in centrifugal pump volute[J].Journal of Drainage and Irrigation Machinery Engineering,2011,29(2):93-98.

        [11]趙威,柯兵,李奇,等.水泵內(nèi)部流動(dòng)誘導(dǎo)噪聲的預(yù)測(cè)[J].中國(guó)水運(yùn),2013,13(2):126-128.

        ZHAO Wei,KE Bing,LI Qi,et al.Flow-induced noise simulation in centrifugal pump[J].China Water Transport,2013,13(2):126-128.

        [12]LIU J M,ZHANG T,ZHANG Y O.Numerical study on flow-induced noise for a steam stop- valve using large eddy simulation[J].J.Marine Sci.Appl.,2013,12(3):351-360.

        [13]張碩,孫江龍,呂續(xù)艦.基于AMESim和Fluent 的載流管道系統(tǒng)降噪研究[J].艦船科學(xué)技術(shù),2013,35(8):91-95.

        ZHANG Shuo,SUN Jiang-long,LV Xu-jian.Research on noise reduction in the pipe system based on AMESim and Fluent[J].Ship Science and Technology,2013,35(8):91-95.

        [14]郭濤.管路的流致振動(dòng)及噪聲研究[D].武漢:華中科技大學(xué),2012.

        GUO Tao.Research on flow induced vibration and sound of pipe[D].Wuhan:Huazhong University of Science and Technology,2012.

        [15]LAFON P,CAILLAUD S,DEVOS J P,et al.Aeroacoustical coupling in a ducted shallow cavity and fluid/structure effects on a steam line[J].J.Fluids Struct.,2003,18:695-713.

        猜你喜歡
        模型
        一半模型
        一種去中心化的域名服務(wù)本地化模型
        適用于BDS-3 PPP的隨機(jī)模型
        提煉模型 突破難點(diǎn)
        函數(shù)模型及應(yīng)用
        p150Glued在帕金森病模型中的表達(dá)及分布
        函數(shù)模型及應(yīng)用
        重要模型『一線(xiàn)三等角』
        重尾非線(xiàn)性自回歸模型自加權(quán)M-估計(jì)的漸近分布
        3D打印中的模型分割與打包
        亚洲av毛片在线网站| 国产va精品免费观看| 日本肥老熟妇在线观看| 日本在线观看一二三区| 亚洲av无码专区亚洲av网站| 小12箩利洗澡无码视频网站| 久久久亚洲精品午夜福利| 一区二区精品天堂亚洲av| 亚洲av中文无码乱人伦在线咪咕 | 国产av无码专区亚洲av| 亚洲AV永久无码精品一区二国 | 久久久亚洲av成人网站| 久久亚洲国产精品成人av秋霞| 一区视频在线观看免费播放.| 国产片在线一区二区三区| 亚洲av无码专区在线观看下载| 乱码午夜-极国产极内射| 日本一区免费喷水| 一区二区三区视频亚洲| 秘书边打电话边被躁bd视频| 在线观看91精品国产免费免费| 国产精品国产三级国av在线观看| 精品人妻中文av一区二区三区| 国产男女无遮挡猛进猛出| 香蕉成人啪国产精品视频综合网 | 久久半精品国产99精品国产 | 国产女厕偷窥系列在线视频| 另类专区欧美在线亚洲免费| 国产三级精品三级在专区中文 | 特级黄色大片性久久久| 精品国产精品国产偷麻豆| 7878成人国产在线观看| av免费看网站在线观看| 国产白浆一区二区三区性色| 人禽伦免费交视频播放| 免费一级欧美大片久久网| 中文字幕亚洲入口久久| 亚洲av美国av产亚洲av图片| 国产成人一区二区三中文| 国产av精品一区二区三区不卡| 国产老熟妇精品观看|