999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

統計能量法計算水下圓柱殼輻射噪聲準確性的驗證與分析

2017-08-05 01:37:04張愷紀剛周其斗劉文璽
中國艦船研究 2017年4期
關鍵詞:模態

張愷,紀剛,周其斗,劉文璽

海軍工程大學艦船工程系,湖北武漢430033

統計能量法計算水下圓柱殼輻射噪聲準確性的驗證與分析

張愷,紀剛,周其斗,劉文璽

海軍工程大學艦船工程系,湖北武漢430033

[目的]統計能量法(SEA)是解決結構高頻振動與聲輻射問題的有效方法,但是該方法通常假定流體為“輕質流體”,在分析水中結構時其計算結果可能不準確。[方法]分別運用SEA方法和有限元耦合邊界元法(FEM/BEM)計算水下圓柱殼模型的輻射聲壓級,以驗證SEA預報水下圓柱殼輻射噪聲的準確性。運用SEA計算不同的子系統劃分方式和不同的內損耗因子誤差時圓柱殼的輻射聲壓級,分析影響SEA計算結果準確性的因素。[結果]在400 Hz以下時SEA和FEM/BEM的計算結果相差很大,在400 Hz以上基本一致;不同子系統的劃分方式造成的誤差在5 dB左右;內損耗因子誤差100%時造成的誤差在2~3 dB。[結論]通過研究發現,在模態密度足夠時可以使用SEA計算水下圓柱殼的輻射噪聲,對于低頻沿周向劃分子系統不可靠,可能導致計算結果不準確;對于高頻沿周向劃分子系統比沿軸向劃分子系統得出的計算結果更準確;對于能量高的子系統其內損耗因子誤差對仿真結果影響更大,應采取更精確的方式確定其內損耗因子。研究結果對于運用SEA研究水下結構振動與噪聲問題有一定參考價值。

統計能量法;聲輻射;子系統;內損耗因子

0 引 言

潛艇的輻射噪聲是潛艇聲學研究的重要內容。在低頻,有限元法(FEM)和邊界元法(BEM)是解決復雜結構振動與聲輻射問題的有效方法,但是在高頻由于模態密度的增大,有限元法難以求解。統計能量法(SEA)是解決結構高頻振動與聲輻射問題的有效方法。驗證與分析SEA計算水中結構輻射噪聲的準確性,對于正確使用SEA研究水中結構高頻噪聲有著重要意義。

20世紀60年代提出SEA方法后,經過50多年的不斷的完善與發展,出現了AutoSEA,SEAM等成熟的商業軟件,并已成功地應用于車輛、建筑和航空航天等領域。對于SEA方法用于水下結構振動與聲輻射問題也有一定的研究。彭臨慧等[1]用SEA方法分析了湍流邊界層脈動壓力激勵下水下結構的自噪聲,根據工程條件進行了簡化,得到相應的自噪聲工程估算關系。童宗鵬等[2]用SEA方法對水下航行器結構在寬頻帶范圍內的水下振動和聲輻射規律進行了分析與研究,其將模型簡化為圓錐殼和圓柱殼,計算了模型的噪聲級。

FEM/BEM方法是目前計算水下結構振動與輻射噪聲比較成熟的方法。本文將以FEM/BEM方法的計算結果作為參照,將SEA的計算結果與FEM/BEM的結果對比,以驗證SEA預報水下圓柱殼輻射噪聲的準確性,同時研究不同的子系統劃分方式和不同的內損耗因子誤差對SEA仿真結果的影響。

1 基本原理

1.1 SEA計算輻射聲壓的原理

根據能量守恒,穩態時輸入子系統的功率等于子系統自身損耗的功率和流入其他子系統的功率之和。假設一個結構有N個子系統,則子系統i的功率平衡方程為

式中:Pi,in為子系統的輸入功率;i為子系統的能量關于時間的導數,對于穩態系統i=0;Pid為子系統的內損耗功率;Pij為子系統i流入子系統j的功率。

將式(1)寫成矩陣的形式

左側矩陣為總損耗因子矩陣,它與結構本身的動力學特性有關,可以通過實驗或半經驗公式得到;左側列向量為子系統的能量;右側列向量為子系統輸入功率,可以由激振力和輸入阻抗得到。求解該矩陣方程得到每個子系統的能量。由子系統的能量得到子系統相關的動力學參數,如位移、速度、加速度、應力和聲壓等。

