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

艦船低頻水下輻射噪聲的聲固耦合數值計算方法

2018-02-27 11:13:58楊德慶
振動與沖擊 2018年3期
關鍵詞:有限元

李 清, 楊德慶, 郁 揚

(上海交通大學 船舶海洋與建筑工程學院 高新船舶與深海開發裝備協同創新中心, 上海 200240)

艦船水下輻射噪聲頻率范圍覆蓋1 Hz~10 kHz。計算水下輻射噪聲需建立不同頻段的分析模型,采用不同的數值計算方法求解。目前計算艦船低中頻域水下輻射噪聲最常用流程是:首先運用大型結構動力學有限元軟件(例如MSC/NASTRAN、ABAQUS或ANSYS等)進行艦船結構流固耦合振動響應的計算;接著將船體濕表面振動速度直接作為聲學計算的外載荷邊界條件,利用聲學有限元或邊界元軟件(例如LMS Virtual.Lab Acoustics、MSC/ACTRAN軟件)計算水下聲學物理量[1-5]。對于艦船水下輻射噪聲的合理建模方法、計算方法、及其計算精度和計算效率等尚未結合流固耦合與聲固耦合的理論本質進行更深入地探討。

本研究認為,目前常用的低頻域艦船水下輻射噪聲計算方法采用的是基于無黏、無旋、不可壓縮流體假設下的流固耦合模式,把船體浸水濕表面振動速度作為邊值問題計算聲場。這種計算方法忽略了水作為可壓縮性聲學介質對船體振動的影響,評估結果在計算精度上尚有提升空間。本論文以某SWATH船為研究對象,聚焦低頻域1 Hz~100 Hz艦船水下輻射噪聲數值計算,提出艦船水下輻射噪聲計算基于聲固耦合動力學方程的兩種標準數值計算方法,該方法是基于無黏、無旋、可壓縮流體假設下的聲固耦合模式。使用聲功率級作為評價標準,探討了聲固耦合有限元數值模型中船體周圍流體域特征尺度的選取,比較分析了上述兩種標準方法與常規的基于流固耦合的兩種方法在計算特性方面的差異,為艦船水下輻射噪聲預報規范化提供參考與建議。

1 聲固耦合動力學基本理論

與流體介質接觸的結構物在受到流體動力激勵或機械載荷等而發生振動時,其周圍流場亦發生變化,流場變化反過來使結構所受流體動力載荷發生變化,形成反饋的流體-結構相互作用,這類問題稱為流固耦合問題。結構在空氣或水等流場中發生振動導致流場的壓力脈動而產生聲音,流場的壓力脈動反過來對結構的振動產生影響,形成反饋的聲-結構相互作用,這類問題稱為聲固耦合問題,屬于廣義的流固耦合問題。二者本質上是分別在力學范疇和聲學范疇內考慮結構振動響應與流體脈動壓力場之間的聯系。艦船水下輻射噪聲的產生機理是上述流固耦合或聲固耦合作用的結果。本文研究艦船聲固耦合振動及聲輻射特性時假定流場絕對靜止,即不考慮流體動力激勵。

1.1 聲固耦合動力學方程有限元列式

均質、無黏和絕熱狀態下流體內縱向波的聲學波動方程為

(1)

式中:p為瞬時聲壓;c0為聲波在介質中的傳播速度。聲固耦合問題的邊界條件為:

(1) 在流固交界面SI上

(2)

(2) 在固定邊界Sb上

(3)

(3) 在自由表面SF上

(4)

(4) 在無限遠邊界Sr上

(5)

對于聲固耦合問題,采用有限元離散模式將結構振動和聲場在同一個耦合環境中計算,結構有限元耦合聲學有限元的聲固耦合動力學方程的位移-壓力格式如下[6-8]

(6)

(7)

(8)

(9)

(10)

(11)

而在聲固耦合問題中,聲波傳播依賴于彈性媒質的可壓縮性,流體內各質點的運動相位是有差異的,運動相位差的存在說明了阻尼的存在,從而耗散能量;將流體視為不可壓縮時流場內各質點運動是同相位的,因而它不產生阻尼[10]。聲場對結構的影響表現為附加質量(式(7)及式(10))和阻尼(式(8))的共同影響。尤其對于機械載荷的高頻激勵,聲波波長遠小于船體振動的特征長度,此時水的壓縮性對船體振動產生較大影響。特別地,對于沖擊波等瞬態載荷作用下流場和聲場的可壓縮性均需充分考慮[11-12]。

