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

渦輪葉片雙層干摩擦阻尼器減振特性研究

2022-11-21 03:39:44滿吉鑫邱欣可
振動與沖擊 2022年21期
關鍵詞:模型系統

滿吉鑫, 高 慶, 曾 武, 邱欣可

(1.中國科學院 工程熱物理研究所先進燃氣輪機實驗室,北京 100190;2.中國科學院 先進能源動力重點實驗室(工程熱物理研究所),北京 100190;3.中國科學院 輕型動力創新研究院,北京 100190;4.中國科學院大學 工程科學學院,北京 100049)

燃氣輪機葉片系統工作狀態下受到氣流激勵的作用,該激勵往往具有很寬的頻率范圍,當葉片自身的固有頻率與其相一致時,就會產生共振。共振的存在往往會導致葉片的高周疲勞失效,要想完全避免共振幾乎是不可能。而系統固有的結構阻尼對于渦輪葉片的減振貢獻很小[1-2],為了減小共振應力及其帶來的危害,發展了一種附加的阻尼器結構,干摩擦阻尼器,依靠接觸面之間的干摩擦作用消耗振動能量,達到減小振動的目的。

干摩擦接觸普遍認為是一個復雜的非線性接觸過程,建立合適的摩擦接觸模型對于研究緣板干摩擦阻尼器的振動行為以及預測葉片的響應具有十分重要的作用。Griffin[3]對渦輪葉片系統中的干摩擦阻尼結構進行了綜述,其中重點介紹了干摩擦阻尼建模以及計算的方法。通過將庫侖摩擦元件與彈簧相串聯的方式,提出了一種具有恒定法向載荷的一維整體滑動接觸模型。許多學者也在其模型的基礎上進行了改進,進而研究了恒定法向載荷下的干摩擦阻尼作用[4-6]。Menq等[7]研究了干摩擦阻尼系統中摩擦接觸面法向載荷可變時對于系統減振特性的影響。但是此種模型僅在接觸面法向正壓力較小時得到較好的應用,當接觸面法向正壓力過高時,無法對接觸面狀態進行詳細描述。為了解決接觸面正壓力過大時整體滑動模型的局限性,準確描述接觸面上應力分布不均的問題,Iwan等[8-10]提出了一種彈簧串聯、并聯模型,并在之后的研究中得到了廣泛的應用;Menq等[11-12]基于此種模型發展了一種接觸面壓力較大時的微滑移模型,利用此種模型對渦輪葉片緣板干摩擦阻尼減振過程進行了研究,結果表明,當接觸面法向正壓力較高時,此種模型對接觸面的描述更加準確,并在此模型的基礎上進行了相關試驗,證明了此種模型在高法向正壓力下的有效性;何尚文等[13]在此基礎上,提出一種改良的微滑移摩擦模型,對簡化的B-G型葉片干摩擦緣板阻尼器減振特性進行了研究,同時分析了不同受力情況下接觸面上摩擦力的分布。漆文凱等[14]提出了一種不考慮摩擦接觸面之間完全黏滯階段的解析微滑移模型。Zhang等[15]采用考慮微滑動過程的分布式摩擦力學模型,得到緣板阻尼關鍵參數對系統振動響應的影響規律。本文在漆文凱等研究的基礎上,利用Masing Rule[16]對微滑移過程進行了解析求解,將局部滑動過程拓展到整體滑動,進而建立了整體-局部統一滑動摩擦模型。依據所建立的摩擦模型,推導了阻尼器中摩擦力與位移的解析關系,對摩擦力-位移遲滯曲線進行了仿真。針對現有某型號燃氣輪機中的雙層阻尼結構,即緣板-阻尼器、凸起-阻尼器兩層阻尼結構進行簡化,建立雙層干摩擦阻尼器有限元模型,并發展帶摩擦阻尼系統的響應分析方法,對不同位置處干摩擦阻尼器作用過程中系統的振動特性進行分析,研究阻尼器設計中的關鍵參數。