對于流體負載下的結構振動,給結構附加一個與頻率和流體屬性有關的“附加質量”和一個與結構在“輕質流體”下振動時向半無限流域輻射聲功率有關的“附加阻尼”,求出流體負載下結構各子系統的能量Ei。

子系統的能量與響應速度的關系為

式中:Mi為子系統的質量;νirms是子系統的響應速度的均方根。

具有相同面積、相同均方速度的剛性平板輻射功率 p′rad為

式中:A為平板面積;ρa為流體密度;Ca為流體聲速。

輻射比的定義為

式中:σrad為結構的輻射比;prad為結構的輻射聲功率。根據式(4)和式(5)已知結構振動的均方速度和輻射比可以求出輻射聲功率 prad。

在遠場,聲波以球面波的形式傳播,根據輻射聲功率計算測點所在波振面的平均聲強,得到測點處的聲壓級。假設該聲壓級為聲場均勻分布時的聲壓級。

1.2 驗證原理

FEM/BEM是計算水下結構振動與輻射噪聲比較成熟的方法。在計算比較簡單的結構時,FEM/BEM方法可以計算到較高頻率。在計算圓柱殼多個測點處的輻射聲壓級時,取平均聲壓級可以近似得到聲場均勻分布時的聲壓級。當FEM/BEM兩種仿真方法分析模型的阻尼相同時,將SEA方法的計算結果與其計算的平均聲壓級進行對比,則可在一定程度上驗證SEA計算水下圓柱殼輻射噪聲的準確性。

本文的FEM/BEM方法,采用從結構到流體的解耦方式,通過結構濕表面的法向速度vn將流體聲壓與結構振動聯系在一起。

通過結構的有限元方程式(6)和流體的Helmholtz方程式(7),求解得到濕表面法向速度vn和任意聲場點聲壓 p(x)。

式中:M,C,K分別為結構的質量、阻尼、剛度矩陣;v為結構的響應速度;f為激振力;fp為流體附加作用力;G(r)為格林函數,對于自由空間G(r)基本解為 e-ikr/4πr;p(x)為聲場點聲壓,p(x′)為濕表面聲壓;pI為入射聲壓。

2 SEA預報水下結構輻射噪聲的準確性

分別使用SEA仿真軟件AutoSEA2和FEM/BEM軟件,來計算水下圓柱殼在距圓柱殼中心6.09 m處的輻射聲壓級。將FEM/BEM方法的計算結果取空間平均后與SEA能量方法的計算結果進行對比。

2.1 SEA計算圓柱殼輻射聲壓級

應用SEA軟件AutoSEA2建立圓柱殼的SEA模型,如圖1所示。表1為模型的結構參數。模型分為7個子系統,包括3個殼體子系統、2個端蓋對應的平板子系統、內部聲腔子系統以及測點處的半無限流域子系統。激振力作用于艙壁,方向沿徑向向內,大小為4.448 2 N。計算頻帶為10~2 000 Hz的1/3倍頻程。模型參考了文獻[3]的實驗模型。

圖1 圓柱殼SEA模型Fig.1 SEA model of cylindrical shell

表1 模型結構參數Table 1 Parameters of model structure

通過仿真軟件AutoSEA2計算測點處的輻射聲壓級。AutoSEA2可以自動求解子系統的模態密度、耦合損耗因子等參數,需要用戶確定的是子系統的內損耗因子。子系統的內損耗因子為ηi=ηis+ηir+ηib,其中,ηis為結構損耗因子,是由結構的材料決定的,結構損耗因子取有限元分析中的結構阻尼系數0.06;ηir為聲輻射損耗因子,反映了子系統由于聲輻射損耗的能量,在AutoSEA2中根據輻射效率可以自行計算出ηir;ηib為子系統的邊界連接阻尼損耗因子,子系統間連接為剛性時可以忽略[3]。

在子系統的Auxiliary選項設置流體對結構的負載,將圓柱殼外域設為水,軟件會自動附加一個與頻率有關的附加質量。

AutoSEA2計算圓柱殼輻射聲壓級的方法為:首先通過結構的功率平衡方程,計算出流體作用下結構各子系統的均方法向速度;之后通過各子系統的輻射比,求出各子系統的輻射聲功率。因為SEA將子系統的各參量取平均,故可以假設各個子系統輻射的聲場不相干,根據近場、遠場條件分別采用平面波、球面波模型計算波振面上的平均聲壓級。

