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

計(jì)及熱效應(yīng)的小展弦比機(jī)翼氣彈穩(wěn)定性分析*

2014-02-19 04:18:24韓曉林費(fèi)慶國(guó)
振動(dòng)、測(cè)試與診斷 2014年3期
關(guān)鍵詞:模態(tài)結(jié)構(gòu)分析

趙 衛(wèi), 韓曉林, 費(fèi)慶國(guó)

(1.東南大學(xué)土木工程學(xué)院 南京,210096) (2.東南大學(xué)江蘇省工程力學(xué)分析重點(diǎn)實(shí)驗(yàn)室 南京,210096)

引 言

隨著高/超聲速飛行器的發(fā)展,高/超聲速氣流引起氣動(dòng)加熱產(chǎn)生的熱問(wèn)題越來(lái)越受到重視。結(jié)構(gòu)在熱環(huán)境下的動(dòng)力學(xué)問(wèn)題是發(fā)展高/超聲速航天飛行器面臨的重要課題。20世紀(jì)50年代末至60年代初,美國(guó)對(duì)熱氣動(dòng)彈性問(wèn)題的研究出現(xiàn)了一個(gè)高潮,這些研究為之后航天飛機(jī)的設(shè)計(jì)打下了堅(jiān)實(shí)的基礎(chǔ)[1]。Garrick[2]綜述了熱氣動(dòng)彈性力學(xué)的發(fā)展,在傳統(tǒng)氣彈力三角形的基礎(chǔ)上補(bǔ)充了由熱效應(yīng)引起的作用力,提出了一種新的描述熱氣動(dòng)彈性力學(xué)各學(xué)科關(guān)系的四面體關(guān)系圖。隨著計(jì)算機(jī)運(yùn)行能力的提高和計(jì)算流體動(dòng)力學(xué)(computational fluid dynamics,簡(jiǎn)稱CFD)的發(fā)展,氣動(dòng)加熱的計(jì)算精度和效率得到很大提高,為結(jié)構(gòu)在熱環(huán)境中的動(dòng)力特性分析提供了準(zhǔn)確的“溫度載荷”。在結(jié)構(gòu)和氣動(dòng)力耦合方面,Krist等[3]基于RANS方程開(kāi)發(fā)了CFL3D求解器,為基于CFD的高/超聲速飛行器的結(jié)構(gòu)氣彈穩(wěn)定性分析提供了有效工具。由于CFL3D能夠準(zhǔn)確模擬非定常氣動(dòng)力,McNamara等[1,4]基于 CFL3D分析了超聲速氣流中氣動(dòng)彈性穩(wěn)定性,并取得了一定的進(jìn)展。20世紀(jì)90年代初,美國(guó)開(kāi)展了NASP(25倍馬赫數(shù))、TAV,Hy-Tech,Hyper-X等飛行器研究,為熱氣彈分析提供了大量的實(shí)驗(yàn)數(shù)據(jù)[5-6]。Heeg等[7]針對(duì) NASP 驗(yàn)證模型提出了熱氣動(dòng)彈性分析的5個(gè)步驟,但其并沒(méi)有考慮熱應(yīng)力引起的附加剛度。針對(duì)熱效應(yīng)對(duì)結(jié)構(gòu)振動(dòng)特性的影響,史曉明等[8-9]研究了變厚度彈翼的熱動(dòng)力學(xué)特性,對(duì)比有限元分析和實(shí)驗(yàn)結(jié)果,發(fā)現(xiàn)結(jié)構(gòu)前3階固有頻率較常溫均有所下降。在熱氣動(dòng)彈性數(shù)值計(jì)算方面,吳志剛等[10-11]提出了分層求解的思想,采用分開(kāi)建立方程的松耦合解決方案,即分別進(jìn)行氣動(dòng)加熱計(jì)算、結(jié)構(gòu)溫度場(chǎng)分析、熱模態(tài)分析、氣動(dòng)彈性分析,在一個(gè)時(shí)間步上順序求解,從而使問(wèn)題簡(jiǎn)化。基于此分層求解思想,文獻(xiàn)[10]分析了高/超聲速穩(wěn)態(tài)熱環(huán)境下翼面的結(jié)構(gòu)、氣動(dòng)和熱的耦合動(dòng)力學(xué)問(wèn)題,采用單向耦合方法進(jìn)行計(jì)算,發(fā)現(xiàn)受熱結(jié)構(gòu)的動(dòng)力特性和顫振特性均可能發(fā)生變化,尤其是對(duì)于根部固支翼面。

