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

熱彈耦合梁結構動力響應的區間數值分析

2016-07-26 02:21:48云永琥陳建軍曹鴻鈞
振動與沖擊 2016年1期

云永琥, 陳建軍, 曹鴻鈞

(西安電子科技大學 電子裝備結構設計教育部重點實驗室,西安 710071)

?

熱彈耦合梁結構動力響應的區間數值分析

云永琥, 陳建軍, 曹鴻鈞

(西安電子科技大學 電子裝備結構設計教育部重點實驗室,西安710071)

摘要:研究了含有區間參數梁結構在溫度載荷和力載荷共同作用下的動力響應問題,考慮材料變形與傳熱的相互影響,建立了梁在熱彈耦合下的動力學有限元模型,并給出對結構瞬態熱傳導方程與動力學方程進行相互交替迭代求解的計算方法。針對結構響應不確定性問題,以不確定參數作為約束變量,通過尋求結構響應函數的區間范圍,將區間問題轉化為優化問題,并利用遺傳算法給出了結構響應函數的區間界限。通過算例結果與概率有限元方法的計算結果比較,表明文中所提出方法的可行性和有效性,并獲得在熱彈耦合作用下梁結構的固有振動頻率有所增加,而振動響應振幅則逐漸減弱的結論。該方法只需已知不確定參數所在范圍的界限,而無需其他統計信息,為解決區間參數熱彈耦合梁問題提供了一種途徑。

關鍵詞:熱彈耦合;區間有限元;全局優化;遺傳算法;動力響應;固有頻率

在航天工程中,各種航天器廣泛采用可展開的大型外懸式柔性結構,如太陽能帆板、天線和裝載各種有效載荷的桅桿等,它們在太空中主要承受隨著航天器在軌進出陰影和日照區時的瞬變熱載荷,其中柔性結構因瞬態溫度梯度的巨變將有可能誘發熱致振動[1]。這種結構變形導致熱傳導邊界條件發生變化,影響了結構溫度變化,反過來由于高溫和溫度梯度而產生的熱應力,則改變了結構的剛度分布,從而影響了結構的動態特性。這種動態過程中溫度場和應力場的相互耦合作用,即為熱彈耦合效應。

Zhao等[2]研究了考慮在輻射下層疊板的熱誘發振動問題。Boly[3]分析了梁在突加熱載荷作用下發生振動的問題,提出了耦合系數的概念,并對梁結構進行了耦合分析。安翔等[4]對大型空間結構的熱彈耦合動力學問題進行了分析,提出了在全區間內考慮耦合而在單個時間步內考慮非耦合的有限元計算方法。Sun等[5]研究了微尺度圓板的熱彈耦合問題,得出了熱阻尼與環境溫度、板的尺寸以及邊界條件的變化關系。李智勇等[6-8]關于熱載荷對結構振動的影響進行了相應的研究。

現有的結構熱彈耦合研究主要集中于確定性的分析模型,然而在許多實際工程結構中,由于加工、裝配和建模中可能存在誤差或不確定因素,使得結構的物理參數、幾何參數、初始和邊界條件也具有不確定性,雖然大多數情況下這些誤差或不確定性可能很小,但由多種誤差或不確定性的共同作用則有可能使結構響應產生較大的偏差,尤其是在復雜的機構系統中[9]。因此研究材料參數、尺寸參數,邊界條件和初始條件為不確定性的熱彈耦合動力學問題,既具有理論意義又具有很現實的工程背景。

工程結構中不確定性問題可采用隨機方法、模糊方法和區間方法求解,但前兩種方法建立的結構分析模型,需提供不確定性結構參數的概率密度函數或隸屬度函數,其計算量一般也比較大。如文獻[10]將不確定性結構參數視為隨機變量,利用隨機因子法推導出結構動力響應數字特征的計算表達式,對隨機參數彈性桿在瞬態溫度場下的結構響應問題進行了分析。但在實際的工程結構分析過程中,有時很難獲得較多的統計數據來描述不確定性參數的概率分布或隸屬函數,此時區間分析方法較之前兩者顯得更加有效。該方法將不確定性的結構參數視為未知變量,使它們在具有已知邊界的區間內取值,這是一種只需要通過較少數據信息就可以描述結構不確定性的分析方法。文獻[11-12]利用區間分析方法對區間不確定性結構的動態特性及動力響應進行了分析。