因此,流體不可壓縮假設下的流固耦合動力學方程與流體可壓縮假設下的聲固耦合動力學方程的有限元列式是不相同的,兩者不能等價,前者是后者的低頻近似處理方式。

1.2 聲固耦合動力學數值解法

對于復雜結構的聲固耦合問題,大多數情況下無法求得解析解而多采用數值解。常用的聲學數值計算方法較多,主要有聲學有限元法、聲學邊界元法、聲線法、統計能量法、能量有限容積法、有限元-統計能量混合法和混合能量有限元-統計能量分析法等。各種聲學數值方法適用于不同頻率范圍,其中聲學有限元法、聲學邊界元法適用于低頻噪聲計算,本文建議采用如下方法求解聲固耦合問題。

1.2.1 聲固耦合有限元及遠場自動匹配層方法

采用聲學有限元進行船體水下聲輻射分析,一方面,需要在船體附連水區域劃分聲學實體網格,該過程較為繁瑣,且數值計算結果在一定程度上依賴于聲場區域大小、網格尺寸和網格質量;另一方面,因為輻射場網格區域不可能無限大,在求解半封閉空間遠場輻射問題上,需要在聲學輻射邊界定義聲輻射邊界條件,利用邊界元方法求出遠場聲輻射。LMS Virtual. Lab Acoustics聲學軟件的聲輻射邊界條件是從無限元技術,到PML(Perfect Matched Layer)完美匹配層吸收邊界條件,最終發展到AML(Automatic Matched Layer)自動匹配層邊界條件[13-14]。作為目前最先進的技術,AML方法能夠根據計算頻率自動生成并調整PML層,完全取代需要手動繪制聲學匹配層的PML方法,很容易就能滿足低頻和高頻計算要求。目前在LMS Virtual. Lab Acoustics操作環境下可采用耦合有限元方法并設定AML吸聲條件求解船舶水下聲輻射問題。

1.2.2 聲固耦合邊界元方法

式(6)已列出了結構有限元耦合聲學有限元的具體格式,同樣地可以將聲場通過邊界元格式離散與結構有限元耦合求解。邊界元法分為直接邊界元法和間接邊界元法,所需網格為面網格。直接邊界元法要求網格封閉,而間接邊界元法的網格可封閉,也可不封閉。由于直接邊界元法的網格是封閉的,所以直接邊界元法可以計算封閉網格內部聲場,或計算封閉網格外部聲場,但是不能同時計算內部聲場和外部聲場;由于間接邊界元法可以不封閉,所以間接邊界元法可以同時計算內聲場和外聲場,艦船濕表面這類非封閉結構宜于采用間接邊界元法進行求解。對于直接邊界元法和間接邊界元法聲振耦合分析,結構振動位移與聲場分布是在考慮耦合邊界處速度連續的基礎上同步求解的。滿足單層勢σ=0的結構有限元耦合流體間接邊界元的聲固耦合系統頻域下動力學方程如下

(12)

式中:Lc為耦合矩陣;D為間接邊界元影響矩陣;μ為節點雙層勢向量,即結構表面聲壓差。盡管聲學有限元在求解內聲場和聲輻射方面有極大優勢,但在求解超大結構外場聲輻射問題上,聲學邊界元法仍然是最優的選擇。對于船舶水下聲輻射,用邊界元方法只需提取船體濕表面的面網格作為邊界元網格即可完成計算。

2 艦船低頻水下輻射噪聲聲固耦合計算方法

艦船水下輻射噪聲計算常見計算方法是低頻域(0 Hz~100 Hz) 使用聲學有限元/邊界元(FEM/BEM)方法建模與計算,中高頻域(125 Hz~10 kHz)采用FEA-SEA混和法及統計能量SEA法建模與計算,最后進行綜合評價。本文對艦船低頻水下輻射噪聲預報數值計算方法進行比較研究,將水下輻射噪聲預報數值計算方法分為三類耦合分析模式。

2.1 艦船水下輻射噪聲數值計算的三種模式