筆者以小展弦比機(jī)翼作為研究對(duì)象,根部完全固支,探索其在高/超聲速氣流的瞬態(tài)熱環(huán)境下的動(dòng)力學(xué)特性和氣彈穩(wěn)定性,考慮了熱應(yīng)力引起的剛度效應(yīng),研究熱應(yīng)力對(duì)結(jié)構(gòu)剛度的影響,利用活塞理論計(jì)算非定常氣動(dòng)力,采用p-k法求解氣彈方程[12-15],分析討論熱效應(yīng)對(duì)結(jié)構(gòu)振動(dòng)特性及顫振邊界的影響,同時(shí)還研究了不同邊界約束條件下熱應(yīng)力對(duì)結(jié)構(gòu)振動(dòng)特性和氣彈穩(wěn)定性的影響。

1 有限元模型

小展弦比機(jī)翼根部弦長(zhǎng)和尖部弦長(zhǎng)分別為1.0m和0.3m,展長(zhǎng)為0.4m,前緣后掠角為33.7°,機(jī)翼根部完全固支。有限元模型[16]如圖1所示。

材料為超硬鋁合金7075-Al,該型號(hào)鋁合金材主要用于制造飛機(jī)結(jié)構(gòu)及要求強(qiáng)度高、抗腐蝕性強(qiáng)的高應(yīng)力構(gòu)件。假定材料密度不隨溫度變化,ρ=2 750kg/m3,泊松比μ=0.3。其他材料結(jié)構(gòu)參數(shù)和熱特性參數(shù)隨溫度變化如表1所示。

表1 不同溫度下7075-Al材料參數(shù)Tab.1 Material properties of 7051-Al at differnet temperature

圖1 小展弦比機(jī)翼有限元模型示意圖Fig.1 Finit element model of low aspectration wing

2 瞬態(tài)溫度場(chǎng)分析

2.1 溫度場(chǎng)分布方程

三維瞬態(tài)溫度場(chǎng)在直角坐標(biāo)系下滿足微分方程

當(dāng)一個(gè)方向上(若為z方向)溫度變化為零時(shí),方程就為二維問(wèn)題的熱傳導(dǎo)方程

若只有z方向有溫度變化,則方程為一維熱傳導(dǎo)方程

其中:kx,ky,kz分別為材料沿x,y,z方向的熱傳導(dǎo)系數(shù);φ為瞬態(tài)溫度場(chǎng)的場(chǎng)變量φ(x,y,z,t);ρ為大氣密度;c為比熱;Q為內(nèi)部熱源。

筆者假設(shè)機(jī)翼沿厚度方向的溫度分布變化為0,即z向無(wú)溫度梯度,同時(shí)機(jī)翼內(nèi)部不產(chǎn)生內(nèi)源,則式(2)可寫(xiě)成

2.2 瞬態(tài)溫度場(chǎng)分析

利用式(4)計(jì)算翼面溫度分布。本研究中“熱源”為熱流密度,假設(shè)機(jī)翼前緣的熱流密度為1.0×105W/m2,考慮空氣對(duì)流的影響,對(duì)流傳熱系數(shù)h=10.0W/m2°C,環(huán)境溫度為20°C,選取50,500,1 000,2 000和3 000s時(shí)段的溫度視為“溫度載荷”進(jìn)行結(jié)構(gòu)模態(tài)分析和顫振分析,分析不同時(shí)刻下的頻率和顫振邊界,并探討其變化規(guī)律。

3 熱模態(tài)分析

熱效應(yīng)對(duì)結(jié)構(gòu)振動(dòng)特性的影響主要是改變結(jié)構(gòu)總體剛度。溫度對(duì)結(jié)構(gòu)剛度的影響主要表現(xiàn)在兩個(gè)方面:a.高溫下結(jié)構(gòu)材料性能發(fā)生退化,直接降低了結(jié)構(gòu)剛度;b.不均勻的溫度場(chǎng)使結(jié)構(gòu)產(chǎn)生熱變形,熱變形在外界約束條件下會(huì)在結(jié)構(gòu)內(nèi)部產(chǎn)生預(yù)拉或預(yù)壓應(yīng)力,或二者同時(shí)存在,使結(jié)構(gòu)發(fā)生硬化或軟化。為此,熱環(huán)境下結(jié)構(gòu)剛度矩陣包含原剛度陣受溫度變化影響后的值KT和附加的應(yīng)力剛度Kσ,即

