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

考慮幾何非線性的氣動彈性模型縮比方法

2019-05-05 02:01:46柴睿譚申剛黃國寧
北京航空航天大學學報 2019年4期
關鍵詞:模態(tài)變形方法

柴睿, 譚申剛, 黃國寧

(航空工業(yè)第一飛機設計研究院, 西安 710089)

近些年來,無人機得到了長足發(fā)展。為了延長航行時間并降低能源消耗,無人機一般采用大展弦比高柔性機翼設計以追求較高的氣動效率[1]。同時為了減輕飛機重量,機翼設計大都選擇輕質材料,所以大展弦比無人機機翼一般都是柔性機翼。然而,機翼的展弦比增大以及重量減小會帶來新的挑戰(zhàn),這就要求對機翼在正常飛行載荷下的大變形問題做進一步研究,因為大變形問題會帶來幾何非線性行為及相關的氣動彈性問題[2-3]。這些大變形會改變機翼的固有頻率,從而使其氣動彈性行為發(fā)生顯著變化[4]。

在進行氣動彈性特性理論計算的同時,通過風洞試驗來確定飛機的氣動彈性特性是必不可少的[5]。縮比模型的氣動彈性試驗是飛機研制過程中最基本的保證,縮比模型應盡可能地反映出全尺寸飛機的氣動彈性響應[6]。使用縮比模型來進行試驗可以更好地理解機翼的物理行為,同時降低制造全尺寸飛機的風險與成本[7]。

氣動彈性模型縮比設計要求同時考慮氣動力和結構的物理特性。氣動力的模擬可以通過氣動外形幾何縮比再進行分析計算得到。而通常情況下,通過結構的幾何尺寸縮比對結構進行模擬是不現實的。根據縮比設計的要求,指定幾何縮比的結構,其制造所用的材料通常是不存在的,而且很大程度上用于制造全尺寸飛機的加工技術在縮比模型上是不可能實現的。唯一可行的方法是對內部結構進行重新設計并優(yōu)化,比如對結構的質量、剛度、阻尼特性進行縮比,并且與全尺寸飛機保持協調[8]。

目前研究考慮非線性的氣動彈性模型縮比的相關文獻較少,只是最近幾年該課題才熱門起來,這些研究主要針對的是具有大展弦比高柔性機翼的高空長航時無人機[9]和聯翼布局飛機[10-11]。

Wan和Cesnik[9]推導得到包含大變形引起的非線性剛度和預應力影響的氣動彈性運動方程的縮比參數,研究表明線性和幾何非線性的氣動彈性運動方程的縮比因子在形式上是相同的。同時Wan和Cesnik[9]對大柔性飛機2種形式的縮比設計進行了分析,其中一種設計有意忽略了Froude數,他們對包括結構模態(tài)、非線性氣動彈性配平、不同載荷下的線性長周期穩(wěn)定性以及非線性瞬態(tài)陣風響應等特性進行了計算。結果表明:包含Froude數的縮比模型響應與全尺寸模型吻合得很好,所以Froude數的相似在驗證的問題上不能被忽略。

Bond等[10]提出了一個非線性氣動彈性模型縮比方法,用來設計具有幾何非線性特征的聯翼飛機縮比模型。其設計的結構可以匹配前3階頻率、前3階模態(tài)振型以及第1階線性屈曲特征解,且縮比模型的氣動彈性頻率和阻尼在整個速度包線上都匹配得很好。同時,結果表明,縮比模型的非線性靜響應,在60%屈曲載荷下都匹配得很好。

Ricciardi等[11]發(fā)展的幾何非線性氣動彈性模型縮比設計方法可以直接匹配縮比的線性和非線性靜響應,同時可以滿足模態(tài)頻率的約束。他們用該方法設計的模型同利用傳統方法并經過優(yōu)化的模型進行比較,發(fā)現用前面提到的方法設計并優(yōu)化的縮比模型其非線性氣動彈性響應同目標響應吻合較好,而利用傳統方法設計的縮比模型其非線性氣動彈性響應同目標響應相比誤差較大。

基于上述觀點,本文提出一個新的非線性氣動彈性模型縮比方法,稱為動力學有限元模型的非線性靜響應-模態(tài)協同優(yōu)化方法,其中非線性靜響應的匹配是基于等效靜態(tài)載荷法來實現的。這是首次將處理非線性靜力學的等效靜態(tài)載荷法與大展弦比機翼結構動力學相結合的優(yōu)化方法。

