• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    高速飛行器熱結構工作時變模態(tài)參數(shù)辨識

    2015-06-24 13:48:37周思達劉莉李昱霖周小陳
    航空學報 2015年1期
    關鍵詞:模態(tài)有限元結構

    周思達, 劉莉, 李昱霖, 周小陳

    北京理工大學 宇航學院, 北京 100081

    高速飛行器熱結構工作時變模態(tài)參數(shù)辨識

    周思達*, 劉莉, 李昱霖, 周小陳

    北京理工大學 宇航學院, 北京 100081

    高速飛行器由于其很高的飛行速度而無可避免地受到氣動加熱作用的影響,進而引起結構特性的時變。采用理論或有限元方法(FEM)進行數(shù)值分析,難以獲取反映結構在飛行(工作)狀態(tài)下的真實模態(tài)參數(shù)。通過辨識獲取高速飛行器熱環(huán)境下的時變結構模態(tài)參數(shù)是一項十分具有挑戰(zhàn)性的任務。針對此問題,引入?yún)?shù)化時頻域的最大似然方法,對氣動加熱作用下的高速飛行器升力面結構的時變模態(tài)參數(shù)進行了辨識。通過模擬真實飛行狀態(tài)的數(shù)值算例研究,說明參數(shù)化時頻域的最大似然方法能夠很好地辨識出低信噪比(SNR)情況下的模態(tài)頻率和模態(tài)振型,驗證了參數(shù)化時頻域最大似然方法適用于具有顯著時變特征的高速飛行器熱結構的時變結構模態(tài)參數(shù)辨識,可為將來相關的工程研究和應用提供良好的理論支持。

    時變; 模態(tài)參數(shù); 參數(shù)辨識; 高速飛行器; 熱結構; 參數(shù)化; 最大似然

    高速飛行器在速度上的顯著優(yōu)勢,因此在應用領域受到了廣泛關注。關于高速飛行器的研究已經(jīng)成為相關領域的熱點。當然,高速飛行器也因其高飛行速度給設計和分析帶來了眾多挑戰(zhàn)。如何準確地獲取高速飛行器在氣動加熱環(huán)境下的熱結構時變模態(tài)參數(shù)就是挑戰(zhàn)之一[1]。首先,高速飛行器在飛行中才會產(chǎn)生相應的氣動加熱效應,因此,相關的模態(tài)分析是工作狀態(tài)下的;第二,由于熱環(huán)境隨時間改變,高速飛行器結構在熱作用下的動力學特性,即模態(tài)參數(shù),是隨時間變化的。另外,雖然有限元方法(FEM)等數(shù)值方法能夠?qū)Ω咚亠w行器熱環(huán)境下的結構動力學特性進行預測,但無論是理論還是有限元等數(shù)值方法,無法最為真實地反映復雜結構在環(huán)境作用下的工作模態(tài)參數(shù)。這就需要基于實驗的工作模態(tài)分析手段,來獲取結構真實的模態(tài)參數(shù)[2]。最近,也有學者進行了高速飛行器熱結構動力學特性的實驗研究。吳振強等[3]利用實驗得到了不同溫度下的飛行器壁板結構動特性;白云鶴[4]采用非參數(shù)化時頻域的方法對熱結構模態(tài)頻率進行了辨識。

    在基于響應的時變結構模態(tài)參數(shù)辨識領域,近年來產(chǎn)生了許多新方法,其中,Poulimenos和Fassois[5]提出的泛函序列時間相關自回歸滑動平均(FS-TARMA)方法、Zhou等[6]提出的參數(shù)化時頻域兩步最小二乘方法、Xu等[7]提出的基于小波變換的非參數(shù)化時頻域方法等具有代表性。

    本文將利用所提出的參數(shù)化時頻域最大似然方法[8]對高速飛行器熱結構在工作狀態(tài)下的結構模態(tài)參數(shù)進行辨識,建立模擬真實情況的數(shù)值算例,并通過算例討論方法在不同測量噪聲情況下的性能,說明方法對高速飛行器熱結構模態(tài)參數(shù)辨識的適用性。為將來高速飛行器熱結構的實驗研究,特別是模態(tài)參數(shù)辨識提供一定的理論參考和技術支持。

    1 時變結構的時頻域參數(shù)化動力學模型

    線性時不變結構的傳遞函數(shù)可利用參數(shù)化的公分母模型來表達[9]。假設作用在結構上的激勵為白噪聲,結構響應的功率譜函數(shù)也可以表示為參數(shù)化的公分母模型。將時不變結構功率譜函數(shù)的公分母模型拓展為時間相關的公分母模型,以表示時變結構的時間相關功率譜函數(shù),如:

    (1)

    式中:tτ和ωf分別為采樣時刻和角頻率,下標τ=1,2,…,Nτ表示采樣時間,f=1,2,…,Nf表示采樣頻率;k=1,2,…,NoNr(即k=(i-1)No+o),No和Nr分別為輸出點和參考點數(shù)量;分子多項式和分母多項式皆定義為

    (2)

    (3)

    (4)

    進一步寫為

    (5)

    為方便以下推導,將式(5)的參數(shù)向量寫為

    (6)

    式中:θ為參數(shù)化的公分母模型的待定參數(shù)向量。

    另外,如式(2)所示,分子多項式和分母多項式都含有時頻基函數(shù),其數(shù)值特性直接影響后續(xù)辨識的數(shù)值條件,某些情況下甚至影響辨識的成敗。本文采用一種正交多項式與z域多項式的混合形式(Orthogonal Polynomial andz-domain Mapping hybrid,OPZM)的時頻基函數(shù)[8],即時間多項式采用正交多項式(如Chebyshev多項式、Jacobi多項式等),頻率多項式采用z域多項式:

    (8)

    式中:ω為角頻率;Ts為采樣長度。

    2 參數(shù)化時頻域的最大似然辨識

    根據(jù)頻域的最大似然費用函數(shù)和時頻域時變結構動力學參數(shù)化的時間相關公分母模型(式(1)),時頻域的最大似然辨識方法的費用函數(shù)為

    (9)

    為得到較大的動態(tài)范圍,采用以下對數(shù)形式的總方程誤差:

    (10)

    矩陣形式的總方程誤差可寫為

    (11)

    與線性最小二乘辨識不同,由于如式(10)所示的總方程誤差是關于待估參數(shù)θ的非線性函數(shù),因此最大似然辨識方法的費用函數(shù)在最小化時需求解非線性代數(shù)方程,可采用Gauss-Newton法。

    根據(jù)Gauss-Newton迭代,在第m個迭代步,待估參數(shù)的增量Δθm可由如下線性方程求得:

    (12)

    式中:J為m個迭代步時的Jacobi矩陣。

    在每一個迭代步后對待估參數(shù)θ進行更新,如:

    θm+1=θm+Δθm

    (13)

    待估參數(shù)θ的初始值由參數(shù)化時頻域的最小二乘辨識給出[10]。利用更新的待估參數(shù)θ,可以得到更新的Jacobi矩陣和更新的總方程誤差(式(10))。對應Jacobi矩陣(如式(12)定義)的分塊矩陣分別為

    (14)

    (15)

    其中,

    (16)

    (17)

    根據(jù)式(14)和式(15)重寫式(12),得

    (18)

    式(18)中的前NoNr個方程可寫為

    (19)

    消去Δβk,由式(18)可得

    (20)

    求解方程式(20)可得當前迭代步中待估參數(shù)α的增量Δα,并利用式(19)回代得到待估參數(shù)β的增量Δβ。按照式(6)將Δα和Δβ寫成待估參數(shù)θ的增量Δθ。每一個迭代步都是一個求解縮減的正則方程的過程,迭代的終止條件為前后兩個迭代步的所有響應的功率譜函數(shù)估計值與測量值之間相對誤差的最大值之差小于ζ,定義為

    (21)

    一旦得到待估參數(shù)向量θ,便可計算出各個時刻的模態(tài)參數(shù),具體方法參考文獻[11]。

    3 氣動加熱下的升力面結構模態(tài)參數(shù)辨識

    本文提出的方法能夠針對工作狀態(tài)下的結構,辨識其時變的結構模態(tài)參數(shù)。下面通過一個數(shù)值仿真算例,說明其針對高速飛行器結構在時變熱環(huán)境下,如何辨識結構隨時間變化的模態(tài)參數(shù)。

    為實現(xiàn)與“凍結”的結構進行比較,從而驗證所提方法,算例將不考慮結構變形對氣動力和熱傳導的影響,結構受到熱的影響是單向的。另外,不考慮熱應力導致的結構微分剛度對結構整體剛度的影響,僅考慮材料受溫度的影響。這些假設在實驗中都是難以滿足的,但是對于方法的驗證至關重要。當然,以上忽略的效應對結構的影響,在實際中可以包含在結構特性隨時間變化的整體效應中,并且本文采用的參數(shù)化時頻域最大似然辨識方法基于結構響應,反映結構與外部因素作用的宏觀模態(tài)特性,而不依賴復雜耦合問題內(nèi)在機理的建模,因此,本文數(shù)值算例的假設仍不失一般性。

    本文算例參考文獻[12],建立如圖1所示的升力面。

    圖1 升力面外形尺寸Fig.1 Profile of lifting surface

    3.1 氣動加熱計算

    對于高速飛行器,氣動加熱是其顯著特點。為模擬真實狀態(tài)下升力面的氣動加熱效應,采用高精度CFD計算方法,獲取結構表面的溫度場邊界,并用于后續(xù)計算。假設飛行器的飛行狀態(tài)曲線如圖2所示(速度和攻角曲線)[13-14]。

    圖2 飛行狀態(tài)曲線[13-14]Fig.2 Flight state curves[13-14]

    在一系列時刻點,利用Fastran氣動加熱分析軟件進行CFD計算仿真,獲得升力面上下表面的溫度場分布,參考文獻[14]。以10 s時刻為例,上表面(迎風面)的溫度分布如圖3所示。

    圖3 t=10 s時刻升力面上表面的溫度分布Fig.3 Temperature distribution on upper face of lifting surface at t=10 s

    3.2 分布溫度邊界下的瞬態(tài)熱傳導計算

    為得到升力面結構在上下表面分布溫度邊界情況下的結構各處溫度,采用有限元方法進行瞬態(tài)熱傳導分析。用于瞬態(tài)熱傳導的有限元模型如圖4所示,單元為四節(jié)點的熱傳導單元(DS4)。

    圖4 升力面有限元模型Fig.4 Finite element model of lifting surface

    由于在表面處,CFD和有限元的網(wǎng)格并不一一對應,因此,需利用CFD的計算結果(以圖3為例),在結構上下表面有限元節(jié)點上分別進行二維插值,從而得到各個時刻(前20 s)上下表面的分布溫度邊界,結構材料屬性如表1所示。

    表1 結構材料屬性

    由于材料的熱傳導系數(shù)和熱容系數(shù)與溫度無關,因此,在材料特性和分布溫度邊界條件已知的情況下,利用線性有限元的瞬態(tài)熱傳導分析,計算得到結構各節(jié)點處的溫度分布。翼稍處某節(jié)點溫度隨時間的變化曲線如圖5所示。

    圖5 翼稍處某節(jié)點溫度隨時間的變化曲線Fig.5 Temperature vs time at a certain point of wing tip

    3.3 狀態(tài)“凍結”下的參考熱模態(tài)計算

    由于材料的彈性模量與溫度相關,導致結構剛度與溫度相關,因此,為驗證本文采用的參數(shù)化時頻域的最大似然辨識方法,利用一系列狀態(tài)“凍結”的結構獲取其模態(tài)參數(shù)與前者比較。狀態(tài)“凍結”指某時刻下給定由熱傳導分析得到的分布溫度場,邊界條件為固定翼根部(固定節(jié)點的平動和轉(zhuǎn)動自由度)。另外,為得到顯著的時變特性,本文將結構的彈性模量設置為一個較大的范圍,如表1所示。

    結構有限元模型網(wǎng)格與熱傳導有限元模型一致,如圖4所示。單元類型為四節(jié)點減縮積分的Mindlin殼單元,厚度方向積分點個數(shù)為5。利用有限元模態(tài)分析(Lanczos特征值提取方法),獲取各時刻下的結構模態(tài)。模態(tài)頻率如圖6所示;10 s時刻第一、二階模態(tài)振型如圖7所示。

    圖6 狀態(tài)“凍結”下的各時刻模態(tài)頻率Fig.6 “Frozen”-configuration modal frequency at different time instants

    圖7 10 s時刻的第一、二階模態(tài)振型Fig.7 Mode shapes for the first and second modes at 10 s

    3.4 溫度變化下的結構動力學響應計算

    為辨識提供數(shù)據(jù),采用基于有限元方法的結構瞬態(tài)響應分析,獲取結構的響應信號。結構邊界條件與有限元模態(tài)分析一致。激勵點位于如圖4所示的三角指示處,激勵方向為沿坐標z軸方向,激勵形式為白噪聲(辨識時假設未知)。利用Newmark-β法計算響應(Newmark-β積分參數(shù)γ=0.5,β=0.25),積分步長為0.000 5s(約為第二階模態(tài)頻率最大值的37倍)。

    為模態(tài)參數(shù)辨識提供響應數(shù)據(jù),分別在結構表面的6個位置(如圖4中的圓點所示)處z方向“測量”速度響應,采樣頻率為200 Hz。為模擬與真實情況較為貼近的測量噪聲污染的情況,對“測量”的響應加入信噪比(SNR)分別為10 dB (10%測量噪聲)和3 dB (50%測量噪聲)的噪聲。信噪比分別為10 dB和3 dB時,響應點1處的響應如圖8所示;信噪比分別為10 dB和3 dB時,響應點1處的時間相關自功率譜函數(shù)(采用平滑偽Wigner-Ville分布(SPWVD)計算[15])如圖9所示。

    圖8 響應點1處的速度響應歷程Fig.8 Velocity history at point 1

    圖9 響應點1處的時間相關自功率譜函數(shù)Fig.9 Time-dependent auto power spectra at point 1

    3.5 工作狀態(tài)熱結構模態(tài)參數(shù)辨識

    3.4節(jié)通過有限元的瞬態(tài)響應分析,得到了“測量”的結構響應。這些響應將為結構模態(tài)參數(shù)辨識提供最基本的數(shù)據(jù)。

    利用第2節(jié)的參數(shù)化時頻域的最大似然方法進行辨識。首先,利用時頻分析(此處用SPWVD)估計各響應的時間相關功率譜函數(shù),并且為了估計振型,使用響應的互功率譜,即以響應點1為參考點,響應點1~6分別與其計算互譜,如圖9所示(自譜為自身的互譜的特例)。

    圖10 時變模態(tài)頻率的辨識結果Fig.10 Estimated modal frequency of time-varying structure

    圖11 辨識得到的10 s時刻的模態(tài)振型Fig.11 Estimated shapes for the first and second modes at 10 s

    如圖10所示,采用參數(shù)化時頻域的最大似然方法,能夠很好地辨識出在氣動加熱環(huán)境下,結構時變的模態(tài)頻率,甚至在信噪比較低(50%噪聲,信噪比為3 dB)的情況下,除很少的間斷外,辨識值和有限元方法得到的理論參考值十分吻合。另外,比較圖7和圖11,發(fā)現(xiàn)辨識得到的模態(tài)振型也十分接近于理論振型,取得較好的結果,同樣可以得到其余任意時刻的模態(tài)振型,在此限于篇幅,僅給出10 s時刻的模態(tài)振型。

    4 結 論

    本文利用參數(shù)化時頻域的最大似然方法對氣動加熱作用下的高速飛行器升力面結構的時變模態(tài)參數(shù)進行了辨識。不同于傳統(tǒng)短時的方法,需要將響應數(shù)據(jù)分段進行一系列時不變模態(tài)參數(shù)辨識,本方法在時頻域上整體地一次辨識出結構隨時間變化的模態(tài)參數(shù),數(shù)值算例說明其不僅能夠很好地辨識出低信噪比情況下的模態(tài)頻率,還能準確地辨識出結構的模態(tài)振型,并驗證了參數(shù)化時頻域最大似然辨識方法適用于具有顯著時變特征的高速飛行器熱結構的時變結構模態(tài)參數(shù)辨識。另外,為實現(xiàn)對比,本文采用的算例僅考慮了氣動熱作用對結構的單向作用,但是辨識方法本身是基于結構響應的,辨識結果可反映結構與外部因素作用的宏觀模態(tài)特性,能夠包括結構與氣動和熱耦合等。因此,本方法仍可為實際的的工程研究和應用提供良好的理論支持。

    [1] Yu Y H, Li S, Wang Y D. Study on dynamic characteristics of metallic thermal protection system with thermal environ[J]. Journal of Vibration, Measurement & Diagnosis, 2013, 33(1): 171-175 (in Chinese). 余艷輝, 李書, 王遠達. 熱環(huán)境下的金屬熱防護系統(tǒng)的動力學特性研究[J]. 振動測試與診斷, 2013, 33(1): 171-175.

    [2] Zhang L, Brincker R, Andersen P. An overview of operational modal analysis: major development and issues[C]∥1st International Operational Modal Analysis Conference (IOMAC), 2005: 179-190.

    [3] Wu Z Q, Cheng H, Zhang W, et al. Effects of thermal environment on dynamic properties of aerospace vehicle[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(2): 334-342 (in Chinese). 吳振強, 程昊, 張偉, 等. 熱環(huán)境對飛行器壁板結構動特性的影響[J]. 航空學報, 2013, 34(2): 334-342.

    [4] Bai Y H. Research of structural modal experiment and parameter identification for high speed vehicles under flying environments[D]. Harbin: Harbin Institute of Technology, 2013 (in Chinese). 白云鶴. 高速飛行器飛行環(huán)境下結構模態(tài)試驗及參數(shù)識別方法研究[D]. 哈爾濱: 哈爾濱工業(yè)大學, 2013.

    [5] Poulimenos A G, Fassois S D. Output-only stochastic identification of a time-varying structure via functional series TARMA models[J]. Mechanical Systems and Signal Processing, 2009, 23(4): 1180-1204.

    [6] Zhou S D, Heylen W, Sas P, et al. Parametric modal identification of time-varying structures and the validation approach of modal parameters[J]. Mechanical Systems and Signal Processing, 2014, 47(1): 94-119.

    [7] Xu X, Shi Z Y, You Q. Identification of linear time-varying systems using a wavelet-based state-space method[J]. Mechanical Systems and Signal Processing, 2012, 26(6): 91-103.

    [8] Zhou S D, Heylen W, Sas P, et al. Maximum likelihood estimator of operational modal analysis for linear time-varying structures in time-frequency domain[J]. Journal of Sound and Vibration, 2014, 333(11): 2339-2358.

    [9] Verboven P. Frequency-domain system identification for modal analysis[D]. Brussels: Vrije Universiteit Brussel, 2002.

    [10] Zhou S D, Liu L, Yang W, et al. Matrix fraction polynomial model-based least square estimation of modal parameters for linear time-varying structures[J]. Journal of Vibration and Shock, 2014, 33(6): 118-123 (in Chinese). 周思達, 劉莉, 楊武, 等. 基于矩陣分式多項式時變結構模態(tài)參數(shù)最小二乘辨識[J]. 振動與沖擊, 2014, 33(6): 118-123.

    [11] Zhou S D. Study on theory and experiment of time-frequency-domain modal parameter estimation for linear time-varying mechanical structures[D]. Beijing: Beijing Institute of Technology, 2012 (in Chinese). 周思達. 線性時變結構時頻域模態(tài)參數(shù)辨識理論及實驗研究 D]. 北京: 北京理工大學, 2012.

    [12] McNamara J J, Friedmann P P, Powell K G, et al. Aeroelastic and aerothermoelastic behavior in hypersonic flow [J]. AIAA Journal, 2008, 46(10): 2591-2610.

    [13] Zhu H, Liu L, Zhou S, et al. Integrated aerodynamic thermal structure design optimization method of lifting surfaces [J]. Journal of Aircraft, 2012, 49(5): 1521-1526.

    [14] Li Y L. Studies on aerodynamic-structural-thermal multidisciplinary design optimization and method of flight vehicles[D]. Beijing: Beijing Institute of Technology, 2014 (in Chinese). 李昱霖. 氣動熱結構多學科分析及高效優(yōu)化策略研究 [D]. 北京: 北京理工大學, 2014.

    [15] Roshan-Ghias A, Shamsollahi M B, Mobed M, et al. Estimation of modal parameters using bilinear joint time-frequency distributions[J]. Mechanical Systems and Signal Processing, 2007, 21(5): 2125-2136.

    Tel: 010-68918752

    E-mail: zhousida@bit.edu.cn

    *Corresponding author. Tel.: 010-68918752 E-mail: zhousida@bit.edu.cn

    Operational identification of time-varying modal parameters for thermal structures of high-speed aerial vehicles

    ZHOU Sida*, LIU Li, LI Yulin, ZHOU Xiaochen

    SchoolofAerospaceEngineering,BeijingInstituteofTechnology,Beijing100081,China

    The high-speed aerial vehicles are unavoidably impacted by the aerodynamic heating effect due to their high speed, which may further result in the time-varying characteristics of their structures. The theoretical and finite element method (FEM)-based numerical approaches hardly acquire the real modal parameters of the in-flight (operational) aerial vehicles. Focusing on this problem, the parametric time-frequency-domain maximum likelihood method is introduced into the application of the high-speed aerial vehicles and the time-varying modal parameters are estimated by this parametric time-frequency-domain maximum likelihood method. Based on a high-fidelity simulation example, it is proved that the introduced method can identify the modal frequency and mode shapes even in the low signal noise ratio (SNR) cases and the method is suitable for the modal parameter estimation of thermal structures of high-speed aerial vehicles with notable time-varying natures, which can theoretically support the relative engineering studies and applications.

    time-varying; modal parameter; parameter identification; high-speed aerial vehicles; thermal structures; parametric; maximum likelihood

    2014-07-01; Revised: 2014-09-10; Accepted: 2014-09-23; Published online: 2014-09-26 10:41

    s: National Natural Science Foundation of China (11372036, 11402022); Beijing Institute of Technology Foundation for Basic Research (20120142009)

    2014-07-01; 退修日期: 2014-09-10; 錄用日期: 2014-09-23; 網(wǎng)絡出版時間: 2014-09-26 10:41

    www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0267.html

    國家自然科學基金(11372036, 11402022); 北京理工大學基礎研究基金(20120142009)

    Zhou S D, Liu L, Li Y L, et al. Operational identification of time-varying modal parameters for thermal structures of high-speed aerial vehicles[J].Acta Aeronautica et Astronautica Sinica, 2015, 36(1): 373-380. 周思達, 劉莉, 李昱霖, 等.高速飛行器熱結構工作時變模態(tài)參數(shù)辨識[J].航空學報, 2015, 36(1): 373-380.

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

    10.7527/S1000-6893.2014.0267

    V214.1; TH113.1

    A

    1000-6893(2015)01-0373-08

    周思達 男, 博士, 講師。主要研究方向:飛行器結構動力學、時變結構動力學建模與辨識。

    *通訊作者.Tel.: 010-68918752 E-mail: zhousida@bit.edu.cn

    URL: www.cnki.net/kcms/detail/10.7527/S1000-6893.2014.0267.html

    猜你喜歡
    模態(tài)有限元結構
    《形而上學》△卷的結構和位置
    哲學評論(2021年2期)2021-08-22 01:53:34
    論結構
    中華詩詞(2019年7期)2019-11-25 01:43:04
    論《日出》的結構
    國內(nèi)多模態(tài)教學研究回顧與展望
    創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
    基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
    磨削淬硬殘余應力的有限元分析
    由單個模態(tài)構造對稱簡支梁的抗彎剛度
    計算物理(2014年2期)2014-03-11 17:01:39
    基于SolidWorks的吸嘴支撐臂有限元分析
    箱形孔軋制的有限元模擬
    上海金屬(2013年4期)2013-12-20 07:57:18
    欧美乱色亚洲激情| 少妇裸体淫交视频免费看高清| 成人三级黄色视频| or卡值多少钱| 欧美黄色淫秽网站| 天堂av国产一区二区熟女人妻| 白带黄色成豆腐渣| 久久精品国产亚洲av香蕉五月| a级毛片a级免费在线| 欧美乱妇无乱码| 精品日产1卡2卡| 久久久久免费精品人妻一区二区| 男人舔奶头视频| 国产91精品成人一区二区三区| 欧美zozozo另类| 国产成人一区二区三区免费视频网站| 在线国产一区二区在线| 欧美日本亚洲视频在线播放| 性色av乱码一区二区三区2| 非洲黑人性xxxx精品又粗又长| 中文字幕久久专区| 国产高清视频在线播放一区| 女生性感内裤真人,穿戴方法视频| 禁无遮挡网站| 人妻丰满熟妇av一区二区三区| 黑人巨大精品欧美一区二区mp4| 国产不卡一卡二| 国产黄色小视频在线观看| 很黄的视频免费| 日日摸夜夜添夜夜添小说| 一级黄色大片毛片| 国产真人三级小视频在线观看| 国产一级毛片七仙女欲春2| 国产成人系列免费观看| 国产一区二区激情短视频| 亚洲精华国产精华精| 免费在线观看亚洲国产| 国产1区2区3区精品| 看片在线看免费视频| 欧美色欧美亚洲另类二区| 51午夜福利影视在线观看| 淫秽高清视频在线观看| 91av网站免费观看| 国产一区二区在线av高清观看| 黄色视频,在线免费观看| 中文字幕熟女人妻在线| 久久久久久国产a免费观看| 亚洲一区高清亚洲精品| 国产精品爽爽va在线观看网站| 极品教师在线免费播放| or卡值多少钱| 我要搜黄色片| 精品免费久久久久久久清纯| 黄频高清免费视频| 男女之事视频高清在线观看| 亚洲欧美一区二区三区黑人| 亚洲黑人精品在线| 亚洲精品久久国产高清桃花| 黄色日韩在线| 一个人看视频在线观看www免费 | 97超视频在线观看视频| 99精品欧美一区二区三区四区| 亚洲18禁久久av| 女生性感内裤真人,穿戴方法视频| 国产成人福利小说| 国产精品女同一区二区软件 | 露出奶头的视频| 日韩有码中文字幕| 免费观看精品视频网站| www.自偷自拍.com| 久久精品影院6| 免费搜索国产男女视频| 日韩欧美国产在线观看| 三级男女做爰猛烈吃奶摸视频| 午夜福利18| 香蕉av资源在线| 成人无遮挡网站| 婷婷丁香在线五月| 欧美日韩瑟瑟在线播放| 色哟哟哟哟哟哟| 国产成年人精品一区二区| 天天躁日日操中文字幕| 欧洲精品卡2卡3卡4卡5卡区| 老司机在亚洲福利影院| 国产精品一及| 又大又爽又粗| 成熟少妇高潮喷水视频| 亚洲午夜精品一区,二区,三区| 欧美又色又爽又黄视频| 一个人免费在线观看的高清视频| 别揉我奶头~嗯~啊~动态视频| 午夜福利成人在线免费观看| 亚洲欧美激情综合另类| 丁香欧美五月| 国产一区二区在线观看日韩 | 日本撒尿小便嘘嘘汇集6| 97超级碰碰碰精品色视频在线观看| 国产成人av教育| 在线看三级毛片| 欧美中文日本在线观看视频| 欧美丝袜亚洲另类 | 一边摸一边抽搐一进一小说| 色综合欧美亚洲国产小说| 一级毛片高清免费大全| 国产三级中文精品| 日韩三级视频一区二区三区| 18禁裸乳无遮挡免费网站照片| 亚洲精品乱码久久久v下载方式 | 国产极品精品免费视频能看的| 国产免费av片在线观看野外av| 一个人观看的视频www高清免费观看 | 黄片大片在线免费观看| 久久热在线av| 欧美日韩综合久久久久久 | 亚洲乱码一区二区免费版| 国产一区二区在线av高清观看| 亚洲欧美日韩高清专用| 中文字幕av在线有码专区| 国产一区二区三区视频了| 久久久久久久久久黄片| 后天国语完整版免费观看| 久久午夜综合久久蜜桃| 亚洲精品一卡2卡三卡4卡5卡| 亚洲片人在线观看| 国产精品 国内视频| 狂野欧美激情性xxxx| 欧美高清成人免费视频www| 日本免费一区二区三区高清不卡| 国产三级黄色录像| 国产真实乱freesex| 一级毛片精品| 国产av在哪里看| or卡值多少钱| 色视频www国产| 午夜日韩欧美国产| 搡老岳熟女国产| 少妇裸体淫交视频免费看高清| 免费在线观看成人毛片| 一个人看视频在线观看www免费 | 欧美一区二区精品小视频在线| 色综合站精品国产| 夜夜夜夜夜久久久久| 国产成人精品久久二区二区免费| 少妇人妻一区二区三区视频| 岛国在线免费视频观看| 久久热在线av| 国产乱人视频| 久久中文字幕一级| 首页视频小说图片口味搜索| 国产精品久久久av美女十八| 亚洲aⅴ乱码一区二区在线播放| 91麻豆av在线| 九九在线视频观看精品| 天堂动漫精品| 又爽又黄无遮挡网站| 国产精品亚洲美女久久久| 日本 av在线| 在线十欧美十亚洲十日本专区| 悠悠久久av| АⅤ资源中文在线天堂| 99在线视频只有这里精品首页| 看黄色毛片网站| 亚洲专区中文字幕在线| 99国产精品一区二区蜜桃av| 久久久久久久久免费视频了| 久久精品人妻少妇| 日日干狠狠操夜夜爽| 午夜福利成人在线免费观看| 精品电影一区二区在线| 91久久精品国产一区二区成人 | 黄色片一级片一级黄色片| 亚洲自偷自拍图片 自拍| 成人三级黄色视频| 国产精品久久视频播放| 成人欧美大片| 美女免费视频网站| 久久久久久久精品吃奶| 观看免费一级毛片| 国产精品一区二区精品视频观看| 岛国视频午夜一区免费看| 午夜福利视频1000在线观看| 丁香六月欧美| 亚洲片人在线观看| 精品一区二区三区视频在线 | 精品福利观看| 九色成人免费人妻av| 狂野欧美白嫩少妇大欣赏| 9191精品国产免费久久| 黄片小视频在线播放| 99精品久久久久人妻精品| 国产1区2区3区精品| 亚洲熟女毛片儿| 国产91精品成人一区二区三区| 欧美成狂野欧美在线观看| 国产精品电影一区二区三区| 少妇丰满av| 一边摸一边抽搐一进一小说| 日日夜夜操网爽| 嫩草影视91久久| 亚洲av第一区精品v没综合| 欧美黄色淫秽网站| 老熟妇乱子伦视频在线观看| 九色成人免费人妻av| 国内精品久久久久久久电影| 国产又黄又爽又无遮挡在线| 欧美日韩国产亚洲二区| 麻豆成人av在线观看| 日韩精品青青久久久久久| 久久婷婷人人爽人人干人人爱| 日本五十路高清| 午夜日韩欧美国产| 成人永久免费在线观看视频| 国产精品一区二区三区四区免费观看 | 久久中文看片网| 国产精品日韩av在线免费观看| 国产伦人伦偷精品视频| 一个人观看的视频www高清免费观看 | 亚洲av第一区精品v没综合| 俺也久久电影网| 精品日产1卡2卡| 一级毛片高清免费大全| 久久国产精品影院| 亚洲国产精品999在线| 亚洲av熟女| 啦啦啦韩国在线观看视频| 美女午夜性视频免费| 中亚洲国语对白在线视频| 精品一区二区三区四区五区乱码| or卡值多少钱| 成人无遮挡网站| 狂野欧美白嫩少妇大欣赏| 国产精品久久电影中文字幕| 欧美日本视频| 麻豆久久精品国产亚洲av| 色综合亚洲欧美另类图片| 国产高清有码在线观看视频| 热99re8久久精品国产| 最近最新中文字幕大全电影3| 婷婷精品国产亚洲av| 色综合亚洲欧美另类图片| 俺也久久电影网| 99热这里只有是精品50| 久久精品91蜜桃| 精品国产三级普通话版| 18禁美女被吸乳视频| 欧美xxxx黑人xx丫x性爽| 狂野欧美白嫩少妇大欣赏| 丰满人妻熟妇乱又伦精品不卡| 天天一区二区日本电影三级| 日韩精品青青久久久久久| 日本 av在线| 又粗又爽又猛毛片免费看| 免费av毛片视频| 国产精品自产拍在线观看55亚洲| 免费看十八禁软件| 亚洲国产精品久久男人天堂| 免费一级毛片在线播放高清视频| 亚洲欧美精品综合一区二区三区| 国产成人系列免费观看| 少妇裸体淫交视频免费看高清| 久久草成人影院| 九九在线视频观看精品| 精品国产美女av久久久久小说| 欧美成人一区二区免费高清观看 | 国产亚洲精品av在线| 全区人妻精品视频| 亚洲乱码一区二区免费版| 99国产综合亚洲精品| 亚洲国产色片| 天堂av国产一区二区熟女人妻| 精品电影一区二区在线| 视频区欧美日本亚洲| bbb黄色大片| 亚洲国产欧美网| 午夜福利高清视频| 久久中文看片网| 俺也久久电影网| 久久午夜综合久久蜜桃| 麻豆成人av在线观看| 天堂动漫精品| 国产成人影院久久av| 色尼玛亚洲综合影院| 久久精品91蜜桃| 女生性感内裤真人,穿戴方法视频| 精品一区二区三区视频在线 | 国产97色在线日韩免费| 小蜜桃在线观看免费完整版高清| 精品一区二区三区视频在线观看免费| 亚洲片人在线观看| 观看美女的网站| 97超级碰碰碰精品色视频在线观看| 黄色片一级片一级黄色片| 观看免费一级毛片| 小说图片视频综合网站| 国产视频一区二区在线看| 大型黄色视频在线免费观看| 变态另类丝袜制服| www.www免费av| 欧美丝袜亚洲另类 | 黄色片一级片一级黄色片| 久久久久久大精品| 身体一侧抽搐| 亚洲人成网站在线播放欧美日韩| 国产精品久久久久久人妻精品电影| 不卡一级毛片| 宅男免费午夜| 久久亚洲真实| 18禁黄网站禁片午夜丰满| 黄色丝袜av网址大全| 亚洲最大成人中文| 国产精品久久久久久久电影 | 欧美日韩国产亚洲二区| 男人舔女人下体高潮全视频| а√天堂www在线а√下载| 国产成人一区二区三区免费视频网站| 在线播放国产精品三级| 久久久久免费精品人妻一区二区| 久久久久亚洲av毛片大全| 成人特级av手机在线观看| 成人无遮挡网站| 法律面前人人平等表现在哪些方面| xxx96com| 又紧又爽又黄一区二区| 国产精品九九99| 亚洲一区高清亚洲精品| 久久久久免费精品人妻一区二区| 亚洲成人中文字幕在线播放| 一个人观看的视频www高清免费观看 | 国产成人aa在线观看| 久久久久免费精品人妻一区二区| 久久亚洲精品不卡| 真人一进一出gif抽搐免费| 99国产精品99久久久久| 色在线成人网| 欧美日韩中文字幕国产精品一区二区三区| 久久精品国产99精品国产亚洲性色| 波多野结衣巨乳人妻| 国产乱人视频| 久久久国产精品麻豆| 99久国产av精品| 欧美三级亚洲精品| 宅男免费午夜| 性色avwww在线观看| 精品99又大又爽又粗少妇毛片 | 国内毛片毛片毛片毛片毛片| 狂野欧美白嫩少妇大欣赏| 成人欧美大片| 99国产极品粉嫩在线观看| 99在线人妻在线中文字幕| 别揉我奶头~嗯~啊~动态视频| 男女那种视频在线观看| 国语自产精品视频在线第100页| 五月伊人婷婷丁香| 亚洲av熟女| 亚洲国产中文字幕在线视频| 成人高潮视频无遮挡免费网站| 9191精品国产免费久久| av中文乱码字幕在线| 国内毛片毛片毛片毛片毛片| 女人被狂操c到高潮| 精品久久久久久久人妻蜜臀av| 亚洲一区二区三区色噜噜| 久久99热这里只有精品18| 嫁个100分男人电影在线观看| 精品一区二区三区av网在线观看| av在线蜜桃| 午夜福利免费观看在线| 18禁黄网站禁片午夜丰满| 日日摸夜夜添夜夜添小说| 色在线成人网| av视频在线观看入口| 一级毛片女人18水好多| 麻豆久久精品国产亚洲av| 久久久水蜜桃国产精品网| 麻豆国产av国片精品| 国产高清三级在线| 97碰自拍视频| 亚洲无线在线观看| 久久久国产成人精品二区| h日本视频在线播放| 亚洲av免费在线观看| 亚洲色图 男人天堂 中文字幕| 黑人操中国人逼视频| 美女大奶头视频| 国产高清有码在线观看视频| 精品国内亚洲2022精品成人| 国产高清视频在线播放一区| 久久精品aⅴ一区二区三区四区| 丰满人妻一区二区三区视频av | 国产成人一区二区三区免费视频网站| 欧洲精品卡2卡3卡4卡5卡区| 国产精品久久久av美女十八| 午夜福利免费观看在线| 熟女人妻精品中文字幕| 国产伦精品一区二区三区四那| ponron亚洲| 亚洲精品色激情综合| 一本久久中文字幕| 丰满人妻熟妇乱又伦精品不卡| 久久精品亚洲精品国产色婷小说| 亚洲国产精品999在线| 视频区欧美日本亚洲| 少妇的丰满在线观看| 又黄又粗又硬又大视频| 国产不卡一卡二| 欧美成人一区二区免费高清观看 | 亚洲一区二区三区色噜噜| 中文字幕人妻丝袜一区二区| 国产亚洲av嫩草精品影院| 国产欧美日韩精品一区二区| 啦啦啦免费观看视频1| 国产成人精品久久二区二区91| 亚洲精品456在线播放app | 一个人免费在线观看电影 | 午夜两性在线视频| 久久精品综合一区二区三区| 精品久久久久久成人av| 午夜福利在线在线| 国产精华一区二区三区| 日韩大尺度精品在线看网址| 18禁观看日本| 亚洲性夜色夜夜综合| 欧美3d第一页| 国产黄片美女视频| 亚洲,欧美精品.| 操出白浆在线播放| 色哟哟哟哟哟哟| 国产免费av片在线观看野外av| 欧美在线一区亚洲| 男女那种视频在线观看| 国产乱人视频| 久久中文字幕一级| www.熟女人妻精品国产| 悠悠久久av| 听说在线观看完整版免费高清| 成人亚洲精品av一区二区| 国产精品99久久久久久久久| 国产伦精品一区二区三区视频9 | 国产精品亚洲一级av第二区| 天堂√8在线中文| av片东京热男人的天堂| 亚洲在线自拍视频| 日本免费一区二区三区高清不卡| 日本免费a在线| 亚洲第一电影网av| 观看美女的网站| 日本三级黄在线观看| 亚洲国产精品合色在线| 床上黄色一级片| av黄色大香蕉| 在线视频色国产色| 麻豆成人午夜福利视频| 国产黄片美女视频| 国产精品精品国产色婷婷| 18禁黄网站禁片午夜丰满| 欧美日韩综合久久久久久 | 亚洲精品粉嫩美女一区| 欧美日韩亚洲国产一区二区在线观看| 人人妻人人澡欧美一区二区| 真人一进一出gif抽搐免费| 白带黄色成豆腐渣| 国产成年人精品一区二区| 一进一出抽搐动态| 免费看光身美女| 色综合亚洲欧美另类图片| 欧美成狂野欧美在线观看| 亚洲av中文字字幕乱码综合| av欧美777| 成人av在线播放网站| 亚洲国产精品合色在线| 国产三级在线视频| 亚洲av日韩精品久久久久久密| 校园春色视频在线观看| 黄色成人免费大全| 一二三四社区在线视频社区8| 少妇的丰满在线观看| 国产精品亚洲美女久久久| 色在线成人网| 中亚洲国语对白在线视频| 在线国产一区二区在线| 国产蜜桃级精品一区二区三区| 久久久色成人| 中文字幕高清在线视频| 亚洲美女视频黄频| 精品一区二区三区视频在线 | 欧美色欧美亚洲另类二区| h日本视频在线播放| 亚洲熟妇中文字幕五十中出| 国产淫片久久久久久久久 | 日韩欧美精品v在线| 国产真实乱freesex| 国产免费男女视频| 国产男靠女视频免费网站| 露出奶头的视频| 国产精品久久电影中文字幕| 国产精品综合久久久久久久免费| 成在线人永久免费视频| 精品国产三级普通话版| 18禁黄网站禁片午夜丰满| 变态另类成人亚洲欧美熟女| or卡值多少钱| 国产精品永久免费网站| 国产亚洲av高清不卡| 午夜福利免费观看在线| 亚洲国产精品久久男人天堂| 一区二区三区高清视频在线| 亚洲黑人精品在线| 久久国产精品人妻蜜桃| 2021天堂中文幕一二区在线观| 午夜福利在线观看吧| 日韩大尺度精品在线看网址| 午夜福利18| 男插女下体视频免费在线播放| 国产探花在线观看一区二区| 全区人妻精品视频| a级毛片a级免费在线| www.精华液| 小说图片视频综合网站| 十八禁网站免费在线| 欧美极品一区二区三区四区| ponron亚洲| cao死你这个sao货| 午夜日韩欧美国产| 999久久久精品免费观看国产| 性色av乱码一区二区三区2| 精品久久久久久成人av| av天堂在线播放| 亚洲国产欧美一区二区综合| 亚洲国产欧美网| 亚洲片人在线观看| 免费在线观看成人毛片| 在线观看免费视频日本深夜| 51午夜福利影视在线观看| 看片在线看免费视频| 国产成年人精品一区二区| 特大巨黑吊av在线直播| 一进一出抽搐gif免费好疼| 久久精品91无色码中文字幕| 亚洲无线观看免费| 色吧在线观看| 老司机福利观看| 国产三级在线视频| 亚洲人成电影免费在线| 老司机在亚洲福利影院| 岛国视频午夜一区免费看| 国产97色在线日韩免费| 动漫黄色视频在线观看| 久久久久久久精品吃奶| 91九色精品人成在线观看| 在线观看一区二区三区| 国产69精品久久久久777片 | 国产精品一区二区三区四区久久| 99久国产av精品| 亚洲 国产 在线| 国模一区二区三区四区视频 | 亚洲真实伦在线观看| 黄片大片在线免费观看| 韩国av一区二区三区四区| 国产亚洲欧美在线一区二区| 欧美大码av| 亚洲国产看品久久| 老司机午夜十八禁免费视频| 国内精品美女久久久久久| 日本a在线网址| 久久久久久久久免费视频了| 国产精品av视频在线免费观看| 亚洲欧美日韩高清专用| 在线观看免费视频日本深夜| 成人午夜高清在线视频| 精品国产乱子伦一区二区三区| 在线观看一区二区三区| 成年女人毛片免费观看观看9| 此物有八面人人有两片| 久9热在线精品视频| 成人国产综合亚洲| 18禁美女被吸乳视频| 国产探花在线观看一区二区| 精品乱码久久久久久99久播| 欧美丝袜亚洲另类 | 亚洲国产色片| 国产精品一区二区精品视频观看| 偷拍熟女少妇极品色| 丁香六月欧美| 长腿黑丝高跟| 99热这里只有精品一区 | 亚洲片人在线观看| 日日干狠狠操夜夜爽| 亚洲午夜理论影院| 久久久水蜜桃国产精品网| 国产v大片淫在线免费观看| 日韩免费av在线播放| 久9热在线精品视频| 搞女人的毛片| 黄色视频,在线免费观看| 亚洲av成人一区二区三| 国产欧美日韩一区二区三| 99久久成人亚洲精品观看| 亚洲国产看品久久| 国产日本99.免费观看| 丁香六月欧美| 久久香蕉国产精品| 桃色一区二区三区在线观看| 欧美在线黄色| xxx96com| 非洲黑人性xxxx精品又粗又长| 99久久久亚洲精品蜜臀av| 久9热在线精品视频| 国产男靠女视频免费网站| 99久国产av精品| 国产精品久久久人人做人人爽| 夜夜躁狠狠躁天天躁| 天天添夜夜摸| 757午夜福利合集在线观看| 禁无遮挡网站| 此物有八面人人有两片|