本文針對Euler-Bernouli梁結構,考慮其物性參數、溫度和外力載荷均為區間變量,建立了其熱彈耦合動力學區間有限元模型,將智能算法引入到區間結構分析之中,并將區間有限元模型的計算轉換成為區間參數結構的全局優化問題來求解,以避免直接采用區間運算法則而使得計算結果易于擴張[13]。通過算例表明了該方法的可行性與有效性,并考察了耦合效應對結構動力學響應的影響。

1熱彈耦合簡支梁有限元方程的建立與求解

以圖1矩形截面簡支梁為分析對象,初始環境溫度為T0,在t=0的初始時刻,梁的上表面同時受到階躍均布力f和熱流q的共同作用,上下表面與外界進行對流換熱,忽略左右兩側面以及梁前后表面的換熱。由于熱流在梁上表面均勻分布,并沿著梁軸向任意截面上的溫度處處相等,因此熱流僅沿梁厚度方向傳導。求解簡支梁熱彈耦合動力響應需要分別構建以下兩種有限元模型。

圖1 簡支梁結構圖Fig.1 Simply supported beam structure diagram

1.1梁動力學有限元模型

將梁沿軸向離散為m個單元、m+1個節點,每一節點具有橫向位移v和截面轉角φ兩個自由度。利用兩節點Hermite插值函數并根據哈密頓原理推導出梁單元的動力學有限元方程,并通過對各個單元矩陣進行組集,忽略阻尼效應,可得梁的動力學有限元方程為[14]:

(1)

K=Ks+Kσ

(2)

式中:u為結構的位移向量,M為總質量矩陣,K為總剛度矩陣。FB為總力載荷向量,FT為總溫度載荷向量。以上矩陣均由文獻[15]可得。其中梁結構在熱環境下總剛度矩陣K由結構自身的彈性剛度矩陣Ks=∫VBTDBdV和熱應力剛度矩陣Kσ=∫VGTΓGdV疊加組成,式中B為幾何矩陣,D為彈性矩陣,G為形函數對空間坐標微分矩陣,Γ為應力矩陣。這是因為需要考慮熱環境引起的熱應力對結構振動特性的影響,在總剛度矩陣中增加一項熱應力矩陣Kσ。

1.2梁熱分析有限元模型

梁熱分析模型如圖2所示,將梁橫截面沿厚度方向離散分為m層、m+1個節點。由自由能密度和熵密度的計算,并利用最小勢能原理,導出一維耦合熱傳導有限元方程為[6]:

(3)

其中各個矩陣和向量的表達式如下:

式中:T為待求的節點瞬態溫度向量;C為熱容矩陣;Kk為熱傳導剛度矩陣;Kh為對流換熱矩陣;P為節點熱載荷列陣;N為單元節點溫度的形函數;H為溫度與結構變形的耦合矩陣項,它表明待定溫度場不僅與熱源、熱力學物性參數及換熱邊界條件有關,還受到結構變形應變率的影響,其在一定程度上改變物體內部的熱量傳遞;k為熱傳導系數,h為換熱系數,ρ為質量密度,c為比熱容,μ為泊松比,q為熱流量,E為彈性模量,α為熱膨脹系數。

圖2 熱分析模型Fig.2 Thermal analysis model

1.3熱彈耦合梁結構固有頻率分析

對于在熱環境下考慮熱彈耦合效應時梁結構固有頻率計算即為求解式(4)的廣義特征值問題:

(K-ω2M)Φ=0

(4)

式中:ω和Φ分別為結構的固有頻率和振型向量。

