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

聲-固耦合系統高頻局部能量響應的FEM-SmEdA分析

2019-10-21 06:17:50王攀陳強李彥斌
振動工程學報 2019年4期
關鍵詞:有限元法

王攀 陳強 李彥斌

摘要: 針對高頻激勵下聲-固耦合系統的局部能量響應預示問題,采用有限元(FEM)和統計模態能量分布分析(SmEdA)相結合的方法預示復雜聲-固耦合系統局部能量響應。首先,以平板/聲腔耦合系統為研究對象,驗證FEM-SmEdA方法的準確性;隨后,為了驗證該方法的適用性,分別將該方法應用于湍流邊界層(TBL)激勵下平板/聲腔耦合系統局部響應預示及加筋艙段/聲腔耦合系統局部響應預示。利用雙模態方程法(DMF)對基于FEM-SmEdA法的TBL激勵作用下結構局部響應預示結果進行驗證,基于FEM-SmEdA方法揭示TBL激勵作用下艙段結構和聲腔的能量分布隨頻率的變化規律,并研究加筋對聲-固耦合系統能量分布的影響。結果表明:FEM-SmEdA方法能夠準確地預示高頻環境下聲-固耦合系統的局部能量響應;該方法適用于TBL激勵下結構局部響應預示;隨著分析頻率的升高,艙段和聲腔的能量分布更加均勻;加筋后結構的能量分布更加均勻,加筋對聲腔能量響應的抑制更加明顯。

關鍵詞: 聲-固耦合系統; 統計模態能量分布分析; 模態能量; 能量分布; 有限元法

中圖分類號: O421+.6; V414.4 文獻標志碼: A 文章編號: 1004-4523(2019)04-0590-12

DOI:10.16385/j.cnki.issn.1004-4523.2019.04.005

引 言

聲-固耦合問題廣泛存在于航空航天結構中,尤其是對于高頻激勵下的薄壁結構,結構和聲場之間極易產生聲-固耦合效應,引起結構的振動并改變聲場的分布,進而影響結構的可靠性和儀器設備的功能。因此,結構的振動水平及聲場的噪聲水平是結構設計階段必須考慮的重要指標,研究薄壁結構在高頻激勵下的聲-固耦合問題具有重要的理論意義和應用價值。

目前,獲取聲-固耦合響應的方法主要有試驗方法、理論方法和數值方法。試驗方法的結果真實可信,其不足之處在于耗費巨大、試驗周期較長且僅能實現有限實驗條件和工況。理論分析方法較難適用于復雜工程結構。近年來,隨著計算機性能的提高,數值方法成為一種有效的分析手段[1]。對于低頻聲-固耦合問題,通常采用有限元法(Finite Element Method, FEM)[2]和邊界元法等離散化方法預示結構響應。在結構單元尺寸滿足同一波長內包含6個單元的條件后[3],離散化方法能較為精確地預示聲-固耦合系統在低頻段的局部響應。但當分析頻率升高時,若要保證FEM有足夠的精度,須對單元進行細化,這將導致計算規模呈指數級增長。對于高頻聲-固耦合問題,通常采用統計能量分析(Statistical Energy Analysis, SEA)[4-5]預示結構的高頻能量響應。SEA將結構劃分為若干子系統,通過引入模態能量均一化和雨流(Rain On the Roof, ROR)載荷等基本假設,構建子系統間的功率流平衡方程,進而預示子系統在空間和分析頻帶內的平均能量響應[6]。但SEA的各項假設在工程應用中較難完全滿足,而且僅能得到空間平均能量響應,難以得到結構局部能量響應[7]。

