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

基于降維算法和等效桿長的可展結構精度分析

2017-11-22 02:05:57祁俊威王春潔丁建中
航空學報 2017年6期

祁俊威, 王春潔,2,*, 丁建中

1.北京航空航天大學 機械工程及自動化學院, 北京 100083 2.北京航空航天大學 北京市數字化設計與制造重點實驗室, 北京 100083

基于降維算法和等效桿長的可展結構精度分析

祁俊威1, 王春潔1,2,*, 丁建中1

1.北京航空航天大學 機械工程及自動化學院, 北京 100083 2.北京航空航天大學 北京市數字化設計與制造重點實驗室, 北京 100083

考慮鉸鏈間隙和桿件尺寸誤差的不確定性并通過概率統計方法對其進行研究,提出了一種基于單變量降維算法(UDRM)和等效桿長模型的可展結構精度分析方法。利用UDRM將可展結構的精度性能函數解耦為多個桿件尺寸誤差的獨立作用形式,建立精度分析模型。引入等效桿長模型,等效桿件替代原桿件進行精度計算。將鉸鏈間隙與原始桿件尺寸誤差合并到等效桿件的尺寸誤差中,同時證明了等效桿長尺寸誤差近似服從正態分布。以某衛星可展開天線為算例,結合高斯求積公式求解展開狀態下精度指標的分布期望和方差。通過與蒙特卡羅模擬(MCS)和一次二階矩(FOSM)法計算結果的對比驗證了本文精度分析方法的正確性和高效性。

鉸鏈間隙; 尺寸誤差; 可展開結構; 精度分析; 降維算法; 等效桿長

可展開結構在空間衛星中應用越來越廣泛,對其工作構型精度的要求越來越高。空間可展開結構通常具有尺寸大等特點,其在軌工作精度和機械誤差緊密相關。可展開結構多設計為鉸接多閉環連桿機構,其機械誤差主要來源于構件間運動副間隙以及構件自身尺寸誤差[1-4]。這些誤差不僅使可展開結構實際工作狀態與理論設計狀態存在一定偏差,嚴重時還可能造成其在軌工作失效。

鉸鏈間隙和桿件尺寸誤差使得可展開結構工作階段呈現出結構不確定性[5-6],衡量不確定性的大小可以借助隨機變量期望和方差的概念。現今已有多種求解功能函數響應期望和方差的方法,如一次二階矩(First Order Second Moment, FOSM)法[7-8]、蒙特卡羅模擬(Monte Carlo Simulation, MCS)[9-10]、響應面法(Response Surface Method, RSM)[11-12]等。上述方法中FOSM法需已知功能函數的梯度信息,但當功能函數為隱式或形式復雜時梯度運算困難;MCS為了獲得較高精度的結果所需計算量巨大,耗時長;RSM面對非線性較高的工程問題需要用到高階多項式進行擬合,但多項式的擬合范圍有限從而影響到結果精度。Rahman和Xu提出了一種單變量降維算法(Univariate Dimension Reduction Method, UDRM)[13],該方法能夠避免對功能函數梯度的求解,可極大降低計算工作量,在結構可靠度應用方面取得了很好的效果[14-16]。

本文提出了單變量降維和等效桿長模型的可展開結構精度分析方法。將可展開結構中的桿件尺寸誤差和鉸鏈間隙視作隨機變量,建立可展開結構精度分析的通用模型。以某可展開天線為研究對象,結合高斯求積公式求解精度指標的期望和方差。最后通過與蒙特卡羅模擬和一次二階矩法對比驗證本文精度分析方法的正確性和高效性。

1 單變量降維算法

假設n維向量S=[s1s2…sn]T為系統輸入變量,g(S)為由S決定的系統輸出。則g(S)降維展開表達式為

g12…n(s1,s2,…,sn)

(1)

式中:g0為對g(S)產生的零階影響;gi(si)為變量si單獨對g(S)產生的一階影響;gi1i2(si1,si2)為變量si1和si2聯合對g(S)產生的二階影響;gi1i2…ik(si1,si2,…,sik)為變量si1,si2,…,sik共同對g(S)產生的k階影響;g12…n(s1,s2,…,sn)為所有變量對g(S)產生的影響。

選取隨機變量空間中任一參考點S0=[s01s02…s0n]T,并將g(S)在參考點處展開可得

(2)

式中:

(3)

可求得輸出響應的期望和方差為

(4)

式中:E(·)為期望算子;D(·)為方差算子;