1 渦輪葉片雙層阻尼結構介紹

圖1為某型號燃氣輪機的雙層阻尼結構模型圖。可以看到:與一般的緣板阻尼結構不同,該結構在緣板凹槽位置①處存在一處干摩擦阻尼裝置,當葉片旋轉時,依靠自身的離心力使得阻尼器摩擦接觸面壓緊,產生摩擦阻尼作用,從而耗散能量;在榫根靠上的位置②處做了凸起處理,施加了另一處干摩擦阻尼裝置,該位置處的阻尼結構不僅起到了摩擦減振的作用,同時能夠依靠后部位置處的封嚴結構進行氣體封嚴。

圖1 雙層阻尼結構模型圖Fig.1 Double friction damper model

阻尼裝置表面的正壓力是影響阻尼器阻尼效果的關鍵參數之一,而針對本文中出現的雙層干摩擦阻尼結構,阻尼裝置正壓力的來源為旋轉時自身離心力,大小如式(1)所示

(1)

2 整體-局部統一滑動模型的建立及仿真

2.1 整體-局部統一滑動模型

為了模擬葉片緣板與阻尼器接觸點之間的摩擦界面,應用了如圖2所示的微滑動摩擦模型。該模型源自Menq等對于Iwan模型的改進,Csaba[17]在此基礎上做了部分簡化,但簡化后的模型仍能夠完整的顯示微滑動摩擦界面的重要特性。摩擦界面是通過將彈性矩形板壓在半無限大的剛性表面上來模擬的,矩形板AB的長度為l,在其上表面施加有均布載荷q,u為平板右端的位移距離,右側F為施加的外激勵,當阻尼塊處于局部滑動階段時,外激勵大小等于摩擦阻尼力的大小。E,A分別為板的彈性模量和AB的橫截面積,接觸區域內摩擦因數μ為常數。假設接觸面上的摩擦作用符合庫侖摩擦定律,當右端外激勵F出現時,右端首先開始滑動。滑動的區域隨著F的增大而逐漸增大,直到滑動區域拓展到整體,發生整體滑動,摩擦力表示為:f(x)=μq(x)。其滑動區域分布如圖2所示。

圖2 局部滑動Fig.2 Partial slip

設滑動區域的長度為Δ,滑動區域上單位長度的摩擦力為τ,外激勵F為簡諧激勵,做周期性變化。在滑動區域內任選一段微元體進行分析,如圖3所示。

圖3 微元體分析Fig.3 Differential analysis

微元體中受力平衡為

(2)

彈性體的應變為

(3)

因此微元體右端受力表示為

(4)

將式(4)代入式(2)可得

(5)

當F從0開始逐漸增大時,B點開始滑動,出現了滑動區域且滑動區域自右向左逐漸拓展,此時摩擦力的大小與外激勵大小相等。摩擦接觸過程滿足以下方程及邊界條件

(6)

由式(6)可以得到位移和力之間的函數關系為

(7)

彈性矩形板末端的力可以通過對摩擦力進行積分求得,從而可以表示為滑動區域長度的函數

(8)

將式(8)代入式(7)中可得

(9)

因此,矩形板末端的力和位移以及他們之間的函數關系可以分別表示為

(10)

(11)

當F繼續增大,Δ=l,F=μql時,滑動區域拓展到左端A點,系統達到臨界滑動,隨后開始整體滑動,這時的滑動摩擦力為f=F=μql。

當矩形塊達到整體滑動階段時,可以對式(11)得到的力和位移關系曲線進行拓展,以允許在滑動過程中出現整體滑動階段。假設矩形板右端受到的力為Fmax,此時矩形板達到最大滑動位移umax,之后力開始卸載,當每次矩形塊的位移方向改變時,都會從微滑動階段開始,直到當Δ=l時,微滑動階段結束,整體滑動開始。界面的整體滑動位移由一個變量u′決定,這個變量的含義是矩形板的整體位移與產生整體滑動臨界位移u0之間的差值