另外該方法與Ricciardi等[11]提出的方法主要區(qū)別是縮比方法的實現過程不同。本文方法中剛度和質量分布設計是在2個優(yōu)化循環(huán)過程中進行的,而Ricciardi等的方法只是在一個優(yōu)化循環(huán)中進行的。而且本文驗證的實例是一架具有超大展弦比高柔性機翼的某型無人機,它的幾何非線性效應比上述文獻中驗證的飛機要強得多。

1 理論方法

1.1 相似準則

縮比因子的選擇是氣動彈性縮比模型設計過程中關鍵的第1步。通常選定氣動彈性縮比模型設計的3個基本縮比因子:kb、kv和kρ。這里,kb為縮比模型與全尺寸模型的長度比,kv為縮比模型與全尺寸模型的速度比,kρ為縮比模型與全尺寸模型所處環(huán)境的空氣密度比。

為了滿足縮比模型和全尺寸模型之間的氣動彈性相似性,在前3個基本縮比因子[12-13]的基礎上,必須滿足下面描述的4個相似準則。當要求在重力作用下氣動彈性響應滿足相似要求時,必須滿足重力相似準則作為附加條件。

1) 幾何相似

幾何相似的縮比因子即前述的kb。幾何相似意味著縮比模型與全尺寸模型之間的所有尺寸比(包括翼展比,弦長比等)都與kb相同。同時還假定縮比模型和全尺寸模型具有相同的翼型。

2) 質量相似

滿足質量相似可以寫為

(1)

這意味著縮比模型的分布質量應該通過質量縮比因子km與全尺寸模型的分布質量成比例。

3) 剛度相似

滿足剛度相似可以描述為

(2)

這意味著縮比模型的分布剛度應該通過剛度縮比因子kK與全尺寸模型的分布剛度成比例。

基于固有頻率與質量和剛度的關系,頻率的縮比因子即為

(3)

4) 氣動相似

當縮比模型和全尺寸模型之間的馬赫數與雷諾數分別滿足一致性要求時,幾何相似便會帶來氣動相似。但是,馬赫數相似和雷諾數相似都是不容易滿足的。對于大展弦比高柔性無人機的典型低速飛行,馬赫數效應可以忽略。但是,縮比模型與全尺寸模型之間的雷諾數應該相同或至少在相同的數量級上。

雷諾數的縮比因子可以定義為

kRe=kρkvkb/kμ

(4)

式中:kμ為空氣動黏度的縮比因子。

基于升力方程,升力的縮比因子為

(5)

5) 重力相似

Froude數決定了重力載荷下撓度與空氣動力和慣性載荷下撓度的比值。當顫振試驗中考慮重力效應,縮比模型和全尺寸模型之間的Froude數應相同,可以通過式(6)來滿足[10]:

(6)

與線性結構相比,考慮幾何非線性的結構具有其特殊的剛度特性,其特殊性取決于當前的結構變形。顯然,當縮比模型和全尺寸模型在未變形狀態(tài)和任一給定的變形狀態(tài)下均滿足幾何相似,質量相似和Froude數相似要求,則縮比模型和全尺寸模型在任何變形狀態(tài)下上述相似關系也是成立的。因此,這取決于考慮幾何非線性的結構隨結構變形而變化的剛度是否可以在縮比模型和全尺寸模型之間保持相似。這將直接影響線性氣動彈性相似準則對考慮幾何非線性的結構的適用性。

Wan和Cesnik[9]的研究表明線性和幾何非線性氣動彈性運動方程的縮比因子形式上是相同的。因此,若縮比模型是基于上述縮比因子設計的,則縮比模型與全尺寸模型之間的氣動彈性行為在小變形狀態(tài)和大變形狀態(tài)下都是相似的。由式(1)~式(6)可知,對于考慮幾何非線性的氣動彈性縮比模型而言,當重力相似準則作為附加條件時,在長度比kb、空氣密度比kρ確定的情況下,作為模型設計目標的模型質量、模型固有頻率、模型顫振速度和顫振頻率便可隨之確定。

1.2 縮比方法

1.2.1 傳統線性縮比方法