2.2 FEM/BEM計算圓柱殼輻射聲壓級

圖2為本文模型的有限元模型。左側有限元模型單元數為1 428個,單元尺度為0.1 m;右側為聲場網格,在6.09 m處沿周向等間隔劃分為72個單元。徑向激振力作用在圓柱殼上,F=4.448 2 cos(2πfot)N(式中,fo為頻率,t為時間),計算頻率為10~2 000 Hz,頻率點間隔10 Hz。當頻率為2 000 Hz時鋼中的結構波長為0.725 m,因此單元尺度滿足1/6波長,有限元的計算結果是可靠的。

圖2 圓柱殼FEM模型Fig.2 FEM model of cylindrical shell

采用基于Nastran開發的FEM/BEM軟件求解式(1)、式(5)和式(6)得到測點處的聲壓級。因為激振力是非對稱的,所以圓柱殼的徑向輻射聲場有指向性。在穿過圓柱殼中心且垂直于軸線的平面內,距離圓柱殼中心6.096 m處平均設置72個聲場點,計算72個聲場點的聲壓級 Li。取72個聲場點輻射聲壓級的平均聲壓級 Lˉ,近似獲得圓柱殼6.096 m處的平均輻射聲壓級。式(8)和式(9)是平均聲壓級 Lˉ的計算方法,其中參考聲壓 pref=1×10-6Pa。

2.3 比較兩種方法的計算結果

FEM/BEM和SEA這2種方法的計算結果如圖3所示。圖3(a)是兩種方法計算的輻射聲壓級(dB),圖3(b)是對FEM/BEM按1/3倍頻程取平均之后的結果。

圖3 圓柱殼的計算輻射聲壓級Fig.3 Radiated sound pressure level of cylindrical shell caculated by SEA and FEM/BEM methods

由圖3(a)可知,FEM/BEM的計算結果起伏較大,有許多波峰、波谷,波峰對應頻率是該模型的“譜峰頻率”;SEA計算的結果比較平滑,因為SEA的計算頻帶是10~2 000 Hz的1/3倍頻程,得到的是每個頻帶上的平均聲壓級。

圖3(b)中,將FEM/BEM的計算結果按照10~2 000 Hz的1/3倍頻程取頻帶上的平均值。400 Hz之前,兩種方法的計算結果相差很大。因為在低頻時,SEA的子系統模態密度過低,其計算結果不可靠。400 Hz之后,兩種方法計算結果吻合很好。通過對比兩種方法的結果,可以發現模態密度足夠大時,SEA方法可以較準確地預報水下圓柱殼的輻射噪聲。

3 子系統的劃分方式對SEA仿真結果的影響

合理劃分子系統對SEA計算結果的準確性有重要影響。同一子系統的模態要滿足模態相似和模態均分的假設,并且保證有足夠高的模態密度。如果不能滿足這些假設,計算結果可能會出現較大誤差。通過研究不同劃分方式下圓柱殼輻射聲壓級的結果,研究子系統劃分方式對仿真結果的影響。

3.1 子系統的劃分原則

子系統是指相似的共振模態組成的振形群,同一振形群的模態具有相同的振動能量和耦合損耗因子。子系統的劃分首先要滿足模態相似原則,其次要滿足模態密度足夠高以保證計算的可靠性。對于復雜模型,通常按照自然幾何邊界條件、動力學邊界條件、材料介質特性等建立子系統模型[4]。AutoSEA軟件可以根據幾何子系統自動劃分振形群子系統,因此只需要將幾何形狀劃分成合適的子系統即可。

3.2 按不同方式劃分子系統

采用第2節中的模型,通常將2個蓋板劃分為2個平板子系統,將內部聲腔劃分為1個聲場子系統。而圓柱殼則有多種劃分方式,可以將圓柱殼沿周向劃分為不同數量的圓柱殼子系統,或者沿周向劃分為不同數量的單曲率板子系統。如圖4所示,將圓柱殼分別沿軸向劃分1,3,5個子系統,沿周向劃分為2,4個子系統。分別計算圓柱殼在6.096 m處的輻射聲壓級。

3.3 對比不同方式劃分子系統時計算的輻射聲壓級