u′=umax-u0

(12)

外力F是簡諧變化的,針對反復的卸載以及再加載過程,為了建立一個振動周期內的摩擦力-位移穩態遲滯曲線,需要對卸載以及再加載過程進行研究。根據Masing Rule的研究,材料的周期性加卸載過程中獲得的遲滯回線路徑與單調加載具有相同的形式,只是將其擴展了兩倍,并且進行映射,從而創建了閉合的對稱形式的遲滯回曲線。該假設的具體數學關系表示為

(13)

式中:un為卸載位移,是外力F的函數,是由臨界位移u0減去函數M的兩倍得出;F0為達到臨界滑動時對應的摩擦力。函數M為單調加載下力與位移之間的函數關系。再加載函數ur如式(14)所示

ur(F)=-un(-F)

(14)

利用式(13)、式(14)兩式可以得到一個完整加載周期中的力和位移關系曲線即摩擦過程的遲滯回曲線。

2.2 摩擦接觸仿真計算

利用之前得到的力和位移關系式進行仿真計算,繪制不同滑動狀態下的遲滯回曲線。給定矩形板的長度=0.015 m,EA=3.6×106N·m,q=50 000 N/m,μ=0.3。μ為滑動摩擦因數,正壓力為N=ql=500 N。仿真得到的摩擦力-位移遲滯曲線如圖4所示。

圖4 摩擦力-位移遲滯曲線仿真Fig.4 Simulation of friction-displacement hysteresis curve

圖4中分別包含了當Δ=0.5l即局部滑動狀態下的遲滯回曲線,當Δ=1.0l時系統達到臨界狀態時的初始加載曲線及遲滯回曲線,當Δ=1.5l時系統處于整體滑動狀態時的遲滯回曲線。從式(10)可知,當Δ=l時,系統達到臨界狀態,此時外力F=μql=225 N。

3 振動響應求解

干摩擦阻尼系統具有非線性特性,為了便于在頻域內計算葉片的響應特性,在本文所應用的摩擦模型前提下,阻尼裝置和緣板之間始終保持接觸,即摩擦接觸面沒有顯著的法向分離特征,高階諧波的成分影響較小,采用一次諧波即可以滿足對于此工程性動力學問題的動力響應特性的近似計算需求[18-19],且其計算效率較高,便于將其引入具有高自由度的復雜結構系統有限元模型。利用白鴻柏等[20]所提到的線性化方法,基于等效黏性阻尼的思想,用等效橢圓代替摩擦過程中力與位移形成的遲滯回曲線。每個循環中阻尼器實際消耗的能量為W=∮udF。

當系統處于局部滑動狀態時,摩擦力視為具有復剛度的彈簧,k=keq+iωceq。由等效線性化方法可以得到等效黏性阻尼功,由式(15)給出

W=πωcequ2

(15)

式中:ω為激振角頻率;u為當前位移值。

求出等效黏性阻尼系數為

(16)

等效剛度為

(17)

當系統處于整體滑動狀態時,利用等效線性化和一次諧波平衡法相結合可以得到

(18)

其中,式(18)中的第一個式子的分子為發生整體滑動時一個周期內的阻尼功,包括遲滯回曲線中平行四邊形的面積與局部滑動時的面積。u0為臨界位移值,kd為發生臨界滑動時的等效剛度值。