(5)

其中:ρ為隨機變量的概率密度函數。

當g為隱式函數或表達式復雜時,若直接對式(5)進行積分運算,計算十分困難。因此本文利用高斯求積公式對誤差函數的期望和方差進行求解,式(5)可轉化為

(6)

2 等效桿長模型分析

當可展開結構構型變得復雜時,桿件和鉸鏈數目增多,單變量降維方法得到的展開項也隨之增多,直接求解計算量增大。為了減少變量數目、降低計算量,引入等效桿長模型。等效桿件由一根桿件和其所連接的一個或兩個鉸鏈間隙組成,將其作為一個整體計算等效桿長,再將等效桿長替代原始桿長進行精度指標響應的求解,簡化計算過程。

2.1 等效桿件尺寸確定

無質量虛擬連桿法是分析含間隙機構運動特性的主要手段之一[17]。將間隙用無質量的連桿代替,分析時只需確定虛擬桿的長度和方位,來表示間隙矢量的大小和方向即可。機構中由鉸鏈連接的兩個構件,如圖1所示,銷軸外徑為ra,襯套內徑為rb,可以得到間隙圓半徑大小為rc=rb-ra。

連桿兩端均有鉸鏈間隙與只有一端有鉸鏈間隙兩種情況下的分析方法相同。下面以連桿兩端均含有鉸鏈間隙為例,等效桿長示意圖如圖2所示,其計算方法為[18]

(7)

式中:li為原始桿長;Li為等效桿長;x1i、y1i和x2i、y2i分別為以襯套中心為局部坐標的兩個銷軸的圓心坐標值,銷軸圓心在襯套內位置隨機分布,應滿足

(8)

圖1 含間隙鉸鏈模型Fig.1 Joint clearance model

圖2 等效桿長模型Fig.2 Effective link length model

2.2 等效桿長誤差分析

現對等效桿長模型進行誤差分析,首先在各參數名義值處進行Taylor級數二階展開,可得

(Δy2-Δy1)2]+O

(9)

式中:F=(x2-x1+l)(Δx2-Δx1+Δl)+(y2-y1)(Δy2-Δy1);Δx1、Δy1、Δx2、Δy2和Δl分別為x1、y1、x2、y2和l的誤差值;O為高階微小量,計算中可以忽略;(1/L)3?(1/L),計算中可以被忽略;第3項為高次項,計算時也可以被忽略。因此在計算等效桿長誤差ΔL時展開至一階項就能夠保證計算精度。

間隙尺寸x1、y1和x2、y2可看做參數名義值均為0的隨機變量,則L的誤差可以化簡為

(10)

因此易得L誤差的期望與方差為

(11)

通常假設構件尺寸的加工誤差服從正態分布,則其尺寸的期望和方差有

(12)

li~N(μli,σli)i=1,2,…,n

在空間無自重無外力的自由漂浮狀態下,可以假設銷軸中心在間隙圓內均勻分布,可得銷軸圓心坐標的概率密度函數為[6]

ρ(xi)=

(13)

式中:rci為第i個鉸鏈的間隙圓半徑;xi為間隙沿x軸的分量。yi與xi具有相同形式的概率密度函數,由此計算可得xi和yi的期望和方差為

(14)

2.3 等效桿長誤差分布確定

Δx1、Δx2和Δl為相互獨立的誤差變量,則ΔL的概率密度函數可通過如下卷積的方式進行計算:

(15)

式中:*為卷積符號。

進而可以得到等效桿長誤差的卷積形式為

ρΔx1,Δx2,Δl(ΔL)=ρΔx1*ρΔx2*ρΔl(ΔL)

(16)

將式(13)代入式(16),發現間隙尺寸部分的卷積是一個瑕點積分的形式,即使采用數值積分也難以進行求解,導致等效桿長誤差ΔL概率密度函數的精確表達式難以求得。因此本節采用Edgeworth級數對其概率密度函數進行估計[19-20]。

概率密度函數Edgeworth級數展開式定義為

(17)

式中:φ(·)為標準正態分布概率密度函數;ΔL*為經標準正態化變換的等效桿長誤差,即ΔL*=(ΔL-μΔL)/σΔL;κn為ΔL第n階中心距,其表達式為

(18)

(19)

通常利用隨機變量前4階矩信息進行擬合,能夠得到較好的近似結果[21-22]。計算可得ΔL的概率分布Edgeworth級數估計的前4階展開式為

[(ΔL*)4-6(ΔL*)2+3]

