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

通用航空飛行器的參數化建模與多目標設計優化

2022-12-29 09:54:38徐小平張天天
航空兵器 2022年6期
關鍵詞:優化方法設計

徐小平,張天天,沈 洋

(1.國防科技大學 北京學科交叉中心,北京 100101;2.航天工程大學 宇航科學與技術系,北京 101400;3.國防科技大學 空天科學學院,長沙 410073)

0 引 言

通用航空飛行器(Common Aero Vehicle,CAV)是一種采用通用導航、制導與控制系統以及通用空氣動力學頭罩的機動再入飛行器,能夠裝載各種類型的戰斗部、傳感器或情報、監視和偵察系統,并且能與多種發射系統匹配[1]。由CAV衍生出的高超聲速技術飛行器(Hypersonic Technology Vehicle, HTV)項目開展了多次飛行試驗,具有廣闊的發展前景[2]。早期,波音公司提出的基于雙錐構型的先進機動再入飛行器(Advanced Maneuvering Reentry Vehicle,AMaRV)是一種能夠直接用于發展CAV項目的飛行演示系統。同時,洛克希德·馬丁公司也發展了兩個與CAV直接相關的項目,其中高性能機動再入飛行器(High Performance Maneuvering Reentry Vehicle,HPMaRV)系統具有詳細的仿真和風洞試驗數據。作為再入飛行器,采用雙錐構型的AMaRV升阻比較小,被稱為低性能通用航空飛行器(Low-Performance Common Aero Vehicle,CAV-L),而采用升力體構型的HPMaRV升阻比較大,被稱為高性能通用航空飛行器(High-Performance Common Aero Vehicle,CAV-H)。盡管如此,CAV-H構型的尖前緣在帶來大升阻比的同時,也對防熱系統提出巨大挑戰。使CAV在高超聲速再入環境下具有良好的空氣動力學性能且不被燒蝕成為一項重要課題[3]。

無論是升阻比還是加熱環境均與CAV的構型密切相關,通過構型的設計和優化有望在任務要求和材料受熱極限的約束下發現可行方案。而對CAV構型的優化設計需要具備高保真度的參數化建模方法、性能評估方法和數值優化方法。其中,參數化建模方法使設計者能夠將復雜的三維外形用簡單的數學參數進行表示,實現對外形的精細控制[4-5];性能評估方法能夠在已知構型的基礎上通過計算或者實驗獲得目標特性,如升阻比、容積率、熱流峰值等;數值優化方法能夠通過尋優,找到使目標特性最符合預期的輸入參數設置。

氣動外形的設計與優化一直是飛行器設計領域的熱點,大量的新式翼型、飛行器構型在優化中產生,使飛行器設計水平顯著提高。其中,三維類別/形狀轉換(Class/Shape Function Transformation,CST)建模[6]、三維B樣條[7-8]和自由網格變形方法(Free-Form Deformation,FFD)[9]作為應用最為廣泛的參數化建模方法在氣動優化過程中得到廣泛應用。Chen等[10]在跨聲速機翼的流固耦合穩定性優化中采用CST參數化建模方法。Xia等[11]采用FFD參數化建模方法,結合改進的粒子群優化算法和徑向基函數對二維翼型和運載火箭進行氣動外形優化。李治宇等[12]將構型參數化、網格自動化及流場快速求解程序進行有機結合,發展了快速有效的氣動布局優化設計方法, 并設計了新型高升阻比融合升力體氣動布局。

考慮到單純提高升阻比可能會給飛行器的結構、容積和防熱等性能帶來新的問題,在氣動優化過程中往往還要同時考慮其他設計目標或約束。Zhang等[13]采用多目標設計優化的方法研究了高超聲速滑翔飛行器升阻比與容積率之間的關系,發現兩者難以同時達到最優。Deng等[14]在進行帶翼高超聲速飛行器的升阻比優化時,考慮了飛行器的容積率和操穩特性等條件帶來的約束。張龍龍等[15]基于CST參數化建模方法對翼身融合升力體構型進行參數化設計和多目標優化。Zuhal等[16]基于代理模型對翼型進行多目標氣動外形優化,比較了不同加點方法對優化效率提升的影響。