振動的響應求解過程利用了多程序協調求解,因目前的有限元無法獨立的處理干摩擦過程中涉及的非線性問題,所以需開展分析程序的二次開發。結合ANSYS的APDL參數化語言以及MATLAB軟件實現對結構附加雙層干摩擦阻尼器后振動響應的求解:基于有限元軟件實現對結構的有限元建模以及線性系統的諧響應求解,并可以輸出結構特定位置處,如雙層阻尼結構作用處的單元位移等結果。基于摩擦滑動模型以及滑動狀態判斷理論表達式,在MATLAB中進行自主編程,求解對應滑動狀態下兩處干摩擦阻尼結構的等效剛度和等效阻尼并傳遞給ANSYS。在ANSYS的諧響應分析模塊中作為邊界條件將等效阻尼和等效剛度附加在有限元模型中,再次進行迭代計算。以特定位置處單元位移前后計算結果誤差小于給定誤差時為迭代計算結束,最后計算結果為真實位移以及兩處干摩擦阻尼結構所提供的等效阻尼以及等效剛度。基本流程如圖5所示。

圖5 多程序協調求解計算流程圖Fig.5 Multi-program coordination calculation process

4 有限元算例分析

為了保留實際葉片中重點研究部位的特征,同時便于進行機理研究,建立了與實際葉片特征相似的有限元模型,如圖6所示。①位置處代表上層阻尼器與緣板位置處的接觸,②位置處代表下層阻尼器與凸起之間的接觸。

圖6 簡化有限元模型Fig.6 Simplified finite element model

利用本文第2章中的響應求解方法對存在①,②兩層緣板阻尼結構的穩態響應進行分析。因干摩擦阻尼裝置實質上提供了系統的等效剛度和等效阻尼,所以在ANSYS軟件中利用通用性較好的Matrix27單元來模擬摩擦過程中的等效剛度和等效阻尼。對線性系統進行諧響應分析可以得到①,②阻尼結構特定位置處單元節點的位移,將得到的位移應用于所建立的摩擦模型,可以進行兩個位置處的等效剛度和等效阻尼的求解。將得到的結果利用Matrix27單元,作為邊界條件再次代入線性系統進行諧響應分析,迭代求解,直到系統響應滿足收斂條件后,得到真實的位移以及等效剛度和等效阻尼。阻尼塊和葉片材料均為結構鋼,其具體參數為:彈性模量E=2×1011Pa,密度ρ=7 850 kg/m3,泊松比ε=0.3,摩擦因數為0.2。忽略角度等影響因素的條件下,將阻尼塊簡化為矩形塊,其結構參數為厚度2 mm,長度40 mm,寬度30 mm。載荷及邊界條件為葉片根部固支,外激勵施加在葉片中部,振動響應拾取點為葉尖。

4.1 等效剛度和等效阻尼

圖7、圖8分別為①,②兩處阻尼器作用時,某一固定頻率下產生的等效剛度和等效阻尼隨接觸面正壓力的變化情況。可以看到,雙層阻尼器不同位置處的等效剛度和等效阻尼均隨著正壓力的增大而增大;其次,在振動過程中①阻尼所提供的等效剛度和等效阻尼均大于②阻尼,由式(10)、式(16)、式(17)可知,這與滑動區域長度的變化有關。

圖7 等效剛度隨正壓力變化曲線Fig.7 Variation curve of equivalent stiffness with normal force

圖8 等效阻尼隨正壓力變化曲線Fig.8 Variation curve of equivalent damping with normal force

4.2 正壓力對系統減振效果的影響

阻尼裝置的阻尼特性與作用在其上的正壓力大小有關,針對本文所涉及的雙層阻尼裝置,對正壓力的探討分為如下兩個方面:第一部分對上下兩層阻尼器施加相同的正壓力,通過計算結果對比分析不同位置處阻尼裝置的阻尼效果;第二部分選取某特定轉速工況,按照與真實模型相同的壓力分布比例進行加載,判斷雙層阻尼裝置存在的合理性以及其共同作用下系統振動特性的變化規律。