流固耦合分析模式:該模式首先求解機械載荷下基于流固耦合動力學方程的船體濕表面速度頻響,之后將其作為聲學激勵邊界條件,導入聲學FEM或BEM模型進行后續分析,這是目前國內預報艦船水下聲輻射問題的普遍做法。盡管該模式的聲學模型較簡單,但在計算中忽略了聲學介質可壓縮性對船體結構振動的影響。

聲固耦合分析模式:該模式首先建立全船結構FEM模型并計算整船結構模態(干模態),之后導入統一的聲學計算環境中,建立耦合聲學FEM或BEM模型,最終進行聲固耦合一體化分析。該模式充分考慮了流體對結構的耦合作用,對流體域采用有限元或邊界元離散建模;同時,船體振動模態計算不必先行進行流固耦合計算,只需計算或導入船體干模態。結構-流體耦合作用的影響是通過耦合動力學方程同時求解而加以考慮,完全實現聲固耦合系統在同一求解器下的共同分析,是嚴格遵循耦合方程推導結果的處理方式。

流固耦合+聲固耦合的關聯(組合)分析模式:該模式先以流固耦合模式求解出船體濕模態,然后將濕模態導入聲固耦合分析模式中,重新對流體和結構進行建模,最終求解聲固耦合動力學方程。本文認為該模式重復考慮了流體耦合作用,將導致水下噪聲計算結果偏小。

2.2 艦船水下輻射噪聲計算的聲固耦合方法

本文提出兩種基于聲固耦合模式的艦船水下聲輻射數值計算方法,分別是:① 耦合聲學有限元及自動匹配層FEM/AML方法;② 耦合聲學間接邊界元IBEM方法。本文結合Virtual. Lab Acoustics實現上述算法。

耦合聲學有限元及自動匹配層(FEM/AML)方法:設定聲學有限元環境,在流體網格對應的自由液面處賦予空氣聲阻抗屬性模擬吸聲邊界,通過導入船體的干模態,直接輸入機械載荷及噪聲源,進行船體結構振動及水下聲輻射的同步計算。

耦合聲學間接邊界元方法:設定聲學邊界元環境,通過構建面內壓力為零的反對稱平面模擬自由液面軟邊界,通過導入船體的干模態,直接輸入機械載荷及噪聲源,進行結構振動及水下聲輻射的同步計算。

2.3 艦船水下輻射噪聲計算評估標準

進行艦船水下輻射聲場計算時,由于輻射聲場具有指向性,近場、遠場的區分界限不明確。本文建議艦船水下輻射噪聲水平評價的物理量選取為聲功率,亦可進一步換算為離聲中心1 m遠處的聲源級。聲功率表示單位時間內通過垂直聲傳播方向面積的平均聲能量,無論近場和遠場,解決了以聲壓作為評價指標對于流體域遠場近場選取的嚴格要求。

3 小水線面雙體船低頻水下輻射噪聲預報算例

為驗證上述方法,以某小水線面雙體船為研究對象,進行低頻域1 Hz~100 Hz艦船水下機械輻射噪聲數值計算,探討本文提出的聲固耦合計算方法與常規的流固耦合算法在建模難度、計算規模、計算精度和求解效率等方面的差異。

在Genuine Inter(R) CPU @ 2.9 GHz(20核)128 G(RAM) 64位操作系統高性能服務器環境下,采用幅值為1 kN的垂向激振力作為虛擬載荷,載荷頻率范圍為1 Hz~80 Hz,取1/3倍頻程的中心頻率,激勵位置在主機機腳處。采用模態疊加法計算聲學響應時,截取船體前15 000階(80.98 Hz)干模態,對應固有頻率大于分析頻段最大頻率,基本可兼顧工程計算的精度與效率。

對于聲學有限元與邊界元模型,在將聲場劃分成網格時,有限元或邊界元的網格大小要劃分一致,且流體模型的計算精度取決于整體網格尺寸,通常要假設在最小波長內有6個單元,本文研究的SWATH低頻聲輻射問題要求網格尺寸小于2.5 m,以下聲學有限元或邊界元計算模型的單元特征長度均符合上述要求。根據本船主尺度顯示聲壓分布圖的平面取長1 000 m、深500 m區域。

3.1 聲固耦合聲學FEM/AML方法

