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

航天器太陽帆板多自由度減振裝置設計

2021-02-06 11:24:02朱仕堯雷勇軍
振動與沖擊 2021年3期
關鍵詞:模態振動系統

朱仕堯,雷勇軍,郭 欣

(國防科技大學 空天科學學院,長沙 410073)

現代航天器上外伸的太陽帆板會引起兩類振動問題。①殘余振動:太陽帆板具有頻率低、阻尼小等特點,導致其在航天器姿態調整和進出陰影產生的振動響應需要較長時間衰減,影響航天器姿態穩定時間和快速響應能力[1];②驅動擾動:為盡可能多獲取能量,太陽帆板需要不斷旋轉以持續面向太陽,由此引起的連續擾動會影響高精度航天器指向精度和敏感設備的動態環境[2]。

一方面,殘余振動屬于自由振動問題,發生在帆板面外方向,主要由彎曲模態振動引起。抑制殘余振動的目的是為了加快衰減速度,縮短振動穩定所需的時間,可以從被動和主動兩個方面開展減振研究。哈勃空間望遠鏡阻尼器[3]是太陽帆板被動減振方法的典型應用,其利用阻尼層剪切變形耗散振動能量。劉超等[4]用聚合物阻尼材料替代阻尼器中的黏彈性阻尼材料,并證明模態阻尼比隨溫度升高而增大。Jia等[5-7]分別研究了平臺式、斜桿支撐式和質量調諧式太陽帆板阻尼器設計。在主動減振方面,Weck等[8-12]以嵌入太陽帆板支架和基板內部形狀記憶合金絲或粘貼于基板表面的壓電陶瓷片為作動器,采用分散化控制、線性最優控制、非線性黃金分割自適應控制或混合靈敏度H∞魯棒控制等理論設計控制率,提高太陽帆板振動衰減速率。

另一方面,驅動擾動屬于強迫振動問題,發生在帆板旋轉運動方向,主要由驅動裝置驅動運動和太陽帆板扭轉模態耦合引起。抑制驅動擾動的目的是通過避免結構共振,降低擾動幅度。目前對太陽帆板驅動擾動的抑制主要集中在改進驅動裝置運動平穩性,減小輸出力矩的波動幅度。Chen等[13-14]討論了太陽帆板步進驅動擾動特征。Zhou等[15]引入輸入成形的驅動指令以改善驅動指向精度。Atlas等[16]采用一種修正的正弦驅動指令,減小系統誤差的影響,提高驅動穩定性。Zhu等[17-19]從力矩平衡方程出發,在參考電流中引入補償項,抵消諧波擾動成分,提高驅動運動穩定性。

綜上,殘余振動和驅動擾動屬于不同的振動類別,發生在不同的自由度,且涉及不同的振動模態,上述研究中將太陽帆板驅動系統(Solar Array Drive System, SADS)的減振問題單獨考慮殘余振動或驅動擾動問題。這無疑簡化了驅動系統減振問題的難度,但較少涉及驅動系統動力特性的系統性改進。為此,本文提出了在太陽帆板與驅動裝置之間設置一種剛度可調的多自由度減振裝置(Multi-degree-of-freedom Vibration Reduction Device, MVRD),以改善結構系統多種振動模式下的動力學特性,同時削弱上述兩種振動干擾(如圖1所示)。首先基于虛功原理和混合坐標描述,對含減振裝置的驅動系統進行動力學建模;然后系統分析減振裝置等效參數對扭轉和彎曲模態動力學特性的影響規律。在此基礎上,設計減振裝置原理樣機,通過驅動擾動特性試驗和模態參數測試試驗,驗證了減振裝置對兩類減振問題的有效性。

圖1 太陽帆板驅動系統組成示意圖

1 驅動系統動力學建模

1.1 基本假設

坐標體系定義如圖1所示?;炯僭O如下:①航天器本體:本體結構為剛性體,且約束平動和轉動位移;②驅動裝置:太陽帆板轉速很慢,短時間內可認為驅動裝置轉角為小量;③減振裝置:相對于展開狀態太陽帆板,減振裝置質量和結構尺寸為小量,僅體現其剛度和阻尼特性;④太陽帆板:小變形情況下,忽略太陽帆板材料參數和展開鎖定裝置連接剛度等非線性因素。