對上下兩層阻尼裝置施加以相同的正壓力,施加于葉片中部的外激勵大小相同,激勵大小為F=1 N。圖9~圖11分別為僅存在①阻尼、②阻尼以及存在①和②雙層阻尼時系統的幅頻曲線。頻率計算范圍在130~180 Hz,可以看到阻尼裝置的存在主要對葉片峰值響應及出現峰值響應時的頻率兩個方面產生影響。3種阻尼情況下都存在這樣的現象:隨著接觸面正壓力的增大,共振峰值呈現先減小后增大的趨勢;系統的共振頻率呈現增大的趨勢,當增大到某值時,阻尼器的存在只起到改變剛度的作用且剛度保持為一恒定值,因此共振頻率不發生變化。

從3種不同阻尼情況下的對比可以看到,阻尼器的位置變化會影響到系統的共振頻率與最優正壓力。當存在②阻尼時,最優正壓力為200 N且共振頻率在140~150 Hz變化;僅存在①阻尼時,最優正壓力為60 N且共振頻率在140~170 Hz變化,調節系統共振頻率的范圍相對更大;同時存在①和②阻尼時,最優正壓力為60 N,共振頻率在140~170 Hz變化,但最優正壓力對應的共振峰值有所下降。證明①和②阻尼同時存在時,在特定的正壓力下阻尼效果優于僅存在單層阻尼。

選取圖9~圖11中的振動峰值及對應頻率繪制圖12,葉片受到1 N激振力下,不同阻尼作用時系統的振幅、頻率隨接觸面正壓力的變化曲線。可以看到,僅存在②阻尼時,最優正壓力的值相對于僅存在①阻尼和存在雙層阻尼時較大且最優正壓力對應的振幅也更高,其單獨存在時阻尼效果最差;壓力-振幅特性曲線中僅存在①阻尼和存在雙層阻尼的趨勢大致相似,但是雙層阻尼對應最優正壓力和振動幅值都更小,體現出更好的阻尼效果;不同阻尼作用狀態時系統的共振頻率隨正壓力的增大而增大到某一定值,僅②阻尼時共振頻率變化范圍相對更小,僅為6 Hz。

圖12 不同阻尼作用時正壓力-頻率/振幅曲線Fig.12 Normal force-frequency/amplitude curve under different damping

在本文的算例中①阻尼位置處的最優正壓力為60 N,小于②阻尼位置處的最優正壓力200 N,意味著在相同的轉速和半徑條件下,①阻尼位置也即相對遠離榫根處使用質量更小的阻尼裝置就可以達到最佳的減振效果。這為雙層阻尼器的設計提供了指導。

針對本文所建立的有限元模型,選取部分能夠突出規律特性的轉速工況,如第1章所述按照與真實模型相同的壓力分布比例對上下兩層阻尼裝置進行正壓力加載。施加于葉片中部的外激勵大小相同,激勵大小為F=1 N。圖13為不同轉速工況下雙層阻尼同時作用時系統的幅頻特性曲線。可以看到,在轉速達到2 500 r/min時,雙層阻尼系統的減振效果達到最佳。

圖13 ①和②阻尼同時存在時幅頻曲線Fig.13 Amplitude-frequency curve both ① and ② damping

類比圖12,將相同轉速下不同位置處的共振峰值及對應頻率繪制如圖14所示。可以看到,雙層阻尼的存在可以調節系統的共振頻率,使得系統在更高的轉速下達到更優的減振效果;僅②阻尼作用時系統趨于穩定的共振頻率以及最優減振效果對應的轉速均最小;從離心力形成的角度分析,同一轉速時②阻尼產生的離心力更大,但因其位置相對更靠下,產生的貢獻更小,因此②阻尼質量變化對于系統的振動特性的影響更小。

圖14 不同阻尼作用時轉速-頻率/振幅曲線Fig.14 Rotational velocity-frequency/amplitude curve under different damping

4.3 外激勵對于系統減振效果的影響

