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

基于伴隨理論的無人機(jī)氣動(dòng)彈性優(yōu)化

2023-10-31 11:39:26黃胤錚王曉喆王柳青
軟件導(dǎo)刊 2023年10期
關(guān)鍵詞:有限元優(yōu)化結(jié)構(gòu)

黃胤錚,王曉喆,王柳青

(1.北京航空航天大學(xué) 飛行學(xué)院;2.北京航空航天大學(xué) 無人系統(tǒng)研究院,北京 100083)

0 引言

無人機(jī)因具有體積小、成本低、對惡劣環(huán)境適應(yīng)性好等優(yōu)點(diǎn),在民用和軍用領(lǐng)域均有巨大的應(yīng)用潛力[1]。但無人機(jī)的作戰(zhàn)能力會(huì)受到飛機(jī)結(jié)構(gòu)、氣動(dòng)布局、任務(wù)載荷、巡航時(shí)間等因素影響,無人機(jī)任務(wù)完成度以及生存概率也與其飛行性能和作戰(zhàn)性能密切相關(guān)。隨著航空技術(shù)的發(fā)展,未來空域的環(huán)境復(fù)雜度都會(huì)有不同程度的上升,對無人機(jī)的整體性能提出了更高要求。

隨著有效載荷不斷增加,無人機(jī)的結(jié)構(gòu)重量系數(shù)不斷下降,結(jié)構(gòu)柔度不斷增加,氣動(dòng)、結(jié)構(gòu)等學(xué)科的耦合關(guān)系進(jìn)一步增強(qiáng),若在優(yōu)化過程中采用串列優(yōu)化方法會(huì)造成學(xué)科間的反復(fù)迭代,在降低設(shè)計(jì)效率的同時(shí),也影響了整體性能的提升。氣動(dòng)彈性優(yōu)化充分考慮了氣動(dòng)與結(jié)構(gòu)之間的耦合效應(yīng),并以此挖掘設(shè)計(jì)潛力[2],實(shí)現(xiàn)飛行器在復(fù)雜環(huán)境下經(jīng)濟(jì)性和使用性能的提升[3]。在大型客機(jī)[4-5]和超音速客機(jī)[6]的設(shè)計(jì)中,氣動(dòng)彈性優(yōu)化已經(jīng)得到了較為廣泛的應(yīng)用。基于敏度信息的優(yōu)化算法在求解大規(guī)模的多學(xué)科優(yōu)化問題時(shí),收斂較快[7-8],而以遺傳算法為代表的進(jìn)化式算法需要大量樣本,不適合直接應(yīng)用于多設(shè)計(jì)變量的氣動(dòng)/結(jié)構(gòu)耦合優(yōu)化問題中[9]。1988 年,Jameson[10]首次將伴隨方程法應(yīng)用于氣動(dòng)目標(biāo)函數(shù)梯度的求解過程中,并對機(jī)翼進(jìn)行優(yōu)化設(shè)計(jì)。由于伴隨理論的計(jì)算消耗與設(shè)計(jì)變量的數(shù)量無關(guān),適合處理具有大規(guī)模設(shè)計(jì)變量的多學(xué)科優(yōu)化問題,因此近年來得到了廣泛關(guān)注[11-14]。Anderson 等[15]采用多學(xué)科建模優(yōu)化方法和伴隨理論,實(shí)現(xiàn)了風(fēng)力渦輪機(jī)葉片的高保真結(jié)構(gòu)優(yōu)化和負(fù)載應(yīng)力最小化;Wang 等[16]采用基于伴隨方法的優(yōu)化設(shè)計(jì)框架對機(jī)翼受動(dòng)載荷進(jìn)行了分析優(yōu)化;Batay 等[17]采用多學(xué)科設(shè)計(jì)優(yōu)化工具,結(jié)合伴隨求解器,對風(fēng)力渦輪機(jī)進(jìn)行了并行氣動(dòng)設(shè)計(jì)優(yōu)化;Wilke[18]采用自主開發(fā)的多目標(biāo)優(yōu)化框架實(shí)現(xiàn)了對直升機(jī)旋翼槳葉的優(yōu)化設(shè)計(jì);李潤澤等[19]在多目標(biāo)優(yōu)化過程中結(jié)合伴隨方法,提高了超臨界機(jī)翼氣動(dòng)優(yōu)化設(shè)計(jì)效率;劉曉東等[20]采用基于伴隨理論的氣動(dòng)優(yōu)化設(shè)計(jì)方法,實(shí)現(xiàn)了飛翼布局飛行器的氣動(dòng)優(yōu)化。

