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

基于傅里葉級數法的開孔板振動固有特性分析

2017-08-05 01:37:08王旻昊李凱邱永康李天勻朱翔
中國艦船研究 2017年4期
關鍵詞:振動分析

王旻昊,李凱,邱永康,李天勻,3,4,朱翔,3

1華中科技大學船舶與海洋工程學院,湖北武漢 430074 2中國艦船研究設計中心,湖北武漢430064 3船舶與海洋水動力湖北省重點實驗室,湖北武漢 430074 4高新船舶與深海開發裝備協同創新中心,上海 200240

基于傅里葉級數法的開孔板振動固有特性分析

王旻昊1,李凱2,邱永康1,李天勻1,3,4,朱翔1,3

1華中科技大學船舶與海洋工程學院,湖北武漢 430074 2中國艦船研究設計中心,湖北武漢430064 3船舶與海洋水動力湖北省重點實驗室,湖北武漢 430074 4高新船舶與深海開發裝備協同創新中心,上海 200240

[目的]開口板結構普遍存在于各類工程結構中,對其振動特性的研究直接關系到整體結構的減振降噪和穩定性分析。為研究針對彈性薄板在任意位置開與板平行的矩形口的自由振動特性研究問題,[方法]通過改進傅里葉級數形式表示開口矩形板的位移容許函數,用區域劃分思想將開口板沿開口延伸線劃分為多個區域板,采用沿邊界均勻分布的線性模擬彈簧模擬經典邊界條件和區域板間連續邊界條件,將邊界表達為彈性勢能的形式,從而將有約束問題轉化為無約束問題,并結合位移連續條件和能量泛函變分方法,對未知傅里葉展開系數一次變分求極值以求解標準特征值方程。然后將得到的開口矩形板的固有頻率值及其對應振型與有限元軟件(ANASYS)計算結果進行對比,最后分析不同邊界條件、開口尺寸和開口位置對開口板自振特性的影響。[結果]結果驗證了方法的有效性和精確性,[結論]所得結果可為相關實際工程應用提供理論參考。

任意開口位置;開口矩形板;改進傅立葉級數;能量法

0 引 言

開口矩形板結構廣泛存在于各個行業中,例如汽車、航空航天、土木建筑等。在船舶領域,為提高隱身性能,方式之一是采用一體化上層建筑,雷達布置在上層建筑結構開口處。結構開口會影響其原有振動特性,存在與開口處動力設備共振的風險。所以開口矩形板的動態性能研究對船舶上層建筑結構的減振設計具有重要的理論指導意義和工程實用價值。

關于開口板結構的動態性能,國內外學者開展了大量的研究。Paramasivam[1]通過擴展網格模型對有限差分算子進行改進,研究了不同邊界條件下開口對板固有頻率的影響;Cho等[2]基于假定模態法,利用拉格朗日運動方程求解多自由度系統矩陣的特征值,分析了不同形狀開口板的自由振動問題;Lam等[3]對含有開口或裂紋的矩形板結構的自由振動性能進行研究分析,對開口板進行分區處理,以正交多項式模擬板位移函數得出中心開口矩形板的振動頻率,其為用板分區處理思想研究開口板問題提供了新的思路;李成等[4]根據非均質各向異性彈性理論,對含橢圓孔的正交各向異性板的孔邊進行應力分析,提出了積分方程法求解方案;李自林等[5]利用三角形薄板廣義協調元分析了復合式多層四邊簡支、四邊固定無孔和開孔矩形薄板的自由振動,求出了前幾階頻率系數并與ANSYS軟件計算結果進行了對比;Reddy[6]根據 Reissner-Mindlin 型剪切變形理論和非線性的von Kármán應變——位移關系理論,對開口板大振幅彎曲振動進行了研究;邱昌林等[7]采用有限元與間接邊界元相結合的方法,以開有圓孔、四邊簡支、無障板的鋼制平板為對象,開展了開孔板水下振動及聲輻射特性研究;Hegarty等[8]通過最小二乘配點法,對中心開圓孔的矩形板的自由振動問題進行了研究,最后得出板自由振動頻率與開口尺寸之間的關系曲線;Aksu等[9]結合有限差分技術的變分原理,對含有1~2個開口矩形板的固有頻率進行了預測;Tan[10]運用一種易于處理和封閉的形式,對含橢圓形開口的各向異性和正交異性板的有限寬度修正因子進行研究進而分析了其性能;張宇力等[11]選取典型船舶板架,采用ANSYS系統,對開口和不開口的板分別進行了特征值屈曲分析和極限承壓屈曲分析,同時,運用便于簡化處理的相當梁系法,將加筋板架近似看作由帶附連翼板的交叉梁系構成,轉換以后給出了一套工程近似計算方法。