(20)

式中:航天領域大型衛星制造公差一般小于0.05 mm,中小衛星制造公差一般小于0.03 mm,本文取rc為0.2 mm進行分析;ΔL*經標準正態化處理后取值主要落在(-3,3)之間。現對2項取值進行估計

(ΔL*)4-6(ΔL*)2+3=

((ΔL*)2-3)2-6∈[-6,30)

計算可得式(20)中的第2項取值在[-5×10-5,2.5×10-4)之間,要遠小于1,在實際應用中可以忽略不計。因此ΔL*的概率密度分布估計可化簡為

f(ΔL*)≈φ(ΔL*)

(21)

一組桿件尺寸誤差和鉸鏈間隙的分布參數如表1所示。

由圖3可知兩種方法的結果吻合很好,證明了本節推導方法的正確性,等效桿長的誤差服從正態分布,因此在求解誤差函數期望和方差時可利用Gauss-Hermite求積公式。等效桿長替代原始桿長和鉸鏈間隙,最終誤差函數的期望和方差表示為

(22)

表1桿件尺寸誤差和鉸鏈間隙分布參數

Table1Distributionparametersoflinklengtherrorandjointclearance

VariableDistributionformExpectation/mmVariance/mm2Jointclearance(0.1mm)Uniformincircle00.0025LengtherrorNormal00.01

圖3 MCS和Edgeworth級數推導結果比較 Fig.3 Comparison of results of MSC and Edgeworth series

式中:L為各桿件等效長度組成的向量;μL為等效桿長均值組成的向量。

本文方法的計算流程如圖4所示。

圖4 精度分析模型計算流程圖Fig.4 Flow chart of calculation with precision analysis model

3 算例分析

3.1 天線展開構型

本文以某型號可展開天線為例進行分析。可展開天線由1套可折疊的空間支撐桁架來保證其剛度和精度。天線在軌展開之后由支撐桁架固定在展開工作狀態。支撐桁架關于展開方向與天線板垂線方向所構成的平面對稱,在進行展開狀態精度分析時可簡化為平面桿系處理。該天線展開構型如圖5所示。

圖5 可展開天線的展開構型Fig.5 Deployment state of deployable antenna

以根鉸鏈b1轉動中心為坐標原點建立如圖5所示的全局坐標系。鉸鏈b1和d1與星體固連,d3和d4均為復合鉸鏈。鉸鏈d2、c1和c2處設置鎖定機構,達到180° 時自動鎖死。

3.2 指向精度評定指標

引入位置偏差和角度偏差兩個指標對天線板展開結束后的指向精度進行量化評判,其原理如圖6所示。

取每塊天線板框架與支撐桁架連接點作為測量的擬合點,即a1、a2、a3和a4,計算其在引入桿件尺寸誤差情況下的坐標值。4個擬合點的橫縱坐標組成向量X和Y,即

式中:xai為ai點橫坐標值;yai為ai點縱坐標值。

圖6中兩塊天線板理想位置表示為直線P0,兩塊天線板實際位置通過線性回歸算法擬合后表示為直線P1,其方程為

α1x+y+α2=0

(23)

式中:x和y為擬合平面上點的坐標值,α1和α2為平面方程參數,其計算方法為

位置偏差定義為直線P1與坐標軸y的截距,其表達式為

epos=-α2

(24)

角度偏差定義為直線P1與平面P0的夾角,其表達式為

圖6 位置偏差及角度偏差測量原理Fig.6 Measurement principles for position and angular errors

eang=arctan(-α1)

(25)

在理想狀態下(不考慮鉸鏈間隙和桿件尺寸誤差),位置偏差和角度偏差值為0 mm和0°;在實際狀態下(考慮鉸鏈間隙和桿件尺寸誤差),位置偏差和角度偏差會偏離理想值,且偏離越多表示誤差越大。

3.3 計算結果及分析

天線各桿件原始尺寸及其誤差分布參數如表2 所示。

間隙對鎖定鉸鏈的影響較小,且當鎖定角度為π時,由間隙導致的鎖定角度偏差對等效桿長的影響也十分微小,因此本文在計算指向精度時將鎖定鉸鏈連接的兩桿視為一個固定構件[6]。按照一定的工藝方法,可假設所有的鉸鏈間隙都具有同樣的尺寸,令所有間隙圓半徑rc=0.05 mm。通過單變量降維算法和Gauss-Hermite積分可以直接求解精度指標的期望和方差。計算結果如表3 所示。