從線彈性的控制方程可以看出,經過對全尺寸飛機適當的縮比,對于不同結構形式的縮比模型,只要具有相同的幾何縮比,以及相同的分布剛度和質量縮比,將帶來相同的模態(tài)特性(固有頻率和模態(tài)振型)[14]。因此可以利用這一理論從全尺寸飛機那里獲得縮比模型的目標值(質量、慣性矩、固有頻率、模態(tài)振型和氣動外形)。通過設計一個與全尺寸模型相似的模態(tài)特性的縮比模型,得到與內部結構無關的分布剛度和質量相似性,即可采用簡化的內部結構。一般情況下,為了匹配目標縮比值,有2種不同的傳統線性縮比方法:剛度質量耦合匹配模態(tài)響應法(方法A);剛度質量解耦匹配模態(tài)響應法(方法B)。圖1和圖2分別是方法A和B的優(yōu)化設計流程圖。

圖1 傳統線性縮比方法(方法A)Fig.1 Traditional linear scaling method (Method A)

圖2 傳統線性縮比方法(方法B)Fig.2 Traditional linear scaling method (Method B)

方法A類似于Bisplinghoff等[12]提出的方法,方法B則是French和Eastep[6]提出的。不過上述2種方法都是針對大部分的線性問題,對于非線性問題,這2個方法的適用性還需要進行進一步的驗證。

1.2.2 考慮幾何非線性的縮比方法

與方法B類似,通過分別匹配解耦的剛度和質量分布來匹配模態(tài)特性。該方法(方法C)的 流程如圖3所示。方法B與方法C的區(qū)別體現在第1階段的剛度優(yōu)化過程中。在方法C中,縮比模型的剛度優(yōu)化引入了等效靜態(tài)載荷法,而且通過非線性靜力分析來獲得非線性廣義位移場。這個優(yōu)化過程可以在MATLAB中實現,靜力分析與模態(tài)分析使用MSC.Nastran軟件,優(yōu)化準則使用粒子群尋優(yōu)算法。當第1階段的剛度優(yōu)化完成后,就可以進行第2階段的質量優(yōu)化,質量優(yōu)化的目的是在匹配了非線性靜響應后匹配模態(tài)特性(固有頻率和模態(tài)振型)。

圖3 考慮幾何非線性的縮比方法(方法C)Fig.3 Scaling method considering geometric nonlinearity (Method C)

研究發(fā)現,對非線性靜力響應的直接優(yōu)化很困難。對于這項任務,基于等效靜態(tài)載荷法來應用等效線性系統的效率會更高[15]。等效靜態(tài)載荷法就是通過不斷迭代更新線性系統中的設計變量來獲得與非線性靜力分析結果相同的廣義位移場。

首先給定當前或初始設計變量g0gggggg和載荷條件{F},然后非線性靜力分析通過求解非線性平衡方程(式(7))得到非線性變形{xNL}。

[KNL(d)]{xNL}={F}

(7)

由線性剛度矩陣[KLin(d)]乘以非線性變形{xNL}得到等效載荷{Feq}:

{Feq}=[KLin(d)]{xNL}

(8)

由等效載荷{Feq}和線性剛度矩陣[KLin(d)]將可得到等效線性靜態(tài)廣義位移{xLin}:

[KLin(d)]{xLin}={Feq}

(9)

接下來在等效線性系統中便可以執(zhí)行高效的優(yōu)化,等效靜態(tài)載荷法需要一個迭代過程來更新線性化系統的設計變量,優(yōu)化收斂時的線性化的最佳響應與非線性系統的響應相匹配。

(10)

式中:i為迭代次數;ε為誤差參數(滿足設計要求的小值)。應用等效靜態(tài)載荷法的過程如圖4所示。

該優(yōu)化過程是在2個域(分析域和設計域)里實施的。首先,在分析域中進行非線性靜力分析,獲得非線性廣義位移場;然后,利用位移場通過等效靜態(tài)載荷法獲得新的載荷條件,把新的載荷條件應用到設計域的線性靜力優(yōu)化中;最后,通過不斷迭代更新設計域中的設計變量來使非線性系統和線性系統的廣義位移場誤差滿足要求。該過程的具體實現如圖5所示。

本文目的是對具有幾何非線性特征的結構進行氣動彈性縮比,設計出的縮比模型要與全尺寸 飛機的模態(tài)特性吻合,所以筆者根據傳統線性縮比方法結合等效靜態(tài)載荷法來實現這一縮比過程。

