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

靜壓作用對吸聲覆蓋層性能的影響與分析

2023-07-03 07:22:58陳文炯周祥超
船舶 2023年3期
關鍵詞:有限元變形結構

陳文炯 盧 辰 周祥超

(大連理工大學 船舶工程學院 大連 116023)

0 引 言

隱身性是艦艇和潛艇等武器裝備的重要特征之一,敷設在結構表面的吸聲覆蓋層是實現水下聲隱身的重要技術途徑之一[1]。自二戰時期Alberich型吸聲覆蓋層[2]問世以來,空腔型覆蓋層憑借多吸聲機理協同作用的優勢,廣泛應用于現代海軍潛艇的吸聲減振設計。空腔型覆蓋層由橡膠基體與內嵌空腔結構所構成,通過橫波與縱波之間的波形轉換、基體材料對于聲能的損耗及空腔諧振等特性,實現高頻段下的聲能吸收特性,從而降低回聲和目標強度,提升聲隱身性能[3]。然而,由于現代水下聲吶探測的工作頻率不斷向低頻移動,航行器下潛深度也不斷增加,導致覆蓋層表面受到高靜水壓力的影響。因此,對靜壓工況下開展吸聲覆蓋層的相關研究有重要意義。

數十年來,國內外學者對空腔型覆蓋層開展了大量研究。HENNION[4]利用有限元法建立了周期性吸聲覆蓋層的分析模型,研究了空腔型覆蓋層的反射和透射特性。姚熊亮等[5]與商超等[6]基于水中空腔結構體有限元方法,計算了含混合橢圓柱形空腔等復合腔型聲學覆蓋層的吸聲特性,并將聲學覆蓋層應用于雙殼動力艙段,通過實驗驗證了所設計聲學覆蓋層在艇體結構中的降噪性能。陳文炯等[7]通過一維聲學模型,計算出空腔結構的等效介質參數,結合傳遞矩陣法和遺傳優化算法對復合聲學覆蓋層的空腔結構進行設計與優化,設計了指數振蕩型空腔形狀,并討論其特有的吸聲特性。YE 等[8]采用等效參數法,討論了典型空腔在常壓下的聲學特性。葉韓峰等[9-10]建立了平面波斜入射的單層和多層含空腔型覆蓋層理論和有限元模型,研究了聲波在斜入射條件下的吸聲覆蓋層結構,以及材料參數對吸聲性能的影響。

目前雖然已有上述這些對于吸聲覆蓋層內空腔形狀的研究工作,但多數研究并沒有考慮到靜壓對于覆蓋層結構與吸聲性能的影響。覆蓋層受到靜水壓力后,基體材料與空腔形狀會發生變形并且結構內存在預應力,導致結構的動態力學性能發生變化,從而影響覆蓋層的吸聲性能。目前對靜壓下空腔型吸聲覆蓋層的研究主要是考慮變形對于吸聲性能的影響。汪慧銘[11]與董文凱等[3]應用有限元法對靜壓下的聲學覆蓋層吸聲性能進行研究,將仿真步驟分為兩步:一是求解靜水壓力下的結構變形量;二是將變形量導出構件新有限元模型進行吸聲性能仿真,并計算靜壓對于3 種不同空腔結構吸聲性能的影響。然而,靜水壓力對于基體材料與空腔的變形是無規律的,對靜壓下的覆蓋層進行二次建模的誤差也難以避免,因此楊立軍等[12]基于COMSOL 中的“移動網格”模塊,直接將靜壓下所計算的變形后模型應用于吸聲性能的計算中。

可以看出:目前基于靜壓對覆蓋層影響的研究主要考慮了空腔和基材的幾何變形,但尚未考慮靜壓作用下(即含有預應力狀態)引起的吸聲性能變化。而預應力的作用會引起結構剛度等參數的變化并導致固有頻率發生改變,從而引起吸聲性能發生改變。因此,模擬靜壓狀態下進行覆蓋層吸聲性能的分析有著重要意義。

本文利用有限元方法,建立了靜壓作用狀態下吸聲覆蓋層性能分析模型。該模型考慮靜壓引起的空腔結構變形和預應力的作用,討論了圓柱形、圓錐形和喇叭形這3 類典型空腔型吸聲覆蓋層的吸聲性能在靜壓狀態下的變化;對比了不同空腔型覆蓋層在靜壓保持狀態下的最大變形量、頻段平均吸聲系數和吸聲系數變化量;獲得了靜壓對不同空腔型覆蓋層及覆蓋層不同頻段的影響結果;比較了靜壓下圓柱形空腔覆蓋層特定頻率處的振動幅值;最后,討論了靜壓變化對空腔型覆蓋層吸聲性能產生影響的原因。