圖5所示為沿軸向劃分為1,3,5個子系統,沿周向劃分為2,4個子系統時得到的輻射聲壓級。圖6所示為不同子系統劃分方式下的圓柱殼子系統的模態密度。

圖4 不同的子系統劃分方式Fig.4 Schemes of different subsystems

圖5 不同的子系統劃分方式計算的輻射聲壓Fig.5 Radiated sound pressure levels caculated by different schemes

由圖5中可知,不同的子系統劃分方式對計算結果的準確性有較大影響,造成的誤差在5 dB左右,在低頻時計算結果相差更大。合理的劃分方式對計算結果的準確性非常重要。

采用SEA方法計算時,子系統的模態數N至少大于5,子系統的大小直接影響子系統的模態數。由圖6中可知,模態數主要與子系統的數量有關,與軸向和周向的劃分方式基本無關。隨著子系統數目增多,模態數減少,N>5對應的頻率更高。但是當頻率到達400 Hz左右時,模態密度迅速增大,遠遠超過5。因此,在計算較低頻率時,應當將子系統劃分得盡量大。在計算較高頻率時,子系統的大小基本無影響。

圖6 不同子系統劃分方式的模態數Fig.6 Mode orders regarding different schemes of subsystem

由圖5可知,沿軸向劃分子系統時,隨著子系統數目的增加計算結果增大,不同的劃分方式造成的誤差在5 dB左右。沿周向劃分子系統時,在2 000 Hz以下,計算結果相差巨大,計算結果不可靠;在2 000 Hz以后的誤差較小,在2 dB以內。因此在計算較低頻率時,沿周向劃分子系統是不合理的,這會導致計算結果不可靠;在高頻時沿周向劃分子系統,誤差相對較小,計算結果更準確。

4 內損耗因子誤差對SEA仿真結果的影響

內損耗因子(ILF)是由系統阻尼特性決定的反映系統能量損耗的參數。確定內損耗因子的方法主要是經驗公式和實驗等方法,這些方法都很難得到精確的內損耗因子真實值,總是存在一定的誤差。內損耗因子的誤差不僅會引起子系統本身能量損耗的誤差,還會導致子系統間流動能量的誤差[5]。通過改變所有子系統的內損耗因子和改變單個子系統的內損耗因子,以研究內損耗因子誤差對仿真結果的影響。計算模型采用第2節中建立的模型。

4.1 所有子系統內損耗因子存在相同誤差

計算所有子系統內損耗因子為0.06以及內損耗因子誤差 (Δη)分別為10%,20%,50%,100%時圓柱殼輻射聲壓級的誤差,計算結果如圖7所示。

圖7 內損耗因子誤差對輻射聲壓的影響Fig.7 The influence of different ILF error on sound pressure level of cylindrical shell

由圖7可知,在50%誤差范圍以內時,內損耗因子的誤差對輻射噪聲的影響基本呈線性規律。每10%的誤差大約導致0.3 dB的誤差,100%誤差時輻射噪聲誤差在2~3 dB之間。

4.2 單個子系統內損耗因子存在誤差

分別計算不同的子系統內損耗因子為100%時的輻射聲壓級誤差,并且計算出內損耗因子為0.01時各個子系統的能量,以研究單個子系統內損耗因子誤差的影響。

圖8所示為受激勵的圓柱殼子系統1(圓柱殼1)、沒有激勵的圓柱殼子系統2(圓柱殼2)以及平板子系統內損耗因子誤差分別為100%時的輻射聲壓級誤差。圖9所示為內損耗因子為0.01時的各子系統的能量。

圖8 Δη=100%時3個子系統輻射聲壓級誤差Fig.8 The errors of radiation sound pressure level of three subsystems when Δη=100%

圖9 η=0.06時3個子系統的能量Fig.9 The energy of three subsystem whenη=0.06

由圖8可以看出,不同子系統的內損耗因子敏感性不同。圓柱殼1內損耗因子誤差對輻射噪聲的影響最大,接近模型整體內損耗因子改變導致的誤差;圓柱殼2和艙壁內損耗因子的誤差對輻射噪聲影響很小。圖9中,3個子系統的能量差距很大。結合之前的分析可以發現,子系統內損耗因子敏感度與子系統的能量是有一定關系的。能量大的子系統,內損耗因子誤差對計算結果的影響更大。