燃氣輪機工作時,隨著工作狀態的變化,葉片所受到氣流帶來的激振力也會發生變化,因此研究了不同激振力下葉片的響應特性。針對本文有限元模型算例,為了突出雙層阻尼器不同位置阻尼作用時的規律特點,選取0.5 N,1 N和2 N 3個值作為外激勵條件進行計算分析。

圖15為不同阻尼作用下,接觸面正壓力保持為一恒定值時,不同外激勵下的葉片穩態響應幅頻曲線。由圖15可知,當系統受到的外激勵增大時,系統的共振頻率有所減小,即出現了剛度軟化現象,反之共振頻率會有所增加。不同形式的阻尼作用時,激振力變化趨勢呈現的規律是一致的。

圖15 不同激振力下3種阻尼狀態時的幅頻曲線Fig.15 Amplitude-frequency curve of three damping states under different exciting forces

圖16為不同阻尼作用狀態,接觸面正壓力保持為一定值時系統的振動峰值隨外激勵的變化情況。可以看到,僅②阻尼作用時,系統的振動峰值在外激勵大于2 N時產生突變。僅①存在時,在6 N之后才發生振幅突變,而①和②阻尼都存在時,振幅發生突變對應的外激勵值更大。

圖16 不同阻尼作用狀態時振幅峰值隨外激勵變化曲線Fig.16 Variation curve of peak amplitude with external excitation under different damping states

對于外激勵變化的系統來說,僅②阻尼存在時系統峰值響應變化更加敏感;①和②雙層阻尼的存在使得系統可以在更大的外激勵變化范圍內保持較小的峰值變化。

4.4 滑動摩擦因數對于系統減振效果的影響

干摩擦阻尼結構在工作時,由于接觸面之間的摩擦,在長時間的工作后,難以避免產生磨損使得表面狀態發生變化從而影響摩擦因數。因此對于摩擦因數對系統影響的研究也是十分有必要的。

取外激勵的大小為1 N,分別做出僅存在①阻尼、僅存在②阻尼以及存在①和②阻尼時系統的峰值響應隨接觸面法向正壓力變化的正壓力-振幅曲線,摩擦因數分別取0.2,0.3,0.4,如圖17所示。

由圖17可知,曲線整體上呈現先減小后增大的趨勢,證明了每一個摩擦因數下,系統都存在一個最優的正壓力,能夠使得系統的響應最低。同時,可以看到摩擦因數的改變對系統響應的影響主要體現在兩個大的方面:一是隨著摩擦因數的逐漸增大,系統的最優正壓力逐漸減小,這是積極的一面,意味著系統在保持一定轉速的條件下,最優正壓力時對應的阻尼塊質量較小;二是隨著摩擦因數的增大,當系統處于最優正壓力狀態時,峰值響應有所上升,即系統響應的最大降幅減小,這是消極的一面。

圖17 不同摩擦因數下3種阻尼狀態時的正壓力-振幅曲線Fig.17 Normal force-amplitude curves of three damping states under different friction coefficients

圖18所示為接觸面正壓力60 N,僅存在②阻尼時系統的幅頻曲線。可以看到,在一定正壓力下摩擦因數的改變對共振頻率的影響幾乎可以忽略;而對系統的響應幅值則有較大的影響,共振以及附近區域,摩擦因數越大,減振效果越好;遠離共振區域時,減振效果不隨摩擦因數的增加而出現明顯的變化。

圖18 不同摩擦因數下幅頻曲線Fig.18 Amplitude-frequency curve under different friction coefficients

5 結 論

通過本文研究,得到以下結論:

(1) 葉片系統中,阻尼裝置的存在主要起到提供等效阻尼和等效剛度的作用,從而影響系統共振頻率和峰值;存在①和②兩層阻尼作用的葉片減振系統中,遠離葉根固定約束也即①阻尼對于系統的減振特性影響程度更大,①位置處使用質量更小的阻尼裝置就能達到更優的減振效果。