1 理論與數值模型

吸聲覆蓋層受靜壓作用產生變形,且處于預應力作用狀態,現有的解析解無法準確預報該情況,因此采用有限元數值模型來進行分析。文中的公式主要描述的是吸聲覆蓋層吸聲的機理以及有限元方程,是由前人所提出,而本文的重點是將這些方法和公式運用于靜壓下吸聲覆蓋層吸聲性能的計算中。

1.1 水下吸聲覆蓋層基本參數

設定空腔型覆蓋層在x和y方向都無限大,并且空腔在xy平面上周期性分布,有限元分析時只取含單個空腔的部分。

吸聲覆蓋層整體結構與單胞結構示意圖見圖1。

圖1 吸聲覆蓋層整體結構與單胞結構示意圖

如圖 1 (a) 所示,在z方向,計算區域劃分為3 個部分:區域1 為半無限空間的水域,區域2 包含空腔型覆蓋層、背襯潛艇鋼層與有限厚度的流體層,區域3 包括覆蓋層背襯殼體內部的空氣層。選擇區域2 中的1 個單胞作為計算單胞,該單胞包含三維有限結構及其聲場區域,其中S-和S+分別對應聲波的入射和透射界面。

如圖 1 (b)所示,以流體外邊界面S-和S+將入射端和透射端半無限空間劃分為兩部分,與結構相鄰的部分分別為水層和空氣層。入射端和透射端聲壓梯度的等效節點載荷分別為Φ-和Φ+。在入射端,Φ-可以分解為入射波場和反射波場這兩部分的累加,見式(1):

令入射聲波為單位幅值的簡諧平面波,設入射角與z軸的夾角為θ、與x軸的夾角為φ,不考慮波形隨時間的變化,則入射波聲壓如式(2)所示:

為簡化計算模型,計算中僅考慮入射聲波為沿z軸的垂直入射波,且幅值為1 Pa,故式(2)可以簡化為式(3):

反射波表示為x和y方向,波數分量分別為knx和kmy的各階平面波的疊加,因此反射波聲壓可表達為式(4)形式:

式中:Rpnm是相對于入射波聲壓幅值為單位值時,第(n,m)階反射波的聲壓幅值,對于垂直入射僅考慮(0, 0)階波即可;knm為z方向的波數;knx=nπ/d1,kmy=mπ/d2;d1、d2如圖1(a)所示。

從能量角度出發,可用各階波的反射或透射系數的均方根,求得遠場反射系數和透射系數[8],即:

式中:N=(2Mx+1)(2My+1);Mx、My分別為M矩陣的行列階數;n=-Mx∶1∶Mx,n=-My∶1∶My。

對應吸聲覆蓋層的吸聲系數為:

通常,當聲波垂直入射時,反射系數和透射系數分別等于(0, 0)階波的反射系數和透射系數。

1.2 有限元方程

針對吸聲覆蓋層的有限元分析,應對其聲-固耦合物理場進行分析,其本質是固體有限元方程與流體有限元方程。

1.2.1 固體有限元方程

固體結構振動的有限元形式為:

式中:KS和MS分別為固體結構的整體剛度矩陣和質量矩陣,N/mm 和kg;δ為節點位移,mm;Fm、Fp分別為節點所受機械激勵的等效節點載荷和流體對結構作用的等效節點載荷,即靜壓和入射波聲壓,N。

KS和MS由單元剛度矩陣和單元質量矩陣集成產生,即:

式中:ρs為材料密度,kg/m3;Nδ為位移插值形函數矩陣;D為彈性矩陣,即彈性體應力與應變關系的矩陣;Bδ為應變矩陣, 表示應變和應力的關系。

1.2.2 流體有限元方程

理想流體介質中,聲波的波動方程矩陣形式為:

式中:c為流體介質中的聲速,m/s;p為聲壓,Pa;t為時間,s。

據此對流體區域及其邊界進行有限單元離散化,引入流體聲單元剛度矩陣、質量矩陣、流體聲單元與結構單元的耦合矩陣,并分別集成為對應的整體矩陣來實現單元求和運算。最終,得到如式(11)所示流體有限元方程:

式中:Kp和Mp分別為流體的剛度矩陣和質量矩陣;R為覆蓋層結構和水的耦合矩陣;p為結點壓力,N;Φ為流體外邊界上聲壓梯度的法向等效節點載荷,N。

將流體有限元方程與結構有限元方程聯立,得到結構與聲場耦合的有限元方程:

式中:ω為頻率,Hz。

1.3 模型有效性驗證

基于上述數值計算模型,本文利用COMSOL軟件實現了靜壓下的空腔型吸聲覆蓋層數值仿真。為便于研究,在保證孔隙率不變的情況下,三維棱柱體單胞可以等效成圓柱體單胞單元,進而使用二維軸對稱模型代替三維模型,降低計算成本。

空腔型覆蓋層有限元模型以及有限元計算流程圖見圖2。

圖2 空腔型覆蓋層有限元模型以及有限元計算流程圖

圖3 常壓與靜壓情況下的計算結果與文獻結果對比

圖4 計算結果收斂性驗證

空腔型覆蓋層的二維單胞模型及相應邊界條件如圖 2 (a)所示。其中水層上邊界的平面波輻射即為入射的平面聲波,模型對稱邊界即對應周期性邊界條件。圖 2 (b) 為該有限元模型數值計算的關鍵流程圖。

本模型利用移動網格模塊進行靜壓下的變形傳遞,將結構的變形情況傳遞到頻域進行計算,避免因二次建模而導致的誤差。為了將靜壓作用下的結構響應傳遞至吸聲計算,在同一個研究中必須包含2 個計算步驟:計算預應力的穩態計算步驟以及計算聲固耦合的頻域擾動計算步驟。此處利用“穩態-頻域擾動”這一研究類型將預應力狀態傳遞至聲固耦合的頻域計算中,對于平面波則采用Linper 函數來模擬線性擾動。

為驗證該靜壓下有限元計算模型的有效性,本文基于此計算方法復現了文獻[12]中含橢球形空腔吸聲覆蓋層在常壓和靜壓下的結果。該覆蓋層的結構尺寸和材料參數均參照原文中設定,其中橡膠材料的參數為根據實驗數據擬合得到的表達式,并且隨頻率變化而變化。常壓與靜壓情況下所得吸聲系數計算結果與文獻結果對比見圖 3。

由圖 3(a)可見,有限元仿真所得的常壓下的計算結果以及考慮了靜壓引起變形的計算結果,均與文獻[12]中的數據有很好的一致性,對比所得最大差值小于0.06,故可證明該建模方式的正確性。此外,本文構建的模型不僅可以考慮靜壓引起結構幾何變形,并且能夠實現考慮預應力影響的吸聲性能計算。

由圖 3(b)可見,靜壓作用下考慮預應力作用對吸聲結果的影響較大,2 000 Hz 頻率點以上的頻段尤為明顯。

由上述驗證及對比結果可知,該有限元計算模型滿足靜壓下的吸聲覆蓋層計算要求。該模型在考慮了變形的基礎上同時考慮到預應力的影響,這對靜壓下空腔型覆蓋層吸聲計算來說是必要的。

1.4 收斂性驗證

通過對圖2(a)建立有限元模型,并分別計算了當最大有限元網格為最小波長1/5、1/6、1/7、1/8、1/9 和1/10 時的吸聲覆蓋層吸聲系數,得到如圖 4 所示的驗證結果。

可見,在不同網格尺寸情況下,計算結果的誤差非常小。因此,本文通過最小波長1/5 對有限元模型進行網格劃分,既可保證計算結果的收斂性,又可減少計算的時間成本。

2 靜壓保持狀態下的吸聲覆蓋層吸聲性能分析

基于上述計算,本章節首先分析靜壓對于基體材料動態力學性能的影響,并進行頻響分析;而后研究靜壓保持狀態下的空腔型覆蓋層吸聲性能;最后討論了3類典型空腔型覆蓋層靜壓作用下的吸聲性能變化。

2.1 靜壓對基體材料的力學性能影響

本節僅考慮不同靜壓對于含有空腔的基體材料力學性能的影響,模型各層結構的基本尺寸和材料參數見表1。

表1 模型各層的基本尺寸和材料參數