隨著現(xiàn)代無人機(jī)對高機(jī)動(dòng)以及長航時(shí)性能的需求不斷提高,采用基于伴隨理論的氣動(dòng)彈性優(yōu)化方法深度挖掘多學(xué)科耦合所帶來的設(shè)計(jì)潛力是必要且迫切的。本文結(jié)合渦格法和有限元分析法,通過松耦合的方式實(shí)現(xiàn)氣動(dòng)彈性分析,并基于序列二次規(guī)劃和伴隨理論開展非線性敏度優(yōu)化。分別針對長航時(shí)和高機(jī)動(dòng)兩類典型無人機(jī),以最小燃油消耗為目標(biāo),對機(jī)翼的氣動(dòng)外形和空間梁結(jié)構(gòu)等參數(shù)進(jìn)行優(yōu)化設(shè)計(jì),并對優(yōu)化設(shè)計(jì)的結(jié)果進(jìn)行對比分析。

1 氣動(dòng)彈性仿真

1.1 氣動(dòng)渦格法

對于飛行器的氣動(dòng)分析本質(zhì)上是對空間N-S 方程的求解。渦格法在位勢流方程的求解過程中較為常見,與空間N-S 方程的求解相比,利用基于求解位勢流的氣動(dòng)計(jì)算方法能夠以更高的計(jì)算效率對啟動(dòng)導(dǎo)數(shù)進(jìn)行比較準(zhǔn)確的求解[21]。

根據(jù)邊界不可穿透定理,可以在控制點(diǎn)處添加馬蹄渦與氣動(dòng)表面相切的條件,如式(1)所示。其中,CC為控制點(diǎn)處誘導(dǎo)速度的系數(shù)矩陣,Γ為每個(gè)渦線渦強(qiáng)組成的列向量,V∞為控制點(diǎn)來流速度矩陣,N為控制點(diǎn)處單位法向量。

通過求解式(1)可得到每個(gè)渦的渦流強(qiáng)度和氣動(dòng)力,如式(2)所示。其中,F(xiàn)n為渦段上的升力,ρ∞為來流密度,V∞為來流速度,Vn為當(dāng)前渦段中點(diǎn)誘導(dǎo)速度,Γn為當(dāng)前渦段上的渦強(qiáng)。

渦格法渦線布置方法如圖1所示。

Fig.1 Vortex lattice method vortex line layout method圖1 渦格法渦線布置方法

1.2 結(jié)構(gòu)有限元法

對于組成機(jī)翼的各部分而言,機(jī)翼上下表面的蒙皮主要承載機(jī)翼上的彎曲載荷,機(jī)翼的翼梁和副梁主要承載剪切載荷,蒙皮和梁共同承載扭轉(zhuǎn)載荷。在分析過程中,將機(jī)翼翼梁簡化為空間梁單元的有限元模型進(jìn)行分析。空間梁單元的有限元模型如圖2所示[22]。

Fig.2 Finite element model of space beam element圖2 空間梁單元有限元模型

通過對空間梁模型的分析,可以獲得結(jié)構(gòu)分析部分的控制方程如下:

其中,f為作用在節(jié)點(diǎn)上的力和力矩,d為節(jié)點(diǎn)位移列陣。通過求解式(3),可以得到各節(jié)點(diǎn)上的靜力學(xué)載荷分布。

1.3 仿真模型

氣動(dòng)彈性的耦合求解采用松耦合方法,氣動(dòng)與結(jié)構(gòu)之間的數(shù)據(jù)傳遞如圖3、式(4)和式(5)所示[23]。

Fig.3 Aerodynamic load transfer圖3 氣動(dòng)載荷傳遞

其中,i=1,2,分別代表氣動(dòng)網(wǎng)格對應(yīng)有限元模型中的左右兩個(gè)結(jié)構(gòu)節(jié)點(diǎn),T為氣動(dòng)網(wǎng)格上通過渦格法計(jì)算得到的合力,rsai為結(jié)構(gòu)節(jié)點(diǎn)指向氣動(dòng)力作用點(diǎn)的方向向量,F(xiàn)i、Mi分別為通過該牽引方案得到的節(jié)點(diǎn)力和節(jié)點(diǎn)扭矩。