為了克服SEA中的模態能量分析均分假設,Maxit等基于雙模態方程(Dual Modal Formulation, DMF)[8-9],推導了子系統模態能量平衡方程,提出統計模態能量分布分析(Statistical Modal Energy Distribution Analysis, SmEdA)[10-11]。SmEdA理論以分析頻帶內子系統每階模態的能量為基本未知量,考慮了同一子系統中不同模態能量差異較大的情形。基于SmEdA理論預示子系統的模態能量,僅需非耦合子系統在耦合邊界的模態信息即可獲取子系統的模態能量。對于聲-固耦合子系統而言,提高了計算效率,尤其是聲腔子系統在高頻段模態數密集的情況下。此外,SmEdA理論以每階模態為基本輸入參數,可考慮同一子系統內模態受到不均勻激勵的情況,從而可考慮局部激勵,克服SEA理論中雨流載荷的假設。Aragonès等[12]基于SmEdA理論對典型聲腔/平板/聲腔耦合系統進行了聲功率傳輸路徑分析,識別出與主要傳遞路徑相關的平板模態,并在此基礎上對平板加筋的降噪效果進行研究。Totaro等[13]基于SmEdA理論預示局部激勵作用下平板/聲腔耦合系統的動能和勢能密度的空間分布。Qiao等[14]利用耦合系統中耦合損耗因子與輻射效率的關系,基于SmEdA理論研究了圓柱殼和圓錐殼的輻射效率。各子系統的模態信息為SmEdA理論的重要輸入條件之一,用于計算各模態間的耦合損耗因子、各模態的阻尼耗散功率及各模態上的載荷輸入功率。對于簡單聲-固耦合系統的模態信息可由解析法給出,而對于復雜聲-固耦合系統,可借助FEM獲取各子系統的模態信息[15],進而結合SmEdA理論,預示各子系統的高頻局部能量響應。目前,統計模態能量分布分析法的應用尚限于ROR載荷和單點力載荷,ROR載荷即空間上處處不相關的面壓載荷。工程上常把噪聲載荷處理成完全不相關的面壓載荷[16],而航天器面臨的實際噪聲環境存在空間相關性。因此,有必要開展考慮載荷空間相關性的統計模態能量分布分析法研究。

本文采用FEM和SmEdA結合的方法,對于復雜聲-固耦合系統局部能量響應進行預示。首先,以平板/聲腔耦合系統為研究對象驗證FEM-SmEdA方法的準確性;隨后,分別將該方法應用于湍流邊界層(Turbulent Boundary Layer,TBL)激勵下平板/聲腔耦合系統局部能量響應預示及艙段/聲腔耦合系統的高頻局部響應預示,利用DMF法對基于FEM-SmEdA法TBL激勵下結構局部響應預示結果進行驗證,并基于FEM-SmEdA法揭示結構和聲腔能量隨著頻率的變化規律,研究加筋對聲-固耦合系統能量分布的影響。

由式(24)和(25)可知,基于SmEdA法預示結構的能量響應,不需要整體結構的剛度矩陣,僅需節點質量及節點剛度。因此,基于SmEdA可以預示結構的局部能量響應,克服了傳統的SEA方法僅能預示結構平均能量響應的不足。此外,由于結構的特定模態信息在分析頻帶內是不變的,根據方程(24)和(25)可知,能量分布的準確性取決于頻帶內基于SmEdA獲得的模態能量。

2 方法驗證

以平板/聲腔耦合系統為對象,基于FEM-SmEdA開展結構的局部能量預示研究,通過與有限元方法的對比驗證本文方法的準確性。平板和聲腔耦合示意圖如圖 2所示,平板和聲腔的幾何尺寸及材料屬性分別如表1和2所示。參考文獻[12,15,18],平板和聲腔的結構阻尼均取0.01。平板四周簡支,聲腔除與平板耦合面,其余各面為固定邊界。

2.1 模態能量分析

為保證有限元計算結果的準確性,以中心頻率為630 Hz的1/3倍頻程頻帶[561, 707] Hz作為分析頻段Δω開展研究。基于商用有限元軟件獲得平板及聲腔的模態信息,頻帶Δω內平板和聲腔分別有11階和20階模態振型。將平板和聲腔的模態計算結果代入式(2),可以求得如圖 3所示的平板與聲腔之間模態耦合損耗因子。