表2天線桿件參數統計特征

Table2Parametricstatisticalcharacteristicsofantennalinks

LinknameOriginalvalue/mmErrorexpectation/mmErrorvariance/mmld1d3467300.01ld3d4 21000.01la1d3480100.01la2d3197500.01la3d4197500.01la4d4480100.01

表3 天線指向精度計算結果Table 3 Calculation results of antenna pointing precision

從表3可以看出,FOSM法和本文方法與MCS相比位置偏差期望的絕對誤差分別為0.18 μm 與0.14 μm,角度偏差期望的絕對誤差分別為(1.90×10-5)°與(3×10-8)°;位置偏差方差的相對誤差分別為4.84%和0.81%,角度偏差方差的相對誤差分別為10.68%和0.17%,在期望估計方面3種方法結果相差不多,但是對方差估計的結果,本文方法與蒙特卡羅法吻合更好,說明了本文精度分析方法的正確性;同時比較3種方法計算所需的迭代次數可以看出,本文方法所需的迭代次數遠小于另兩種方法,說明了在保證結果精度的前提下該方法具有高效性。

4 結 論

1) 提出了一種基于單變量降維算法和等效桿長模型的可展開結構精度分析方法,該方法無需精度指標函數的梯度信息。

2) 引入等效桿長模型,將鉸鏈間隙和桿件尺寸誤差轉化到桿件等效桿長誤差上,等效桿長作為中間變量與可展結構誤差函數建立聯系,可減少計算變量數目、降低迭代次數。

3) 以某衛星可展開天線為算例,通過與MCS和FOSM法的計算結果對比,驗證了本文方法的正確性與高效性。

4) 本文方法具有流程簡單、易于編程、求解速度快等優點,對于可展開結構的精度分析具有普遍適用性,可為后續分析與設計提供依據。

[1] 劉榮強, 田大可, 鄧宗全. 空間可展開天線結構的研究現狀與展望[J]. 機械設計, 2010, 27(9): 1-10.

LIU R Q, TIAN D K, DENG Z Q. Research actuality and prospect of structure for space deployable antenna[J]. Journal of Machine Design, 2010, 27(9): 1-10 (in Chinese).

[2] 劉明治, 高桂芳. 空間可展開天線結構研究進展[J]. 宇航學報, 2003, 24(1): 82-87.

LIU M Z, GAO G F. Advances in the study on structure for space deployable antenna[J]. Journal of Astronautics, 2003, 24(1): 82-87 (in Chinese).

[3] CHEN G L, WANG H, LIN Z Q. A unified approach to the accuracy analysis of planar parallel manipulators both with input uncertainties and joint clearance[J]. Mechanism & Machine Theory, 2013, 64(6): 1-17.

[4] MERLET J P. Computing the worst case accuracy of a PKM over a workspace or a trajectory[C]//5th Chemnitzer Parallel Kinematic Seminar. Zwickau: Verlag Wissenschaftliche Scripten, 2006.

[5] SUN D Y, CHEN G P. Kinematic accuracy analysis of planar mechanisms with clearance involving random and epistemic uncertainty[J]. European Journal of Mechanics—A/Solids, 2016, 58: 256-261.

[6] LI X, DING X L, CHIRIKJIAN G S. Analysis of angular-error uncertainty in planar multiple-loop structures with joint clearances[J]. Mechanism & Machine Theory, 2015, 91: 69-85.

[7] 姚澤良, 李寶平, 周雪峰. 結構可靠度分析的一次二階矩方法與二次二階矩方法[J]. 西北水力發電, 2005, 21(3): 20-23.

YAO Z L, LI B P, ZHOU X F. Simple and quadratic of the two ranks quadrature of the structure reliability[J]. Journal of Northwest Hydroelectric Power, 2005, 21(3): 20-23 (in Chinese).

[8] 陳勝軍, 賈方. 曲柄滑塊機構運動精度的概率分析與計算[J]. 機械設計與制造, 2013(7): 203-206.

CHEN S J, JIA F. Probability analysis and calculation of kinematic accuracy for slider-crank mechanism[J]. Machinery Design & Manufacture, 2013(7): 203-206 (in Chinese).

[9] PADMANABHAN D, AGARWAL H, RENAUD J E, et al. A study using Monte Carlo simulation for failure probability calculation in reliability-based optimization[J]. Optimization & Engineering, 2006, 7(3): 297-316.