忽略結(jié)構(gòu)阻尼的影響,考慮熱效應(yīng)的結(jié)構(gòu)振動(dòng)方程為

其中:M為結(jié)構(gòu)的質(zhì)量矩陣;ω為結(jié)構(gòu)的固有頻率;φ為結(jié)構(gòu)的模態(tài)振型。

對(duì)兩種工況進(jìn)行分析:工況1為同時(shí)考慮熱應(yīng)力引起的附加剛度及材料特性退化引起結(jié)構(gòu)剛度的下降;工況2為只考慮材料特性的退化引起結(jié)構(gòu)剛度的下降。

熱模態(tài)分析分為兩步:a.加載給定加熱時(shí)刻下的“溫度載荷”,更新剛度矩陣,選中存儲(chǔ)每一步的剛度矩陣;b.調(diào)用修正后的剛度矩陣進(jìn)行模態(tài)計(jì)算。利用SOL 106求解器分析不同時(shí)刻溫度載荷下的應(yīng)力,選中“正則化模態(tài)分析”選項(xiàng)進(jìn)行模態(tài)分析,其結(jié)構(gòu)固有頻率如表2所示。對(duì)應(yīng)的前4階常溫模態(tài)如圖2所示,常溫模態(tài)頻率與文獻(xiàn)[16]中的數(shù)據(jù)相符。由表2可知,加熱到3ks時(shí),前4階固有頻率較常溫下降了9.4%,17.3%,12.9%和4.5%,其中扭轉(zhuǎn)模態(tài)(第2階和第3階模態(tài))頻率下降比較大。

表2 不同溫度載荷下機(jī)翼的固有頻率(工況1)Tab.2 Model frequencies at different temperature loads(case 1) Hz

圖2 常溫下機(jī)翼前4階模態(tài)振型圖Fig.2 Mode shapes of the first 4modes

表3 不同溫度載荷下機(jī)翼的固有頻率(工況2)Tab.3 Natural frequencies of the first 4modes at different temperature loads(case 2) Hz

忽略熱應(yīng)力引起的附加剛度效應(yīng)時(shí),固有頻率在不同時(shí)刻溫度載荷下的數(shù)據(jù)如表3所示。由表3可知,該工況下翼面加熱到3ks時(shí),前4階固有頻率較常溫時(shí)下降了4.4%,4.9%,5.5%和4.3%。

對(duì)比兩種工況下的頻率變化曲線見(jiàn)圖3。只考慮材料退化對(duì)結(jié)構(gòu)剛度的影響時(shí),結(jié)構(gòu)各階固有頻率隨著氣動(dòng)加熱的進(jìn)行一直降低。當(dāng)考慮熱應(yīng)力的變化產(chǎn)生附加剛度時(shí),熱效應(yīng)在氣動(dòng)加熱初期提高了結(jié)構(gòu)固有頻率,表明此時(shí)熱應(yīng)力的變化增大了結(jié)構(gòu)的總體剛度;但隨著氣動(dòng)加熱的進(jìn)行,熱應(yīng)力的引入在很大程度上降低了結(jié)構(gòu)總體剛度,特別是對(duì)扭轉(zhuǎn)剛度的影響很大,如圖3中的第2、第3支模態(tài)。

圖3 熱效應(yīng)對(duì)結(jié)構(gòu)固有頻率的影響Fig.3 Influences of thermal effect to the natural frequencies

4 顫振計(jì)算

4.1 氣彈運(yùn)動(dòng)方程

在超聲速、高超聲速非定常氣動(dòng)力的計(jì)算中,活塞理論能很好地滿足工程精度的要求[14]。活塞理論是一種無(wú)粘非定常氣動(dòng)理論,同時(shí)也是一種簡(jiǎn)化的氣動(dòng)力理論,只適用于在超聲速氣流下對(duì)機(jī)翼進(jìn)行顫振分析。實(shí)驗(yàn)證明,馬赫數(shù)在2~5之間,用這種理論對(duì)超聲速翼面計(jì)算都能得到滿足精度要求的結(jié)果[12-14]。

在等熵條件下由動(dòng)量方程可得