在平板坐標為(0.7, 0.2)的一點施加垂直于板面方向的單位力激勵F,由式(6)得到平板模態上載荷的輸入功率,由式(7)得到平板和聲腔在分析頻帶Δω內的模態能量,如圖 4所示。圖中dB表示兩個量的比值大小,對于輸入功率dB_P=10lg(P/ref_P),P表示輸入功率,ref_P表示參考值,且ref_P=1 W。對于振動能量dB_E=10lg(E/ref_E),E表示振動能量,ref_E表示參考值,且ref_E=1 J。由圖 4可知,在集中力載荷作用下,即子系統受到局部激勵的情況下,各階模態上的載荷輸入功率不同;平板和聲腔的各階模態能量均具有明顯的差異,平板的各階模態能量差值最大為20 dB,聲腔差值最大為25 dB,不滿足傳統SEA中的模態能量均一化假設;平板的各階模態能量變化趨勢與模態上的載荷輸入功率一致。

圖5為分別基于FEM-SmEdA和FEM預示得到分析頻帶Δω內平板的動能分布,圖 6為分別基于FEM-SmEdA和FEM預示得到分析頻帶Δω內平板的勢能分布。

從圖5及6對比可知,兩者存在一定差異的。分析主要存在以下兩點原因:(1)SmEdA理論基于僅考慮了共振傳輸,由于本文研究對象結構阻尼為0.01,對于小阻尼系統,子系統間的非共振傳輸可忽略,因此本文研究中只考慮分析頻帶內的共振模態。(2)該方法強調“系統高頻局部能量響應預示”,適用于子系統有較多模態數的分析頻段,對于足夠大的頻帶,滿足模態數足夠多,可以忽略同一子系統中不同模態間的能量,即忽略Tinp(Q,ω)和Vinp(Q,ω)的影響。隨著分析頻段中心頻率的升高,分析頻段的模態數增多,忽略同一子系統中不同模態間的能量所造成的誤差減小,基于FEM-SmEdA預示的局部響應結果更準確。

3 湍流邊界層激勵下結構局部能量響應預示 ?正如已經在文獻[19]指出的,SmEdA理論以每階模態為基本輸入參數,可考慮同一子系統內模態受到不均勻激勵的情況,克服SEA理論中雨流載荷的假設。為了證明FEM-SmEdA也適用于預測復雜環境激勵下結構的局部響應預示,在平板/聲腔耦合系統中考慮TBL激勵。如1.2節已經指出的,平板和聲腔的特定模態信息是不變的,因此能量分布的精度由基于FEM-SmEdA法在頻帶中獲得的模態能量決定。由于TBL激勵下系統的局部響應難以用有限元方法預測,因此在第3.2節中用DMF方法驗證了TBL激勵下子系統模態能量預示的準確性。

3.1 湍流邊界層激勵

3.2 模態能量分析

在第2節所示的平板/聲腔耦合系統中施加TBL激勵,給出了U=200 m/s的來流速度,分別基于FEM-SmEdA和DMF預示中心頻率為630 Hz的1/3倍頻程頻帶內板和聲腔的模態能量,如圖 7和8所示。平板和聲腔的模態能量具有明顯的離散性,平板模態能量隨模態變化的能量差值可達20 dB,聲腔模態能量隨模態變化的能量差值可達30 dB,表明在考慮TBL激勵的情況下,能量均分的假設是不合理的。結果表明,基于FEM-SmEdA和DMF方法預示的子系統各階模態的能量值基本一致。如圖 8所示,聲腔的模態能量略有差異,且兩個結果之間的差值小于1 dB。因此,基于FEM-SmEdA得到的TBL激勵下子系統的模態能量是準確合理的。

3.3 結構局部能量響應預示

基于FEM-SmEdA預示TBL激勵下平板的動能和勢能分布,如圖9和10所示。通過與圖5和6對比可知,在中心頻率均為630 Hz的1/3倍頻程頻帶內,TBL激勵下平板的局部響應與單點力作用不同。