在整個分析過程中,可認為質量矩陣M不受溫度場的影響,由于結構溫度的不均勻變化以及結構不能自由的熱變形而使體內產生了熱應力,由式(2)可知,結構熱應力矩陣Kσ的變化,造成總體剛度矩陣K的改變,最終在式(4)中將導致結構固有頻率的變化。為此在求解熱環境下梁的固有頻率時,需綜合考慮熱應力對剛度矩陣的影響。

1.4熱彈耦合梁結構動力響應求解

由于結構溫度場的變化會影響結構的變形,反之結構變形也會影響溫度場的分布。因此在計算初始邊界時,先計算結構瞬態溫度場分布情況,然后將溫度值作為載荷施加到結構的動力分析當中,從而得到結構在受到力和溫度載荷共同作用下的響應。對于本文熱彈耦合動力學有限元方程的求解,首先對瞬態熱傳導方程式(3)進行求解,再對動力學方程式(1)進行求解,將兩者聯立并相互交替迭代計算。

對式(3),采用Galerkin迭代格式進行求解,取θ=2/3的無條件穩定迭代格式[16]如下:

(2KT+3C/Δt)Tt+Δt=

(5)

而式(1)則利用Newmark法將其轉化為求解如下擬靜力方程:

(6)

式中:

(7)

(8)

2區間參數熱彈耦合響應求解全局優化模型

考慮到結構制造誤差和外界環境等多種不確定性因素的影響,本文將結構的物性參數:ρ、c、k、h、α、μ和E等均視為區間參數,同時,將結構所受外力f、溫度載荷q以及結構初始溫度T0亦視為區間參數,將它們統一以區間參數向量形式表為β=(β1,β2,…,βm)T,其所在的范圍為:

(9)

將式(9)代入式(4)、(1)和(2)中,則得到區間參數梁結構的模態、動力學以及瞬態溫度場方程:

(K(β)-ω2M(β))Φ=0

(10)

FB(β)+FT(β)

(11)

(12)

由式(10)~(12)可知,物性參數和載荷的區間性將會引起梁結構固有頻率、動力響應以及瞬態溫度場的區間變化,有效的估計梁結構的動力學響應以及溫度場響應的區間界值具有十分重要的意義。

目前求解類似于式(10)~(12)的區間有限元方法主要有攝動法、單調性法和端點組合法等。這些方法的共同特點是需要做一些前提假設,例如區間攝動法假設區間參數具有很小的變化范圍,單調性法和端點組合法則是基于位移響應函數為區間變量的單調函數的假設,且對于區間參數較多情況下計算量較大,理論推導及編程實施復雜,故局限于低維簡單問題的分析。然而在實際的工程問題中,由于區間參數數目多、參數區間離散程度大,以及參數之間相關性的影響,使得上述方法計算的結果區間范圍與實際相比過于保守。

對于含有區間參數的結構響應結果,工程中最為關心的是其最大值與最小值,如在結構的可靠性設計中往往需要考慮其響應的極值情況。故在本文中,直接以區間參數和區間載荷作為設計變量,以結構固有頻率ω(β)、位移響應u(β,t)和瞬態溫度場T(β,t)作為目標函數,利用全局優化方法求解結構響應極值。這不僅可以解決此類問題,還可得到結構響應范圍的準確估計值[13]。對以上目標函數的區間范圍可通過求解如下約束優化模型獲得:

max/minω(β)

s.t.[K(β)-ω2M(β)]Φ=0

β∈βI

(13)

max/minu(β,t)

K(β)u(β,t)=FB(β)+FT(β)

β∈βI

(14)

max/minT(β,t)

(15)

β∈βI

以下將給出求解區間參數梁結構熱彈耦合動力方程響應上下界的計算方法。

3考慮區間不確定性的優化算法

本文采用遺傳算法對優化模型式(13)、(14)和(15)進行求解。