由于活塞前進(jìn)速度|v|?c∞,故有|v/c∞|?1,上式展開(kāi)后則可略去高階微分項(xiàng)。當(dāng)只保留一階項(xiàng)時(shí),稱為一階活塞理論,即為

當(dāng)保留二階項(xiàng)時(shí),可得到二階活塞理論,即

作用在機(jī)翼表面上的壓力分布可得到上、下表面的壓力差為

對(duì)于本研究中的均勻厚度機(jī)翼,此時(shí)厚度效應(yīng)為零,則上、下翼面壓力差

在熱環(huán)境下,熱氣動(dòng)彈性顫振方程[10]為

Kj包括熱應(yīng)力引起的附加剛度,其結(jié)果可由式(5)計(jì)算可得。針對(duì)方程(12),采用p-k法求解,即給定一系列速度,反復(fù)迭代求解顫振數(shù)據(jù)。

4.2 飛行環(huán)境

氣流馬赫數(shù)為3.0,空氣密度ρ=1.226kg/m3,飛行高度為海平面,空氣密度比為1.0,機(jī)翼攻角為0°。結(jié)構(gòu)單元為實(shí)體單元,氣動(dòng)網(wǎng)格點(diǎn)的位移和力通過(guò)無(wú)限板樣條方法插值到結(jié)構(gòu)單元上,在NASTRAN中采用p-k法進(jìn)行氣彈分析[15]。

4.3 結(jié)果分析

沿用熱模態(tài)分析中的兩種工況,工況1下,結(jié)構(gòu)在超聲速氣流中的顫振邊界如表4所示。與常溫下的顫振邊界相比,氣動(dòng)加熱到3ks時(shí),顫振速度和顫振頻率較常溫時(shí)下降了33.67%和26.49%。工況2下,顫振邊界如表5所示。相比較而言,忽略熱應(yīng)力的影響時(shí),3ks時(shí)段的顫振速度和顫振頻率較常溫時(shí)下降了5.11%和6.40%,與表3中固有頻率的下降幅度非常接近。

圖4反應(yīng)出熱應(yīng)力的引入對(duì)結(jié)構(gòu)顫振邊界的影響是很大的,盡管加熱初期熱應(yīng)力提高了固有頻率,增大了顫振速度和顫振頻率,但隨著氣動(dòng)加熱的進(jìn)行,熱應(yīng)力卻在很大程度上降低了顫振邊界,相比只有材料性能退化對(duì)顫振邊界的影響下降很多。

表4 不同溫度載荷下的顫振邊界(工況1)Tab.4 Flutter boundaries of the wing at different temperature loads(case 1)

表5 不同溫度載荷下的顫振邊界(工況2)Tab.5 Flutter boundaries of the wing at different temperature loads(case 2)

為了進(jìn)一步探討熱應(yīng)力的變化對(duì)顫振邊界的影響,筆者考慮了另外一種邊界約束條件下的顫振邊界,在根部弦處(0.2,0,0)和(0.8,0,0)兩點(diǎn)處固支,其對(duì)應(yīng)的顫振邊界如表6所示。

表6 不同溫度載荷下的顫振邊界(兩點(diǎn)約束)Tab.6 Flutter boundaries of the wing under two-node fixed boundary conditions

圖4 熱效應(yīng)對(duì)顫振特性的影響Fig.4 Influence of thermal effect to the flutter boundary

分析表6中的顫振數(shù)據(jù),在邊界約束較少的情況下,其顫振邊界在氣動(dòng)加熱初期也得到了提高。隨著氣動(dòng)加熱的進(jìn)行,顫振速度和顫振頻率較常溫時(shí)下降了11.72%和12.76%。對(duì)比機(jī)翼兩種邊界約束條件下分別在兩種工況下的顫振邊界曲線,如圖5所示。從圖中可以看出,邊界約束越少,熱應(yīng)力的變化對(duì)結(jié)構(gòu)顫振邊界的影響也越小。

5 結(jié) 論

1)高溫導(dǎo)致結(jié)構(gòu)材料特性的退化直接減小原始結(jié)構(gòu)剛度,降低結(jié)構(gòu)固有頻率。

2)熱應(yīng)力的引入在氣動(dòng)加熱初期提高了結(jié)構(gòu)各階固有頻率,然而隨著氣動(dòng)加熱的進(jìn)行,熱應(yīng)力卻在很大程度上降低了其各階固有頻率,其中扭轉(zhuǎn)模態(tài)頻率下降最為嚴(yán)重。