4 艙段/聲腔耦合系統局部能量響應預示

4.1 艙段/聲腔耦合系統能量分布預示 ?以廣泛用于航空航天器的艙段/聲腔耦合系統為研究對象,開展湍流邊界層激勵下耦合結構的局部能量響應預示。艙段和聲腔的有限元模型如圖11所示,其幾何尺寸及材料屬性分別如表3和4所示。

基于FEM獲得艙段及聲腔的模態結果,由式(2)可得分析頻段Δω1(561.2-2828.4 Hz)內艙段與聲腔之間的模態耦合損耗因子,如圖12所示。由圖12可知,較大的模態耦合損耗因子多存在于共振模態間,非共振模態間的耦合損耗因子相對較小。由文獻[18]可知:對于小阻尼系統,子系統間的非共振傳輸可忽略,因此式(7)中可只考慮分析頻帶內的共振模態。本文研究對象屬于小阻尼系統,為提高計算效率,分析中僅考慮共振模態間的功率傳輸。

基于FEM-SmEdA計算艙段和聲腔在中心頻率分別為630,1000,2500,4000 Hz的1/3倍頻程頻帶內的能量空間分布,計算結果如表 5所示。由表 5可知:(1)在不同頻段內,結構及聲腔的動能和勢能分布基本一致,艙段動能的最大值高于勢能的最大值,而聲腔勢能的最大值高于動能的最大值。(2)結構的能量響應遠大于聲腔。(3)隨著分析頻率的升高,艙段和聲腔的能量分布更加均勻。當分析頻帶中心頻率分別為630,1000,2500 Hz時,艙段和聲腔的能量分布不均勻程度較大;當分析頻帶中心頻率為4000 Hz時,艙段和聲腔的能量分布較為均勻。隨著分析頻率的升高,艙段和聲腔的能量分布更接近SEA中子系統能量平均假設,此時使用SEA方法計算艙段的平均能量即可較為準確地表征子系統的能量大小。

4.2 加筋艙段/聲腔耦合系統能量分布預示

為了驗證FEM-SmEdA的適用性,以復雜加筋艙段/聲腔耦合系統為研究對象,基于FEM-SmEdA預示聲-固耦合系統的能量分布。加筋艙段的幾何模型如圖13所示。其幾何尺寸及材料屬性如表3所示。筋條包括8根縱向筋,3根環向筋。筋條材料屬性與艙段一致,筋條高度為0.02 m、厚度為4 mm。

基于FEM獲得加筋艙段及聲腔的模態結果,在中心頻率為1000 Hz的1/3倍頻程頻段Δω2(890.9-1122.5 Hz)內,加筋艙段和聲腔分別有14階和12階模態振型。將加筋艙段和聲腔的模態結果代入式(2),可得如圖 14所示的Δω2頻段內加筋艙段與聲腔之間模態耦合損耗因子。由圖 14可知,Δω2頻段內的共振模態間的模態耦合損耗因子只有少數模態耦合損耗因子較大,絕大多數模態耦合損耗因子較小,不滿足均一化假設。

在加筋艙段表面施加湍流邊界層激勵,基于FEM-SmEdA分別計算加筋艙段和聲腔在Δω2內的能量分布,如表 6所示。

通過對比表 5和6可知,加筋后艙段結構和聲腔的能量分布更加均勻;加筋后艙段和聲腔能量分布最大值變小。由4.1節圖 12可得分析頻段Δω2(890.9-1122.5 Hz)內艙段與聲腔之間的模態耦合損耗因子,如圖 15所示。對比圖 14和15可知,加筋后艙段與聲腔的耦合強度減弱。這主要是由于筋條改變了艙段的固有頻率,并同時抑制了部分模態振型,使得艙段結構向聲腔內傳遞的能量減少。