熱彈耦合動力響應問題與時間有關,對目標函數是在整個時間歷程上進行求解,需要在時間域上離散。當通過遺傳算法計算得到結構的動力響應在當前時間迭代步的最大與最小值,之后在下一個時間迭代步內,重復利用遺傳算法進行計算,進而求出在整個時間域上的動力響應曲線的最大和最小包絡線(即為目標函數的區間范圍)。具體求解步驟如下:

步驟1確定結構中區間變量個數和區間有限元方程參數的區間變化范圍,給定遺傳算法參數,并在給定的區間內創建初始種群;設置遺傳最大進化代數M;

步驟2在當前時間迭代步內,根據給定的參數值,對熱彈耦合動力學有限元方程進行求解得到結構的響應值,從而計算出群體中各個個體的適應度;

步驟3對終止條件進行判斷,若不滿足循環停止條件,則M=M+1,轉到步驟(4);否則將當前迭代步進化過程中所得到的最大適應度個體作為最優解輸出,然后轉到步驟(2)繼續在下一時間迭代步上進行上述遺傳進化計算;

步驟4重復進行遺傳進化操作,經過選擇、交叉、變異運算之后得到下一代種群,并計算新種群個體適應度,轉到步驟(3);

步驟5當在給定時間域上逐步完成計算之后,則得到整個時間域上動力響應的最大和最小包絡線,計算結束。

4算例

求圖1所示矩形截面簡支梁在溫度和力載荷共同作用下動力學響應的區間。梁長L=600 mm,寬b=100 mm,高a=50 mm,材料為鋁,結構的區間參數見表1。在遺傳計算中,取種群規模為50,交叉概率為0.95,變異概率為0.05,最大進化代數為100。

求解該梁熱彈耦合動力學響應需分別采用兩種有限元模型。圖1為梁動力學分析模型,沿其長度方向被離散為4個單元、5個節點,利用結構動力學方程計算得出節點的位移響應。圖2為梁熱分析模型,沿梁橫截面的厚度方向亦被離散為4個單元、5個節點,利用熱傳導有限元方程可求得截面溫度場,從而算出梁截面溫度載荷向量。

表1 梁結構參數取值區間

表2 熱流密度為1.63×106/(W·m-2 )時各階固有頻率比較

由表2可知,在考慮熱彈耦合效應時,梁結構內產生了附加熱應力剛度矩陣,導致總體剛度的增加,因此結構的前三階固有頻率要比不考慮熱彈耦合效應的固有頻率大,各階固有頻率增大的比例幅度基本相同。

圖3給出了梁的前3階固有頻率隨熱流載荷的變化規律。由圖可知,隨著外加熱流密度的增大,梁結構不均勻膨脹而產生的熱應力導致了結構的剛度變化,從而使得梁前三階固有頻率均單調升高。

圖4為梁上下表面溫差均值分別在耦合與非耦合(即耦合項為0時,結構變形對溫度場未產生影響)情況下的時間歷程計算結果。由于熱彈耦合與結構的變形應變率相關,在未考慮耦合效應時,梁上下表面溫差在初始階段升高很快,隨著時間逐漸收斂,趨于恒定值。而在耦合效應的作用下,梁上下表面溫差的變化與非耦合規律一致,但是上下表面溫差則隨著位移響應產生相應的小幅度波動。

圖3 梁固有頻率隨熱載變化曲線Fig.3 The natural frequency of beam with the heat load

圖4 簡支梁上下表面溫差均值時間歷程Fig.4 Time history of the midpoint temperature at the simply supportedbeam mid-span

圖5為梁中點位移響應均值隨時間變化歷程的計算結果。當梁上表面受熱載荷作用時,較高的溫度使得線性膨脹量增大,但均布力載荷也同樣施加在梁的上表面,限制了梁在y方向的熱膨脹,由于均布力載荷的影響大于熱載荷的影響,因此中點的位移則沿作用力的方向。從圖中可以看出,在力載荷的沖擊作用下,梁中點迅速向下運動達到最大值,并圍繞該位置產生微幅振動。但由于熱彈耦合效應,施加給梁結構的機械能將轉化為熱流并耗散,從而抑制了結構向下運動的幅度,且還使得位移振幅隨時間逐漸減弱。若不考慮熱彈耦合效應,位移響應在各個周期的振動狀態都相同。這與文獻[6,17-18]中的結論一致。