[10] 吳建云, 王春潔, 汪瀚. 基于蒙特卡洛法的衛星天線板展開精度分析[J]. 航天返回與遙感, 2013, 34(6): 89-94.

WU J Y, WANG C J, WANG H. Accuracy analysis of satellite antenna plate deployment based on Monte Carlo method[J]. Spacecraft Recovery & Remote Sensing, 2013, 34(6): 89-94 (in Chinese).

[11] 彭茂林, 楊自春, 曹躍云, 等. 基于響應面法的可靠性穩健設計優化[J]. 航空動力學報, 2013, 28(8): 1784-1790.

PENG M L, YANG Z C, CAO Y Y, et al. Reliability robust design optimization based on response surface method[J]. Journal of Aerospace Power, 2013, 28(8): 1784-1790 (in Chinese).

[12] 劉成立, 呂震宙. 結構可靠性分析中考慮高次項修正的組合響應面法[J]. 航空學報, 2006, 27(4): 594-599.

LIU C L, LV Z Z. Response surface combination improved by high order term for structure reliability analysis[J]. Acta Aeronautica et Astronautica Sinica, 2006, 27(4): 594-599 (in Chinese).

[13] RAHMAN S, XU H. A univariate dimension-reduction method for multi-dimensional integration in stochastic mechanics[J]. Probabilistic Engineering Mechanics, 2004, 19(4): 393-408.

[14] ZHANG X F, PANDEY M D, ZHANG Y M. A numerical method for structural uncertainty response computation[J]. Science China Technological Sciences, 2011, 54(12): 3347-3357.

[15] WANG J G, ZHANG J F, DU X P. Hybrid dimension reduction for mechanism reliability analysis with random joint clearances[J]. Mechanism & Machine Theory, 2011, 46(10): 1396-1410.

[16] 孟廣偉, 馮昕宇, 李鋒,等. 基于降維算法和Edgeworth級數的結構可靠性分析[J]. 北京航空航天大學學報, 2016, 42(3): 421-425.