為更加深入地說明加筋對艙段與聲腔間的耦合強度影響,對加筋前后的艙段在分析頻帶內施加相同的輸入功率,對比加筋前后艙段和聲腔的總能量,如表 7所示。由表 7可知,加筋后艙段的能量減少了21.4%,聲腔能量減少了82.0%,加筋對聲腔響應抑制更加明顯。

5 結 論

本文結合FEM和SmEdA理論,形成了適用于聲-固耦合系統的高頻局部能量響應預示方法。首先,以平板/聲腔耦合系統為研究對象,通過與FEM的對比驗證了FEM-SmEdA方法的準確性;隨后,分別將該方法應用于TBL激勵下結構局部響應預示及艙段/聲腔耦合系統的高頻局部能量響應預示,利用DMF法對基于FEM-SmEdA法的TBL激勵下結構局部響應預示結果進行驗證,并將該方法應用于艙段/聲腔耦合系統,揭示艙段結構和聲腔的能量分布隨頻率的變化規律,并研究加筋對聲-固耦合系統能量分布的影響。結果表明:

1)FEM-SmEdA預示得到的能量響應大小和空間分布與FEM預示結果基本一致,這說明FEM-SmEdA方法能夠準確地預示聲-固耦合系統的局部能量響應;

2)基于FEM-SmEdA和DMF方法預示的子系統各階模態能量結果的基本一致。因此,基于FEM-SmEdA得到的TBL激勵下子系統的模態能量是準確合理的,該方法適于TBL激勵下結構的局部響應預示。

3)對于TBL激勵下的艙段/聲腔耦合系統,艙段及聲腔的動能和勢能的分布基本一致,艙段的能量響應遠大于聲腔;隨著分析頻率的升高,艙段和聲腔的能量分布更加均勻;

4)對于TBL激勵下的艙段/聲腔耦合系統,加筋后艙段結構和聲腔的能量分布更加均勻;加筋后艙段和聲腔能量分布最大值變小;與加筋后艙段的總能量減少值相比,聲腔總能量減少較大,故加筋對聲腔響應抑制更加明顯。

參考文獻:

[1] 李彥斌, 張 鵬, 吳邵慶,等. 復合材料加筋板計及熱效應的聲-固耦合分析[J]. 振動工程學報, 2015,28(4):531-540.

Li Yanbin, Zhang Peng, Wu Shaoqing, et al. Structural- acoustic coupling analysis of a composite stiffened panel in a thermal environment[J]. Journal of Vibration Engineering, 2015,28(4):531-540.

[2] Yao L Y, Yu D J, Cui X Y, et al. Numerical treatment of acoustic problems with the smoothed finite element method[J]. Applied Acoustics, 2010, 71(8):743-753.

[3] 李彥斌, 張 鵬, 鄒元杰,等. 隨機基礎激勵下承力筒-蒙皮結構的聲-固耦合分析[J]. 宇航學報, 2015, 36(2):236-242.

Li Yanbin, Zhang Peng, Zou Yuanjie, et al. Structural- acoustic coupling analysis of a bearing cylinder-skin structure under random base excitation[J]. Journal of Astronautics, 2015, 36(2):236-242.

[4] Lyon R H. Theory and Application of Statistical Energy Analysis[M]. Elsevier, 2014.

[5] Chen Qiang, Fei Qingguo, Wu Shaoqing, et al. An efficient transient analysis method for time-varying structures based on statistical energy analysis[J]. Mechanics Research Communications, 2018, 91:93-99.

[6] Chen Qiang, Fei Qingguo, Li Yanbin, et al. Prediction of the transient energy response for complex vibro-acoustic systems[J]. Journal of Mechanical Science and Technology, 2019, 33(2): 495-504.

[7] Mace B R. Statistical energy analysis: Coupling loss factors, indirect coupling and system modes[J]. Journal of Sound & Vibration, 2005, 279(1):141-170.

[8] Karnopp D. Coupled vibratory-system analysis, using the dual formulation[J]. Journal of the Acoustical Society of America, 1966, 40 (2):380-384.