盡管氣動優化理論已經得到廣泛關注和蓬勃發展,然而類似CAV這樣的三維曲面構型,其參數化建模的精度與參數數量之間的矛盾依然存在,高超聲速條件下其升阻比與熱流峰值之間的矛盾仍然難以調和,優化設計精度和效率還需進一步提升。本文提出了針對不同類型CAV構型需靈活選擇三維參數化方法的理念,為CAV的高效參數化建模提供參考;通過多目標設計優化,獲得以CAV構型升阻比、容積率和峰值熱流為目標的三維非劣解前緣,并對其進行深入分析,提高了CAV構型設計水平。

1 CAV構型參數化建模方法

1.1 傳統參數建模

已知CAV-L為雙錐構型,其長度處于107~144 in之間,尾部直徑為36 in。此類旋轉體構型只需對二維輪廓進行參數化設計,并結合雙錐構型的特點,即可實現用少量參數實現三維構型的參數化建模,如圖1所示。

圖1 CAV-L構型設計參數示意圖

可見,對于此類簡單的雙錐構型,只需5個參數即可確定其旋轉體外形,分別為頭部鈍化半徑R,一級頭錐角度θ,一級頭錐長度Ln,總長度L和尾部直徑D,其他特征參數如二級圓臺長度、錐角等均可由上述參數的確定而得到。由于尾部直徑D給定,鈍化半徑可根據防熱需求另行確定且對構型整體形狀和性能影響較小,因此,該構型的實際有效控制變量分別為L,Ln和θ。依據公開資料,設定L=144 in,Ln=50 in,θ=12°。

同時,給定鈍化半徑R=2 in,則可得CAV-L的三維構型。波音公司設計的實驗構型是在旋轉體構型的基礎上通過安裝氣動控制舵面以及穩定裝置來滿足飛行需求。本節則通過削平底面并切割的方式得到控制舵,并在構型兩側對稱安裝固定控制面,得到的三維概念構型如圖2所示。該構型與公開資料中的CAV-L構型高度相似。

圖2 CAV-L構型3D模型概念圖

1.2 三維CST曲面參數化建模方法

CAV-H為升力體構型,雖然具有面對稱的特點,但其上下表面均為曲面,難以用傳統的參數化方法進行表示,因此,需要采用新型三維曲面參數化方法進行建模,并為構型優化打下基礎。

類別/形狀轉換(CST)參數化方法是一種兩層曲線參數化方法。該方法的出發點是建立一個具有目標曲線基礎形狀的類函數,然后建立一個形狀函數對類函數進行進一步修飾。CST方法的表達式為

(1)

在二維CST方法的基礎上,首先建立橫、縱坐標之間的CST函數關系,然后引入第三維坐標,建立其與橫、縱坐標之間的三維關系,從而能夠對三維曲面進行參數化建模。三維參數化方法的表達式為

(2)

ΔεM, N(φ,η)

(3)

式(3)建立了Z方向坐標與X和Y方向坐標之間的關系,其中X方向坐標歸一化后用φ表示,Y方向坐標歸一化后用η表示,兩者的類函數和形狀函數組合構成了曲面表達式;nx和ny分別為兩個Bernstein方程的階數。由類函數和形狀函數表達式可知,在X和Y分別取極值時類函數和形狀函數均等于0,即曲面邊線均處于Z=0的平面上,為此,三維方程引入ΔεM,N(φ,η)來定義曲面邊緣所處的空間位置。ΔεM,N(φ,η)定義了曲面邊緣在X-Z平面和Y-Z平面上投影的形狀,具有多種形式,仍然可以采用標準二維CST方法表示,即

(4)

用該三維CST方法設計的類CAV-H構型如圖3所示。該方法的詳細介紹可參見文獻[13]。

圖3 采用三維CST參數化方法設計的CAV-H構型

1.3 FFD自由網格變形參數化建模方法

FFD自由網格變形方法通過Bernstein函數,將構型表面所有點的空間坐標映射到控制網格的節點坐標,從而實現目標的參數化。Bernstein函數在[0,1]空間內的分布如圖4所示,其中黑色虛線代表各項和的分布,為一條恒等于1的直線,即Bernstein系數均取1時所有輸入點得到的輸出均為1;此外,當各項系數均取對應項編號與階數的比值(例如第i項系數取i/l,l為Bernstein函數的階數)時,其所得輸出等于輸入值,如圖4中的紅色斜線所示。通過調整各系數的取值,即可對應調整函數輸出,若將歸一化的控制網格節點坐標直接作為Bernstein系數,則當其均勻分布時,構型形狀不變,而當控制節點坐標改變時,則會對應改變Bernstein系數,從而改變構型形狀。