上述文獻中對開口板邊界條件、開口尺寸和開口位置對其整體振動性能的研究較少。文獻[3]中以正交多項式形式表示開口板的位移容許函數,使得最后得出結果的精度完全依賴于所選取板的假定振型與實際情況的吻合程度,不僅影響計算效率,而且其只能針對中心開口的對稱結構,限制了開口位置的多樣性。

Li[12]提出了任意支撐梁振動分析的改進傅里葉級數方法,并隨后被拓展應用到矩形板[13]和圓柱殼[14]等結構的振動分析之中;文獻[13]表明任意邊界條件下板的假定振型函數可以不變地用一種改進傅里葉級數形式表示。本文將引入改進傅里葉級數方法建立任意邊界條件下開口矩形板的振動分析模型,并應用區域劃分思想將開口在任意位置的開口板劃分為8塊區域板,采用沿邊界均勻分布的位移約束彈簧和轉角約束彈簧模擬板的邊界條件,然后運用能量泛函變分方法對結構振動問題進行求解,并與有限元仿真運算結果進行對比分析以說明文中方法的準確性,最后討論邊界條件、開口尺寸和開口位置對開口板振動特性的影響,以便為工程應用提供理論指導。

1 理論分析

1.1 開口矩形板模型描述

文中研究對象為任意位置矩形開口的矩形板,如圖1所示。開口矩形板的長度為a,寬度為b,矩形開口位于圖中所示位置,其角點坐標分別為(a1,b1)和(a2,b2),板厚為h。考慮到開口板結構和邊界的任意性,將開口板劃分為8塊區域板進行研究。

圖1 開口矩形板示意圖Fig.1 Rectangular plate with a rectangular opening

在圖1中(1)~(12)邊界及(21)~(24)邊界均為對地邊界,(13)~(20)為區域板間的連接邊界。以上邊界均可采用兩類沿邊界均勻分布的線性彈簧來模擬,兩類線性彈簧分別為位移約束彈簧和旋轉約束彈簧,通過設置這兩類彈簧的剛度來模擬開口板任意邊界。例如,對于對地邊界,在自由邊界條件下,邊界位移和轉角均無約束,將兩類彈簧剛度均設為0,則可獲得自由邊界;在固支邊界條件下,邊界位移和轉角均為0,將兩類彈簧剛度都設為∞,則可模擬固支邊界;在簡支邊界條件下,邊界位移為0,轉角無約束,將位移約束彈簧和轉角約束彈簧的剛度值分別取0和∞,則可模擬經典簡支邊界;當其取0~∞時,可模擬任意彈性邊界等。對于區域板間的連接邊界,其位移和轉角均連續,區域板間位移和轉角均相對為0,故兩類模擬彈簧值均設置為∞。

當模擬彈簧剛度設置為∞時,本文的取值為1×1012,由后文算例驗證,當剛度取值為 1×1012,一般彈性邊界可完全退化為經典邊界。根據以上彈簧模擬方法,可得到開口板的物理模型如圖2所示。

圖2 開口板物理模型描述圖Fig.2 Physical model of rectangular plate with an opening

根據上述物理模型針對每個單獨的區域板進行單獨研究,每個區域板均取獨立的坐標系,區域矩形板的位移容許函數可表示為[3]

式中:Amn為其未知傅里葉級數展開系數;簡諧時間因子eiωt表示矩形板在不同時刻的位移函數;j=1,2,3,…,8,分別表示8塊板的位移函數;M,N為截斷項取值;?m(x)為x方向的特性梁容許函數,ψn()

y為y方向的特性梁容許函數。所假定的位移容許函數的優劣將直接影響到計算結果的準確性,因而選取合適的位移容許函數非常重要。文中引入改進傅里葉級數形式來表示開口矩形板的位移容許函數,該級數形式具有收斂性好、精度高等特點,可以滿足任意邊界條件,改進傅里葉級數方法可使板的位移容許函數在整個求解域內三階導數連續且四階導數各點均存在,可以有效克服邊界處可能出現的不連續現象。根據改進傅里葉級數方法,式(1)中?m(x)與ψn(y )可分別表示為[13]:

式中:x*,y*為無因次坐標;m=1,2,3,…,M;n=1,2,3,…,N。

1.2 含開孔板振動固有特征的能量分析