此外,還對覆蓋層表面分別施加5 MPa 和10 MPa靜壓,以進行頻響分析,從而計算出不同靜壓下,空腔壁面Z方向位移幅值分量,見圖 5。

由圖5 可見:隨著靜壓增加,空腔壁面的位移也顯著增加。并且對于該結構而言,不同靜壓下,在3 800 Hz 頻率點處,結構空腔壁面的位移均達到峰值,而此時結構內部的應力也達到最大值。

圖5 不同靜壓下的空腔壁面Z 方向位移幅值分量

圖6 給出了5 MPa 和10 MPa 靜壓下,3 800 Hz頻率點處的結構內部應力云圖。

圖6 3 800 Hz 頻率點處覆蓋層應力云圖

由圖6 可見:隨著靜壓增加,覆蓋層內部應力也逐漸增加,并集中于空腔壁面處。這使得在進行覆蓋層吸聲性能計算前,結構內部就存在很大的預應力,不僅使空腔結構發生變形,而且由于空腔內部梯度內應力的存在,使得基體材料分配不均并導致材料屬性發生改變。例如基材的密度在應力較大處增大,應力較小處減小;并且結構固有頻率發生改變,使空腔的諧振特性也產生變化,從而影響吸聲性能。

2.2 3類典型空腔型覆蓋層靜水作用狀態下的吸聲性能

本節根據不同的計算工況(即0 MPa、5 MPa和10 MPa 靜壓工況),對3 類典型空腔型覆蓋層構型進行計算及相關定量分析。

2.2.1 3 類典型空腔構建

圖7 為3 類常見典型空腔型覆蓋層結構示意圖。

圖7 典型空腔型覆蓋層結構示意圖

圖8 3 類空腔型覆蓋層吸聲系數曲線對比圖

圖9 3 類空腔類型分別在不同靜壓下的吸聲系數曲線對比圖

圖10 常壓下第一吸聲峰處的振動幅值

圖11 不同靜壓下的圓柱形空腔吸聲系數與振動幅值

將這3 類空腔分別按照相應的參數化公式進行設計與建模,各參數見表2。

表2 不同空腔類型的參數表

由表2 可知,只要確定p和q的值,就能分別確定3 類覆蓋層的空腔幾何尺寸。本文假設3 類空腔的體積相等(即模型的孔隙率相同),調整p和q參數值便可得到3 類空腔結構,對應尺寸設定見表3。

表3 不同空腔類型的尺寸表mm

2.2.2 3 類典型空腔型覆蓋層靜壓下的吸聲性能對比

首先比較3 類覆蓋層在相同靜壓下的性能差異,這里靜壓取5 MPa 和10 MPa,計算得到3 類空腔型覆蓋層在靜壓下的最大變形量,見表4。

表4 3類空腔型覆蓋層在靜壓下最大變形量對比mm

最大變形量即覆蓋層結構受到靜壓的垂向恒力作用后產生的最大變形,反映了空腔型覆蓋層的耐壓性能。由表4 對比結果可知:相同靜壓下,圓柱形空腔覆蓋層的最大變形量最大,其次是喇叭形空腔覆蓋層,圓錐形空腔覆蓋層的最大變形量最小。

這是因為空腔形狀不同,導致空腔表面應力分布存在差異;不同位置的材料變形量不同,使得材料的屬性不均勻。此外,從圖 6 可知,靜壓對于圓柱形空腔所產生的應力主要集中于空腔下半部分,而截面梯度變化的空腔類型(圓錐形空腔、喇叭形空腔)可以減緩應力集中效應。

下頁圖 8 為3 類空腔型吸聲覆蓋層在相同靜壓下的吸聲系數曲線對比。由該圖可見:在頻率2 000 Hz 以下時,靜壓對于這3 類空腔型覆蓋層吸聲性能的影響較小。在常壓下,這3 類空腔型覆蓋層的吸聲系數曲線頻率在2 000 Hz 以上時有明顯差異,然而隨著靜壓的增加,差異逐漸減小,其中圓柱形空腔與圓錐形空腔所對應的吸聲系數曲線在靜壓10 MPa 時幾乎重合。在靜壓10 MPa 時,圓柱形空腔的變形量最大并且空腔表面存在壓力梯度,使圓柱形空腔被擠壓成了圓錐形空腔,因此,此時的圓柱形空腔與圓錐形空腔有著相似的吸聲性能。靜壓的增加使3 類覆蓋層吸聲系數曲線變得更為陡峭,曲線在頻率6 000 Hz 處有明顯提升。