圖4 Bernstein函數各項及其組合項在[0, 1]空間的分布

綜上所述,可引出FFD方法中控制節點坐標Qi, j, k對設計對象任一點坐標(x,y,z)的控制函數Pffd:

(5)

(6)

式中:i=0, 1, …,l;j=0, 1, …,m;k=0, 1, …,n;l,m,n為x,y,z三個方向上的Bernstein函數階數,為便于計算,將控制點坐標和控制體坐標均認為是歸一化后的坐標。

2 性能計算方法

2.1 高超聲速氣動性能工程估算方法

高超聲速工程估算方法是一種快速估算構型高超聲速氣動性能的方法。該方法認為在高超聲速條件下,物體迎風表面的壓力系數Cp與表面法線和來流之間夾角正弦值的二次方正相關,計算方程如下:

Cp=Cpmaxsin2δ

(7)

式中:Cpmax為相關系數,該值可采用正激波駐點后的壓力系數值。

(8)

式中:p0,p∞,ρ∞和V∞分別為自由來流的總壓、靜壓、密度和速度。由于

(9)

(10)

從而得到

(11)

對于飛行器背風面則可采用Prandtl-Meyer膨脹波理論進行求解。已知當來流馬赫數大于1時,經過鈍角物面的來流速度增加、壓力降低,形成膨脹波,Prandtl-Meyer膨脹波理論則定義了偏轉角度與來流馬赫數之間的函數關系,即

(12)

根據壓力系數與流動轉折角之間的關系:

(13)

則可計算得到背風面物面的壓力系數。

由于飛行器為連續曲面構型,而高超聲速工程估算方法反映的是局部來流迎角與當地壓力系數之間的關系,因此,為計算構型整體受力,需要采用面元法將構型表面分成多個微小面元,分別計算壓力系數分布后求和得到構型整體受力。

面元的劃分有多種方式,最常見的是三角面元和四邊形面元劃分方法,如圖5所示。

圖5 兩種面元劃分方式及其與來流之間的關系

面元劃分相當于將物體表面離散化,對每個面元的氣動受力求解前需要先求出當前面元的單位法向量ni和面元面積si。面元法向量用于求解面元迎風角,面積用于計算面元所受來流壓力。

對于四邊形面元,面元法向量求解方式為

(14)

四邊形面元的面積有多種求法,其中一種為

(15)

四邊形相鄰面元之間的拓撲關系相比于三角面元更加具有規律性,因此,采用四邊形面元劃分方法,如圖6所示。

圖6 CAV-L表面的四邊形網格劃分

2.2 容積率計算方法

容積率(ηv)是反映飛行器有效載荷裝填能力的參數,定義如下:

(16)

式中:V和S分別代表飛行器的體積和表面積。容積率越大,表明同樣容積的飛行器更適合用于裝填有效載荷。換言之,同等表面積的飛行器,容積率越大,其內部容積也越大。例如,長方體的容積率小于同等體積的正方體,同樣表面積的細長體容納有效載荷的能力不如球體。

2.3 熱流估算方法

高超聲速條件下的氣動加熱現象不可忽略,為避免飛行器被燒蝕,需要通過構型設計使飛行器表面加熱約束在材料耐受范圍以內。已知飛行器的駐點處于弓形激波正后方,其氣動加熱現象最為顯著,因此,采用工程算法估算駐點熱流量qs,即

(17)

式中:Rs為駐點處的曲率半徑;ρ0為海平面大氣密度;C為固定常數,取值為2.146×10-4J·m-4.5·s2[17]。可見在來流條件一定的前提下,駐點曲率越大,熱流量越大。

3 優化方法

3.1 基于遺傳算法的單目標優化方法

遺傳算法通過模擬自然界中“優勝劣汰”的生存法則來使優化問題收斂到最優結果附近。算法首先在參數設計空間內隨機抽取一部分樣本點,稱為種群;然后計算種群內每個個體的響應結果(即適應度),并根據適應度大小對種群進行排序,在此基礎上產生下一組種群并按此過程迭代循環直到滿足迭代終止條件。在此過程中用于產生下一組種群(子代)的當前種群被稱為親代;在親代中選擇優勝個體并淘汰劣勢個體的過程稱為選擇;由優勝個體產生子代的過程稱為遺傳。遺傳過程有三種形式:第一種是部分優勝個體(精英)直接進入子代;第二種是優勝個體間通過交換步驟將彼此的信息切斷后互相再組合生成新的個體(又稱交叉);第三種則是在親代信息基礎上隨機引入變異。迭代終止條件包括最大迭代代數、迭代時間、適應度函數值迭代殘差值等。可見,遺傳算法中的種群規模、交叉系數、變異概率和收斂條件等均會對收斂結果產生影響。