設定聲學有限元環境,導入船體結構有限元網格與聲學有限元網格,對應的自由液面處為吸聲邊界,賦空氣聲阻抗特性Zp=416.5 kg/(m2·s),遠場特征邊界賦AML屬性。導入船體干模態,直接輸入機械載荷進行船體結構振動及水下聲輻射的同步計算。計算模型見圖1,輻射聲壓級分布見圖2。

圖1 耦合聲學FEM/AML方法計算模型

圖2 耦合聲學FEM/AML方法鉛垂聲壓分布(25 Hz)

3.2 聲固耦合聲學IBEM方法

設定聲學邊界元環境,導入船體結構有限元網格,并以船體濕表面網格作為聲學邊界元網格,構建面內壓力為零的反對稱平面模擬自由液面絕對軟邊界,導入船體干模態,直接輸入機械載荷進行船體結構振動及水下聲輻射的同步計算。計算模型見圖3,輻射聲壓級分布見圖4。

圖3 耦合聲學IBEM方法計算模型

圖4 耦合聲學IBEM方法鉛垂聲壓分布(25 Hz)

3.3 聲場區域特征尺度選取

在流固耦合有限元分析中,需要區分流體的遠場和近場,流體域特征尺度通常可取分析頻段最大波長λ的兩倍以上,或取5倍的船寬B。為揭示聲場流體域尺寸大小對聲輻射計算結果的影響,本文選取尺度為1.5B、2B、4B、5B、6B、8B和10B的聲場區域,建立多個基于聲固耦合的聲學FEM/AML法模型如圖5,對比聲場區域不同尺度對計算結果的影響。參考圖6和表1,計算結果在一定程度上依賴于聲場區域尺度,但當聲場區域尺度不斷增大并超過4倍船寬時,計算結果基本趨于收斂。在既能保證聲學計算精度又不犧牲計算效率的前提下,本文建議采用聲固耦合聲學FEM方法時,聲場取特征尺度為5倍船寬的區域進行離散,可基本滿足計算精度要求。

(a) 2倍船寬聲場區域

(b) 4倍船寬聲場區域

(c) 8倍船寬聲場區域

圖6 聲場區域尺度對耦合聲學FEM/AML方法計算影響

Fig.6 The effects of different acoustic field sizes on computing results in acoustic coupling FEM/AML model

表1 不同聲場區域尺度下合成聲功率級

3.4 流固和聲固耦合計算方法特性比較

為對比本文聲固耦合方法和常規流固耦合計算模式下無耦合聲學有限元與遠場自動匹配層(FEM/AML)方法,以及無耦合聲學間接邊界元方法的計算精度,圖7一并給出了四種標準算法的SWATH船水下輻射聲功率級曲線。對于同一種離散方法下的兩種模式,由于聲固耦合模式考慮流體的阻尼效應,耗散了船體振動能量,聲固耦合計算方法較流固耦合計算聲功率級偏小;而同一種模式下的兩種離散方法計算結果的差異,主要是由于有限元法和邊界元法在模型邊界條件簡化上的區別而產生的模型誤差。前者通過設定空氣聲阻抗吸聲邊界模擬有限大的自由液面,并在遠場設定自動匹配層條件,而后者直接采用無限大的反對稱邊界模擬自由液面,后者更為精確。

圖7 SWATH水下輻射四種標準方法聲功率級對比

Fig.7 The comparison of SWATH underwater radiation sound power level in 4 standard methods

圖8 SWATH水下輻射各計算方法合成總聲功率級對比

Fig.8 The comparison of SWATH Underwater radiation resultant sound power level in all computing methods

將各中心頻率下對應的聲功率級合成為總聲功率級如圖8,同時也計算了導入濕模態的流固耦合+聲固耦合組合分析模式。對于本文SWATH船,流固耦合模式較聲固耦合模式的聲功率級計算值高約1 dB~3 dB,評估結果偏大,設計上偏于保守;而流固+聲固耦合組合分析模式由于重復考慮耦合作用,計算結果比四種標準算法偏小約5 dB~8 dB,評估結果偏小,設計上偏于冒進,容易造成實際噪聲評估結果超標。