圖5 簡支梁中點位移均值的時間歷程Fig.5 Time history of displacement at the midpoint of the simply supported beam

圖6為通過遺傳算法獲得的區間參數梁結構在熱彈耦合下的區間范圍,并與區間參數在其變化范圍內呈均勻分布下的隨機分析方法(Monte Carlo方法)以及確定性方法計算的結果進行了對比。從圖中可知,相對于確定性的分析方法以及Monte Carlo方法的分析結果,采用遺傳算法得到的位移響應結果是一個區間范圍,其變化規律與前兩者計算結果一致,由于區間分析依賴不確定參數信息較少,使得分析結果稍大于Monte Carlo法的分析結果范圍,然而區間分析方法的計算代價要明顯小于Monte Carlo法的計算代價。

圖6 簡支梁中點位移響應區間的時間歷程Fig.6 Time history of interval response of displacement at the midpoint of the simply supported beam

5結論

本文研究了區間參數梁結構在熱彈耦合下動力響應的區間分析方法,主要結論如下:

(1) 文中所提出的方法可以避開由于直接采用區間運算使得計算結果易于擴張的缺陷,且無需引入更多的假設,因此該方法是可靠的。通過數值算例,表明該方法是可行的。與確定性的計算結果相比,則說明了本文方法能夠較準確地求解區間參數梁結構同時受到熱載荷和力載荷作用時的動力響應區間。

(2) 由于熱彈耦合效應,梁結構的固有頻率相比非耦合狀態下的頻率有所增加,梁的固有頻率隨著外加熱載的增大而升高。

(3) 在突加熱載的情況下,梁的上下表面溫差快速增加,但到一定的時間后溫差趨于平穩。由于熱彈耦合效應與位移變化相關,使得溫差產生了波動。

(4) 熱彈耦合效應對結構位移振動有抑制作用,使結構的振幅隨時間不斷減小,并趨于穩態的平衡位置。

參 考 文 獻

[1] Thornton E A, Kim Y A. Thermally induced bending vibrations of a flexible rolled-up solar array[J]. Journal of Spacecraft and Rockets, 1993, 30(4): 438-448.

[2] Zhao S G, Wang J T, Li K, et al. Thermally induced vibration analysis of laminated plate considering radiation by finite element method[J]. Journal of Mechanics, 2011,27(4): 33-37.

[3] Boly B A. Thermally induced vibrations of beams[J]. Journal of the Aeroautical Science, 1956, 23(2): 179-181.

[4] 安翔, 馮剛. 某空間站太陽電池陣中央桁架熱-結構耦合動力學分析[J]. 強度與環境, 2005, 32(3): 8-13.

AN Xiang, FENG Gang. Thermally induced vibration of the main mast of the space station’s solar arrays[J]. Structure & Environment Engineering, 2005, 32(3): 8-13.

[5] Sun Yu-xin, Saka M. Thermoelastic damping in micro-scale circular plate resonators[J]. Journal of Sound and Vibration, 2010, 329(3): 328-337.

[6] 李智勇, 劉錦陽, 洪嘉振. 作平面運動的二維平面板的熱耦合動力學問題[J]. 動力學與控制學報, 2006, 4(2): 114-121.

LI Zhi-yong, LIU Jin-yang, HONG Jia-zhen. Coupled thermoelastic dynamics of a two-dimensional plate undergoing planar motion[J]. Journal of Dynamics and Control, 2006,4(2): 114-121.

[7] 樹學鋒, 蘭姣霞, 武勇忠. 大撓度圓柱殼在溫度場中的熱彈耦合振動分析[J]. 應用數學和力學, 2004, 25(9): 910-916.

SHU Xue-feng, LAN Jiao-xia, WU Yong-zhong. Analysis of thermal elastic coupling vibration of large deflection cylindrical shell[J]. Applied Mathematics and Mechanics, 2004, 25(9): 910-916.