圖4 分析域和設計域優(yōu)化流程圖Fig.4 Flowchart of analysis and design domain optimization

圖5 等效靜態(tài)載荷法的實現過程圖Fig.5 Implementation procedure of equivalent static loads method

2 模 型

2.1 全尺寸模型

本文縮比的全尺寸模型是一架具有大展弦比柔性機翼的某型機。該機機翼的展弦比為31.25,屬于超大展弦比機翼。該機的巡航馬赫數Ma=0.1,巡航高度為20 km。機翼的主要參數在表1中列出。

整個機翼包括1段10 m長的中央翼,2段10 m長的中外翼,2段10 m長的外翼。整個機翼包括梁、肋和前后緣,翼梁為盒形梁,梁截面寬度為240 mm,高度為170 mm。翼肋采用夾芯肋,前后緣為殼形蜂窩夾芯結構。全尺寸機翼的平面布局如圖6所示,該機結構布局如圖7所示。本文使用MSC.Nastran軟件對全尺寸模型進行求解計算。

本文首先對全尺寸模型進行線性和非線性靜力分析,以評估2種分析之間的差異。圖8顯示了隨著過載系數的增加,翼尖線性與非線性垂向變形量變化趨勢的差異。根據定義,通過增加升力來實現不同的過載系數:

(11)

式中:n為過載系數;L和W分別為飛機的總升力和總重力。分析發(fā)現,隨著過載系數的增加,線性和非線性變形量變化曲線逐漸分離。當過載系數是4時,翼尖非線性變形量與線性變形量相差4.07 m。這是因為線性靜力分析時不考慮升力跟隨后的軸向分力,所以翼尖不會向翼根方向移動。從圖8還可以看出,過載系數為1.5時,非線性效應還不明顯,這是因為機翼剛度稍大,載荷較小時,結構還沒有到達非線性狀態(tài)造成的。當過載系數為3時,截面轉角θ和扭轉角φ沿展向的非線性和線性變化趨勢如圖9所示, 其中b為翼展長度。由圖9可知,當進行非線性靜力分析時,結構會經歷一個剛度硬化效應。在圖9中線性和非線性變化趨勢分離的地方出現在翼展20%站位的位置,這是由剛度變化引起的。

表1 全尺寸機翼的主要參數Table 1 Main parameters of full-scale wing

圖7 全機結構及內部結構Fig.7 Whole aircraft structure and internal structure

圖8 翼尖變形趨勢Fig.8 Variation of wing tip deflection

圖6 機翼平面布局
Fig.6 External configuration of wing

針對本文的研究,結構幾何非線性主要歸因于大位移狀態(tài),處理大位移問題時,利用跟隨力的非線性求解器可以大幅改變系統的載荷和剛度矩陣,從而改變變形狀態(tài)。圖8和圖9說明了線性和非線性靜力分析之間的區(qū)別以及隨著結構的變形增大,不斷地引入附加剛度造成的這些差異。

圖9 截面轉角和扭轉角的變化趨勢Fig.9 Variation of cross-section bend angle and twist angle

2.2 縮比模型

本文縮比模型選取基本縮比因子kb=0.1,kρ=13.764,kv=0.316,進而目標縮比參數可以通過這些縮比因子確定。為了增大優(yōu)化過程中可行解區(qū)域,設計變量選取截面的“剛度底層控制信息”,即截面積A、截面慣性矩Ia和Ib以及截面扭轉常數J。因此在縮比優(yōu)化之前不需要對縮比模型結構形式進行設計,可在優(yōu)化結束后對截面“剛度底層控制信息”進行逆運算來獲取截面幾何信息。為了使縮比模型正確描述全尺寸模型的氣動彈性特性,確定匹配兩者前10階固有頻率和相應的模態(tài)振型,其相似度不低于90%。全尺寸模型顫振主要由第4階模態(tài)(水平二彎)和第5階模態(tài)(一扭)參與,所以主要看第4階與第5階模態(tài)匹配相似度。同時,縮比模型和全尺寸模型的總質量匹配相似度不低于90%。

就外部幾何布局而言,縮比模型是全尺寸模型的精確幾何縮比。在縮比的過程中,還要尋找最佳材料來模擬縮比模型的結構部件。全尺寸模型所采用的材料為7075鋁合金,然而精確縮比后的材料是不存在的。所以在縮比之前要對縮比模型的材料進行研究并選擇。最終決定使用較靠近材料縮比目標值的鑄鋁青銅作為縮比模型的材料。