通過上述研究可以發現,對于本文模型,內損耗因子誤差對計算的輻射聲壓結果影響較小;能量大的子系統內損耗因子誤差對結果的影響更大一些。因此,在確定子系統內損耗因子時,對于激勵比較大,振動比較劇烈的子系統,應當采取更加精確的方式。

5 結 論

本文分別采用SEA和FEM/BEM兩種方法計算了水中圓柱殼的輻射聲壓級,驗證統計能量法計算水中結構輻射噪聲的準確性。研究了子系統劃分方式以及內損耗因子誤差對仿真結果造成的影響。針對本文圓柱殼模型,通過研究得到如下結論:

1)低頻時,SEA方法計算結果與FEM/BEM方法的結果相差較大。高頻時,當模態密度足夠時,兩種方法計算的結果吻合良好。在高頻用SEA方法研究水下圓柱殼的輻射噪聲是可靠的。

2)低頻時,劃分過多的子系統會導致計算結果不可靠,高頻時子系統的數目對結果影響不大;對于圓柱殼,低頻時沿周向劃分子系統計算的結果不可靠,高頻時沿周向劃分子系統誤差相對沿軸向劃分子系統較小。

3)內損耗因子誤差主要引起輻射聲壓級曲線的上、下平移。誤差在50%以內時對計算結果的影響基本呈線性,本算例中內損耗因子每10%誤差造成輻射聲壓級0.3 dB的誤差。

4)不同子系統的內損耗因子誤差敏感度不同,能量大的子系統,內損耗因子誤差對仿真結果的影響更大。因此,在確定內損耗因子時,對激勵較大,振動更劇烈的子系統應采用更加精確的方式。

[1]彭臨慧,陸建輝.水下結構物自噪聲統計能量分析[J],中國造船,2001,42(3):58-63.PENG L H,LU J H.Statistical energy analysis for self noise of underwater structure[J].Shipbuilding of China,2001,42(3):58-63(in Chinese).

[2]童宗鵬,王國治,張志誼,等.水下航行器聲振特性的統計能量法研究[J].噪聲與振動控制,2005,25(1):29-32.TONG Z P,WANG G Z,ZHANG Z Y,et al.Re?search on vibration and acoustic radiation of the subma?rinestructure with SEA method[J].Noise and Vibra?tion Control,2005,25(1):29-32(in Chinese).

[3]CHENL H,SCHWEIKERTD G.Sound radiation from an arbitrary body[J].The Journal of the Acoustical Society of America,1963,35(10):1626-1632.

[4]姚德源,王其政.統計能量分析原理及其應用[M].北京:北京理工大學出版社,1995.

[5]寧瑋,張景繪,王珺.統計能量分析法中參數靈敏度分析[J].系統仿真學報,2009,21(17):5366-5370.NING W,ZHANG J H,WANG J.Sensitivity analysis of parameters in statistical energy analysis method[J].JournalofSystem Simulation,2009, 21(17) :5366-5370(in Chinese).

Accuracy verification and analysis of SEA method for calculating radiation noise pressure of submerged cylindrical shell

ZHANG Kai,JI Gang,ZHOU Qidou,LIU Wenxi
Department of Naval Architecture Engineering,Naval University of Engineering,Wuhan 430033,China

Statistical Energy Analysis(SEA)is an effective method for solving high frequency structural vibration and acoustic radiation problems.When we use it to analyze submerged structures,it is necessary to consider the actions of fluid as'heavy fluid'relative to structures,which differs from when it is used in the air.The simple model of a submerged cylindrical shell is used to calculate at a higher frequency using FEM/BEM.The SEA and FEM method are then used to calculate the radiation sound pressure level,verifying the accuracy of the SEA prediction for submerged structures.The classified method of subsystems and the effect of the error of the internal loss factor on the accuracy of the results are explored.The calculated results of SEA and FEM/BEM are very different below 400 Hz,and basically the same above 400 Hz.The error caused by the division of different subsystems is about 5 dB.The error in the calculation results caused by the error of the internal loss factor is 2-3 dB.It is possible to use SEA to calculate the radiated noise of an underwater cylindrical shell when the modal density is high enough.For the cylindrical shell,dividing the subsystems along the circumference is not reliable at a low frequency,as it may lead to inaccurate calculation results.At a high frequency,it is more accurate to divide the subsystems along the circumference than the axle.For subsystems with high energy,the internal loss factor has a greater effect on the simulation results,so a more accurate way should be taken to determine the internal loss factor of subsystems with high energy.