在完成氣動(dòng)載荷傳遞后,通過有限元分析可以得到對應(yīng)的結(jié)構(gòu)位移,并傳遞回氣動(dòng)網(wǎng)格。

其中,i=1,2,分別代表氣動(dòng)網(wǎng)格對應(yīng)有限元模型中的左右兩個(gè)結(jié)構(gòu)節(jié)點(diǎn),di為平動(dòng)位移,θi為轉(zhuǎn)動(dòng)位移,rsai同式(5),為結(jié)構(gòu)節(jié)點(diǎn)指向氣動(dòng)力作用點(diǎn)的方向向量。在此基礎(chǔ)上,結(jié)合開源框架OpenMDAO 建立無人機(jī)的氣動(dòng)、結(jié)構(gòu)以及氣動(dòng)彈性耦合分析模型[24]。

2 優(yōu)化算法及框架

2.1 伴隨理論

傳統(tǒng)的導(dǎo)數(shù)計(jì)算方法以有限差分法和復(fù)變量差分法為主,在梯度優(yōu)化中,都需要對控制方程反復(fù)求解,計(jì)算量大。而伴隨理論算法可以在保持精度的前提下,一次性求解出目標(biāo)函數(shù)對所有設(shè)計(jì)變量的導(dǎo)數(shù)[3]。其計(jì)算量僅與目標(biāo)函數(shù)和約束條件的數(shù)量有關(guān),與設(shè)計(jì)變量的數(shù)量無關(guān),適合多設(shè)計(jì)變量的優(yōu)化設(shè)計(jì)問題。

本文的無人機(jī)氣動(dòng)彈性優(yōu)化以氣動(dòng)結(jié)構(gòu)耦合系統(tǒng)作為分析系統(tǒng),以氣動(dòng)參數(shù)和結(jié)構(gòu)參數(shù)作為設(shè)計(jì)變量。優(yōu)化目標(biāo)函數(shù)I和控制方程殘差R可設(shè)定為:

其中,x表示設(shè)計(jì)變量;w、u表示氣動(dòng)/結(jié)構(gòu)耦合系統(tǒng)狀態(tài)變量,其中w具體表現(xiàn)為渦格法中網(wǎng)格單元的氣動(dòng)負(fù)載,u具體表現(xiàn)為有限元分析中結(jié)構(gòu)有限元分析得到的結(jié)構(gòu)位移;Ra對應(yīng)氣動(dòng)分析的殘差,Rs對應(yīng)結(jié)構(gòu)分析的殘差。

令優(yōu)化目標(biāo)函數(shù)I和控制方程殘差R分別對設(shè)計(jì)變量x進(jìn)行求導(dǎo),并利用恒等變換和伴隨算子ψ得到伴隨方程。通過求解伴隨方程,可以一次性求解出目標(biāo)函數(shù)針對所有設(shè)計(jì)變量的導(dǎo)數(shù),如式(8)所示:

2.2 序列二次規(guī)劃

常用的梯度優(yōu)化方法包括BFGS 擬牛頓算法、共軛梯度法、序列二次規(guī)劃算法(Sequential Quadratic Programming,SQP)等。其中,序列二次規(guī)劃法對大規(guī)模數(shù)據(jù)的計(jì)算效率較高,在工程中得到了廣泛應(yīng)用。具體形式如式(9)所示:

其中,f(x)為目標(biāo)函數(shù),h(x)為等式約束函數(shù),g(x)為不等式約束函數(shù)。

序列二次規(guī)劃法是一種用于求解非線性最優(yōu)化問題的方法,其基本思路是將原問題轉(zhuǎn)化為一系列二次規(guī)劃(Quadratic Programming,QP)問題進(jìn)行迭代求解。在每次迭代中,通過求解二次規(guī)劃問題得到下一次迭代的值,直到序列{xk}收斂于極值點(diǎn)。具體來說,對于某一次迭代值xk,NLP 會(huì)被近似為該點(diǎn)處的二次規(guī)劃問題。通過求解該二次規(guī)劃問題,可以得到下一次迭代的值xk+1。隨著迭代次數(shù)k的增加,序列{xk}將收斂于極值點(diǎn)[25]。本文采用SLSQP 算法(Sequential Least SQuares Programming),該算法是由Kraft[26]于1988 年提出的一種非線性約束優(yōu)化算法,通過將優(yōu)化問題轉(zhuǎn)化為一系列線性或二次規(guī)劃子問題,用于求解無約束或約束的非線性優(yōu)化問題。具體來說,SLSQP 算法采用牛頓法求解每個(gè)子問題,并利用輔助函數(shù)法處理約束條件。