表2匯總對比了四種標準算法的建模難度、計算規模和求解效率,其中計算時間僅針對聲學計算,不包括結構模態及振動響應求解。通過權衡比較,耦合聲學IBEM方法不需劃分聲場體網格,遠場采用邊界元面網格進行遞推計算,建模難度較低;耦合聲學IBEM方法嚴格遵循聲固耦合動力學理論,模型簡化也更為精確,盡管求解效率有所降低,但保證了相對更高的計算精度。

表2艦船水下聲輻射數值方法計算特性對比

Tab.2Thecomparisonofcomputingpropertiesindifferentshipunderwaterradiationnumericalmethods

聲輻射數值方法建模難度計算規模求解效率(時間)聲學FEM/AML高小高(1h)聲學IBEM低小高(1h)耦合聲學FEM/AML高大低(55h)耦合聲學IBEM低大中(16h)

4 結 論

基于無黏、無旋、可壓縮流體的線性小擾動假設,本文提出艦船水下輻射噪聲數值計算的聲固耦合動力學的兩種標準方法。指出目前普遍采用的基于勢流理論的流固耦合計算方法尚有改進余地,宜采用更為嚴格的聲固耦合分析模式。

由于彈性媒質可壓縮性產生的阻尼效應,聲固耦合模式計算方法的振動聲輻射響應較流固耦合模式偏小,計算結果更為精確,對于本文SWATH船,水下輻射噪聲合成總聲功率級偏差可達到1 dB~3 dB,應當予以一定的關注。

通過權衡比較各算法的計算特性,本文認為基于聲固耦合模式的耦合聲學IBEM方法具備建模難度低、模型精度高等優勢,是進行艦船水下聲輻射預報數值計算的首選算法。

[1] 楊德慶,鄭靖明,王德禹,等. 基于SYSNOISE軟件的船舶振動聲學數值計算[J]. 中國造船,2002,43(4):32-38.

YANG Deqing, ZHENG Jingming, WANG Deyu, et al. Numerical analysis of vibro-acoustic characters of ship with SYSNOISE software[J]. Ship Building of China, 2002, 43(4): 32-38.

[2] 楊德慶,王德禹,劉洪林,等. 某型艇近場噪聲和自噪聲數值計算[J]. 聲學學報,2003(5):421-424.

YANG Deqing, WANG Deyu, LIU Honglin, et al. Numerical analysis of acoustic characters in near field and self-noise of

ship[J]. Acta Acustica, 2003(5): 421-424.

[3] 鄒春平,陳端石,華宏星. 船舶水下輻射噪聲特性研究[J]. 船舶力學,2004,8(1):113-124.

ZOU Chunping, CHEN Duanshi, HUA Hongxing. Study on characteristics of ship underwater radiation noise[J]. Journal of Ship Mechanic, 2004, 8(1): 113-124.

[4] 黃毅. 小水線面雙體船噪聲環境特性研究[D].哈爾濱:哈爾濱工程大學,2011.

[5] 付建, 王永生, 丁科, 等. 螺旋槳激振力作用下船體振動及水下輻射噪聲研究[J]. 船舶力學,2015,19(4):470-476.

FU Jian, WANG Yongsheng, DING Ke, et al. Research on vibration and underwater radiated noise of ship by propeller excitations[J]. Journal of Ship Mechanics, 2015, 19(4):470-476.

[6] ATALLA N. Review of numerical solutions for low-frequency structural-acoustic problems[J]. Applied Acoustics, 1994, 43(3): 271-294.

[7] 王勖成. 有限單元法[M]. 北京:清華大學出版社,2003.

[8] LUIS R, JOSé A G, ANTONIO C. Partitioned solution strategies for coupled BEM-FEM acoustic fluid-structure interaction problems[J]. Computers and Structures, 2015, 152: 45-58.

[9] 鄒明松. 船舶三維聲彈性理論[D]. 北京:中國艦船研究院,2014.

[10] 張升明. 流體的可壓縮性對彈性結構振動的影響[J]. 水動力學研究與進展(A輯),1994,9(4):429-436.

ZHANG Shengming. The influence of fluid compressibility to structure vibration[J]. Journal of Hydrodynamics (Ser. A),1994,9(4):429-436.

[11] 王崢,洪明,劉城. 基于FEM/BEM的浸水結構振動及聲輻射特性國內研究綜述[J]. 船舶力學, 2014,18(11):1397-1414.

