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

        ?

        FAST饋源艙傾覆的理論分析及試驗研究*

        2015-03-22 11:30:45李銘哲
        天文研究與技術(shù) 2015年4期
        關(guān)鍵詞:鋼索饋源索力

        李銘哲,李 輝

        (1. 貴州大學理學院,貴州 貴陽 550025;2. 中國科學院國家天文臺,北京 100012; 3. 中國科學院射電天文重點實驗室,北京 100012)

        FAST饋源艙傾覆的理論分析及試驗研究*

        李銘哲1,2,3,李 輝2,3

        (1. 貴州大學理學院,貴州 貴陽 550025;2. 中國科學院國家天文臺,北京 100012; 3. 中國科學院射電天文重點實驗室,北京 100012)

        FAST望遠鏡饋源艙通過6根鋼索懸浮于空中,在運行過程中其傾角連續(xù)變化,存在發(fā)生傾覆的風險。通過艙索系統(tǒng)的靜力學理論分析和模型試驗對FAST饋源艙的最大傾覆角進行了研究?;谂搩A角最大的優(yōu)化原則和艙-索系統(tǒng)靜力平衡,對索力和姿態(tài)角的限定等約束條件,建立了求解艙最大傾角(傾覆角)的目標優(yōu)化函數(shù),利用牛頓迭代法解得了艙在其運行軌跡面(焦面)的傾覆角,同時研究了重心改變對傾覆角的影響。通過艙索系統(tǒng)3 m縮尺模型對理論分析結(jié)果進行了驗證,得到了模型艙在給定位置的傾覆角,試驗結(jié)果與理論計算結(jié)果比較吻合。最后分析了饋源艙原型在整個焦面上的傾覆角,結(jié)果表明饋源艙在實際工作過程中是安全的,不會發(fā)生傾覆。

        FAST望遠鏡;饋源艙;傾覆;靜力學分析;模型試驗

        從1994年中國天文學家提出建造500 m口徑球面射電望遠鏡(Five-hundred-meter Aperture Spherical radio Telescope, FAST)以來,國家天文臺一直在全力推進其建設,F(xiàn)AST立項建議書在2007年7月10日得到國家發(fā)展改革委員會批準,現(xiàn)在于貴州省平塘縣大窩凼地區(qū)全面建設中,建成后將成為世界上最大口徑的單天線射電望遠鏡,在未來20年-30年保持領(lǐng)先地位。FAST擁有3項中國的自主創(chuàng)新技術(shù)——洼地臺址、主動變形反射面和輕型饋源柔索拖動系統(tǒng),突破了在地面建造全可動天線的百米工程極限[1]。

        FAST的輕型饋源柔索拖動系統(tǒng)又稱為饋源支撐鋼索驅(qū)動系統(tǒng),是一個空間尺度巨大的柔索牽引并聯(lián)機器人,其機構(gòu)主體是6根直徑約50 mm的鋼索及其相應的卷揚機,每根鋼索和卷揚機依附于一座鋼塔作為支撐架,6座塔均高達百余米,并按軸對稱位置分布在直徑600 m的圓周上。6根鋼索共同支撐并牽引作為并聯(lián)機器人終端控制平臺的饋源艙在空中運動,使得饋源艙能夠攜帶艙內(nèi)設備(接收機)跟蹤反射面的瞬時焦點,實現(xiàn)望遠鏡對目標星體的觀測。饋源艙在空中的全部運動軌跡是一個曲率半徑160 m,口徑約206 m且與反射面同球心的球冠面,稱為饋源焦面,如圖1右下角。FAST饋源艙重約30 t,外接圓直徑約13 m,是一個集結(jié)構(gòu)、機構(gòu)、測量、控制等相關(guān)技術(shù)于一體的光機電一體化的復雜系統(tǒng),是FAST的核心部件,其主要功能是克服風擾和系統(tǒng)的其他擾動,通過饋源艙內(nèi)的二次調(diào)整裝置,實現(xiàn)饋源接收機的精調(diào)定位,如圖1右上角。

        饋源艙由鋼索驅(qū)動在饋源焦面上運行時,姿態(tài)會發(fā)生變化,特別是從饋源焦面中心向邊緣移動的過程中,傾角不斷增大,因此運行過程中可能出現(xiàn)傾覆的危險,在運行過程特別是極限位置下是否會傾覆對觀測及設備安全影響甚大,因此有必要對運行過程中饋源艙的傾斜狀態(tài)進行研究。

        FAST饋源艙的姿態(tài)與索力相關(guān),文[2-5]對索牽引并聯(lián)機構(gòu)的靜力學問題進行了分析。文[2]提出了鋼索虛牽的判斷標準及索力均勻的原則,并提出了確定Stewart 平臺工作空間的方法;文[3]通過四繩索系統(tǒng)研究了塔分布圓半徑對柔索牽引并聯(lián)機器人工作空間及索力的影響,提出了新的判斷準則解決繩索虛牽的問題;文[5]研究了進艙纜線對柔索牽引并聯(lián)機器人工作空間的影響以及塔高度及饋源艙重心的改變對優(yōu)化結(jié)果的影響。上述文獻均基于索力優(yōu)化進行求解,未曾考慮饋源艙的傾覆問題,亦未考慮饋源艙重心的改變對饋源艙最大傾覆角的影響。

        圖1 FAST全景圖及剖視圖

        本文試圖通過理論建模的方法求解處于給定位置的饋源艙臨界傾覆角,然后通過縮尺模型試驗驗證理論分析的結(jié)果?;谏鲜龇椒ㄗ詈笄蠼庠椭叙佋磁摰呐R界傾覆角,分析其是否處于安全范圍。

        1 艙臨界傾覆角的理論分析

        以饋源艙和6條支撐鋼索構(gòu)成的艙-索系統(tǒng)為研究對象,分析饋源艙在不同姿態(tài)角時的靜力平衡狀態(tài)及6索的張力,求解饋源艙保持靜力平衡的最大傾角,即饋源艙的臨界傾覆角。

        理論分析基于以下假設:

        (1)饋源艙等效為6自由度剛體;

        (2)支撐索簡化為柔軟不可伸長、具有一定質(zhì)量且只能單向受拉的軟繩;

        (3)鋼索曲線近似為拋物線。

        如圖2,在艙體質(zhì)心位置按右手螺旋分別建立局部坐標系ox0y0z0和oxyz,全局坐標系為OXYZ,全局坐標系以主動反射面最低點為原點,Z軸方向向上。局部坐標系以模型艙的質(zhì)心作為原點,z軸正向豎直向上,坐標系ox0y0z0固連于艙體并隨艙姿態(tài)的改變而旋轉(zhuǎn),oxyz不隨艙姿態(tài)的改變而變化,則可以建立艙體靜力平衡方程組。

        圖2 FAST索力分析坐標系示意圖

        (1)式為空間力系沿X軸、Y軸、Z軸方向的平衡方程,(2)式為繞X軸、Y軸、Z軸的空間力矩平衡方程。Hi為第i根索張力的水平分量;αi為第i根索在XOY平面內(nèi)的投影與X軸正方向的夾角;βi為第i個索艙接點處索張力矢量與XOY平面的夾角;G為饋源艙重量,如圖(3),xm、ym分別為饋源艙重心在全局坐標下的坐標;由此可以建立第i根懸索的靜力微分方程:

        (3)

        以支撐索左端為坐標原點,水平方向為x軸,如圖4。

        圖3 饋源艙靜力平衡分析

        圖4 索模型靜力分析

        得(3)式的拋物曲線近似解為

        (4)

        因此可以建立夾角βi與Hi的函數(shù)關(guān)系式為

        (5)

        其中符號ai、hi分別為從第i座支撐塔出索點到第i個索艙接點的相對水平距離和垂直高度;q為單位長度索自重;Hi為第i根繩索的拉力,如圖4,Tani表示第i根繩索對應的高度差與跨度的比值,即Tani=hi/ai。

        (1)式與(2)式中xi、yi和zi為第i個索艙接點的坐標,它們與局部坐標x0i、y0i、z0i的關(guān)系為[6]

        (6)

        式中的局部坐標x0i、y0i和z0i為已知量,僅取決于索艙接點在錨接平面的相對位置;θ、φ、φ表示饋源艙旋轉(zhuǎn)的3個獨立的空間姿態(tài)角,分別表示艙體的俯仰角、方位角和自旋角。如圖3,θ定義為Z軸與z0軸之間的夾角,即本文需要求解的最大傾角;φ定義為z0軸在xoy平面內(nèi)的投影與x軸正方向的夾角;φ定義為星形桁架繞z軸的旋轉(zhuǎn)角;Sφ表示角度φ的正弦函數(shù);Cφ表示角度φ的余弦函數(shù),其它的表示方法與文[6]類似。

        2 基于目標函數(shù)優(yōu)化方法求解艙臨界傾覆角

        結(jié)合(1)、(2)、(5)、(6)式可以得到包含9個未知量(Hi,θ,φ,φ,i=1, 2 … 6)和6個方程的非線性方程組,在這種情況下解一般不是唯一的,需要通過優(yōu)化分析的方法從中選擇一組最佳解,因此需要定義優(yōu)化的準則,建立優(yōu)化分析的目標函數(shù)。

        本文的目的是尋找在指定位置處滿足艙靜力平衡條件的艙最大傾覆角,即:

        θ=θmax,

        (7)

        為了保證函數(shù)的連續(xù)性及良好的可導性,目標函數(shù)用下列式子表示:

        θ.

        (8)

        鋼索是柔性的,不能承載壓力,于是有一個最小力,同時其拉力也不能無限大,值的選取根據(jù)實際設計方案而定,設定索力的上下限為

        Hmin≤Hi≤Hmax.

        (9)

        當饋源艙的傾角超過90°即可認為饋源艙已傾覆,方位角的范圍為0到360°,根據(jù)文[6],饋源艙自旋角基本在0°,在試驗中發(fā)現(xiàn)模型艙的自旋角很小,不超過15°,因此本文中的約束范圍為0到15°。

        當采用目標函數(shù)優(yōu)化方法求解方程組時,原來的9個未知量就變成了優(yōu)化變量,(8)式就變成了對優(yōu)化變量的強制約束,同時考慮(1)、(2)式以及3個姿態(tài)角的范圍,總約束如下:

        (10)

        通過罰函數(shù)法,構(gòu)建用于求解的無約束目標函數(shù)如下:

        對于(11)式的最小值優(yōu)化可通過牛頓迭代法求解。在迭代求解過程中需要不斷地計算一階導數(shù)(Jacobi矩陣)及二階導數(shù)(Hessian矩陣),為了判斷程序是否結(jié)束,在每次計算中需要比較函數(shù)值的差值是否小于給定的精度σ,用下式表示:

        (12)

        其中下標k表示第k步迭代。對應有9個未知量的方程,對應的圖形很復雜,要尋找全局的最小值,初始值的選擇至關(guān)重要,首先利用網(wǎng)格法尋找滿足條件的解[7],選擇傾角最大的那一組解作為牛頓迭代法初始值,選定初始值后為了更快地收斂,在步長的選擇方面通過linesearch方法確定[8-9],判斷是否滿足步長選擇條件的公式:

        (13)

        上式是為了保證方程f朝著減小的方向選擇步長,因此式中c很小,在本試驗中令c=10-4,為第k步待求的步長,計算步長的具體步驟為:

        (2)計算θ(α0),如果滿足(13)式,αk=α0,結(jié)束,否則進入步驟(3);

        (3)在步驟(2)的基礎(chǔ)上利用已知的3個信息θ(α0)、θ(0)、θ'(0),構(gòu)建一個平方插值多項式來擬合θ(αk),其表達式如下:

        '(0)a1+θ(0),

        (14)

        令該方程式取最小值的解為

        (15)

        計算θ(αk1)=f(xk+αk1pk),判斷是否滿足(13)式,如果滿足則令αk=αk1,程序結(jié)束,否則轉(zhuǎn)到步驟(4);

        (4)構(gòu)建一個3次插值多項式,使其最小化的值計算方式如下:

        (16)

        (17)

        (18)

        計算θ(aki+1)=f(xk+aki+1pk),判斷是否滿足(13)式,如果滿足則令αk=aki+1,程序結(jié)束,否則一直進行步驟(4),直到滿足(13)式。

        在計算過程中為了保證每次計算的αk與上次計算的αk-1相比變化不大,因此做了如下的規(guī)定:如果αk?0.5αk-1,令αk=0.5αk-1;如果αk?0.1αk-1,令αk=0.1αk-1。在得到步長的基礎(chǔ)上對函數(shù)f的最優(yōu)化求解,步驟如下:

        (1)輸入系統(tǒng)基本參數(shù);

        (2)判斷重心ZZ是否小于給定值,例如10mm,如果成立進行步驟(3),否則進行步驟(6);

        (3)輸入初值x0,通過計算旋轉(zhuǎn)矩陣R, a, h, Tan,tanβ等,確定fk,通過linesearch法計算步長ak,如果有合適的ak,轉(zhuǎn)到步驟(4),否則轉(zhuǎn)到步驟(5);

        (4)判斷解向量增量是否滿足判斷條件,即‖Δ{x}‖,滿足進入步驟(5),否則用計算得到的xk代替xk-1,繼續(xù)步驟(3);

        (5)重心ZZ加1,轉(zhuǎn)到步驟(2);

        (6)差值計算重心連續(xù)變化過程中對應傾角的值;

        (7)結(jié)束。

        流程圖如圖5。

        3 艙索系統(tǒng)模型優(yōu)化結(jié)果及試驗

        如圖6,在室內(nèi)搭建一個3m尺度的實驗模型,主要模擬FAST原型的艙-索系統(tǒng)。試驗模型的主要參數(shù)按照相似原則進行選定,如表1。

        3m模型與FAST原型之間遵循幾何相似比:1/200,質(zhì)量相似比1∶178570;其中質(zhì)量相似比是根據(jù)幾何相似比與單位長度索自重比得到,模型主要參數(shù)均由此兩個相似比參數(shù)導出。

        為了充分體現(xiàn)饋源艙在整個饋源焦面上最大傾角,本文選擇了如圖2的3個白圈位置,分別為饋源焦面的最低點或中心點①,邊緣點③和中間點位置②,根據(jù)相似原則,3m模型的索力范圍為0到2.1N,模型饋源艙重心移動到鉸接點平面以上10mm時,圖2(a)中δ,相當于原型饋源艙重心在鉸接點平面以上2m,已基本達到原型饋源艙在鉸接點平面以上的最大高度,因此只考慮重心最大增加到10mm。重心位置從0開始,理論計算結(jié)果如圖7。

        圖5 優(yōu)化求解流程圖

        Fig.5Aflowchartoftheroutinetofindtheoptimalsolutionofacriticaltiltforoverturningofthefocuscabin

        表1 3 m模型及FAST原型分析中所采用的主要幾何和力學參數(shù)

        圖6 3 m模型全景圖及模型艙示意圖

        圖7 選定位置處最大傾角隨重心的變化圖

        圖7表示選定的位置處,隨著重心的改變,在索力允許的范圍內(nèi),模型艙所能達到的最大角度,由圖可知,最大傾角從位置①到③呈現(xiàn)逐漸增大的趨勢,其最大傾角在位置①處收斂到85.04°左右,位置②處收斂到85.68°,位置③處收斂到85.78°,且重心在0到10 mm范圍內(nèi)的改變對所處位置的最大傾角幾乎沒有影響。

        圖8中3張圖片均采用雙縱坐標軸表示,橫坐標表示重心的改變,上圖縱軸表示索力,中圖的縱軸表示角度,圖8顯示的重心改變過程中索力及姿態(tài)角相對于重心為0時對應量的改變量,通過對索力及姿態(tài)角隨重心變化的分析可知,當重心變化時,索力、方位角及自旋角變化明顯,饋源艙的最大傾角不隨重心的改變而變化是由于鋼索的索力及另外兩個姿態(tài)角的變化綜合導致的。

        通過計算選定位置處傾角最大的情況下的6根索長,調(diào)節(jié)3 m模型中所對應的鋼絲繩的長度,通過測量模型艙的高度差的方式,得到了3 m模型在選定的位置處,如圖2,模型艙重心在鉸接點以上10 mm時的傾角,結(jié)果如表2,試驗測得在3個位置處的饋源艙傾角(重心位置10 mm)實測值與理論值的對比。

        3 m模型試驗所得數(shù)據(jù)接近理論計算的最大傾角,且索力均在給定的范圍內(nèi),說明理論分析結(jié)果是可靠的。誤差的出現(xiàn)可能因為模型索艙錨節(jié)點處為固性鏈接以及鋼絲繩本身的剛度較大引起,如圖9。

        最后分析了FAST原型中饋源艙在整個饋源焦面上的臨界傾覆角,索力上下限分別為0及40 t,計算中使用的參數(shù)如表1,得到的結(jié)果如圖10(a),其中10(b)是FAST原型基于索力優(yōu)化的結(jié)果[5]。

        圖8 選定位置處索力及姿態(tài)角隨重心的變化圖

        Fig.8 Variations of cable tensions and attitude angles with the center of gravity (COG) at the chosen locations

        表2 3個試驗點處的饋源艙傾角

        從圖10(a)可知,在整個饋源焦面上,饋源艙在索力允許的范圍內(nèi)的最大傾角呈現(xiàn)從最低點到饋源焦面邊緣逐漸增大的趨勢,其計算結(jié)果最小值出現(xiàn)在最低點,為84.95°,最大值出現(xiàn)在焦面邊緣,為86.84°,與饋源艙傾角在焦面最低點處最小,邊緣處最大的趨勢相符。

        從圖10(b)可知,F(xiàn)AST在執(zhí)行天文觀測任務的過程中饋源艙傾角在0~17.5°范圍,遠小于艙傾覆角,因而饋源艙在整個饋源焦面上是安全的。

        圖9 位置③處重心為10 mm時3 m模型艙的傾角圖

        圖10 FAST原型中饋源艙臨界傾覆角在饋源焦面的分布

        4 結(jié) 論

        本文通過理論分析和模型試驗對FAST饋源艙在饋源焦面上運行過程中的臨界傾覆角及是否存在傾覆的可能性進行了研究。通過模型試驗驗證了理論分析結(jié)果的可靠性。

        理論分析結(jié)果表明,饋源艙臨界傾覆角的最小值出現(xiàn)在焦面最低點,為84.95°,最大值出現(xiàn)在焦面邊緣為86.84°,饋源艙在實際運行過程中最大傾角約為17.5°,遠遠小于最大傾角的最小值,因此饋源艙在整個饋源焦面上運行過程中是安全的。饋源艙重心在一定范圍內(nèi)的改變對艙臨界傾覆角幾乎沒有影響。

        [1] Nan Rendong. Five hundred meter aperture spherical radio telescope (FAST) [J]. Science in China: Series G Physics, Mechanics & Astronomy, 2006, 49(2): 129-148.

        [2] 孫欣, 段寶巖. 巨型柔性Stewart平臺解空間、工作空間的研究及懸索張力的優(yōu)化分析[J]. 機械工程學報, 2002, 38(2): 16-21. Sun Xin, DuanBaoyan. Study on solution space, working space and cable tension’s optimized analysis for huge flexible stewart platform[J]. Chinese Journal of Mechanical Engineering, 2002, 38(2): 16-21.

        [3] 姚蕊, 唐曉強, 李鐵民, 等. 大型射電望遠鏡饋源定位3T索牽引并聯(lián)機構(gòu)分析與設計[J]. 機械工程學報, 2007, 43(11): 105-109. Yao Rui, Tang Xiaoqiang, Li Tiemin, et al. Analysis and design of 3T cable-driven parallel manipulator for the feedback’s orientation of the large radio telescope [J]. Chinese Journal of Mechanical Engineering, 2007, 43(11): 105-109.

        [4] 黃亮, 朱文白, 唐曉強, 等. 大射電望遠鏡索牽引并聯(lián)機構(gòu)索力優(yōu)化分析[J]. 天文研究與技術(shù)——國家天文臺臺刊, 2010, 7(3): 268-276. Huang Liang, Zhu Wenbai, Tang Xiaoqiang, et al. Optimization of cable tensions of a cable-driven parallel manipulator for a large radio telescope[J]. Astronomical Research & Technology——Publications of National Astronomical Observatories of China, 2010, 7(3): 268-276.

        [5] Li Hui, Nan Rendong, K?rcher Hans. Working space analysis and optimization of the main positioning system of FAST cabin suspension, in ground-based and airborne telescopes II[C]// Proceedings of SPIE. 2008.

        [6] 李輝, 朱文白, 潘高峰. 基于索力優(yōu)化的FAST柔索牽引并聯(lián)機構(gòu)的靜力學分析[J]. 工程力學, 2011, 28(4): 185-207. Li Hui, ZhuWenbai, Pan Gaofeng. Equilibrium analysis of FAST rope-drive parallel manipulator based on rope force optimization[J]. Engineering Mechanics, 2011, 28(4): 185-207.

        [7] 朱永生. 實驗物理中的概率和統(tǒng)計[M]. 北京: 科學出版社, 2006: 512-513.

        [8] Dennis J E, Schnabel R B. Numerical methods for unconstrained optimization and nonlinear equations[M]. Englewood Cliffs: Prentice-Hall, 1983.

        [9] Nocedal J, Wright S J. Numerical Optimization[M]. New York: Springer-Verlag, 1999: 35-61.

        A Theoretical Analysis and an Experimental Study of Overturning of the Focus Cabin of the FAST

        Li Mingzhe1,2,3, Li Hui2,3

        (1. College of Science, Guizhou University, Guiyang 550025, China; 2. National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China; 3. Key Laboratory of Radio Astronomy, Chinese Academy of Sciences, Beijing 100012, China, Email: limingzhe@nao.cas.cn)

        The focus cabin of the FAST (Five-hundred meter Aperture Spherical radio Telescope) is supported and driven by six parallel flexible steel cables, so that it is suspended in the air and traces a spherical focal surface. As a result, the orientation of the focus cabin is adaptive to the variable spatial position of the focus cabin. The adaptation nevertheless raises the risk of overturning of the focus cabin during its motion. In this paper we first theoretically analyze the mechanics of overturning of the focus cabin and present verification of the analysis through experiments on a 3m scaled model of the FAST focus cabin-cable system. We then establish an object function of optimization of the critical tilt of overturning of the focus cabin by taking the static-equilibrium equations and the limits on cable tensions/cabin orientations as the constraint conditions. We employ the Newton Iteration Method to obtain optimal solutions of the critical tilts using the function. We have also evaluated and tested the effect of the location of the center of gravity (COG) of the focus cabin on the critical tilts. Our test results show accurate consistency with our theoretical analysis. We finally present the calculation results of the critical tilts of the actual focus cabin at all its accessible locations on the focal surface. The results show that a cabin tilt under a normal circumstance is much less than the critical tilt of overturning of the focus cabin. This means the focus cabin is normally safe.

        FAST telescope; Focus cabin; Overturning; Static force analysis; Model test

        國家自然科學基金 (10973023, 11103046, 211311001) 資助.

        2015-03-06;修定日期:2015-04-13

        李銘哲,男,碩士. 研究方向:索牽引并聯(lián)機構(gòu)動力學. Email: limingzhe@nao cas cn

        O312; TH751; TP242

        A

        1672-7673(2015)04-0424-09

        CN 53-1189/P ISSN 1672-7673

        猜你喜歡
        鋼索饋源索力
        江蘇索力得新材料集團有限公司
        “動中通”衛(wèi)星天線的饋源優(yōu)化設計
        科技傳播(2019年22期)2020-01-14 03:06:28
        她用兩年給“天眼”減重
        科教新報(2019年16期)2019-09-10 01:50:38
        她用兩年給“天眼”減重
        科學導報(2019年24期)2019-09-03 04:33:02
        剩余強度廣義應力與橋梁鋼索的破斷——論檢測索力不能評定鋼索的服役安全性
        上海公路(2017年1期)2017-07-21 13:38:33
        漢字走鋼索
        FAST饋源艙大尺寸同軸度測量方法
        鋼索式液壓提升裝置在電力建設工程中的應用
        預應力鋼絞線網(wǎng)加固混凝土橋梁的索力分布試驗研究
        基于拉索振動特征的索力檢測與評估方法
        综合人妻久久一区二区精品| 五月婷婷俺也去开心| 二区三区视频| 久久视频在线视频精品| 手机在线观看av资源| 亚洲欧美一区二区成人片| 内射后入在线观看一区| 99免费视频精品| 国产成人亚洲系列毛片| 无码国产精品久久一区免费| 免费无码国产v片在线观看| 红杏性无码免费专区| 蜜桃夜夜爽天天爽三区麻豆av| 东京热人妻系列无码专区 | 情色视频在线观看一区二区三区| 国产精品日韩av一区二区三区| 日韩精品极品视频在线观看免费| 色偷偷88888欧美精品久久久| 国产美女黄性色av网站| 午夜国产精品视频在线观看| 欧美人与动人物牲交免费观看久久| 亚洲综合性色一区| 自拍偷拍一区二区三区四区| 人妻精品在线手机观看| 夜夜躁狠狠躁2021| 久久一日本道色综合久久大香| 美利坚合众国亚洲视频| 99在线精品视频在线观看| 亚洲天堂2017无码中文| 亚洲伊人成综合人影院| 在线精品国产亚洲av蜜桃| 50岁熟妇的呻吟声对白| 欧美日韩国产另类在线观看| 日韩在线一区二区三区中文字幕 | 精品免费久久久久久久| 四虎欧美国产精品| 青青青草视频手机在线| 国产精品麻豆va在线播放| 国产av国片精品| 日本一区二区三区四区在线看| 在线视频国产91自拍|