3)只考慮材料退化對(duì)結(jié)構(gòu)剛度的影響時(shí),顫振邊界的下降程度與固有頻率的下降程度接近。

圖5 不同邊界兩種工況顫振邊界曲線Fig.5 Flutter boundaries under two different boundary conditions

4)不同的邊界約束條件下熱應(yīng)力的分布是不同的,邊界約束越小,熱應(yīng)力的變化也會(huì)越小,其對(duì)結(jié)構(gòu)剛度的影響程度也越小;反之,影響會(huì)變大。

5)在翼面溫度分布趨向于穩(wěn)態(tài)的過(guò)程中,顫振邊界伴隨著氣動(dòng)加熱的進(jìn)行不斷發(fā)生變化,不同結(jié)構(gòu)邊界約束條件下,熱應(yīng)力在氣動(dòng)加熱初期均提高了顫振速度和顫振頻率,但最終都在很大程度上降低了顫振速度和顫振頻率。因此在高/超聲速飛行器設(shè)計(jì)過(guò)程中,必須要考慮由于氣動(dòng)加熱引起的熱效應(yīng)對(duì)結(jié)構(gòu)振動(dòng)特性和氣動(dòng)彈性穩(wěn)定性的影響,確保飛行器的飛行安全。

[1] McNamara J J,Thuruthimattam B J,F(xiàn)riedmann P P,et al.Hypersonic aerothermoelastic studies for reusable launch vehicles[C]∥45thAIAA/ASME/ASCE/AHS/ASC Structures,Structural Dynamic & Materials Conference.Palm Springs,California:AIAA,2004:1-35.

[2] Garrick I E.A survey of aerothermoelasticity[J].Aerospace Engineering,1963,22(1):140-147.

[3] Krist S L,Biedron R T,Rumsey C L.CFL3Duser′s manual(version 5.0)[R].[S.l.]:NASA,1998.

[4] Friedmann P P,McNamara J J,Thuruthimattam B J,et al.Aeroelastic analysis of hypersonic vehicles[J].Journal of Fluids and Structures,2004,19:681-712.

[5] Jack J M,Peretz P F,Kenneth G P,et al.Three-dimensional aeroelasticity and aerothermoelastic behavior in hypersonic flow[C]∥Proceedings of 46thAIAA/ASME/ASCE/AHS/ASC Structures,Structural Dynamics & Materials Conference. Austin,Texas:AIAA,2005:1-47.

[6] McNamara J J.Aeroelastic and aerothermoelastic behavior of two and three dimensional lifting surface in hypersonic flow[D].PH.D Thesis:The University of Michigan,2005.

[7] Heeg J,Zeiler T A,Pototzky A S,et al.Aerothermoelastic analysis of a NASP demonstrator model[R].[S.l.]:NASA,1993.

[8] 史曉明,楊炳淵.瞬態(tài)加熱環(huán)境下變厚度板溫度場(chǎng)及熱模態(tài)分析[J].計(jì)算機(jī)輔助工程,2006,15(S):15-18.

Shi Xiaoming,Yang Bingyuan.Temperature field and mode analysis of flat plate with thermal environment of transient heating[J].Computer Aided Engineering,2006,15(S):15-18.(in Chinese)

[9] 王宏宏,陳懷海,崔旭利,等.熱效應(yīng)對(duì)導(dǎo)彈翼面固有振動(dòng)特性的影響[J].振動(dòng)、測(cè)試與診斷,2010,30(3):275-279.

Wang Honghong,Chen Huaihai,Cui Xuli,et al,Thermal effect on dynamic characteristics of a missle wing[J].Journal of Vibration,Measurement & Vibration,2010,30(3):275-279.(in Chinese)

[10]吳志剛,惠俊鵬,楊超.高超聲速下翼面的熱顫振工程分析[J].北京航空航天大學(xué)學(xué)報(bào),2005,31(3):270-273.

Wu Zhigang,Hui Junpeng,Yang Chao.Hypersonic aerothermoelastic analysis of wings[J].Journal of Beijing University of Aeronautics and Astronautics,2005,31(3):270-273.(in Chinese)

[11]張偉偉.超聲速、高超聲速非線性氣動(dòng)彈性問(wèn)題研究[D].西安:西北工業(yè)大學(xué),2004.

[12]陳桂彬,鄒從青,楊超.氣動(dòng)彈性設(shè)計(jì)基礎(chǔ)[M].北京:北京航空航天大學(xué)出版社,2004:107-109.

