亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        大型空間電站在軌展開與組裝動力學與控制

        2021-03-31 02:51:42榮吉利石文靜李偉杰辛鵬飛
        宇航學報 2021年3期
        關鍵詞:變形結構

        榮吉利,崔 碩,石文靜,李偉杰,李 瀟,辛鵬飛,3

        (1. 北京理工大學宇航學院,北京100081;2. 北京空間飛行器總體設計部,北京100094;3. 空間智能機器人系統(tǒng)技術與應用北京市重點實驗室,北京100094)

        0 引 言

        近年來中國航天發(fā)展迅速,對大型空間結構的研究越來越迫切,例如深空探測航天器大型太陽帆、空間站的太陽能陣列、大型在軌組裝電站、大型空間模塊化天線、大型載人空間站等。一系列大型空間站的研發(fā)逐漸提上日程。大型在軌組裝電站由于結構巨大、尺度巨大、構建規(guī)模龐大,需要結構單元在折疊狀態(tài)發(fā)射,后在空間軌道進行結構展開與組裝,工況復雜,地面試驗很難滿足需求,因此采用準確有效的數(shù)值模擬方法預測空間結構動力學行為是設計的重要一環(huán)。文獻[1]提出了一種熱-力學耦合分析有限元方法,以解決大型空間結構的熱變形問題。文獻[2-3]對空間結構及在軌操作領域進行了多層次的深入研究。文獻[4]對網(wǎng)架式星載天線非線性動力學與控制進行研究。文獻[5]以太陽帆塔等細長結構中空間柔性梁為研究對象進行了穩(wěn)定性分析。文獻[6]提出了一種空間太陽能電站追蹤太陽運動的準對日定向姿態(tài)方案。文獻[7]提出了基于混合坐標法的太陽電池翼在軌載荷分析方法。美日等國早在多年前便開始研究空間太陽能發(fā)電技術,國內(nèi)針對大型空間電站的研究依然存在大量空白。

        受火箭運輸能力制約,大型空間電站以可展開、可組裝結構為結構單元,通常由大量梁、桿、薄板等柔性構件組成,且由于其尺度巨大,表現(xiàn)出大柔性、大變形的特點,是典型的剛柔多體系統(tǒng)。大多數(shù)商業(yè)軟件無法有效處理大變形問題,因此自編軟件成為解決這類問題的重要思路。采用精確的能解決大變形問題的剛柔耦合模型對研究空間電站展開與組裝過程中動力學特性具有十分重要的意義。

        此外,由于大型空間電站規(guī)模龐大,仿真計算量巨大,商業(yè)軟件很難滿足計算效率上的需求,在自編軟件中嵌入高效的并行算法能顯著提高計算效率,但本文暫不涉及并行算法相關內(nèi)容。

        多體系統(tǒng)動力學包含多剛體系統(tǒng)動力學和柔性多體系統(tǒng)動力學兩類。文獻[8-9]提出多剛體系統(tǒng)建模的自然坐標方法(Natural coordinate formulation, NCF),又稱為完全笛卡爾坐標法,該方法具有質(zhì)量矩陣為常數(shù)矩陣,并且約束表示簡單清晰等優(yōu)勢。

        針對柔性多體系統(tǒng)動力學,文獻[10]提出了浮動坐標法,又稱為混合坐標法,是目前應用最廣泛的柔性多體系統(tǒng)動力學分析方法。但是該方法僅僅適用于小變形問題,不能精確地描述柔性構件大變形。文獻[11-12]提出了能夠精確描述柔性構件大變形的絕對節(jié)點坐標方法(Absolute nodal coordinate formulation, ANCF),采用斜率矢量代替了傳統(tǒng)有限元中的轉角坐標,將單元的節(jié)點坐標全部定義在全局坐標系下,不存在小變形、小轉動假設,所推導的多體系統(tǒng)微分代數(shù)方程具有常數(shù)質(zhì)量矩陣,不存在科氏力和離心力項,能夠精確描述大變形問題。由NCF和ANCF共同搭建的剛柔耦合多體系統(tǒng)動力學方法稱為絕對坐標[13-15](Absolute coordinate-based method, ACB)方法。

        本文針對某典型百米尺度高剛度桁架類大型空間電站在軌展開與組裝過程中的動力學特性進行研究,采用NCF與ANCF建立可適用于模型大變形的剛柔耦合模型,開發(fā)了一套動力學仿真軟件,并采用商業(yè)軟件MSC.ADAMS與ANSYS聯(lián)合仿真[16-17]驗證其在小變形下仿真模型的準確性,且由于所采用建模方法中不存在小變形、小轉動假設,優(yōu)于大多數(shù)商業(yè)軟件,可適用于結構大變形、大轉動動力學預測;針對線形桁架多級展開過程的結構振動問題選定了運動規(guī)劃方案,能有效抑制結構振動,增強展開過程穩(wěn)定性;探究了組裝過程中組裝速度、組裝機構勁度系數(shù)與阻尼系數(shù)對組裝過程穩(wěn)定性的影響,并提出了對參數(shù)設計的優(yōu)化思路;集成了一套針對大型空間電站的動力學仿真平臺搭建方法,可為進一步開展此類大型空間結構的設計與研發(fā)工作提供有力的技術支持與方法指導。

        1 大型空間電站及其建模

        1.1 大型空間電站構建方案

        針對大型空間結構的構建方法[18]大致分為三類:可展開結構構建、太空成形結構構建以及可直立結構構建。三類方法各有其優(yōu)缺點。

        可展開結構構建將折疊結構在地面組裝,后經(jīng)火箭等發(fā)射進入軌道,在軌道中逐級展開,這種方案減少了在軌道中操作的工作量,但實現(xiàn)風險較大,容易出現(xiàn)結構無法完全展開的問題,可靠性較低。太空成形結構構建方法,主要是將未加工的材料,在軌道成形裝配??芍绷⒔Y構構建將大型空間結構的部件在地面搭建,逐個發(fā)射后在空間中進行組裝。

        本文作為研究對象的典型大型空間電站,其構建方案結合了可展開結構構建與可直立結構構建,將部件在地面折疊后集中發(fā)射,進入軌道后進行部件展開與組裝,能有效提高構建效率。但是大型空間電站在軌展開與組裝過程過于復雜,將大大增加構建難度。因此這種方案更適用于由同類或少數(shù)幾類部件構建的大型空間結構,有利于減小在軌操作的風險和難度。

        電站結構的部件構型參考專利:“一種可重復展收的桁架結構及其胞元”[19]。根據(jù)專利設計,整個大型空間電站可以通過兩類基礎桁架胞元構建而成。構建過程大致可分為:1)單個胞元結構折疊并集中發(fā)射;2)在軌道中單個胞元的展開;3)胞元組裝成為空間電站中心體結構;4)胞元多級展開;5)帶太陽翼的線形桁架向中心體的組裝;6)輻射桁架的展開與組裝;7)環(huán)形桁架的展開組裝;8)多級帶太陽翼桁架向空間電站的組裝等。

        1.2 建模方法

        本文所研究的大型空間電站包含大量柔性與剛性構件,是典型的剛柔多體系統(tǒng),選擇采用自然坐標方法與絕對節(jié)點坐標方法對其建模,可以精確描述大型空間結構在構建過程中的大變形、大轉動。

        1.2.1自然坐標法

        自然坐標法中選取剛體上2個固定點的位置矢量與2個單位矢量作為廣義坐標,見式(1),包含12個自由度:

        (1)

        式中:ri,rj為剛體上固定點i,j的位置矢量,u和v為剛體上固定的單位矢量。以點i為局部坐標系原點,以(rj-ri),u,v為局部坐標系的基矢量構建剛體局部坐標系ξ-η-ζ如圖1所示:

        圖1 NCF描述的空間剛性體Fig.1 Spatial rigid body described by NCF

        剛體上任意一點N在局部坐標系ξ-η-ζ下的位置矢量為:

        (2)

        (3)

        剛體上任意一點N在全局坐標下的位置矢量可表示為:

        r=(1-c1)ri+c1rj+c2u+c3v

        (4)

        將式(4)改寫成如下矩陣形式:

        (5)

        式中:C=[(1-c1)I3,c1I3,c2I3,c3I3],稱為形函數(shù),I3是3階單位矩陣。

        1.2.2絕對節(jié)點坐標方法

        采用ANCF全參數(shù)梁單元對結構中的柔性桿建模,如圖2所示。

        圖2 ANCF描述的全參數(shù)梁單元Fig.2 Full-parameter beam element described by ANCF

        該單元中任意一點N的全局位置矢量為:

        r=[r1,r2,r3]T

        (6)

        e為單元的廣義坐標,定義如下:

        (7)

        通過有限元理論形函數(shù)的一般推導方法,可得到該單元的形函數(shù)S為:

        S=

        [s1I3,ls2I3,ls3I3,ls4I3,s5I3,ls6I3,ls7I3,ls8I3]

        (8)

        式中:I3為3×3單位矩陣;s1=1-3ξ2+2ξ3,s2=ξ-2ξ2+ξ3,s3=η-ξη,s4=ζ-ξζ,s5=3ξ2-2ξ3,s6=-ξ2+ξ3,s7=ξη,s8=ξζ,其中,ξ=x/l,η=y/l,ζ=z/l,l為單元長度。

        ANCF全參數(shù)梁單元上任意一點N的全局坐標均可表示為:

        r=Se

        (9)

        將式(9)對時間求導,可得到:

        (10)

        進而可得到單元的動能為:

        (11)

        單元的質(zhì)量矩陣為:

        (12)

        式中:ρ表示密度,V為體積積分變量,Ve為體積積分域,M為常數(shù)質(zhì)量矩陣。

        1.3 滑移鉸建模與仿真校驗

        桁架類空間結構展開過程中伴隨滑移鉸的滑移,在動力學仿真中涉及節(jié)點沿單元的移動,會導致約束方程的改變,因此對滑移鉸的精確建模十分關鍵。

        1.3.1滑移鉸建模方法

        文獻[20]提出了一種滑動關節(jié)建模方法,采用4個約束方程公式化滑移鉸。

        圖3 滑移鉸Fig.3 Sliding joint

        如圖1所示,滑動關節(jié)在主體上的P,以及在滑動桿上的未變形的局部相對位置Q,ξt為t時刻點P相對于MN的物質(zhì)坐標,滿足ξt=lMP/lMN。由式(13):

        rOQ=S(ξ,η,ζ)eMN

        (13)

        滑移鉸所需要的滿足的約束方程如下:

        (14)

        在式(14)中,第一行的等式包含3個約束,第二行包含1個約束。前3個約束用于將點P與Q的運動限制在相同的全局位置,這導致對應生成三個拉格朗日乘子λ=[λ1,λ2,λ3]T。滑移鉸在沒有摩擦力的情況下,約束力僅存在于垂直滑動方向上,在滑動方向上無約束力,即為第4個約束的含義。

        在ANCF全參數(shù)梁單元中,有:

        (15)

        式中:sx1=-6ξ+6ξ2,sx2=l(1-4ξ+3ξ2),sx3=-lη,sx4=-lζ,sx5=6ξ-6ξ2,sx6=l(-2ξ+3ξ2),sx7=lξ,sx8=lζ,其中,ξ=ξt,η=y/l,ζ=z/l,l為單元長度。

        1.3.2帶滑移鉸的展開動力學仿真校驗

        建立帶滑移鉸的單個胞元展開動力學模型,結構構型參考專利:“一種可重復展收的桁架結構及其胞元”[19]。

        根據(jù)專利內(nèi)容,單胞元結構為單自由度結構。采用ANCF全參數(shù)梁單元對柔性桿進行離散;鎖定機構與其位于的剛性框之間建立滑移鉸約束;各桿之間建立旋轉鉸約束。

        圖4 單個胞元結構示意圖Fig.4 Single cell structure diagram

        各項參數(shù)見表1:

        表1 結構參數(shù)Table 1 Structural parameters

        α0=5°為折疊狀態(tài)初始角度,αf=90°為完全展開時目標角度。

        固定胞元一端的剛性框,連桿間交接方式為旋轉鉸接;在胞元上施加沿Y正方向大小為0.005g的重力加速度。

        由于模型中單個胞元展開過程所施加載荷較小,因此柔性與變形較小,商業(yè)軟件能夠有效處理此類問題,選擇采用MSC.ADAMS與ANSYS聯(lián)合仿真方法[16-17]。首先在MSC.ADAMS中建立單個胞元模型,在ANSYS中生成模態(tài)中性文件(.mnf文件)后導入ADAMS/FLEX模塊,可實現(xiàn)基于模態(tài)頻率對柔性件的建模與仿真。

        繪制A點Y坐標變化曲線圖如圖5所示:

        圖5 A點Y坐標變化曲線圖Fig.5 The Y-position of point A

        本文建模方法的結果與商業(yè)軟件MSC.ADAMS與ANSYS聯(lián)合仿真結果基本一致,驗證了所建立約束尤其是滑移鉸的準確性,同時驗證了基于NCF與ANCF建立的動力學模型的準確性。

        相較于商業(yè)軟件MSC.ADAMS與ANSYS聯(lián)合仿真方法,本文建模方法在處理柔性體尤其是大變形問題時更加精確。對比重力加速度為0.05g、0.5g、1g、10g的情況下兩種仿真方法的結果如圖6所示,分別以“G0.05”、“G0.5”、“G1”、“G10”代指對應載荷工況。

        圖6 A點Y坐標變化曲線圖Fig.6 The Y-position of point A

        表2選取0.05 s、0.10 s、0.15 s、0.20 s對比不同載荷情況下商業(yè)軟件相較于本文方法結果之間的誤差。隨載荷增大,兩種方法之間誤差逐漸增大。當載荷較大時,結構出現(xiàn)較大變形,商業(yè)軟件MSC.ADAMS與ANSYS聯(lián)合仿真結果出現(xiàn)較大誤差。本文建模方法不存在小變形、小轉動假設,在求解柔性體變形問題中更加精確,這一大優(yōu)勢使之更適用于有大柔性、大變形特點的大型空間電站的動力學行為論證。

        表2 誤差對比Table 2 Error comparison

        2 大型空間電站構建過程中典型工況

        大型空間電站構建過程復雜,構建方案與構建方式多樣,且涉及不同的環(huán)境擾動等影響因素,自編軟件可以適用于整個大型空間電站的構建過程,本文以某典型桁架類大型空間電站系統(tǒng)為研究對象,選取結構展開與組裝過程中的兩類典型工況進行分析。

        2.1 線形桁架展開規(guī)劃

        由于火箭單次運輸空間有限,空間電站的結構胞元需要在折疊狀態(tài)下運輸至太空軌道,后在軌道中進行展開與組裝操作,因此針對結構在空間中展開過程進行的動力學仿真尤為關鍵。

        胞元展開控制方式對結構穩(wěn)定性影響很大,良好的展開控制方案能有效抑制結構振動。

        針對圖4中單個胞元的展開過程,選用5種角度控制函數(shù)如下,分別為線性函數(shù)、三次多項式插值函數(shù)、五次多項式插值函數(shù)、擺線運動插值函數(shù)以及余弦插值函數(shù):

        (16)

        (17)

        (18)

        (19)

        (20)

        式中:α如圖4所示。α0=5°為折疊狀態(tài)初始角度,αf=90°為完全展開角度,tf=10 s為單個胞元展開時間。

        由于不同的工況中所選取的樣點初始空間坐標會有區(qū)別,以下統(tǒng)一選擇位移曲線而非坐標曲線進行對比分析。

        圖7~9中,SX,SY,SZ分別表示X,Y,Z方向的位移。

        由圖7可知,線性函數(shù)的結果有明顯的振動,這是由于線性函數(shù)的初始增量太大,結構產(chǎn)生較大的彈性變形,導致展開過程不平滑。其它4種控制函數(shù)在Y方向上位移曲線光滑。由于初始狀態(tài)單個胞元存在α0=5°,A點Y坐標初始值為0.6972 m,故圖7中Y方向位移平衡值為7.3028 m。由圖8和圖9可知,余弦插值函數(shù)有最好的展開穩(wěn)定性與最小的振動殘留,但是位移峰值大于擺線運動插值函數(shù)。擺線運動插值函數(shù)穩(wěn)定性與振動殘留僅次于余弦插值函數(shù)。

        圖7 不同的角度控制函數(shù)下A點Y方向位移曲線圖Fig.7 The Y-displacement of point A from different angle control functions

        圖8 不同的角度控制函數(shù)下A點X方向位移曲線圖Fig.8 The X-displacement of point A from different angle control functions

        圖9 不同的角度控制函數(shù)下A點Z方向位移曲線圖Fig.9 The Z-displacement of point A from different angle control functions

        針對包含8個胞元的折疊結構多級展開過程,分別選用擺線運動插值函數(shù)與余弦插值函數(shù)作為角度控制函數(shù)。

        圖10 8個胞元折疊狀態(tài)圖Fig.10 8 cells folded state diagram

        從遠離固定端P的胞元C1逐級沿Y方向展開,參數(shù)與單胞元一致,相鄰胞元的展開開始時間間隔12 s。在胞元C1上對應圖4中A點。

        圖11 不同的角度控制函數(shù)下A點Y方向位移曲線圖Fig.11 The Y-displacement of point A from different angle control functions

        圖12 不同的角度控制函數(shù)下A點X方向位移曲線圖Fig.12 The X-displacement of point A from different angle control functions

        圖13 不同的角度控制函數(shù)下A點Z方向位移曲線圖Fig.13 The Z-displacement of point A from different angle control functions

        結構展開沿Y方向。由圖11可知Y坐標變化有明顯的階段性,擺線運動插值函數(shù)與余弦插值函數(shù)都具有良好的穩(wěn)定性。

        由圖12和圖13可知,隨著結構多級展開,采用余弦插值函數(shù)作為角度控制函數(shù)時結構振動更明顯,在第7個胞元展開后X方向最大位移為4.402 m,Z方向最大位移達到8.253 m,而采用擺線運動插值函數(shù)在X方向最大位移為0.363 m,Z方向最大位移為4.577 m。由此得出結論,在桁架結構多級展開中,擺線運動插值函數(shù)明顯優(yōu)于余弦插值函數(shù)。

        2.2 線形桁架向中心體的組裝

        結構在軌展開后,需通過組裝的方式逐步拼接構建成完整結構。

        大型空間電站采用先組裝中心體結構,后組裝外圍結構的方式逐步搭建。對帶太陽翼的線形桁架向中心體結構的組裝過程進行動力學仿真。

        圖14 線形桁架向中心體組裝Fig.14 Linear truss assembly to center body

        模型如圖14所示。X正軸與負軸方向上兩個被組裝的線形桁架長度為96 m,各包含12個線形胞元。組裝過程中太陽翼處于折疊收攏狀態(tài),折疊狀態(tài)剛度較大,在空間電站組裝過程中無需考慮其柔性,故采用NCF方法建模。單個太陽翼質(zhì)量為2 t。

        圖15 組裝局部圖Fig.15 Partial assembly diagram

        選取一組組裝參數(shù)見表3,作為對照參數(shù),建立3個對比算例,分別改變組裝速度v、組裝機構勁度系數(shù)K以及組裝機構阻尼C,對比不同組裝參數(shù)對組裝過程穩(wěn)定性的影響。

        表3 參數(shù)Table 3 Parameters

        由于在Y,Z方向上位移變化小于10-5m,組裝機構內(nèi)力小于1 N,相較于X方向可以忽略,因此僅對比X方向上位移與內(nèi)力變化。繪制組裝過程中A點在X方向(為組裝方向)上的位移曲線與組裝機構內(nèi)力變化曲線圖如圖16~17所示。

        圖16 組裝過程中A點X方向位移曲線圖Fig.16 The X-displacement of point A during structure assembly

        圖17 組裝過程中組裝機構內(nèi)力曲線圖Fig.17 Internal force curve of assembly mechanism during structure assembly

        由圖16~17,采用表3參數(shù)時,在t=6.1705 s時點A距離平衡位置最遠,為0.10303 m,t=5.6305 s時組裝機構內(nèi)力出現(xiàn)最大值,為246 N;當v=0.5 m/s時,在t=3.161 s點A距離平衡位置最遠,為0.27628 m,在2.611 s時內(nèi)力峰值為660.5 N;當組裝機構阻尼系數(shù)C=500 N·m-1·s-1,在t=6.2505 s時,點A位移最大值0.13453 m,在t=5.9905 s時,內(nèi)力峰值280 N,達到平衡狀態(tài)所需要的時間大大增長;組裝機構勁度系數(shù)K=1500 N/m時,在t=6.3205 s,點A距離平衡位置位移最大值0.11178 m,當t=5.5905 s時內(nèi)力最大值213.6 N。當組裝速度增大時,組裝過程發(fā)生得更早,但組裝過程位移更大,內(nèi)力也明顯增大;當減小C時,位移峰值與內(nèi)力峰值會增大,但峰值發(fā)生時間會延后;當減小K時,位移峰值增大,內(nèi)力峰值略微減小。

        通過以上結果給出優(yōu)化方案,為使組裝過程更加穩(wěn)定,應該減小組裝速度v,增大阻尼系數(shù)C,適當增大機構勁度系數(shù)K。

        3 結 論

        本文以某典型大型空間電站為研究對象,研究空間電站的在軌展開與組裝過程,集成了一套綜合的針對大型空間電站的動力學仿真平臺搭建方法,采用自編軟件的方式建立了一套動力學仿真平臺,得到如下結論:

        1)經(jīng)過與商業(yè)軟件MSC.ADAMS與ANSYS聯(lián)合仿真結果的對比,本文選用的建模方法以及所建立的空間桁架結構模型的動力學仿真平臺十分精確;與大多數(shù)商業(yè)軟件不同,本文所采用的建模方法可以精確描述結構的大變形、大轉動,更適用于大型空間電站的動力學行為評估與論證。

        2)單個桁架胞元展開過程中,選擇余弦插值函數(shù)作為角度控制函數(shù),結構擁有最好的穩(wěn)定性與最小的振動殘差。桁架結構多級展開過程中,擺線運動插值函數(shù)作為角度控制更有利于提高結構的穩(wěn)定性。

        3)線形桁架向中心體組裝過程中,組裝速度對結構穩(wěn)定性的影響很大;適當增大組裝機構阻尼系數(shù),可以有效減少結構達到穩(wěn)定所需要的時間;增大組裝機構勁度系數(shù)可以減小組裝過程的最大位移。為增強結構組裝過程的穩(wěn)定性,可在工程允許的范圍內(nèi)減小組裝速度,增大空間結構系統(tǒng)阻尼與勁度系數(shù)。

        4)本文涉及的針對大型空間電站的仿真平臺搭建方法與所編寫的軟件平臺適用于同類大型空間平臺搭建過程的動力學預測與控制方案評估,軟件平臺具有良好的擴展性,能夠為進一步開展此類空間結構的設計與研發(fā)提供有效的技術支持。

        猜你喜歡
        變形結構
        《形而上學》△卷的結構和位置
        哲學評論(2021年2期)2021-08-22 01:53:34
        談詩的變形
        中華詩詞(2020年1期)2020-09-21 09:24:52
        論結構
        中華詩詞(2019年7期)2019-11-25 01:43:04
        新型平衡塊結構的應用
        模具制造(2019年3期)2019-06-06 02:10:54
        “我”的變形計
        變形巧算
        例談拼圖與整式變形
        會變形的餅
        論《日出》的結構
        創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
        久久亚洲av成人无码软件 | 无码欧美毛片一区二区三| 国产人妻熟女呻吟在线观看| 欧美最猛性xxxx| 亚洲 欧美 国产 制服 动漫| 国产真人无码作爱视频免费| 任你躁欧美一级在线精品免费| 久久最黄性生活又爽又黄特级片| 熟女人妻中文字幕av| 中文字幕在线观看亚洲日韩 | 国产精品视频白浆免费视频| 国产在线一区二区三区四区不卡| 国产亚洲2021成人乱码| 久久精品国波多野结衣| 日本最新在线一区二区| 久久熟女少妇一区二区三区| 国产黑丝美女办公室激情啪啪| 国产av无码专区亚洲av麻豆| 欧美狠狠入鲁的视频777色| 麻豆精产国品| 免费福利视频二区三区| 久久精品国产色蜜蜜麻豆国语版| 性大毛片视频| 免费国产黄线在线播放| 日本成人中文字幕亚洲一区| 手机在线免费av资源网| 欧美性巨大╳╳╳╳╳高跟鞋| 国产成人综合久久久久久| 青青草成人原视频在线播放视频| 国产自产二区三区精品| 四虎影视久久久免费观看| 高潮毛片无遮挡高清免费| 久久久久久人妻一区二区无码Av| 国模91九色精品二三四| 五月av综合av国产av| 久久精品亚洲乱码伦伦中文| 九九九影院| 亚洲精品国产成人久久av盗摄| 色综合久久中文娱乐网| 91精品一区国产高清在线gif| 一区二区三区婷婷在线|