Statistical Energy Analysis(SEA);acoustic radiation;subsystem;internal loss factor

TB53;U661.44

A

10.3969/j.issn.1673-3185.2017.04.014

http://kns.cnki.net/kcms/detail/42.1755.TJ.20170727.1028.030.html期刊網址:www.ship-research.com

張愷,紀剛,周其斗,等.統計能量法計算水下圓柱殼輻射噪聲準確性的驗證與分析[J].中國艦船研究,2017,12(4):89-94.

ZHANG K,JI G,ZHOU Q D,et al.Accuracy verification and analysis of SEA method for calculating radiation noise pressure of submerged cylindrical shel[lJ].Chinese Journal of Ship Research,2017,12(4):89-94.

2016-12-26< class="emphasis_bold">網絡出版時間:

時間:2017-7-27 10:28

張愷,男,1992年生,碩士生。研究方向:潛艇聲隱身技術。E-mail:694021510@qq.com

紀剛(通信作者),男,1975年生,博士,副教授。研究方向:潛艇結構聲學特性。

E-mail:jigang_-_@sohu.com

猜你喜歡
模態
基于BERT-VGG16的多模態情感分析模型
跨模態通信理論及關鍵技術初探
一種新的基于模態信息的梁結構損傷識別方法
工程與建設(2019年1期)2019-09-03 01:12:12
多跨彈性支撐Timoshenko梁的模態分析
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
利用源強聲輻射模態識別噪聲源
日版《午夜兇鈴》多模態隱喻的認知研究
電影新作(2014年1期)2014-02-27 09:07:36
主站蜘蛛池模板: 91av成人日本不卡三区| 欧美α片免费观看| 久久黄色一级视频| 人妻精品久久无码区| 欧美自慰一级看片免费| 免费国产好深啊好涨好硬视频| 最新午夜男女福利片视频| 亚洲国产精品人久久电影| 欧美日韩91| 亚洲综合狠狠| 亚洲h视频在线| 亚洲国产精品一区二区高清无码久久| 午夜福利免费视频| 国产导航在线| 日韩少妇激情一区二区| 国产亚洲欧美日韩在线一区二区三区| 欧美精品不卡| 成人综合在线观看| 中文字幕调教一区二区视频| 国产美女精品人人做人人爽| 亚洲中文字幕国产av| 国产精品女在线观看| 国产精品免费电影| 青青草久久伊人| 成人国产精品2021| 91成人在线观看| www.精品视频| 视频一区视频二区中文精品| 国产一级毛片yw| 国产成人91精品| 99无码熟妇丰满人妻啪啪| 成人在线第一页| 久久无码av三级| 青青久久91| 亚洲无码精彩视频在线观看| 国产精品第一区| 秋霞国产在线| 亚洲精品无码av中文字幕| 欧洲av毛片| 免费xxxxx在线观看网站| 免费国产小视频在线观看| 91人人妻人人做人人爽男同| 好紧好深好大乳无码中文字幕| 亚洲国产日韩一区| av手机版在线播放| 日本人妻丰满熟妇区| 亚洲国产中文欧美在线人成大黄瓜| 久久中文电影| 精品少妇人妻无码久久| 亚洲成人黄色网址| 国产一区二区福利| a国产精品| 尤物精品视频一区二区三区| 不卡网亚洲无码| 一级毛片在线直接观看| 国内精品九九久久久精品| 伊人激情久久综合中文字幕| 国产精品私拍在线爆乳| 欧美一级大片在线观看| 91九色国产porny| 久久99久久无码毛片一区二区| 激情在线网| 中文字幕第4页| 亚洲成人www| 99re视频在线| 美女一级免费毛片| 亚洲人人视频| 国产真实二区一区在线亚洲| 亚洲最大综合网| 国产午夜精品鲁丝片| 国产99视频在线| 老司国产精品视频| 91www在线观看| 激情乱人伦| 91精品人妻互换| 国产对白刺激真实精品91| 在线观看免费人成视频色快速| 久久久久国产一级毛片高清板| 尤物精品国产福利网站| 97超级碰碰碰碰精品| 欧洲极品无码一区二区三区| 67194成是人免费无码|