結合能量泛函變分法對開口矩形板結構進行研究,以確定式(1)中的未知系數。僅考慮彎曲變形,在無外激勵力作用的情況下,開口板的能量泛函為8塊區域板的總和,表示為

式(4)中區域板的彎曲應變能可表示為

在任意邊界條件下,8塊區域板的動能可表示為

式中:ρ為板材料的密度;j=1,2,3,…,8。

儲存在邊界模擬約束彈簧中的彈性勢能為(此處以區域板1為例,其余板同理可得):

式中,ki和 Ki分別表示圖1中開口板第i個邊界的位移約束彈簧剛度(單位:N/m)和旋轉約束彈簧剛度(單位:(N·m)/rad)。

針對式(1)中區域板的總能量泛函進行變分求極值可得:

式中,j=2,3,…,7。

將式(5)~式(7)代入式(8)中,可寫成矩陣形式為

式中:矩陣K表示8個區域板的總勢能剛度矩陣,包括板的彎曲應變能剛度和儲存在邊界中的彈簧勢能剛度;M表示其質量矩陣;ω為圓頻率。

式中,j=1,2,3,…,8。

通過求解式(9)的特征值和特征向量,即可得到開口矩形板的各階固有頻率及其對應的模態振型,至此,將開口板的自由振動問題通過能量泛函變分法處理成求解標準特征值方程。

2 數值計算與分析

本文對任意邊界條件下,開口在任意位置的開口矩形板的自由振動特性進行分析計算,對比不同邊界、不同開口大小和不同開口位置對開口板固有頻率的影響,并與有限元軟件(ANSYS)計算結果進行了對比分析,以說明本文方法的精確性。板材料參數的取值為:ρ=7 850 kg/m3,μ=0.3,E=2.1×1011Pa。

2.1 收斂性及有效性驗證

首先進行收斂性分析,以驗證隨著式(1)中截斷項M,N的取值和邊界模擬彈簧剛度取值的不斷增大,計算結果的穩定性。在對截斷項M,N的取值進行對比分析中,開口板的邊界條件設置為C-F邊界條件,其中C代表其外邊界的邊界條件為固支,F代表其內孔的邊界條件為自由。板的幾何參數設置如下:a=8 m,b=6 m,h=0.02 m,內開口參數 a1=2 m,b1=1 m,a2=5 m,b2=3 m,計算其前6階的固有頻率值(單位:Hz),并與有限元軟件計算結果進行對比分析,為便于分析,這里取M=N,結果如表1所示。

表1 C-F邊界下開口板固有頻率Table 1 Natural frequencies of rectangular plate with a rectangular opening in C-F boundary

從表1中的對比分析可知:隨著截斷項數M,N不斷增加,采用文中方法所求得的開口板固有頻率逐漸趨于穩定,當M=N=12時,固有頻率已不再發生變化,即認為文中方法已收斂,從而證明了改進傅里葉級數方法在計算開口板振動特性上的收斂性和穩定性。且文中方法計算所得開口板固有頻率值與有限元仿真軟件(ANSYS)計算結果吻合很好,后面的計算中均取截斷項M=N=13。

提取前4階固有頻率對應的模態振型,并與有限元仿真軟件(ANSYS)進行對比,如圖3所示。從圖3中可以看出,計算所得頻率和振型與有限元方法所得結果的吻合度很高,從而證明了文中方法的正確性。

圖3 開口板前4階振型對比圖Fig.3 Comparison of the first four order modes of rectangular plate with a rectangular opening

然后,針對邊界模擬彈簧剛度取值對計算結果的影響進行分析。算例參數不變,邊界條件內邊為自由邊界,外邊位移約束彈簧剛度為K′,轉角約束彈簧剛度為0,計算隨著K′的增大,采用文中方法計算所得開口矩形板的固有頻率的穩定性,并與有限元軟件(ANSYS)計算結果進行對比分析,表2給出前6階固有頻率(單位:Hz)結果的對比。

表2 不同位移約束彈簧剛度K′的開口板固有頻率Table 2 Natural frequencies of rectangular plate with a rectangular opening in differentK′

從表2中的數據對比可以得出,位移約束彈簧剛度從0逐漸增大到1×1012N/m時,開口板的外邊邊界條件從自由(F)邊界演化為簡支(SS)邊界,證明了采用文中方法邊界模擬彈簧剛度取值的收斂性。故而當邊界模擬彈簧剛度為∞時,文中算例可取值為1×1012N/m。