1.2 動力學描述

太陽帆板驅動系統虛功主要由慣性力、外力、內力、約束力以及阻尼力引起。慣性力虛功δWI由驅動裝置和太陽帆板兩部分引入,具體可以表示為

(1)

式中:Mm和Jm為驅動裝置轉動部件的質量和轉動慣量矩陣;r為驅動裝置對本體系原點的矢徑;θd和ωd為驅動裝置轉角和轉速;mk、ak和rk為帆板k節點的質量、加速度及矢徑;M、S、J、P和H分別為太陽帆板的質量、靜矩、慣性矩、模態動量和模態角動量矩陣;A和Ad為太陽帆板和驅動裝置到本體坐標系的轉換矩陣;Ads和Asd為減振裝置和太陽帆板隨體坐標系之間的轉換矩陣;η為太陽帆板模態坐標;Δd和φd為減振裝置平動和轉動方向變形向量。

驅動系統的彈性勢能包括太陽帆板模態變形勢能以及減振裝置的平動和轉動變形勢能。根據功能轉換原理,內力虛功δWN等于虛應變能的相反數

δWN=-δU=-(δΔdKrΔd+δφdKθφd)-δηTΛη

(2)

式中:Kr和Kθ為減振裝置的平動和轉動等效剛度矩陣。類似地,阻尼力虛功δWD具體可以表示為

(3)

式中:Dr和Dθ為減振裝置的平動和轉動等效阻尼矩陣。

在不考慮空間環境干擾和航天器本體激擾情況下,驅動力矩Tl是作用在太陽帆板上的唯一外力,對應的虛位移只發生在驅動方向。因此,外力虛功δWE表示為

δWE=T·δθ=δθTT

(4)

(5)

式中:θ|1和θ|3為θ的第一和三分量,ρo為帆板系原點到本體系原點距離。取虛位移為δq,則式(5)可整理為

(6)

式中:Cq是C(q,t)的雅可比矩陣,其中各參數矩陣為

(7)

引入拉格朗日乘子λ后,約束力虛功δWC表示為

(8)

1.3 模型組集

對太陽帆板驅動系統應用虛功原理

δWI+δWE+δWN+δWC+δWD=0

(9)

將式(1)~式(4)和式(8)代入式(9)并整理得到

由于短時間內的驅動裝置轉角θd和減振裝置扭轉變形φd均為小量,則坐標轉換矩陣取值接近單位陣。因此驅動系統動力學方程表示為如下形式

(10)

式中:左邊三項分別為廣義質量、阻尼和剛度矩陣。根據式(5)中的運動約束關系,式(10)可整理為

(11)

式中:θl和ωl為太陽帆板轉角和轉速;Jm2和J2為轉動慣量在驅動方向分量;|2h和|2l表示第二行或列分量。

2 減振裝置影響分析

本節在對所建立模型進行合理等效基礎上,重點分析減振裝置結構參數對驅動系統扭轉和彎曲模態的影響規律,為減振裝置結構設計提供設計依據。

2.1 動力學特性等效分析模型

根據文獻[20]可知,驅動力矩Tl可近似為

(12)

將式(12)代入式(11),并整理為

(13)

式中各項表達式為

式(13)中減振裝置對系統動力學特性影響體現為六自由度彈性阻尼單元,而驅動裝置對系統動力學特性的影響表現為扭轉彈性邊界條件,如圖2所示。

以某型衛星太陽帆板結構為例進行影響分析,其質量為21.6 kg,幾何尺寸和材料如圖2所示。固支邊界下前三階彎曲固有頻率分別為0.396 Hz,2.176 Hz和5.148 Hz;前三階扭轉模態固有頻率分別為1.492 Hz,5.237 Hz和10.004 Hz;振型見圖3,其中不同填充顏色代表振動方向相反。

圖2 驅動系統動力學特性等效模型

(a)一階彎曲

2.2 扭轉模態動力特性影響分析

將驅動系統扭轉模態動力學方程從式(13)中分離出來進行求解,然后分析減振裝置結構系數對驅動系統模態動力學特性的影響。

2.2.1 特征值問題求解

扭轉模態動力學方程可從式(13)中分離出來:

進行拉普拉斯變換后得到

(MTs2+DTs+KT)qT(s)=0

(14)

由此將扭轉模態特征值問題轉換為關于式(14)系數矩陣的代數方程求解問題:

det[MTs2+DTs+KT]=[(J2s2+Dds+Kd)×

(15)

(16)

2.2.2 扭轉模態參數影響分析

從扭轉模態動力學方程不難發現:減振裝置繞y軸剛度和阻尼系數對系統扭轉模態動力學特性具有重要影響。圖4為fT和ξT隨減振裝置繞y軸等效剛度的變化曲線,其中減振裝置等效阻尼系數取1.5 Nms/rad。fT僅在特定剛度范圍內隨減振裝置剛度系數升高而增大,分別由0.039 Hz、2.280 Hz和6.810 Hz提高到0.748 Hz、4.645 Hz和12.990 Hz,而在其余區間內基本保持不變。ξT1隨減振裝置剛度增大而逐漸減小,并最終穩定在0.004;ξT2和ξT3在等效剛度較小和較大的區間內均保持不變,僅在中間特定剛度范圍內隨剛度增大而逐漸減小。圖5為fT和ξT隨減振裝置繞y軸等效阻尼的變化曲線,其中減振裝置等效剛度為100 Nm/rad。fT對減振裝置阻尼參數并不敏感,基本保持不變。ξT隨減振裝置阻尼參數增大而逐漸增大,ξT2和ξT3的變化幅度明顯大于ξT1,分別可達到0.085和0.160。

圖4 減振裝置繞y軸剛度與扭轉模態參數關系曲線

圖5 減振裝置繞y軸阻尼與扭轉模態參數關系曲線

2.3 彎曲模態動力特性影響分析

2.3.1 特征值問題求解

將彎曲模態動力學方程從式(13)中分離出來

通過拉普拉斯變換可得頻域方程

(MBs2+DBs+KB)qB(s)=0

(17)

類似地,將彎曲模態特征值問題轉換為關于式(17)系數矩陣的代數方程問題:

det[MBs2+DBs+KB]=

(18)

2.3.2 彎曲模態參數影響分析

由動力學方程不難發現:減振裝置沿x軸和繞z軸方向上的剛度和阻尼分量會對系統彎曲模態產生重要影響。圖6給出繞z軸剛度分量對彎曲模態動力學參數的影響關系。在不同剛度范圍內,前三階固有頻率分別從0.005 Hz、0.194 Hz和0.836 Hz明顯增大到0.099 Hz、0.268 Hz和1.313 Hz,其余剛度范圍內均達到穩定。隨著減振裝置剛度參數的增大,ξB1逐漸減小。ξB1在剛度系數為0.1 Nm/rad時取得最大值為0.26,隨著減振裝置剛度參數的增大而逐漸減小,并最終穩定在0.039;剛度系數較小時,ξB2和ξB3平緩變化,并在特定剛度范圍內逐漸減小,在高剛度范圍內達到穩定。圖7給出繞z軸阻尼分量對彎曲模態動力學參數的影響關系。前三階固有頻率在所研究阻尼范圍內基本保持不變,說明其對減振裝置的阻尼系數并不敏感。但前三階模態阻尼比隨減振裝置阻尼系數的增大而增大,且增大速率也在逐漸增大。ξB3增長最明顯,而ξB2增長最慢。

圖6 減振裝置繞z軸剛度與彎曲模態參數關系曲線

圖7 減振裝置繞z軸阻尼與彎曲模態參數關系曲線

圖8和圖9分別給出了沿x軸剛度和阻尼分量對彎曲模態動力學參數的影響關系。前三階固有頻率隨x軸剛度的變化趨勢與繞z軸剛度的變化趨勢基本相同,僅變化幅度更為明顯。在減振裝置剛度為0.1 Nm/rad時,ξB1取得最大值0.49,然后隨剛度增加而逐漸減小,直到40 Nm/rad剛度系數趨于穩定。剛度系數相對較小時,ξB2和ξB3平緩變化,而后隨著剛度增大而明顯減小,并最終趨于穩定。固有頻率對減振裝置x軸阻尼分量不敏感,在研究范圍內的相對變化小于2%。前三階模態阻尼比隨減振裝置沿x軸阻尼分量的增大而增大,且增大率也具有相同的變化趨勢;其中ξB2最大(可達0.865),ξB1最小(僅為0.052)。

圖8 減振裝置沿x軸剛度與彎曲模態參數關系曲線

圖9 減振裝置沿x軸阻尼與彎曲模態參數關系曲線

3 減振裝置結構設計

通過合理的減振裝置結構設計,可以全面抑制驅動擾動和剩余振動。為便于驗證減振裝置對太陽帆板一階彎曲和一階扭轉振動的抑制效果,且考慮到試驗驗證的可行性,選用包含1個支架和1塊基板構成的太陽帆板結構為減振對象,其中部件尺寸和材料如圖2所示。結構總質量為14 kg,扭轉轉動慣量11.74 kg·m2,固支邊界一階彎曲和扭轉頻率分別為0.787 Hz和0.635 Hz。一方面,根據表2所列驅動裝置參數,諧波激勵頻率為0.21 Hz的整數倍[20];而實際驅動系統的第一扭轉固有頻率為0.635 Hz,非常接近三次諧波激勵頻率0.63 Hz。因此,調整第一扭轉固有頻率有利于降低共振響應。另一方面,彎曲模態的衰減時間與模態阻尼比固有頻率均成負相關關系,因此應增大模態阻尼,并盡量避免減小相關固有頻率,可以加速殘余振動的衰減。

為便于調整等效結構參數,減振裝置結構形式相對簡單。如圖10(a)所示,減振裝置主要由高剛度中軸和高損耗因子阻尼環同軸安裝而成。圖10(b)為采用鋁合金材料的中軸結構,其通過上下連接界面分別與頂蓋和底蓋螺栓連接。中軸為可以更換的部件,通過選用不同中軸截面形狀改變減振裝置整體剛度特性,例如十字形截面,因其設計參數更多,便于靈活調整各方向剛度分量。圖10(c)為采用硅橡膠高分子阻尼材料的阻尼環結構,通過環向均勻布局,盡量增大其分布面積,提高減振裝置的阻尼系數。該構型可以保證太陽帆板在發生軸向扭轉振動時,阻尼環各部分均發生明顯的結構變形,有利于提高阻尼器的等效阻尼參數。頂蓋和底蓋結構與中軸和阻尼環連接為整體,并在中心開設通孔,便于中軸替換;內表面平整,便于與阻尼環粘接或共固化連接;外表面按安裝界面設螺紋孔,用于與驅動裝置和太陽帆板連接。

(a)實物圖

經反復調整后,確定中軸長度為30 mm;十字界面各邊形狀相同,單邊長為8 mm,單邊厚為2 mm;阻尼環內外直徑分別為150 mm和110 mm,整體質量為1.169 kg,外包絡為Φ156 mm×56 mm。分別在安裝減振裝置前后進行復模態分析和擾動特性分析,結果見表1和圖11。安裝減振裝置后,一階扭轉固有頻率為0.555 Hz,與擾動頻率明顯錯開,對應的模態阻尼比提高491.04%;擾動力矩波動幅度由0.456 Nm降至0.293 Nm,減少35.74%;一階彎曲模態阻尼比可提高十倍以上,但固有頻率僅能降低19.70%,可以計算一階彎曲模態99%振幅衰減時間從225.9 s大幅降低至24.2 s。

表1 驅動系統復模態分析結果

圖11 驅動系統擾動力矩仿真分析結果

4 減振效果驗證

在試驗室環境中進行了擾動特性和模態特性測試試驗,以評價減振裝置對太陽帆板軸向驅動擾動和彎曲殘余振動的抑制效果。

4.1 擾動特性試驗

4.1.1 試驗系統

所搭建的模擬太陽帆板驅動系統包括:驅動裝置、柔性附件以及減振裝置(如圖12所示)。驅動裝置為低軌直接驅動型SADA(參數見表2);柔性附件由兩塊構型相同(尺寸為1 630 mm×400 mm×5 mm)的鋁合金薄板組成,其與太陽帆板滿足轉動慣量等效(11.8 kg·m2)和扭轉振動基頻等效(0.8 Hz)。為抵消柔性附件對驅動裝置軸承的壓力,在轉動軸線方向配備龍門吊和吊裝組件以卸載重力。

圖12 驅動系統擾動特性測試系統

表2 太陽帆板驅動裝置參數

測試系統由Kistler六自由度測力平臺、Polytec激光測振儀(包括前置端和后置端)、信號放大器、數據采集器、SADA外置驅動器、計算機及分析軟件等組成,如圖12所示。測力平臺與SADA通過轉接工裝固定連接,以測量驅動擾動;激光測振儀直接測量激光柔性附件運動速度變化。測量結果分別經過信號放大器和測振儀后置端進行處理,并通過數據采集器傳遞給計算機。SADA外置驅動器用于產生指定的脈沖信號,控制驅動裝置以0.063 5°/s勻速轉動。

4.1.2 結果分析

圖13給出了勻速驅動指令下達后,驅動系統穩定驅動段擾動力矩時間歷程和頻譜分布曲線。擾動力矩具有周期性波動特點,具體表現為主要低頻波動頻率和若干小幅高頻波動頻率的疊加;頻譜分布曲線中均包含一個明顯的低頻峰值,安裝前后頻譜峰值對應頻率分別為0.625 Hz和0.563 Hz;安裝前的波動幅度和頻譜峰值均明顯高于減振裝置安裝后。

(a)時間歷程

進一步通過圖11和圖13(a)對比不難發現:分析結果與試驗數據一致性較好,擾動力矩主要波動頻率基本一致,擾動幅度偏差在17%以內(仿真0.456 Nm,試驗0.552 Nm),其主要由試驗測得的高頻擾動成分引入。表3給出了穩態擾動力矩波動幅度和均方差、平均轉速和系統基頻數據。對比相同測試狀態下的①~③數據可發現:試驗結果的重復性良好。安裝減振裝置前后的擾動力矩波動幅度平均值相對減小41.61%,力矩均方差降低53.48%。這是由于減振方案的實施,一方面調整了系統一階固有頻率,使其降低至0.563 Hz,偏離了激擾頻率,減小了系統共振響應;另一方面系統模態阻尼比增加后,在一定程度上減小了柔性附件振幅,從而實現減小系統擾動的目的。

表3 驅動系統擾動特性試驗結果

4.2 彎曲模態試驗

4.2.1 試驗系統

由于驅動裝置在非驅動方向可以近似為剛體,因此試驗對象僅由減振裝置和太陽帆板組成。太陽帆板結構和材料特性如第3章所述,其由重力卸荷裝置懸吊,消除重力的影響。減振裝置安裝在支架和墻面之間。由于本試驗是在大氣環境中進行,因此大氣阻尼效應影響不可避免。測試系統由電磁激勵器、力傳感器、加速度傳感器、電荷放大器、數據采集系統、計算機以及模態測試分析軟件等組成,如圖14所示。

圖14 太陽帆板模態特性測試系統

模態試驗采用單點力激勵和多點加速度測量方式,其中激勵點位于支架橫梁上,測點分布在支架和基板表面。模態分析軟件通過電荷放大器將0~12.5 Hz頻段隨機激勵信號輸出至電磁激勵器,數據采集器以400 Hz采樣頻率采集加速度測點響應,傳輸給計算機后,分析軟件后采用頻域辨識方法進行模態辨識,頻率辨識分辨率為62.5 mHz。

4.2.2 結果分析

獲取的一階彎曲模態的固有頻率和模態阻尼比如表4所示。與表1給出的分析結果對比不難發現:分析與試驗測得的固有頻率非常接近(分析0.787 Hz,試驗0.765 Hz),相對偏差小于3%;由于測試過程中存在空氣阻尼,因此測試模態阻尼數據與分析結果的取值存在偏差,但是模態阻尼的變化趨勢相同。安裝減振裝置后,模態阻尼比由3.202%(含空氣阻尼)顯著增加到8.878%(含空氣阻尼),其中減振裝置的作用使模態阻尼增加了177%,說明減振裝置對太陽帆板的阻尼貢獻很大;固有頻率從0.765 Hz降至0.63 Hz,僅降低約18%。對比同種測試狀態下的①~③數據可以發現:模態試驗結果具有良好的重復性。安裝減振裝置前后的太陽帆板99%振幅衰減時間分別為29.893 s和13.103 s,相對減少56.17%。測試所得模態形狀如圖15所示,其中模態位移均沿太陽帆板的展開方向逐漸增大,最大位移出現在帆板外緣。

表4 驅動系統模態試驗結果

(a)無減振裝置

5 結 論

太陽帆板驅動系統的驅動擾動和殘余振動屬于不同的振動類別,發生在不同的振動方向,涉及不同的振動模態。為了系統抑制這兩種振動,本文合理設計了一種安裝于太陽帆板和驅動裝置之間的多自由度減振裝置,通過理論建模、影響分析、結構設計和性能驗證,得到主要結論如下:

(1)該減振裝置由可更換的中心軸和固定的阻尼環組成,可以在調節固有頻率的同時增加系統阻尼。

(2)減振裝置的等效剛度僅在特性剛度區間內對驅動系統固有頻率影響顯著,因此應慎重選擇,以使扭轉模態固有頻率錯開驅動擾動激擾頻率,避免共振。

(3)減振裝置的等效阻尼與系統模態阻尼比呈正相關關系,因此在減振裝置的結構設計中建議采用損耗系數高的阻尼材料以加速帆板殘余振動衰減。

(4)通過理論和試驗方式驗證,減振裝置對兩類振動問題均取得了較好的抑制效果,其中驅動擾動幅值衰減40%以上,99%殘余振動衰減時間減少了56%以上。

猜你喜歡
模態振動系統
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
中立型Emden-Fowler微分方程的振動性
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 亚洲视频色图| 精品日韩亚洲欧美高清a| 国产女人水多毛片18| 国产欧美日韩综合一区在线播放| julia中文字幕久久亚洲| av在线无码浏览| 国产区成人精品视频| 国产美女精品人人做人人爽| 国产精品成人一区二区| 中文字幕在线免费看| 欧美色图第一页| 免费aa毛片| 欧美日韩第三页| 免费播放毛片| 欧美在线网| 亚洲免费黄色网| 欧美性天天| 2022国产无码在线| 免费在线国产一区二区三区精品| 免费A∨中文乱码专区| 欧美自慰一级看片免费| 日韩 欧美 国产 精品 综合| 亚洲精品视频免费| 久久精品只有这里有| 国产一区二区三区精品久久呦| 欧美日本视频在线观看| 午夜福利视频一区| 亚洲天堂视频在线播放| 久久久久人妻精品一区三寸蜜桃| 亚洲国产欧美目韩成人综合| 91娇喘视频| 成人免费午间影院在线观看| 久久综合亚洲鲁鲁九月天| 香蕉99国内自产自拍视频| 国产欧美精品午夜在线播放| 欧美日韩动态图| 国产超薄肉色丝袜网站| 中文字幕在线看| 日韩不卡免费视频| 亚洲天堂免费观看| 国产成人免费手机在线观看视频 | 国产亚洲现在一区二区中文| 白丝美女办公室高潮喷水视频| 粉嫩国产白浆在线观看| 国产精品无码一区二区桃花视频| 白浆免费视频国产精品视频| 99久久精品免费看国产免费软件 | 国产美女一级毛片| 国产精品色婷婷在线观看| 71pao成人国产永久免费视频 | 久久婷婷六月| 毛片视频网| 91精品情国产情侣高潮对白蜜| 激情无码字幕综合| 无码丝袜人妻| 白浆视频在线观看| 91精品国产自产在线老师啪l| 国产精品lululu在线观看| 国产AV无码专区亚洲精品网站| 五月综合色婷婷| a欧美在线| 四虎免费视频网站| 欧美一级色视频| 国产精品成人不卡在线观看| 日本精品影院| 国产办公室秘书无码精品| 亚洲国产成人精品无码区性色| 亚洲国产第一区二区香蕉| 在线观看亚洲人成网站| 天天综合色天天综合网| 三上悠亚在线精品二区| 日本中文字幕久久网站| 91精品啪在线观看国产91九色| 久久99国产综合精品1| 亚洲无线国产观看| 91麻豆国产在线| 久久夜色精品| 精品小视频在线观看| 久久无码av一区二区三区| 亚洲三级视频在线观看| 欧美午夜小视频| 日韩毛片视频|