而后,本文將繼續研究相同空腔在不同靜壓下的吸聲性能。下頁圖 9 為常壓0 MPa、靜壓5 MPa 和10 MPa 下,3 類吸聲覆蓋層在100 ~6 000 Hz 頻段內的吸聲曲線。由該圖分析靜壓對覆蓋層各頻段的吸聲系數曲線的影響,規定吸聲系數變化量為覆蓋層在靜壓下和常壓下對某頻率聲波的吸聲系數差值,具體如式(13)所示。

式中:αF和α0分別為靜壓下和常壓下的吸聲系數。

3 類空腔型覆蓋層的吸聲系數曲線受靜壓影響的變化趨勢一致,靜壓增加使吸聲系數曲線第一吸聲峰向低頻方向移動。在100 ~ 2 000 Hz 頻段內,3 類覆蓋層的吸聲系數均隨靜壓增大而增加;之后一段頻率內,均存在吸聲系數隨靜壓增大而減少的趨勢;但在更高頻率處,吸聲系數仍隨靜壓增大而增加。

3 類空腔型覆蓋層在特定頻段處吸聲系數變化量較大。其中靜壓10 MPa 下,圓柱形空腔覆蓋層700 Hz 和6 000 Hz 頻率處的吸聲系數較常壓時顯著提高,變化量分別為0.05 和0.11;在3 700 Hz 頻率處的吸聲系數大幅下降,變化量約為0.06。對于圓錐形空腔覆蓋層,靜壓10 MPa 下的吸聲系數在700 Hz 和3 700 Hz 頻率處的變化量分別約為0.05和0.06。喇叭形空腔覆蓋層的吸聲系數曲線同樣反映了該趨勢,靜壓10 MPa 下在上述兩頻率時的變化量分別約為0.07 和0.08。由此可見,靜壓對特定頻段的吸聲系數具有較大影響。

針對特定頻段分析靜壓影響吸聲效率的原因。以圓柱形空腔覆蓋層為例,其常壓下吸聲峰值(2 300 Hz 頻率)處的振動幅值如下頁圖 10 所示。

由該圖可見,在此頻率處,覆蓋層基體部分的位移從下至上逐漸增大,覆蓋層發生了整體的拉伸或壓縮,并且基體位移大于空腔位移,由此證明該峰值與空腔共振無關,而是由于覆蓋層結構和背襯結構耦合共振所致。

最后,分析靜壓下振動幅值的變化。振動幅值即結構在受到聲波擾動后所產生的微小振動幅度,圖 11 反映了在不同靜壓下圓柱形空腔的吸聲系數曲線及對應振動幅值曲線。

由圖可知,靜壓的增加使該頻段(100 ~3 000 Hz)的“基體-背襯”系統的位移幅值發生改變,即改變了系統的固有振動幅值和對應頻率,因而引起了該頻段吸聲性能的變化,其中2 000 Hz 頻率以內的振動幅值隨靜壓增加而增大,由此也解釋了靜壓下該頻段吸聲性能略有提升的原因。

3 結 論

本文基于靜壓作用狀態構建了吸聲覆蓋層有限元計算模型,靜壓作用會使結構內部產生變形與預應力。由于現有研究僅考慮了變形對吸聲覆蓋層吸聲性能的影響,因此通過與相關文獻結果進行對比,發現在僅考慮變形以及考慮變形-預應力這2 種不同情況下,覆蓋層吸聲性能有明顯差異,從而證明了考慮預應力的必要性。

本文計算了圓柱形空腔覆蓋層在不同靜壓下空腔表面的位移幅值與結構內部應力狀態,給出了3 種不同的計算工況(即0 MPa、5 MPa 和10 MPa靜壓工況),并進行了定量分析。發現隨著靜壓的增加,不僅位移增加,而且結構內部點的應力差異更加明顯,這會使基體材料分配不均,繼而導致材料參數與結構固有頻率發生改變,從而影響吸聲性能。文中以3 類典型空腔型覆蓋層為例,比較其在靜壓下的變形以及所選頻段吸聲系數的變化情況;以覆蓋層最大變形量、各頻段平均吸聲系數以及吸聲系數變化量為指標進行分析,量化3 類典型空腔型覆蓋層在靜壓下的性能變化;分析了圓柱形空腔覆蓋層受靜壓作用后,特定頻段的結構振動幅值和頻率的改變直接影響特定頻段的吸聲性能,并最終得出以下結論:

(1)在相同的孔隙率下,靜壓使得圓柱形空腔覆蓋層的最大變形量最大,其次是喇叭形空腔覆蓋層,而圓錐形空腔覆蓋層的最大變形量則最小;靜壓對于3 類覆蓋層在2 000 Hz 以下頻率時的影響較小;隨著靜壓的增加,圓柱形空腔和圓錐形空腔的吸聲系數曲線逐漸靠近,在靜壓10 MPa 時幾乎重合。因為在靜壓10 MPa 時的圓柱形空腔變形量最大,并且空腔表面存在壓力梯度,使得圓柱形空腔被擠壓成了圓錐形空腔,故此時圓柱形空腔與圓錐形空腔有著相似的吸聲性能。

(2)對于不同空腔型覆蓋層,靜壓增加使得吸聲系數曲線第一吸聲峰向低頻方向移動。吸聲系數在100 ~ 2 000 Hz 頻段處,吸聲性能略有提高;在2 000 ~ 5000 Hz 頻段處,吸聲性能隨著靜壓增加而逐步降低,但在更高頻段仍隨壓力增加而提高。這是由于靜壓使結構內部發生了整體的拉伸或壓縮,覆蓋層結構與背襯結構發生共振,改變了系統固有振動幅值和對應頻率,從而引起吸聲性能的變化。

以上結論主要適用于考慮靜壓作用(包括覆蓋層的形變和內應力)時,計算吸聲覆蓋層的吸聲性能。

猜你喜歡
有限元變形結構
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
“我”的變形計
例談拼圖與整式變形
會變形的餅
論《日出》的結構
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 亚洲无线一二三四区男男| 久久久91人妻无码精品蜜桃HD| 999精品在线视频| 免费国产无遮挡又黄又爽| 精品久久蜜桃| 国产一级小视频| 制服丝袜 91视频| 日韩不卡高清视频| 91久久国产成人免费观看| 激情综合婷婷丁香五月尤物| 久久国产黑丝袜视频| 2024av在线无码中文最新| 国产精品成人一区二区不卡| 蝌蚪国产精品视频第一页| 国产 日韩 欧美 第二页| 国产喷水视频| 国产白浆一区二区三区视频在线| 欧美高清三区| 婷婷丁香色| 2021最新国产精品网站| 欧美区日韩区| 国产精品短篇二区| 中文字幕亚洲无线码一区女同| 亚洲人成电影在线播放| 成人毛片免费观看| 992tv国产人成在线观看| 中文字幕佐山爱一区二区免费| …亚洲 欧洲 另类 春色| 色成人综合| 国产91蝌蚪窝| 精品偷拍一区二区| 欧美日一级片| 久热这里只有精品6| 片在线无码观看| 精品人妻无码区在线视频| 亚洲区一区| 99精品国产自在现线观看| 国产精品一区在线麻豆| 久久免费视频播放| 亚洲欧洲日产国码无码av喷潮| 欧美成人看片一区二区三区| 国产无遮挡裸体免费视频| 无码一区中文字幕| 99视频在线看| 波多野结衣亚洲一区| 无码在线激情片| 一本大道无码日韩精品影视| 综合社区亚洲熟妇p| 中文字幕色站| 欧美在线国产| www.狠狠| 国产一级视频久久| 在线免费看片a| 久久黄色一级片| 日本高清免费不卡视频| 久久网欧美| 中文字幕日韩久久综合影院| 婷婷久久综合九色综合88| 亚洲精品另类| 在线观看精品国产入口| 亚洲香蕉伊综合在人在线| 爆乳熟妇一区二区三区| 国模视频一区二区| 国产精品主播| 国产成人精品无码一区二| 国产极品美女在线观看| 免费国产小视频在线观看| 中文字幕久久亚洲一区| 色爽网免费视频| 国产精品对白刺激| 2020国产在线视精品在| 性色一区| 亚洲欧美另类视频| 性色一区| 全部免费特黄特色大片视频| 久久大香香蕉国产免费网站| 日本国产精品一区久久久| 思思99思思久久最新精品| 国产特级毛片aaaaaa| 狼友av永久网站免费观看| 国产色婷婷视频在线观看| 精品视频一区二区观看|