2.3 氣動(dòng)彈性優(yōu)化框架

本文所采用的氣動(dòng)彈性優(yōu)化框架如圖4 所示。具體流程如下:①設(shè)定優(yōu)化目標(biāo)并設(shè)置目標(biāo)參數(shù),明確目標(biāo)函數(shù)、約束函數(shù)、設(shè)計(jì)變量與狀態(tài)變量等相關(guān)變量及數(shù)據(jù);②根據(jù)初始數(shù)據(jù)建立氣動(dòng)和結(jié)構(gòu)的物理模型;③對目標(biāo)函數(shù)和約束函數(shù)關(guān)于所有設(shè)計(jì)變量的總梯度進(jìn)行計(jì)算;④利用序列二次規(guī)劃算法進(jìn)行尋優(yōu);⑤通過不斷更新氣動(dòng)和結(jié)構(gòu)的物理模型,直至得到收斂解,完成優(yōu)化設(shè)計(jì)。

Fig.4 Basic process of optimization design based on gradient optimization algorithm圖4 基于梯度優(yōu)化算法的優(yōu)化設(shè)計(jì)基本流程

3 無人機(jī)氣動(dòng)彈性優(yōu)化

3.1 大過載無人機(jī)

在有人機(jī)上,8g是受過專業(yè)訓(xùn)練的航天員一般可以承受的最大過載[27]。而無人機(jī)沒有相關(guān)限制,可以做出更多大過載的機(jī)動(dòng)動(dòng)作來規(guī)避潛在危險(xiǎn)。

3.1.1 工況計(jì)算

大過載無人機(jī)具體的飛行工況如表1 所示。選取NACA0015 翼型,機(jī)翼初始參數(shù)如表2 所示,并選取7075 鋁合金作為翼梁結(jié)構(gòu)材料。

Table 1 High overload UAV flight state parameters表1 大過載無人機(jī)飛行狀態(tài)參數(shù)

Table 2 Initial shape parameters of high overload UAV wing表2 大過載無人機(jī)機(jī)翼初始外形參數(shù)

目標(biāo)函數(shù)為巡航狀態(tài)下的燃油消耗量,對于飛機(jī)燃油消耗,采用Breguet 距離方程進(jìn)行計(jì)算。具體方程形式如下:

其中,Wf為燃油重量,W0為飛機(jī)空重,Ws為結(jié)構(gòu)重量,R為最大航程,V為巡航速度,SFC為燃油消耗率。設(shè)計(jì)變量為機(jī)翼梢根比、機(jī)翼后掠角、機(jī)翼扭轉(zhuǎn)角、機(jī)翼翼梁厚度、機(jī)動(dòng)狀態(tài)下的機(jī)翼迎角。設(shè)計(jì)變量具體數(shù)量及數(shù)值范圍如表3 所示。約束條件為:機(jī)動(dòng)狀態(tài)下機(jī)翼結(jié)構(gòu)應(yīng)力接近但不超過應(yīng)力極限、機(jī)動(dòng)狀態(tài)下升力等于重力。

Table 3 Parameter range and number of design variables for high overload UAV表3 大過載無人機(jī)設(shè)計(jì)變量參數(shù)范圍及數(shù)量

3.1.2 結(jié)果優(yōu)化

經(jīng)過優(yōu)化設(shè)計(jì)后,可得到如表4 所示的3 種優(yōu)化設(shè)計(jì)結(jié)果。無人機(jī)機(jī)翼的外形經(jīng)過優(yōu)化設(shè)計(jì)后的變化如圖5所示。

Table 4 Optimal design results of high overload UAV表4 大過載無人機(jī)優(yōu)化設(shè)計(jì)結(jié)果

Fig.5 Optimal design results of UAV with different wing root chord lengths圖5 不同機(jī)翼翼根弦長無人機(jī)優(yōu)化設(shè)計(jì)結(jié)果