最優(yōu)設計變量的目標是使縮比模型和全尺寸模型的一些特性(質量、固有頻率、模態(tài)振型以及在考慮幾何非線性的方法中的非線性靜響應)匹配。殘差平方和(SSR)是量化響應差異的一個基本參考值。

3 計算結果

本文采用粒子群尋優(yōu)算法對縮比模型的設計變量進行優(yōu)化,來使得縮比模型的各階固有頻率和模態(tài)振型與全尺寸模型對應階的固有頻率和模態(tài)振型相似度達到90%以上。1.2節(jié)提到了2種傳統線性縮比方法(方法A與方法B)以及考慮幾何非線性的縮比方法(方法C)。

表2為傳統線性縮比方法與考慮集合非線性縮比方法結果參數比較,其中方法B和方法C分別給出了階段1和階段2的參數值。從表2可以看出,方法A和方法B相比于方法C耗時較少,所以當結構非線性效應不明顯時,推薦使用傳統線性縮比方法。然而,當考慮了結構的非線性效應,用傳統線性縮比方法設計的模型計算結果誤差較大,無法處理考慮幾何非線性的模型設計。

為了高效地捕捉結構的非線性行為,本文提出使用等效靜態(tài)載荷法進行非線性靜響應匹配,來對動力學有限元模型的非線性靜響應-模態(tài)進行協同優(yōu)化。對使用考慮幾何非線性的方法獲得的非線性翼尖變形量和傳統線性縮比方法設計的模型翼尖變形量以及變形量目標縮比值進行了比較。分析發(fā)現,方法C設計的模型非線性靜變形量相比于方法B,更為靠近目標值(圖10)。

圖11顯示了傳統線性縮比方法(方法A、方法B)和考慮幾何非線性的縮比方法(方法C)在前10階的固有頻率、模態(tài)振型上匹配相似度的差異。目標縮比值是從全尺寸模型的相關計算中獲得到。

表2 傳統線性縮比方法與考慮幾何非線性的 縮比方法結果參數比較Table 2 Comparison of result parameters between traditional linear scaling method and scaling method considering geometric nonlinearity

圖10 非線性翼尖垂向變形Fig.10 Nonlinear wing tip vertical deformation

本文的顫振計算使用偶極子格網法求解非定 常氣動力模型。為了匹配全尺寸模型的氣動彈性響應,分別采用方法B和方法C對模型進行設計,表3是2種模型顫振求解結果與目標縮比值的比較,通過考慮幾何非線性的氣動彈性模型縮比方法設計的模型,其顫振結果誤差得到了顯著降低。圖12(a)為第4階和第5階模態(tài)的V-g(速度-阻尼系數)圖,圖12(b)為第4階和第5階模態(tài)的V-f(速度-頻率)圖。為了便于觀察,圖中是把縮比模型的顫振計算結果進行逆縮比,其逆縮比的值可直接與全尺寸機翼的顫振計算結果進行對比。因為顫振主要是第4階模態(tài)(面內二彎)和第5階模態(tài)(一扭)引起,所以在圖12中,只顯 示了這2階模態(tài)。當采用幾何非線性縮比方法時,氣動彈性響應匹配相似度有所改進,其精度明顯優(yōu)于傳統線性縮比模型得到的結果。

表3 兩種方法結果與目標顫振速度及頻率誤差Table 3 Flutter speed and frequency differences between Method B and C in relation to associated target scaled values %

圖12 V-g 和V-f圖Fig.12 V-g and V-f diagram

4 結 論

本文研究了氣動彈性模型的傳統線性縮比方法和考慮幾何非線性的縮比方法,得到以下結論:

1) 本文提出的考慮幾何非線性的氣動彈性模型縮比方法能夠針對幾何非線性效應明顯的氣動彈性模型進行縮比,縮比模型與全尺寸模型的高階模態(tài)匹配度較高,說明該方法可以捕捉到傳統線性縮比方法無法捕捉到的高階信息,進而縮比模型能夠更好地再現全尺寸飛機的非線性氣動彈性行為。