WANG Zheng, HONG Ming, LIU Cheng. Domestic review of the submerged structure vibration and acoustic radiation characteristics based on FEM/BEM[J]. Journal of Ship Mechanic, 2014, 18(11): 1397-1414.

[12] 湯渭霖,范軍. 水中彈性結構聲散射和聲輻射機理—結構和水的聲-振耦合作用[J]. 聲學學報,2004(5):385-392.

TANG Weilin, FAN Jun. Mechanisms of sound scattering and radiation of submerged elastic structure-vibro-acoustic coupling of structure and water[J]. Acta Acustica, 2004(5): 385-392.

[13] 張冠軍,朱翔,李天勻,等. 水中雙層加筋板結構的聲振耦合特性[C]∥中國造船工程學會船舶力學學術委員會第八次全體會議.大連:中國造船工程學會船舶力學學術委員會,2014.

[14] 李增剛,詹福良. Virtual.Lab Acoustics聲學仿真計算高級應用實例[M]. 北京:國防工業出版社,2011.

猜你喜歡
有限元
基于擴展有限元的疲勞裂紋擴展分析
非線性感應加熱問題的全離散有限元方法
TDDH型停車器制動過程有限元分析
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于I-DEAS的履帶起重機主機有限元計算
基于有限元模型對踝模擬扭傷機制的探討
10MN快鍛液壓機有限元分析
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 免费女人18毛片a级毛片视频| 激情爆乳一区二区| 久久精品国产精品青草app| 好吊色妇女免费视频免费| 国产精品手机在线观看你懂的| 中文字幕无线码一区| 中文字幕伦视频| 午夜精品区| 日韩视频免费| 午夜视频在线观看免费网站| 97se亚洲综合在线韩国专区福利| 国产精品美乳| 精品无码一区二区三区在线视频| 亚洲色成人www在线观看| 国产精品流白浆在线观看| 国产乱子伦精品视频| 国产又爽又黄无遮挡免费观看| 日本不卡在线播放| 精品黑人一区二区三区| 亚洲精品无码av中文字幕| 综合色婷婷| 欧美成人看片一区二区三区| 亚洲欧美日韩另类| 91精品国产麻豆国产自产在线| 国产一级毛片yw| 欧美精品三级在线| 伊人激情久久综合中文字幕| 在线高清亚洲精品二区| 成年女人a毛片免费视频| 日韩在线永久免费播放| 欧美精品啪啪一区二区三区| 一级全免费视频播放| 国产呦视频免费视频在线观看| 国产精品视频系列专区| 国产视频只有无码精品| 国产成人一二三| 国产中文一区a级毛片视频| 992tv国产人成在线观看| 日本黄网在线观看| 亚洲精品日产精品乱码不卡| 亚洲日本在线免费观看| 黑人巨大精品欧美一区二区区| 亚洲成肉网| 国产美女视频黄a视频全免费网站| 免费无码一区二区| 国产丝袜一区二区三区视频免下载| 国产h视频在线观看视频| 中文无码精品a∨在线观看| 在线另类稀缺国产呦| 久久亚洲国产最新网站| 日韩国产亚洲一区二区在线观看| 色播五月婷婷| 成人精品免费视频| 免费A级毛片无码无遮挡| 国产欧美性爱网| 亚洲h视频在线| 精品一区二区三区水蜜桃| 久久无码免费束人妻| 亚洲精品中文字幕无乱码| 久久这里只有精品8| 亚洲精品日产AⅤ| 亚洲无码熟妇人妻AV在线| 国产精品自在在线午夜| 亚洲天堂首页| 色妞永久免费视频| 九九热视频精品在线| 无码又爽又刺激的高潮视频| 91最新精品视频发布页| 久久精品欧美一区二区| 国产在线自在拍91精品黑人| 免费看久久精品99| 国产欧美日韩18| 亚洲欧美在线综合图区| 91精品国产91久久久久久三级| 国产精品福利在线观看无码卡| 欧美黄色网站在线看| 亚洲综合18p| 成年午夜精品久久精品| 97av视频在线观看| 精品久久久久久久久久久| 亚洲中久无码永久在线观看软件| 国产丰满大乳无码免费播放|