通過對不同機(jī)翼翼根弦長情況下機(jī)翼優(yōu)化設(shè)計(jì)結(jié)果的綜合分析,可以看出,在機(jī)翼經(jīng)過優(yōu)化設(shè)計(jì)后,大過載無人機(jī)的燃油消耗量相較于初始狀態(tài)均降低了20%~30%。以機(jī)翼翼根弦長5.5 m 為例,如圖6 所示,在經(jīng)過優(yōu)化設(shè)計(jì)后,可以看到馮米塞斯應(yīng)力在巡航狀態(tài)下遠(yuǎn)低于應(yīng)力極限,符合基本邏輯和設(shè)計(jì)要求。在扭轉(zhuǎn)角方面,機(jī)翼翼尖和翼根的相對扭轉(zhuǎn)角為2.89°,符合機(jī)翼外形設(shè)計(jì)邏輯,改善了升力分布,有利于提升飛機(jī)的巡航性能[28]。設(shè)計(jì)方案總體上較為合理,符合飛行器機(jī)翼外形的設(shè)計(jì)邏輯。

Fig.6 Optimal design results of wing with root chord length of 5.5m圖6 5.5m機(jī)翼翼根弦長優(yōu)化設(shè)計(jì)結(jié)果

3.2 長航時(shí)無人機(jī)

3.2.1 工況計(jì)算

對于長航時(shí)無人機(jī),具體的飛行工況如表5 所示。選取NACA0015 翼型,采用如圖7 所示的矩形機(jī)翼,初始參數(shù)如表6所示,并選取7075鋁合金作為翼梁結(jié)構(gòu)材料。

Table 5 Flight state parameters of long-endurance UAV表5 長航時(shí)無人機(jī)飛行狀態(tài)參數(shù)

Table 6 Initial parameters of long-endurance UAV wing表6 長航時(shí)無人機(jī)機(jī)翼初始參數(shù)

目標(biāo)函數(shù)為燃油消耗量,設(shè)計(jì)變量為機(jī)翼梢根比、機(jī)翼后掠角、機(jī)翼扭轉(zhuǎn)角、機(jī)翼翼梁厚度、巡航狀態(tài)下的機(jī)翼迎角。設(shè)計(jì)變量的具體數(shù)量及數(shù)值范圍如表7 所示。約束條件為:機(jī)翼結(jié)構(gòu)應(yīng)力不超過應(yīng)力極限、巡航狀態(tài)下升力等于重力。

Table 7 Parameter range and number of high overload UAV flight status表7 高過載無人機(jī)飛行狀態(tài)參數(shù)范圍及數(shù)量

3.2.2 結(jié)果優(yōu)化

長航時(shí)無人機(jī)優(yōu)化后的機(jī)翼幾何外形如圖8所示。

Fig.8 Optimized wing geometry of long-endurance UAV圖8 長航時(shí)無人機(jī)優(yōu)化后的機(jī)翼幾何外形

初始狀態(tài)下的燃油消耗量為3 898 kg,長航時(shí)優(yōu)化后的燃油消耗量為3 004 kg,優(yōu)化后的結(jié)構(gòu)重量為1 761 kg。與初始狀態(tài)相比,燃油消耗量降低了22.93%。優(yōu)化設(shè)計(jì)后的具體參數(shù)如表8 所示。長航時(shí)情況下的優(yōu)化設(shè)計(jì)結(jié)果如圖9 所示。在扭轉(zhuǎn)角方面,機(jī)翼翼尖和翼根的相對扭轉(zhuǎn)角為5.15°,處于合理范圍內(nèi)。在馮米塞斯應(yīng)力方面,由于需要獲得最佳的優(yōu)化設(shè)計(jì)結(jié)果,因此需要在設(shè)計(jì)過程中盡可能逼近結(jié)構(gòu)極限。本文的優(yōu)化設(shè)計(jì)方案未超出應(yīng)力極限,符合飛行器外形設(shè)計(jì)邏輯。

Table 8 Optimized design results of long endurance UAV表8 長航時(shí)無人機(jī)優(yōu)化設(shè)計(jì)結(jié)果

Fig.9 Optimized design results of long endurance UAV圖9 長航時(shí)情況下的優(yōu)化設(shè)計(jì)結(jié)果

4 結(jié)語

本文提出一種現(xiàn)代無人機(jī)的高效氣動(dòng)彈性優(yōu)化方法,并通過大過載和長航時(shí)兩類典型無人機(jī)進(jìn)行了驗(yàn)證,具體如下:

(1)提出基于伴隨理論的氣動(dòng)彈性優(yōu)化方法,通過氣動(dòng)渦格法與有限元分析法進(jìn)行氣動(dòng)彈性松耦合分析,并基于伴隨理論和序列二次規(guī)劃法進(jìn)行非線性敏度優(yōu)化。

(2)針對大過載無人機(jī)進(jìn)行結(jié)構(gòu)優(yōu)化設(shè)計(jì),以燃油消耗量為目標(biāo)函數(shù),以機(jī)翼外形參數(shù)和機(jī)翼翼梁結(jié)構(gòu)參數(shù)為設(shè)計(jì)變量,結(jié)合機(jī)動(dòng)狀態(tài)下氣動(dòng)、結(jié)構(gòu)兩方面的約束函數(shù),使無人機(jī)的燃油消耗量分別減少了21.66%、27.93%、32.05%。

(3)針對長航時(shí)無人機(jī)進(jìn)行結(jié)構(gòu)優(yōu)化設(shè)計(jì),以燃油消耗量為目標(biāo)函數(shù),以機(jī)翼外形參數(shù)和機(jī)翼翼梁結(jié)構(gòu)參數(shù)為設(shè)計(jì)變量,結(jié)合巡航狀態(tài)下氣動(dòng)、結(jié)構(gòu)兩方面的約束函數(shù),使無人機(jī)的燃油消耗量減少了22.93%。

綜上所述,本文提出的方法高效、穩(wěn)健、可靠,時(shí)間成本和學(xué)習(xí)成本低,可為無人機(jī)氣動(dòng)彈性優(yōu)化提供理論和方法參考。

猜你喜歡
有限元優(yōu)化結(jié)構(gòu)
超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
論結(jié)構(gòu)
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結(jié)構(gòu)
創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 伊伊人成亚洲综合人网7777| 欧美日本在线一区二区三区| 搞黄网站免费观看| 成人午夜天| 欧美日韩成人| 91免费国产高清观看| 欧美精品在线看| 成人在线观看不卡| 国产噜噜噜| 亚洲精品午夜天堂网页| 人妻无码中文字幕第一区| 国禁国产you女视频网站| 中文一级毛片| 欧美性精品| 国产成人亚洲精品色欲AV| 欧美全免费aaaaaa特黄在线| 国产在线视频自拍| 亚洲色图在线观看| 夜夜操天天摸| 国产噜噜噜视频在线观看 | 午夜精品福利影院| 亚洲成人手机在线| 色综合五月婷婷| 中文字幕亚洲专区第19页| 国产一级在线播放| 亚洲一级毛片在线播放| 日韩免费成人| 天天干天天色综合网| 久久综合色播五月男人的天堂| 欧美视频在线播放观看免费福利资源 | 亚洲综合国产一区二区三区| 91色爱欧美精品www| 在线精品亚洲国产| 婷婷六月天激情| 免费观看亚洲人成网站| 国产成人91精品免费网址在线| 国产免费精彩视频| 九色在线观看视频| 日韩精品一区二区三区swag| 国产精品3p视频| 制服丝袜国产精品| jizz国产在线| 狠狠ⅴ日韩v欧美v天堂| 91尤物国产尤物福利在线| AV不卡在线永久免费观看| a色毛片免费视频| 成人免费一级片| 呦女精品网站| 五月天久久综合| 午夜爽爽视频| 亚洲国语自产一区第二页| 日本日韩欧美| 国产精品偷伦在线观看| 人人91人人澡人人妻人人爽| 无码日韩人妻精品久久蜜桃| 国产成人高清精品免费5388| 国产精品视频系列专区| 日本一区二区不卡视频| 71pao成人国产永久免费视频| 亚洲一级色| 亚洲日本中文字幕乱码中文| 婷婷亚洲最大| 国产亚洲欧美日韩在线一区二区三区| 欧美性久久久久| 免费毛片网站在线观看| 国产在线视频二区| 国产成人调教在线视频| 国产欧美高清| 国产精品自在在线午夜区app| 无码精品一区二区久久久| 国产网站免费观看| av天堂最新版在线| 亚洲男人天堂网址| 亚洲黄色视频在线观看一区| 69av免费视频| 高清国产在线| 欧美色丁香| 国产女人综合久久精品视| 精品国产成人国产在线| 久久精品人人做人人爽97| 在线亚洲小视频| 天天色综合4|