(2) 隨著阻尼裝置接觸面正壓力的增大,不同阻尼狀態時系統共振頻率均逐漸增大且共振峰值呈現先減小后增大的趨勢,存在一個最優正壓力使得系統減振效果最好。僅存在②阻尼時正壓力的變化使得共振頻率變化范圍更小,但系統最優正壓力以及相應振幅均較高,且②處阻尼質量變化對于系統振動特性影響程度更小。

(3) 阻尼裝置接觸面正壓力保持為一定值,隨著外激勵的增加,不同阻尼作用狀態時,系統出現峰值頻率減小也即剛度軟化的現象;①和②雙層阻尼的存在使得外激勵在較大的范圍內變化時,系統振幅峰值都處于較小范圍內。

(4) 摩擦因數增大,不同阻尼作用狀態時系統的最優正壓力減小,但最優正壓力對應的峰值響應有所上升;正壓力保持為一恒定值時,摩擦因數增大,共振區域內減振效果增加,遠離共振區域的減振效果變化較小。

猜你喜歡
模型系統
一半模型
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
基于PowerPC+FPGA顯示系統
半沸制皂系統(下)
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
3D打印中的模型分割與打包
主站蜘蛛池模板: 91亚瑟视频| 欧美成人怡春院在线激情| 日韩在线成年视频人网站观看| 国产乱子伦手机在线| 国产精品久久久久无码网站| 欧美www在线观看| 国产成人高精品免费视频| 婷婷丁香在线观看| 欧美精品一区在线看| 少妇露出福利视频| 精品一区二区三区视频免费观看| 亚洲精品欧美重口| 黄色污网站在线观看| 亚洲国产精品无码久久一线| 成人免费视频一区| 人妻少妇乱子伦精品无码专区毛片| 69国产精品视频免费| 国产精品人成在线播放| 手机在线国产精品| 制服丝袜一区| 毛片网站在线播放| 中文字幕一区二区视频| 欧美不卡视频在线| 91小视频版在线观看www| 国产午夜精品一区二区三| 一本大道香蕉久中文在线播放| 国产99视频精品免费视频7| 99久久婷婷国产综合精| 亚洲视频免| 中文无码毛片又爽又刺激| 免费国产小视频在线观看| 亚洲三级片在线看| 午夜精品久久久久久久无码软件| 99视频精品在线观看| 在线亚洲精品福利网址导航| 午夜高清国产拍精品| 免费国产一级 片内射老| 日韩欧美国产成人| 91麻豆国产在线| 99色亚洲国产精品11p| 国产成人超碰无码| 99热这里都是国产精品| 亚洲熟妇AV日韩熟妇在线| …亚洲 欧洲 另类 春色| 黄色网站在线观看无码| 71pao成人国产永久免费视频| 国产极品美女在线播放| 成人一区在线| 亚洲日韩高清无码| 天天操精品| 亚洲国产精品VA在线看黑人| 日韩欧美国产中文| 午夜无码一区二区三区| 国产日韩欧美一区二区三区在线| 国产精品亚洲一区二区在线观看| 老司机精品久久| 国产草草影院18成年视频| 久久这里只有精品23| 国产成人精品亚洲日本对白优播| 亚洲欧美综合精品久久成人网| 99精品国产自在现线观看| 精品一区二区三区中文字幕| 国产aaaaa一级毛片| 欧美另类图片视频无弹跳第一页 | 国产不卡在线看| 亚洲欧美h| 国产爽爽视频| 女同久久精品国产99国| 午夜a级毛片| 中文字幕色在线| 亚洲Va中文字幕久久一区| 成人在线亚洲| 国产麻豆精品在线观看| 伊人成人在线视频| 中文字幕 日韩 欧美| 97在线公开视频| 日韩东京热无码人妻| 亚洲欧美综合在线观看| 亚洲午夜福利精品无码不卡 | 99久久国产综合精品2020| 日韩国产一区二区三区无码| 亚洲AⅤ永久无码精品毛片|