2) 傳統線性縮比方法求解平均時間為1.2 h,平均誤差為9.03%;幾何非線性氣動彈性模型縮比方法求解時間為2.9 h,誤差為2.80%。可以看出,雖然傳統線性縮比方法求解時間少于幾何非線性氣動彈性模型縮比方法求解時間,但是前者誤差已經無法滿足模型設計要求。所以對于考慮幾何非線性的結構,從縮比設計相似精度要求來看,傳統線性縮比方法是不適用的,推薦使用本文的新方法。

3) 方法A(耦合)求解時間高于方法B(解耦)求解時間,但方法A誤差略低于方法B誤差,雖然兩者對于考慮幾何非線性的結構縮比問題都無法滿足設計要求,但對于幾何非線性效應不明顯的結構,這2種方法都是可以處理相關縮比問題的。

4) 等效靜態(tài)載荷法的優(yōu)化邏輯是把原先非線性求解過程中的剛度無規(guī)律低效更新轉變?yōu)閯偠扔蟹较蚋咝Ц拢褂玫刃ъo態(tài)載荷法來進行非線性靜響應的匹配在保證匹配精度的同時能大大降低時間成本。

本文的后續(xù)工作是要對具有大展弦比機翼的飛機按照文中提出的考慮幾何非線性的縮比方法縮比后進行風洞試驗。從長遠來看,本文所得結論也能應用到設計飛行試驗的模型上,驗證飛行過程中的幾何非線性氣動彈性響應。同時,還需要將本文所提的考慮幾何非線性的縮比方法應用于更復雜的結構以及更高精度的氣動模型,并進一步驗證其有效性。

猜你喜歡
模態(tài)變形方法
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
“我”的變形計
例談拼圖與整式變形
會變形的餅
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
國內多模態(tài)教學研究回顧與展望
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
基于HHT和Prony算法的電力系統低頻振蕩模態(tài)識別
由單個模態(tài)構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 免费又黄又爽又猛大片午夜| 国产91线观看| 欧美视频在线播放观看免费福利资源| 久久国产乱子| 精品少妇人妻av无码久久| 国产香蕉97碰碰视频VA碰碰看| 欧美精品在线视频观看| 凹凸精品免费精品视频| 丰满人妻被猛烈进入无码| 黄色不卡视频| 国产精品国产三级国产专业不| a级毛片一区二区免费视频| 在线国产毛片手机小视频| 免费毛片视频| 国产成人精品无码一区二| 国产一区二区福利| 激情乱人伦| 精品国产福利在线| 亚洲另类第一页| 色网在线视频| 亚洲人网站| 欧美五月婷婷| 一级香蕉人体视频| 亚洲无码视频一区二区三区| 美女被狂躁www在线观看| 中文字幕 日韩 欧美| 国产欧美精品一区二区| 蜜桃视频一区| 在线观看亚洲成人| 日韩资源站| 欧美成人h精品网站| 一区二区三区四区精品视频| 国产欧美视频在线| 国产激情无码一区二区APP | 亚洲高清无码精品| 久久无码高潮喷水| 国产成人一二三| 色丁丁毛片在线观看| 91在线高清视频| 中文字幕 欧美日韩| 日本午夜精品一本在线观看 | 色婷婷在线播放| 国产成人综合亚洲欧洲色就色| 视频二区国产精品职场同事| 亚洲国产成熟视频在线多多| 少妇精品久久久一区二区三区| 亚洲欧洲综合| 97视频免费看| 国产网站免费观看| 国产综合另类小说色区色噜噜| 亚洲中字无码AV电影在线观看| 亚洲精品无码抽插日韩| 精品久久久久久久久久久| 91免费国产高清观看| 国内精品九九久久久精品| 久久夜色撩人精品国产| 国产黄在线观看| 国产97色在线| 欧美国产日韩另类| 免费精品一区二区h| 亚洲精品中文字幕午夜| 亚洲男人的天堂久久香蕉| 国禁国产you女视频网站| 国产亚洲欧美在线专区| 亚洲视频黄| 一本色道久久88综合日韩精品| 就去吻亚洲精品国产欧美| 国产在线观看91精品亚瑟| 91精品国产一区| 试看120秒男女啪啪免费| 视频一区视频二区日韩专区| 国产精品黄色片| 无码高潮喷水在线观看| 在线观看国产精美视频| 免费人成视频在线观看网站| 国产精品自在在线午夜区app| 97se亚洲综合在线天天| 人妖无码第一页| 欧美日本在线观看| 国产精品99r8在线观看| 亚洲国模精品一区| 亚洲国产欧美目韩成人综合|