MENG G W, FENG X Y, LI F, et al. Structural reliability analysis based on dimensionality reduction and Edgeworth series[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(3): 421-425 (in Chinese).

[17] EARLES S W E, WU C L S. Motion analysis of rigid-link mechanism with clearance at a bearing using Lagrangian mechanics and digital computation[C]//Conference on Mechanisms. London: Institution of Mechanical Engineers, 1972.

[18] LEE S J, GILMORE B J. The determination of the probabilistic properties of velocities and accelerations in kinematic chains with uncertainty[J]. Journal of Mechanical Design, 1991,113(1): 84-90.

[19] BLOZNELIS M. An Edgeworth expansion for studentized finite population statistics[J]. Acta Applicandae Mathematicae, 2003, 78(1): 51-60.

[20] 張義民, 賈敬存, 黃賢振. 基于Edgeworth級數法和數據包絡分析法的數控車床可靠性分配[J]. 機械強度, 2016(1): 69-73.

ZHANG Y M, JIA J C, HUANG X Z. Reliability allocation of CNC lathe based on Edgeworth series method and date envelopment analysis[J]. Journal of Mechanical Strength, 2016(1): 69-73 (in Chinese).

[21] 葉江水, 王仲剛, 陳友良, 等. 基于前四階矩的非高斯響應概率密度函數逼近方法研究[J]. 后勤工程學院學報, 2010, 26(1): 12-16.

YE J S, WANG Z G, CHEN Y L, et al. Research on approximate methods of non-Gaussian probability density function based on the first four moments of the response[J]. Journal of Logistical Engineering University, 2010, 26(1): 12-16 (in Chinese).

[22] 劉彥明. 基于四階矩的機械可靠性相關研究[D]. 太原: 太原科技大學, 2011: 23-38.

LIU Y M. Research on mechanical reliability based on the four order moment[D]. Taiyuan: Taiyuan University of Science & Technology, 2011: 23-38 (in Chinese).

(責任編輯: 徐曉)

Precision analysis of deployable structures based on dimensionreduction method and effective link length

QIJunwei1,WANGChunjie1,2,*,DINGJianzhong1

1.SchoolofMechanicalEngineeringandAutomation,BeihangUniversity,Beijing100083,China2.BeijingKeyLaboratoryofDigitalDesignandManufacturing,BeihangUniversity,Beijing100083,China

The uncertainties of joint clearances and link length errors are studied by the method of probability and statistics. A precision analysis method for deployable structure is proposed based on Univariate Dimension Reduction Method (UDRM) and effective link length model. Using the UDRM, the precision function for the deployable structure is decoupled into a combination of independent effects of multiple link length errors to establish the precision analysis model for the structure. The effective link length model is applied to replace the original link length for precision calculation. The effective model converts the joint clearances and link length errors into effective link length errors, which are proved to follow normal distributions. An example of deployable antenna is given to calculate the means and variances in the deployable state with the Gauss quadrature based on the error distributions of link lengths and joint clearances. The correctness and effectiveness of the precision analysis method is verified by comparing the results of Monte Carlo Simulation (MCS) and First Order Second Moment (FOSM) method.

joint clearance; length error; deployable structure; precision analysis; dimension reduction method; effective link length

2016-07-06;Revised2016-08-29;Accepted2016-10-26;Publishedonline2016-11-101418

URL:www.cnki.net/kcms/detail/11.1929.V.20161110.1418.006.html

NationalNaturalScienceFoundationofChina(51635002)

2016-07-06;退修日期2016-08-29;錄用日期2016-10-26; < class="emphasis_bold">網絡出版時間

時間:2016-11-101418

www.cnki.net/kcms/detail/11.1929.V.20161110.1418.006.html

國家自然科學基金 (51635002)

*

.E-mailwangcj@buaa.edu.cn

祁俊威, 王春潔, 丁建中. 基于降維算法和等效桿長的可展開結構精度分析J. 航空學報,2017,38(6):220590.QIJW,WANGCJ,DINGJZ.PrecisionanalysisofdeployablestructuresbasedondimensionreductionmethodandeffectivelinklengthJ.ActaAeronauticaetAstronauticaSinica,2017,38(6):220590.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0273

V11; TH115

A

1000-6893(2017)06-220590-08

*Correspondingauthor.E-mailwangcj@buaa.edu.cn

主站蜘蛛池模板: 天天色天天操综合网| 精品91在线| 日韩高清在线观看不卡一区二区| 国产91蝌蚪窝| 日本三级精品| 2020最新国产精品视频| 国产精品无码影视久久久久久久 | 久久久国产精品免费视频| 国产香蕉在线| 中文国产成人精品久久| 国产视频自拍一区| 国产啪在线91| 久久综合婷婷| 国产九九精品视频| 夜精品a一区二区三区| 免费人成在线观看视频色| 四虎在线高清无码| 国产XXXX做受性欧美88| 一级全免费视频播放| 久久国产精品娇妻素人| 日韩美女福利视频| 欧美视频在线观看第一页| 亚洲人在线| 国产乱人视频免费观看| 免费视频在线2021入口| 69国产精品视频免费| 色135综合网| 欧美不卡视频在线观看| 波多野结衣的av一区二区三区| 手机在线免费不卡一区二| P尤物久久99国产综合精品| 欧亚日韩Av| 四虎成人在线视频| 91久久国产成人免费观看| 久久99国产视频| 无码精品一区二区久久久| 日韩经典精品无码一区二区| 国产成人精品无码一区二| 久久综合丝袜日本网| 成人午夜亚洲影视在线观看| 一边摸一边做爽的视频17国产| 国产一区二区三区日韩精品| 亚洲第七页| 久久成人国产精品免费软件 | 久久久久九九精品影院| 福利视频一区| 欧美午夜视频在线| 中文字幕天无码久久精品视频免费 | 国产精品夜夜嗨视频免费视频| 国产成人三级| 91精品福利自产拍在线观看| 成人免费视频一区| 无码精品福利一区二区三区| 国产精品久久久久久久伊一| 国产在线视频二区| 91网址在线播放| 国产欧美日韩综合在线第一| 国产无码精品在线播放| 亚洲狼网站狼狼鲁亚洲下载| 久久国产精品无码hdav| 免费日韩在线视频| 男女精品视频| 91小视频在线观看免费版高清| 日韩欧美国产精品| 自拍亚洲欧美精品| 国内精品自在自线视频香蕉 | 超碰免费91| 人妻丰满熟妇啪啪| 91精品国产情侣高潮露脸| 国产手机在线ΑⅤ片无码观看| 美美女高清毛片视频免费观看| 亚洲综合狠狠| 久青草免费在线视频| www.亚洲一区| 亚洲最新在线| 国产亚洲欧美在线中文bt天堂| 国产乱子伦手机在线| 欧美日韩在线第一页| 美女被操91视频| 日本一区二区三区精品AⅤ| 亚洲香蕉久久| 欧美在线国产|