2.2 不同邊界條件、開口尺寸和開口位置對開口板振動性能的影響分析

不同邊界條件、開口尺寸和開口位置對開口板振動性能具有直接影響,對比分析其固有頻率值,總結影響規律,可以為工程實際應用提供重要的理論依據。

首先,對比研究不同邊界條件對開口板振動特性的影響。所取邊界均為經典邊界組合,包括F-F,SS-F,C-F,SS-SS和C-C這5種邊界。開口板的幾何參數設置如下:a=8 m,b=6 m,h=0.02 m,內開口參數 a1=1 m,b1=2 m,a2=3 m,b2=4 m,計算其前6階的固有頻率值(單位:Hz),對比結果如圖4及表3所示。

從圖4中可以看出,隨著邊界約束的增大,開口板的同階固有頻率值也逐漸增大,這是因為邊界約束的增大導致開口板整體剛度的上升,從而引起其固有頻率值的增大。因為C-F邊界與SS-SS邊界約束自由度的不同,導致了SS-SS邊界的首階固有頻率小于C-F邊界,其余值均大于C-F邊界這一現象的出現。

研究不同開口尺寸對開口板振動特性的影響。在SS-SS和F-F邊界條件下,設置a=8 m,b=6 m,h=0.02 m,開口中心與板的中心重合,開口邊與外邊的長度比對比分析當 ξ=0.1,0.3,0.5,0.7,0.9時,開口板前6階的固有頻率值,SS-SS和F-F邊界條件下的對比結果分別如圖5和圖6所示。

圖4 不同邊界下開口板前6階固有頻率對比圖Fig.4 Natural frequencies of rectangular plate with square opening in different boundaries

表3 SS-SS邊界條件下不同開口尺寸開口板固有頻率表Table 3 Natural frequencies of rectangular plate with different rectangular openings in SS-SS boundary

圖5 SS-SS邊界下開口板前6階固有頻率對比圖Fig.5 Natural frequencies of rectangular plate with different rectangular openings in SS-SS boundary

圖6 F-F邊界下開口板前6階固有頻率對比圖Fig.6 Natural frequencies of rectangular plate with different rectangular openings in F-F boundary

由圖5和圖6可知,在SS-SS邊界下,開口板的前6階固有頻率值是隨開口尺寸的增加而不斷增大,而在F-F邊界下,則是隨開口尺寸的增大不斷減小。開口尺寸的增加降低了結構的質量,同時也減小了結構的剛度,但因在SS-SS邊界時結構固有頻率以剛度影響為主,質量影響為輔,在F-F邊界下,開口板邊界約束小,結構剛度較小,結構剛度影響弱于質量影響,故而引起上述現象。

最后,討論不同開口位置對開口板自由振動特性的影響。在C-F邊界下,a=5 m,b=5 m,h=0.02 m,內開口大小為1 m×1 m,考慮到矩形板的對稱性,文中研究開口中心沿圖7所示虛對角線移動,對比位置參數 η=1,1.3,1.6,1.9,2.2,2.5這6種情況下開口板固有頻率的變化,提取前3階固有頻率計算結果對比,如圖8所示。

由圖8可以看出,隨著開口位置向開口板的中心的靠近,開口板首階固有頻率值不斷增大,第2階逐漸減小,第3階頻率值呈拋物線型。

圖7 開口位置分布描述圖Fig.7 The opening position of rectangular plate

圖8 C-F邊界下開口板前3階固有頻率對比圖Fig.8 Natural frequencies of rectangular plate with different opening positions in C-F boundary

3 結 論

本文通過引入改進傅里葉級數形式表示開口板的假定振型函數,采用區域板劃分的思想將開口矩形板劃分為8塊區域,通過邊界模擬彈簧解決邊界問題以及區域板之間的位移連續性問題,最后結合能量泛函變分法將開口板自由振動問題轉化為一個求解標準特征值問題。求解任意邊界條件下任意開口尺寸和開口位置開口板的固有頻率值及其相應振型,通過算例對比,分析不同邊界條件、開口尺寸和開口位置對開口板自振特性的影響并總結規律。得出以下結論:

1)隨著傅里葉級數截斷項以及模擬彈簧剛度值的增加,計算結果收斂性良好,數值穩定性很好。

2)文中方法對于處理任意邊界條件下開口在任意位置的開口矩形板的自由振動問題具有較高的精度。