[9] Tabarrok B. Dual formulations for acoustic-structural vibrations[J]. International Journal for Numerical Methods in Engineering, 1978, 13(1):197-201.

[10] Maxit L, Guyader J L. Estimation of SEA coupling loss factors using a dual formulation and FEM modal information, part I: Theory[J]. Journal of Sound and Vibration, 2001, 238 (5):907-930.

[11] Maxit L, Guyader J L. Estimation of SEA coupling loss factors using a dual formulation and FEM modal information, part II: Numerical applications[J]. Journal of Sound and Vibration, 2001, 238 (5):931-948.

[12] ngels Aragonès, Maxit L, Guasch O. A graph theory approach to identify resonant and non-resonant transmission paths in statistical modal energy distribution analysis[J]. Journal of Sound & Vibration, 2015, 350(8):91-110.

[13] Totaro N, Guyader J L. Extension of the statistical modal energy distribution analysis for estimating energy density in coupled subsystems[J]. Journal of Sound & Vibration, 2012,331(13):3114-3129.

[14] Qiao Y, Chen H B, Luo J L. Estimation of shell radiation efficiency using a FEM-SmEdA algorithm[J]. Journal of Vibroengineering, 2013,15(3):1130-1146.

[15] 于士甲, 張 鵬, 李彥斌,等. 基于FEM-MODENA的加筋板聲-固耦合分析[J]. 振動工程學報, 2017,30(2):262-269.

Yu Shijia, Zhang Peng, Li Yanbin, et al. Structural-acoustic coupling analysis of a stiffened plate based on FEM-MODENA[J]. Journal of Vibration Engineering, 2017,30(2):262-269.

[16] Ichchou M N, Hiverniau B, Troclet B. Equivalent ‘rain on the roof loads for random spatially correlated excitations in the mid-high frequency range[J]. Journal of Sound and Vibration, 2009,322(4):926-940.

[17] Karnopp D. Coupled vibratory-system analysis, using the dual formulation[J]. Journal of the Acoustical Society of America, 2005,40(2):380-384.

[18] Maxit L, Ege K, Totaro N, et al. Non resonant transmission modelling with statistical modal energy distribution analysis[J]. Journal of Sound & Vibration, 2013, 333(2):499-519.

[19] Maxit L, Guyader J L. Extension of SEA model to subsystems with non-uniform modal energy distribution[J]. Journal of Sound and Vibration, 2003, 265(2): 337-358.

[20] Corcos G M. The structure of the turbulent pressure field in boundary-layer flows[J]. Journal of Fluid Mechanics, 1964, 18(3):353-378.

[21] Graham W R. A comparison of models for the wavenumber frequency spectrum of turbulent boundary layer pressures[J]. Journal of Sound and Vibration, 1997, 206(4): 541-565.

Abstract: By combining the finite element method (FEM) and the statistical modal energy distribution analysis (SmEdA), a FEM-SmEdA approach is presented to predict the local energy response of structural-acoustic coupling system under high frequency excitation. First, the accuracy of the FEM-SmEdA is verified by a plate/cavity coupling system. Then the approach is applied to a plate/cavity coupling system and a cabin/cavity coupling system under Turbulent Boundary Layer (TBL) excitation to verify its applicability respectively. The energy responses predicted by FEM-SmEdA are verified by Dual Modal Formulation (DMF) method. The influence of analysis frequency band and stiffener on the energy distribution of cabin/cavity coupling system under TBL excitation are investigated by the FEM-SmEdA. Results show that the FEM-SmEdA is capable of predicting the local energy response of structural-acoustic coupling system under the high frequency excitation with sufficient accuracy. This approach is well adapted for predicting the local response of systems under TBL excitation. With the increase of analysis frequency band, the energy distributions of cabin and cavity become more uniform. The energy distribution of stiffened cabin/cavity coupling system is more uniform than that of cabin/cavity coupling system. The suppression of stiffener to response of cavity is more obvious than structure.