3.2 基于NSGA-II的多目標優化方法

當優化問題出現多個目標時則需要采用多目標優化算法,例如本文中不僅希望CAV構型的升阻比大,還希望其容積率大、最大熱流量小。由于不同目標之間可能相互排斥,即一個目標向好的方向發展時,另外的目標反而變得更壞,此時需要引入非劣解的概念。非劣解表示不存在所有目標均優于該解的任何其他解。多目標優化的目的則是在約束范圍內收斂到非劣解集。NSGA-II算法是基于遺傳算法發展而來的多目標優化算法,該算法中親代生成子代的過程仍然按照遺傳算法所需的選擇、交叉、變異等操作。在產生子代后要將親代和子代的結果按照是否為非支配解進行排序,將所有當前的非劣解集篩選出記為Rank0,然后將剔除Rank0后的剩余結果再次排序,篩選出剩余解中的非劣解集并記為Rank1,以此類推。可見Rank值越小表明當前解集中解的質量越高。在標記結束后按Rank值從小到大選擇子代,當達到某一Rank值時出現添加前子代數量小于預定值而添加后子代數量多于預定值,則在該Rank解集中進行擁擠度排序,并按照擁擠度從小到大的順序補足子代所需的個體數。按照此過程循環迭代,直到滿足收斂條件,最終可得到多輪迭代后的非劣解集。在該解集的基礎上,設計者可以按需求任意選擇可行解,所選的解盡管不能保證滿足所有目標最優,但不會存在全面優于該解的其他解。

4 結果與分析

4.1 優化問題定義

根據航程公式,升阻比越大的武器打擊范圍也越大,因此,本文首先對CAV-L構型和CAV-H構型開展單目標優化設計,優化目標是攻角[-5°, 30°]范圍內的最大高超聲速升阻比。其中,CAV-L構型采用傳統參數化建模方法,優化變量包括L,Ln和θ,優化問題表述為

maxL/Dmax

s.t. 107 in≤L≤144 in

15 in≤Ln≤60 in

5°≤θ≤20°

(18)

CAV-H構型采用三維CST參數化方法,優化變量包括Hl,Hu,N1u,M1u,M1l,T1u,T1l,其初始值的集合用V0表示,V0=[0.1,0.25,0.8,1,1,1.5,1],優化過程中優化變量取值上限為1.5V0,取值下限為0.5V0,本次優化不考慮其他約束。

CAV-L和CAV-H構型的氣動性能計算均在馬赫數15、高度25 km的條件下,優化算法均采用遺傳算法,最大迭代代數為100,每代交叉概率為0.8,變異概率為0.2,收斂精度為10-6。由于CAV-L的優化變量較少,故種群規模為40;CAV-H的優化變量較多,種群規模為50。

由于構型容積率和駐點區熱流密度同樣是衡量飛行器性能的重要指標,因此,進一步在CAV-L構型優化的基礎上,開展以高超聲速條件下升阻比最大、構型容積率最大、駐點區熱流密度最小為優化目標的多目標設計優化。采用FFD參數化建模方法,以FFD控制點位置作為設計變量,將所有優化目標轉化為成本型指標,產生的優化問題描述如下:

min-L/Dmax,-ηv,qs

s.t.DVmin≤DV≤DVmax

(19)

多目標優化采用NSGA-II優化算法,設置算法中的交叉概率為0.7、變異概率為0.4、種群規模為300、最大迭代次數為40次。

4.2 以升阻比為目標的單目標優化設計結果

4.2.1 CAV-L構型優化結果

CAV-L原始構型最大升阻比為2.36,處于攻角為10°的飛行條件下。

優化過程中種群最優升阻比和平均升阻比的變化曲線如圖7所示,迭代100代后優化過程收斂,種群升阻比收斂到約2.69。

圖7 CAV-L升阻比優化收斂曲線

圖8所示為優化前后構型的升阻比曲線對比,優化后的最大升阻比出現在更小的攻角條件下,升阻比大小提高了12.7%。優化后的設計參數取值為L=144 in,Ln=40 in,θ=8.66°。

圖8 CAV-L構型優化前后高超聲速升阻比曲線對比