3)隨著邊界約束剛度的增加,開口板的同階固有頻率不斷增大;在邊界約束剛度較大時,振動特性以剛度影響為主,開口板的同階固有頻率隨開口尺寸的增加而增大,在邊界約束剛度較小時,振動特性以質量影響為主,開口板的同階固有頻率隨開口尺寸的增加而減小;在C-F邊界下,開口板的首階固有頻率隨著開口位置向板中心的靠近而增大。

[1]PARAMASIVAM P.Free vibration of square plates with square openings[J].Journal of Sound and Vibra?tion,1973,30(2):173-178.

[2]CHO D S,VLADIMIR N,CHOI T M.Approximatenatural vibration analysis of rectangular plates with openings using assumed mode method[J].Internation?al Journal of Naval Architecture and Ocean Engineer?ing,2013,5(3):478-491.

[3]LAM K Y,HUNG K C.Vibration study on plates with stiffened openings using orthogonal polynomials and partitioning method [J].Computers& Structures,1990,37(3):295-301.

[4]李成,鄭艷萍,王迎佳.積分方程求解復合材料開口板的應力分布[J].玻璃鋼/復合材料,2007(1):9-12.LI C,ZHENG Y P,WANG Y J.Stress analysis of com?posite plate with holes by integral function[J].Fiber Reinforced Plastics/Composites,2007(1):9-12(in Chinese).

[5]李自林,王榮霞,劉興業.用廣義協調元分析復合式多層無孔和開孔矩形薄板的振動[J].地震工程與工程振動,2004,24(2):64-68.LI Z L,WANG R X,LIU X Y.Vibration analysis of complex layered rectangular thin plate with or without holes by generalized conforming element[J].Earth?quake Engineering and Engineering Vibration,2004,24(2):64-68(in Chinese).

[6]REDDY J N.Large amplitude flexural vibration of lay?ered composite plates with cutouts[J].Journal of Sound and Vibration,1982,83(1):1-10.

[7]邱昌林,陳志剛,鄧軼等.開孔平板水下振動及聲輻射特性[J].中國艦船研究,2013,8(6):75-80.QIU C L,CHEN Z G,DENG Y,et al.The characteris?tics of vibration and sound radiation of underwater per?forated plates[J].Chinese Journal of Ship Research,2013,8(6):75-80(in Chinese).

[8]HEGARTY R F,ARIMAN T.Elasto-dynamic analy?sis of rectangular plates with circular holes[J].Interna?tional Journal of Solids and Structures, 1975,11(7/8):895-906.

[9]AKSU G,ALI R.Determination of dynamic character?istics of rectangular plates with cutouts using a finite difference formulation[J].Journal of Sound and Vibra?tion,1976,44(1):147-158.

[10]TAN S C.Finite-width correction factors for anisotro?pic plate containing a central opening[J].Journal of Composite Materials,1988,22(11):1080-1097.

[11]張宇力,曾廣武,段洪,等.大開口對船舶板架穩定性和極限承載力的影響[J].華中科技大學學報(自然科學版),2002,30(5):56-58,64.ZHANG Y L,ZENG G W,DUAN H,et al.The ef?fect of big placket on stability and terminal loading in marine structure[J].Journal of Huazhong University of Science and Technology (Nature Science Edi?tion),2002,30(5):56-58,64(in Chinese).

[12]LI W L.Free vibrations of beams with general bound?ary conditions[J].Journal of Sound and Vibration,2000,237(4):709-725.

[13]LI W L.Vibration analysis of rectangular plates with general elastic boundary supports[J].Journal of Sound and Vibration,2004,273(3):619-635.

[14]DAI L,YANG T J,DU J T,et al.An exact series so?lution for the vibration analysis of cylindrical shells with arbitrary boundary conditions[J].Applied Acoustics,2013,74(3):440-449.

Free vibration characteristics analysis of rectangular plate with rectangular opening based on Fourier series method

WANG Minhao1,LI Kai2,QIU Yongkang1,LI Tianyun1,3,4,ZHU Xiang1,3
1 School of Naval Architecture and Ocean Engineering,Huazhong University of Science and Technology,Wuhan 430074,China 2 China Ship Development and Design Center,Wuhan 430064,China 3 Hubei Key Laboratory of Naval Architecture and Ocean Engineering Hydrodynamics,Wuhan 430074,China 4 Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration,Shanghai 200240,China