[8] Nowruzpour Mehrian S M, Naei M H, Mehrian S Z. Dynamic response for a functionally graded rectangular plate subjected to thermal shock based on LS theory[J]. Applied Mechanics and Materials, 2013, 332(1): 381-395.

[9] 王光遠. 論不確定性結構力學的發展[J]. 力學進展, 2002, 32(2): 205-211.

WANG Guang-yuan. On the development of uncertain structural mechanics[J]. Advance in Mechanics, 2002,32(2):205-211.

[10] 閻彬, 陳建軍, 馬洪波. 隨機彈性桿在隨機瞬態溫度場下的熱響應分析[J]. 電子科技大學學報, 2012, 41(4): 631-636.

YAN Bin, CHEN Jian-jun, MA Hong-bo. Thermal response analysis of stochastic pole structures under random transient temperature field[J]. Journal of University of Electronic Science and Technology of China, 2012, 41(4): 631-636.

[11] 王敏娟, 陳建軍, 林立廣, 等. 區間參數智能梁結構閉環系統動態特性分析[J]. 電子科技大學學報, 2011, 40(1): 152-156.

WANG Min-juan, CHEN Jian-jun, LIN Li-guang, et al. Dynamic characteristic analysis of closed loop systemsfor the intelligent beam with interval parameters[J]. Journal of University of Electronic Science and Technology of China, 2011, 40(1): 152-156.

[12] 李金平, 陳建軍, 朱增青,等. 結構區間有限元方程組的一種解法[J]. 工程力學, 2010, 27(4): 79-83.

LI Jin-ping, CHEN Jian-jun, ZHU Zeng-qing, et al. A method for solving the structural interval finite element equations[J]. Engineering Mechanics, 2010, 27(4): 79-83.

[13] 王登剛, 李杰. 計算具有區間參數結構特征值范圍的一種新方法[J]. 計算力學學報, 2004, 21(1): 56-61.

WANG Deng-gang, LI Jie. A new method for computing the eigenvalues bounds of structures with interval parameters[J]. Chinese Journal of Computational Mechanics, 2004, 21(1): 56-61.

[14] Narasimha M, Appu Kuttan K K, Ravikiran K. Thermally induced vibration of a simply supported beam using finite element method[J]. International Journal of Engineering Science, 2010, 2(12): 7874-7879.

[15] 范緒箕. 高速飛行器熱結構分析與應用[M]. 北京: 國防工業出版社, 2008.

[16] 王勖成. 有限單元法[M]. 北京: 清華大學出版社, 2003.

[17] 尹益輝, 郝志明, 陳裕澤,等. 不同材料參數薄板振動中的熱力耦合效應[J]. 強激光與粒子束, 2001, 13(2): 142-146.

YIN Yi-hui, HAO Zhi-ming, CHEN Yu-ze, et al. Thermo-mechanical coupling effects in vibrations of plates with different properties[J]. High Power Laser and Particle Beams, 2001, 13(2): 142-146.

[18] 王琪, 吉庭武, 謝公南,等. 輕質熱防護系統波紋夾芯結構熱力耦合分析[J]. 應用數學和力學, 2013, 34(2): 172-182.

WANG Qi, JI Ting-wu, XIE Gong-nan, et al.Structural analysis of corrugated-core sandwich panels for lightweight thermal protection system[J]. Applied Mathematics and Mechanics, 2013, 34(2): 172-182.

基金項目:國家自然科學基金項目(51175398)資助

收稿日期:2014-12-01修改稿收到日期:2015-03-31

通信作者陳建軍 男,碩士,教授,博士生導師,1951年生

中圖分類號:O324

文獻標志碼:A

DOI:10.13465/j.cnki.jvs.2016.01.034

Interval numerical analysis for dynamic response of a thermoelastic coupled beam structure

YUN Yong-hu, CHEN Jian-jun, CAO Hong-jun