優化結果表明,總長度L的增加是造成最大升阻比增大的主要原因。該參數組合控制下的CAV-L構型如圖9所示,相比于圖2,優化后的構型更加細長。

圖9 優化后的CAV-L構型

4.2.2 CAV-H構型優化結果

CAV-H構型的優化過程在達到最大迭代次數時停止,升阻比最優值與平均值趨于相等且收斂,升阻比最終收斂值為9.824,如圖10所示。

圖10 CAV-H升阻比優化收斂曲線

根據優化結果,升阻比達到最優值時的設計變量Vopt=[0.05,0.125,0.4,1.088 7,1.09,2.245,1.193 9],前3個設計變量均趨于設計范圍的最小值,第6個設計變量趨于設計范圍最大值,其余3個設計變量處于最大值和最小值之間,即最大升阻比的提高主要是通過Hl,Hu,N1u的減小以及T1l的增大實現的。根據以上設計變量所得類CAV-H構型如圖11所示。

圖11 優化后的CAV-H構型

圖11所示的飛行器構型具有較大的升阻比,其整體呈現扁平狀,相比于圖3展示的參考構型,其厚度更薄,頭部更寬。兩者的升阻比曲線對比如圖12所示。從圖11~12可以看出,優化后的CAV-H構型過于扁平,很難在實際再入過程中承受氣動加熱,且其扁平的形狀也對內部部件的安裝提出巨大挑戰。同時,盡管其最大升阻比較大,但隨著攻角偏離最大升阻比出現的條件(2°),其升阻比值將快速下降,使得實際飛行過程的魯棒性較差,抗干擾能力較弱。

圖12 CAV-H構型優化前后高超聲速升阻比曲線對比

綜上所述,僅以升阻比為構型優化目標難以得到符合工程實際需求的構型,需要采用多目標優化的方法對多個設計目標變量進行優化。

4.3 基于NSGA-II的CAV多目標優化結果

基于CAV-L單目標優化構型所構建的幾何外形FFD參數化設計如圖13所示,考慮幾何約束設計要求的面對稱性和底部平整性,在變形過程中令控制點保持xoz平面對稱、底部位置的控制點在x方向上不移動。

圖13 初始FFD參數化設計

多目標優化所得非劣解集(又稱Pareto前緣)如圖14所示。優化結果表明,相比于統計得到的Pareto前緣解集,初始外形(2.692 8,0.581 1,6.748 5)所處位置是明顯的離群點,且存在較多解方案在3個指標均優于初始外形,說明算法驅動設計空間尋找到了綜合性能良好的優化解集。根據Pareto前緣解集的散點分布情況來看,3個優化目標之中升阻比指標在數值上的提升較為明顯,較強的負相關性存在于升阻比與容積率之間、升阻比與熱流密度(效益)之間,說明在當前設計方法下,位于Pareto前緣的上述兩組指標對存在著較強的矛盾。由于最大升阻比大的構型一般為細長體,而細長體構型的容積率往往會比較小,另外熱流密度的減小往往需要增大構型最小曲率半徑,從而增大阻力,造成3個優化目標無法同時達到最優解。

圖14 CAV構型多目標優化所得Pareto前緣解集

各優化目標上的優秀外形、指標性能及FFD設計布局如圖15所示。Pareto解集中的最大升阻比能夠達到2.973,分析升阻比最優的外形,FFD布局直觀展示了相對于初始外形的變形特征,根據FFD控制點的位移,其優化邏輯是令頭部附近靠近下方的控制點前移,通過拉伸飛行器機身下前部以擴大升力面,同時令頭部附近側面控制點向內收縮,以減小外形阻力。

圖15 各優化目標上的優秀外形及其FFD設計

容積率最優的外形是初始外形的支配解(即3個指標均更優),其優化邏輯是令頭部附近控制點向后收縮,減小外形較細(橫截面面積較小)幾何區域的表面積,從而提高體積效率;另外,令機身區域的控制點變寬,從而直接增加飛行器機身區域的體積。但上述現象與升阻比的優化邏輯顯然存在矛盾。

分析針對熱流密度目標的優化方式,該目標主要與局部外形特征相關。FFD結構顯示其頭部附近控制點形成“凹陷狀”布局,使外形頭部幾何沿流線方向鈍化,該行為與升阻比優化邏輯存在一定的沖突。但若對頭部外形進行更加細致的FFD參數化設計,分別增加對駐點區與升力面的變形自由度,則有可能緩解當前優化水平的升阻比-熱流密度指標對的沖突現象。