Plate structures with openings are common in many engineering structures.The study of the vibration characteristics of such structures is directly related to the vibration reduction,noise reduction and stability analysis of an overall structure.This paper conducts research into the free vibration characteristics of a thin elastic plate with a rectangular opening parallel to the plate in an arbitrary position.We use the improved Fourier series to represent the displacement tolerance function of the rectangular plate with an opening.We can divide the plate into an eight zone plate to simplify the calculation.We then use linear springs,which are uniformly distributed along the boundary,to simulate the classical boundary conditions and the boundary conditions of the boundaries between the regions.According to the energy functional and variational method,we can obtain the overall energy functional.We can also obtain the generalized eigenvalue matrix equation by studying the extremum of the unknown improved Fourier series expansion coefficients.We can then obtain the natural frequencies and corresponding vibration modes of the rectangular plate with an opening by solving the equation.We then compare the calculated results with the finite element method to verify the accuracy and effectiveness of the method proposed in this paper.Finally,we research the influence of the boundary condition,opening size and opening position on the vibration characteristics of a plate with an opening.This provides a theoretical reference for practical engineering application.

arbitrary opening position;rectangular plate with an opening;improved Fourier series method;energy variational method

U661.44

A

10.3969/j.issn.1673-3185.2017.04.016

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

王旻昊,李凱,邱永康,等.基于傅里葉級數法的開孔板振動固有特性分析[J].中國艦船研究,2017,12(4):102-109.

WANG M H,LI K,QIU Y K,et al.Free vibration characteristics analysis of rectangular plate with rectangular opening based on Fourier series method[J].Chinese Journal of Ship Research,2017,12(4):102-109.

2017-01-03< class="emphasis_bold">網絡出版時間:

時間:2017-7-27 10:30

國家自然科學基金資助項目(51379083,51579109)

王旻昊,男,1992年生,碩士生。研究方向:結構振動噪聲,機械結構靜、動力學仿真分析。

E-mail:710959552@qq.com

李天勻(通信作者),男,1969年生,教授,博士生導師。研究方向:結構振動噪聲分析。

E-mail:ltyz801@hust.edu.cn

猜你喜歡
振動分析
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
隱蔽失效適航要求符合性驗證分析
This “Singing Highway”plays music
電力系統不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
電力系統及其自動化發展趨勢分析
中西醫結合治療抑郁癥100例分析
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 中国国产A一级毛片| 99草精品视频| 亚洲成AV人手机在线观看网站| 亚洲成人黄色在线| 亚洲91在线精品| 日本免费精品| 久久中文无码精品| 国产免费福利网站| 亚洲性日韩精品一区二区| 日本精品视频| 国产美女无遮挡免费视频| 久久成人免费| 亚洲欧美自拍视频| 国产激情无码一区二区免费| 99这里精品| 久久综合结合久久狠狠狠97色| 波多野结衣久久高清免费| 亚洲丝袜中文字幕| 亚洲第一成年人网站| 喷潮白浆直流在线播放| 欧美国产中文| 91精品日韩人妻无码久久| 久久免费精品琪琪| 亚洲国产精品无码AV| 3344在线观看无码| 四虎AV麻豆| 国产在线自乱拍播放| 国产天天色| 国产情侣一区| 99精品国产自在现线观看| 亚洲国产日韩在线成人蜜芽| 找国产毛片看| 国产精品开放后亚洲| 中字无码精油按摩中出视频| 欧美中出一区二区| 91国内视频在线观看| 99在线视频免费| 国产日产欧美精品| 午夜a级毛片| 精品国产成人国产在线| 国产一级毛片yw| 看你懂的巨臀中文字幕一区二区| 国产拍在线| 尤物在线观看乱码| 精品久久777| 国产视频久久久久| 香蕉网久久| 成人免费视频一区二区三区| 国产精品成人观看视频国产 | 国产亚洲精品无码专| 国产一二三区在线| 九九九精品视频| 无码专区国产精品一区| 97国产在线观看| 国产精品成| 一级成人a毛片免费播放| 欧美性精品| 亚洲第一黄色网| 四虎永久在线精品影院| 国产91色| 亚洲aⅴ天堂| 国产男女免费完整版视频| 全部免费特黄特色大片视频| 国产永久在线视频| 国产精品自拍合集| 成年人久久黄色网站| 亚洲AV无码精品无码久久蜜桃| 免费看av在线网站网址| 免费国产不卡午夜福在线观看| 久久精品中文字幕免费| 亚洲男人在线| 午夜视频www| 97国产在线视频| 国产午夜不卡| 99精品影院| 日本欧美午夜| 国产白浆在线观看| 中文字幕无码中文字幕有码在线| 色综合色国产热无码一| 超清无码熟妇人妻AV在线绿巨人| 永久天堂网Av| 色爽网免费视频|