[13]葉正寅,張偉偉,史愛(ài)明,等.流固耦合力學(xué)基礎(chǔ)及其應(yīng)用[M].哈爾濱:哈爾濱工業(yè)大學(xué)出版社,2010:107-108.

[14]趙永輝.氣動(dòng)彈性力學(xué)與控制[M].北京:科學(xué)出版社,2007:280-282.

[15]Rodden W P,Johnson E H.MSC/NASTRAN version 68,aeroelastic analysis user′s guide[R].Santa Ana:MSC.Software Corporation,2002.

[16]葉獻(xiàn)輝,楊翊仁.氣動(dòng)加熱下三角機(jī)翼顫振[J].西南交通大學(xué)學(xué)報(bào),2008:43(1):62-66.

Ye Xianhui,Yang Yiren.Flutter of delta wing under aerodynamic heating[J].Journal of Southwest Jiaotong University,2008,43(1):62-66.(in Chinese)

猜你喜歡
模態(tài)結(jié)構(gòu)分析
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
隱蔽失效適航要求符合性驗(yàn)證分析
論結(jié)構(gòu)
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動(dòng)化發(fā)展趨勢(shì)分析
論《日出》的結(jié)構(gòu)
國(guó)內(nèi)多模態(tài)教學(xué)研究回顧與展望
創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長(zhǎng)
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
由單個(gè)模態(tài)構(gòu)造對(duì)稱簡(jiǎn)支梁的抗彎剛度
主站蜘蛛池模板: 香蕉精品在线| 日韩欧美中文| 美女扒开下面流白浆在线试听| 国产精品白浆在线播放| av在线手机播放| 99精品伊人久久久大香线蕉 | 人妻无码一区二区视频| 日本在线视频免费| 一级成人a毛片免费播放| 日韩av电影一区二区三区四区| 亚洲精品在线影院| 98超碰在线观看| 免费A级毛片无码无遮挡| 亚洲第一在线播放| 在线播放精品一区二区啪视频| 中文字幕2区| 青青青视频91在线 | 精品视频一区二区观看| 97精品伊人久久大香线蕉| 久久香蕉国产线看观看式| 九色综合伊人久久富二代| 亚洲第一区精品日韩在线播放| 77777亚洲午夜久久多人| 无码啪啪精品天堂浪潮av| 免费高清毛片| 亚洲开心婷婷中文字幕| 久久国产精品电影| 国产福利小视频在线播放观看| 国产乱子伦视频在线播放| 久久亚洲精少妇毛片午夜无码 | 青青草原国产精品啪啪视频| 亚洲综合国产一区二区三区| 最新亚洲人成无码网站欣赏网| 精品国产成人高清在线| 亚洲色精品国产一区二区三区| 欧美丝袜高跟鞋一区二区| 99尹人香蕉国产免费天天拍| 亚洲区一区| 美美女高清毛片视频免费观看| 日韩精品无码不卡无码| 亚洲中文久久精品无玛| 谁有在线观看日韩亚洲最新视频| 啊嗯不日本网站| 国产99视频精品免费视频7| 国产人免费人成免费视频| 高清国产va日韩亚洲免费午夜电影| 欧美19综合中文字幕| 国产在线精彩视频二区| 国产成人高精品免费视频| 伊人久久精品亚洲午夜| 久久精品视频一| 国产精品视频系列专区| 激情亚洲天堂| 国产精品美女免费视频大全| 国产激情无码一区二区三区免费| 久久人与动人物A级毛片| 国产精品毛片一区| 国产精品嫩草影院av| 久久精品丝袜| 欧美三级不卡在线观看视频| 亚洲va欧美ⅴa国产va影院| 亚洲色成人www在线观看| 国产精品久久久久久久久久久久| 午夜限制老子影院888| 国产91精选在线观看| 91青青草视频| 国产地址二永久伊甸园| 欧美综合区自拍亚洲综合天堂| 视频二区国产精品职场同事| 欧美有码在线| 久久综合丝袜日本网| 天天视频在线91频| 美女一区二区在线观看| 国产亚洲精品91| 天天综合网站| 国产免费羞羞视频| 久久久亚洲国产美女国产盗摄| 99久久国产精品无码| 色综合久久久久8天国| 三级国产在线观看| 操操操综合网| 成人福利在线视频免费观看|