綜上所述,在設計范圍內,相較于初始外形,得到的Pareto解集升阻比指標中位數提升5.2%、容積率中位數提升1.4%、熱流密度中位數減小6.6%,該優化過程達到了預期目的,Pareto前緣上的點可為工程中具體構型的選擇提供參考依據。在后續優化中,基于FFD參數化設計,主要考慮進一步細化對外形頭部變形自由度的構造,局部弱化升阻比與其他目標之間的相關性。

5 結 論

本文針對兩類CAV構型研究了三維參數化建模方法在航空飛行器曲面建模上的應用前景,指出CAV-L和CAV-H的參數化建模流程,開展了以高超聲速條件下升阻比最大為目標的構型優化,探索了以升阻比最大、容積率最大以及熱流量最小為目標的構型多目標設計優化方法,并得出以下結論:

(1)需根據構型設計要求和典型特征選擇合適的三維參數化建模方法,為提高優化效率和結果可靠性,應盡量減少建模參數數量;

(2)最大升阻比大的CAV構型往往比較細長,其升阻比相對攻角的魯棒性變差;

(3)最大升阻比、容積率和最大熱流量三個指標無法同時達到最優,僅以升阻比為目標的優化設計必然犧牲其他性能指標。

本文存在的不足之處包括:尚未進一步研究如何剔除對設計目標影響較小的設計變量,以進一步提高優化效率;采用的工程估算方法相對粗糙,仍需進一步完善;多目標優化結果尚待進一步挖掘分析。

猜你喜歡
優化方法設計
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
主站蜘蛛池模板: 综合社区亚洲熟妇p| 国产亚洲美日韩AV中文字幕无码成人| 国产无码高清视频不卡| 婷婷五月在线视频| 免费无码AV片在线观看中文| 97精品伊人久久大香线蕉| 韩国福利一区| 欧美日韩中文字幕二区三区| 超碰91免费人妻| 91精品国产自产在线观看| 亚洲av成人无码网站在线观看| 国产一区二区三区精品久久呦| 色AV色 综合网站| 欧美特级AAAAAA视频免费观看| 久草热视频在线| 亚洲三级色| 欧美h在线观看| 亚洲天堂2014| 国产浮力第一页永久地址| 91丝袜美腿高跟国产极品老师| 伊人久久大香线蕉aⅴ色| 爆操波多野结衣| 久久国产免费观看| 午夜国产理论| 国产精品高清国产三级囯产AV| 97在线免费| 成人福利在线看| 亚洲日韩精品欧美中文字幕 | 欧美α片免费观看| 欧美三级不卡在线观看视频| 99久久精品久久久久久婷婷| 99久久精品国产自免费| 欧美高清国产| 久久久久国产精品熟女影院| 欧美精品高清| 四虎成人在线视频| 亚洲精品无码av中文字幕| 午夜性爽视频男人的天堂| 国产精品网曝门免费视频| 伦精品一区二区三区视频| 国产午夜精品鲁丝片| 亚洲最大福利视频网| 四虎永久免费网站| 麻豆精品国产自产在线| 福利在线一区| 日韩无码黄色网站| 亚洲福利片无码最新在线播放| 日韩精品久久久久久久电影蜜臀| 538国产在线| 亚洲第一成人在线| 尤物视频一区| 日韩精品资源| 国产国产人免费视频成18| 久久精品日日躁夜夜躁欧美| 香蕉综合在线视频91| 最新亚洲人成无码网站欣赏网| 久久久久国产精品免费免费不卡| 亚洲综合婷婷激情| 欧美亚洲第一页| 亚洲国产精品人久久电影| 欧美成人一区午夜福利在线| 国产理论精品| 亚洲欧洲日韩综合| 亚洲一区二区三区麻豆| 国产18在线播放| 国产菊爆视频在线观看| 人妻夜夜爽天天爽| 国产尤物视频网址导航| 性网站在线观看| 亚洲成人免费看| 久久香蕉国产线看观| 国产麻豆va精品视频| 国产成a人片在线播放| 激情视频综合网| 亚洲国产成人综合精品2020| 国产自在自线午夜精品视频| 国产成人一级| 免费福利视频网站| 亚洲AV无码乱码在线观看裸奔| 色爽网免费视频| 日韩AV无码一区| 在线观看免费AV网|