(MOE Key Laboratory of Electronic Equipment Structure Design, Xidian University, Xi’an 710071, China)

Abstract:Dynamic response of a thermoelastic coupled beam structure with interval parameters was studied under both thermal load and force load. Considering interaction of material deformation and heat conduction, the dynamic model of the beam structure was set up using the finite element method. The calculation method was proposed for solving the transient temperature field and dynamical response with iterative technqique. For structural response uncertainty, taking uncertain parameters as a constraint variables,the interval bounds of structural response function were solved through solving the corresponding optimization problems, and the genetic algorithm was used to solve the global optimization model. Compared with the probabilistic analysis method,the numerical example indicated the feasibility and validity of the proposed method. It was shown that the natural frequencies of the beam increase and the amplitudes of its vibration gradually decrease due to thermoelastic coupled effect. The proposed method only needed to know the limits of the range of the uncertain parameters, without needing other statistical information, and provided a way to solve such complex problems of thermoelastic coupled beam structures with uncertainty.

Key words:thermoelastic coupling; interval finite element method; global optimization; genetic algorithm; dynamic response; natural frequency

第一作者 云永琥 男,博士生,1986年1月生

主站蜘蛛池模板: 女高中生自慰污污网站| 一级毛片在线播放免费观看| 色男人的天堂久久综合| 免费高清a毛片| 亚洲精选高清无码| lhav亚洲精品| 国产69囗曝护士吞精在线视频| 欧美日韩国产成人高清视频| 午夜激情福利视频| 亚洲香蕉久久| 婷婷六月激情综合一区| 好紧好深好大乳无码中文字幕| 综合色区亚洲熟妇在线| 日本少妇又色又爽又高潮| 国产人成网线在线播放va| 69免费在线视频| 免费国产不卡午夜福在线观看| 试看120秒男女啪啪免费| 无码又爽又刺激的高潮视频| 国产新AV天堂| 野花国产精品入口| 国产黄在线免费观看| 999国产精品永久免费视频精品久久| 2020最新国产精品视频| 四虎精品黑人视频| 毛片网站在线看| 久久频这里精品99香蕉久网址| 日韩毛片免费| 国产欧美专区在线观看| 五月婷婷精品| 国产高清自拍视频| 成人福利免费在线观看| 午夜丁香婷婷| 五月天综合婷婷| 99国产精品国产高清一区二区| 天天做天天爱夜夜爽毛片毛片| 一本大道无码日韩精品影视| 欧美午夜视频| 亚洲性影院| 亚洲国产一区在线观看| 亚洲精品国产精品乱码不卞 | 国产精品嫩草影院av| 国产中文在线亚洲精品官网| 亚洲国产综合精品中文第一| 呦系列视频一区二区三区| 国产一区二区三区免费观看| 凹凸国产熟女精品视频| 97精品久久久大香线焦| 一级毛片免费的| 久久久成年黄色视频| 免费在线观看av| 国产精品三级专区| 亚洲国产成人超福利久久精品| 毛片网站在线看| 亚洲最大看欧美片网站地址| 19国产精品麻豆免费观看| 乱人伦视频中文字幕在线| 五月婷婷伊人网| 亚洲精品自拍区在线观看| 动漫精品中文字幕无码| 亚洲一区二区三区麻豆| h视频在线观看网站| 亚洲国产清纯| 伊人国产无码高清视频| 亚洲视频a| 国产午夜小视频| 在线亚洲精品福利网址导航| 亚洲精品国产日韩无码AV永久免费网| 亚洲全网成人资源在线观看| 国产丝袜无码精品| 国产精品久久自在自线观看| 无码一区18禁| 园内精品自拍视频在线播放| 99热这里只有成人精品国产| 91午夜福利在线观看精品| 国产成人麻豆精品| 欧美a在线看| 国产永久免费视频m3u8| 国产香蕉国产精品偷在线观看| 人妻无码AⅤ中文字| 午夜无码一区二区三区| 亚洲成a人片77777在线播放|