Key words: structural-acoustic coupling system; SmEdA; modal energy; energy distribution; FEM

作者簡介: 王 攀(1993-),男,碩士研究生。電話: (025)83790168; E-mail: 220160935@seu.edu.cn

通訊作者: 費慶國(1977-),男,博士,教授,博士生導師。電話: (025)83790168; E-mail: qgfei@seu.edu.cn

猜你喜歡
有限元法
正交各向異性材料裂紋疲勞擴展的擴展有限元法研究
基于有限元法的高頻變壓器繞組損耗研究
基于有限元法副發動機托架輕量化設計
專用汽車(2016年8期)2016-03-01 04:16:43
傳遞矩陣法與有限元法計算電機轉子臨界轉速的對比分析
Sine-Gordon方程H1-Galerkin非協調混合有限元法的誤差分析
三維有限元法在口腔正畸生物力學研究中發揮的作用
RKDG有限元法求解一維拉格朗日形式的Euler方程
計算物理(2014年1期)2014-03-11 17:00:14
集成對稱模糊數及有限元法的切削力預測
有限元法在機械設計方向中的教學實踐
基于HCSR和CSR-OT的油船疲勞有限元法對比分析
船海工程(2013年6期)2013-03-11 18:57:25
主站蜘蛛池模板: 亚洲香蕉久久| 亚洲日本中文字幕天堂网| 欧美不卡在线视频| 国产福利不卡视频| 久久婷婷人人澡人人爱91| 中文字幕资源站| 免费日韩在线视频| 日本道中文字幕久久一区| 国产日韩欧美中文| 国产99在线| 国产女人在线视频| 四虎影视库国产精品一区| 丁香婷婷在线视频| 成人在线亚洲| 亚洲午夜天堂| 精品亚洲欧美中文字幕在线看| 波多野结衣在线se| 国产自在线播放| 9久久伊人精品综合| 2021无码专区人妻系列日韩| 国产成人高精品免费视频| 国产精品久久自在自2021| 欧美亚洲一区二区三区在线| 国产精品偷伦在线观看| 亚洲精品欧美日本中文字幕| 国内熟女少妇一线天| 国产精品视频免费网站| 欧美α片免费观看| 黑人巨大精品欧美一区二区区| 亚洲色图欧美视频| 精品人妻无码中字系列| 免费观看国产小粉嫩喷水| 无码福利日韩神码福利片| 十八禁美女裸体网站| 久久免费视频播放| 国产网站在线看| 四虎永久在线视频| 国产小视频在线高清播放| 极品私人尤物在线精品首页| 伊人蕉久影院| 青青久在线视频免费观看| 亚洲AV无码乱码在线观看裸奔 | 亚洲视频a| 欧美高清三区| 亚洲综合在线最大成人| 伊人国产无码高清视频| 免费看a级毛片| 国产日本视频91| 另类欧美日韩| 国产成人高精品免费视频| 最新国产精品第1页| 国产成人一二三| 亚洲精品视频网| 国产嫖妓91东北老熟女久久一| 欧美综合成人| 毛片在线看网站| 成人中文在线| 人妻一区二区三区无码精品一区| 欧美伦理一区| 高清视频一区| 亚洲国产成人超福利久久精品| 婷婷激情亚洲| 国外欧美一区另类中文字幕| 亚洲欧美人成人让影院| 亚洲国产精品人久久电影| 在线永久免费观看的毛片| 国产精品视频a| 亚洲永久精品ww47国产| 国产精品片在线观看手机版| 福利国产微拍广场一区视频在线| 中国黄色一级视频| 亚洲午夜天堂| 在线观看免费AV网| 国产一区二区精品高清在线观看| 欧美高清三区| 精品国产美女福到在线不卡f| 亚洲欧美激情另类| 亚洲视频四区| 国产18在线| 亚洲国产成人久久77| 无码 在线 在线| 国产永久在线视频|