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

高速列車(chē)受電弓艙氣動(dòng)外形多目標(biāo)優(yōu)化設(shè)計(jì)

2024-03-13 01:53:54朱隴輝武振鋒霍艷忠
關(guān)鍵詞:優(yōu)化模型

朱隴輝,武振鋒,周 琪,霍艷忠

(1.蘭州交通大學(xué)機(jī)電工程學(xué)院,蘭州 730070; 2.中國(guó)鐵路呼和浩特局集團(tuán)有限公司集寧機(jī)務(wù)段,內(nèi)蒙古烏蘭察布 012000)

當(dāng)列車(chē)高速運(yùn)行時(shí),其所引發(fā)的空氣動(dòng)力學(xué)問(wèn)題與低速時(shí)相比復(fù)雜程度提高[1]。對(duì)于明線運(yùn)行的列車(chē)而言,氣動(dòng)阻力和氣動(dòng)升力是影響其氣動(dòng)性能的兩個(gè)主要指標(biāo)[2]。高速列車(chē)的氣動(dòng)性能主要取決于其外部形狀,但列車(chē)的關(guān)鍵部位,如轉(zhuǎn)向架、受電弓導(dǎo)流艙以及風(fēng)擋等也會(huì)對(duì)列車(chē)的氣動(dòng)性能產(chǎn)生較大的影響[3-5]。對(duì)受電弓而言,其整體氣動(dòng)阻力過(guò)大會(huì)導(dǎo)致列車(chē)運(yùn)行總阻力增加[6]。若受電弓正向氣動(dòng)升力過(guò)大,則導(dǎo)致受電弓滑板與接觸線磨損嚴(yán)重,負(fù)向氣動(dòng)升力過(guò)大則會(huì)增加列車(chē)軸重,進(jìn)而使輪軌之間的磨損加劇,所以受電弓處的氣動(dòng)升力最理想的狀態(tài)是使其保持在0 N附近波動(dòng)[7-10]。

對(duì)于高速列車(chē)受電弓氣動(dòng)外形優(yōu)化設(shè)計(jì),已有大量學(xué)者進(jìn)行了相應(yīng)的研究。張亮[11]采用多目標(biāo)優(yōu)化理論,對(duì)受電弓安裝位置和受電弓導(dǎo)流罩進(jìn)行優(yōu)化分析,綜合考慮多個(gè)目標(biāo)的優(yōu)化程度對(duì)受電弓進(jìn)行外形優(yōu)化研究。肖友剛等[12]采用大渦模擬法和FW-H方程對(duì)受電弓絕緣子進(jìn)行外形優(yōu)化,結(jié)果表明橢圓形狀的絕緣子降低噪聲程度最大。燕永釗等[13]提出了一種新型的受電弓流線型滑板弓頭方案,得出該型方案在350 km/h和450 km/h速度下的阻力值,并于原模型在350 km/h的阻力值進(jìn)行比較,新型受電弓模型整體氣動(dòng)阻力在開(kāi)口和閉口狀態(tài)下分別降低25.13%與24.19%。肖守訥等[14]提出兩種受電弓導(dǎo)流罩形式,對(duì)其進(jìn)行噪聲研究,分別分析了低頻段和高頻段各導(dǎo)流罩噪聲分布,得出最優(yōu)結(jié)構(gòu)且得到了對(duì)未來(lái)導(dǎo)流罩的工程化運(yùn)用有一定指導(dǎo)意義的研究方法。CARNEVALE等[15]研究發(fā)現(xiàn),高速列車(chē)運(yùn)行時(shí)受電弓氣動(dòng)升力與列車(chē)運(yùn)行速度有關(guān),氣動(dòng)阻力占列車(chē)總氣動(dòng)阻力的8%~14%。HARA等[16]借鑒機(jī)翼設(shè)計(jì)理論,對(duì)受電弓滑板進(jìn)行空氣動(dòng)力外形優(yōu)化設(shè)計(jì),并利用風(fēng)洞試驗(yàn)驗(yàn)證了優(yōu)化后滑板外形的合理性。IKEDA等[17]將傳統(tǒng)雙受電弓改為單受電弓,降低滑板質(zhì)量來(lái)減小弓網(wǎng)磨損,利用此方法來(lái)提高受流質(zhì)量且降低了受電弓整體的噪聲。

但是,目前研究人員對(duì)于受電弓優(yōu)化研究多集中于提高其1個(gè)氣動(dòng)指標(biāo),對(duì)多個(gè)氣動(dòng)指標(biāo)綜合起來(lái)進(jìn)行優(yōu)化的研究較少,且關(guān)注部位多為受電弓本身?xiàng)U件外形優(yōu)化。因此本文從受電弓安裝艙形狀入手,提出了矩形、橢圓形、膠囊形和六邊形4種內(nèi)置式受電弓艙進(jìn)行CFD計(jì)算,通過(guò)優(yōu)選法得出明線狀態(tài)下改善受電弓整體氣動(dòng)阻力與升力的最優(yōu)方案。采用目前研究中應(yīng)用廣泛的多目標(biāo)優(yōu)化理論,對(duì)優(yōu)選法下4種受電弓內(nèi)置式導(dǎo)流艙中的最優(yōu)方案進(jìn)行多目標(biāo)優(yōu)化,進(jìn)一步提高受電弓整體的氣動(dòng)性能。

1 氣動(dòng)特性計(jì)算結(jié)果與分析

1.1 幾何模型構(gòu)建

CFD仿真計(jì)算采用常見(jiàn)的3編組模型,由于重點(diǎn)關(guān)注受電弓及其附近部位的氣動(dòng)特性變化,故為了提高后期仿真計(jì)算的效率,該模型去掉轉(zhuǎn)向架及一些凸出部位。其中,高速列車(chē)頭車(chē)(尾車(chē))長(zhǎng)度為27.5 m,中間車(chē)長(zhǎng)度為25.0 m,列車(chē)頭部流線形長(zhǎng)度為10 m,車(chē)輛寬3.3 m,高3.5 m。建立高速列車(chē)模型如圖1所示。

圖1 高速列車(chē)模型(單位:m)

在實(shí)際運(yùn)行中,受電弓有升弓和臥弓兩種狀態(tài)。由于受電弓氣動(dòng)升力對(duì)于弓網(wǎng)接觸的影響在本文的考慮范圍之內(nèi),故只研究受電弓升弓狀態(tài)。受電弓模型如圖2所示。

圖2 受電弓模型

所述4種內(nèi)置式受電弓導(dǎo)流艙模型如圖3所示。從圖3可以看出,4種內(nèi)置式受電弓導(dǎo)流艙模型橫向和縱向最大長(zhǎng)度一致,深度也相同,但是其輪廓形狀有所不同。

圖3 4種內(nèi)置式受電弓導(dǎo)流艙

1.2 計(jì)算方法選擇

現(xiàn)以列車(chē)明線單車(chē)運(yùn)行為仿真模擬工況,速度分別設(shè)定為250 km/h和350 km/h。列車(chē)最高運(yùn)行速度為350 km/h,馬赫數(shù)<0.3,故可認(rèn)為運(yùn)行時(shí)列車(chē)周?chē)諝獠豢蓧嚎s,因此將列車(chē)周?chē)鲌?chǎng)按不可壓縮粘性流場(chǎng)處理。數(shù)值模擬計(jì)算時(shí)采用三維不可壓縮Navier-Stokes方程對(duì)列車(chē)周?chē)鲌?chǎng)進(jìn)行求解,湍流模型采用k-ωSST模型。流場(chǎng)的控制方程為

(1)

式中,ρ為空氣密度;φ為流場(chǎng)通量;t為時(shí)間;u為流場(chǎng)速度矢量;Г為擴(kuò)散系數(shù);S為源項(xiàng)。

1.3 計(jì)算域及邊界條件設(shè)定

CFD計(jì)算的計(jì)算域設(shè)定如圖4所示,其縱向長(zhǎng)度為200 m,橫向長(zhǎng)度為80 m,高40 m。邊界abdc設(shè)為速度進(jìn)口條件,大小及方向與來(lái)流速度一致。其中,列車(chē)靜止不動(dòng),以風(fēng)與列車(chē)的相對(duì)運(yùn)動(dòng)來(lái)模擬列車(chē)的實(shí)際運(yùn)行;邊界efhg設(shè)為壓力出口條件,為1個(gè)標(biāo)準(zhǔn)大氣壓;邊界aegc、abfe和bdhf設(shè)為無(wú)滑移壁面條件;為更真實(shí)模擬列車(chē)運(yùn)行實(shí)際情況,將地面cdhg設(shè)置為無(wú)滑移運(yùn)動(dòng)壁面邊界條件,速度與來(lái)流速度一致。

圖4 計(jì)算域尺寸(單位:m)

1.4 網(wǎng)格無(wú)關(guān)性驗(yàn)證

為更好地節(jié)約計(jì)算資源,在對(duì)計(jì)算模型進(jìn)行網(wǎng)格劃分之前,需要對(duì)計(jì)算區(qū)域確定合理的網(wǎng)格數(shù)量。現(xiàn)對(duì)矩形受電弓艙計(jì)算模型劃分了3套網(wǎng)格,設(shè)定進(jìn)口速度v=350 km/h,模擬計(jì)算整車(chē)的氣動(dòng)阻力,三套網(wǎng)格的數(shù)量與相關(guān)計(jì)算結(jié)果見(jiàn)表1。

表1 網(wǎng)格無(wú)關(guān)性數(shù)據(jù)

A套網(wǎng)格數(shù)量設(shè)為680萬(wàn),C套網(wǎng)格數(shù)量設(shè)為1 050萬(wàn),A與B之間的誤差超過(guò)3%,但B與C之間的誤差小于3%,說(shuō)明計(jì)算區(qū)域的網(wǎng)格數(shù)量設(shè)定在825萬(wàn)個(gè)左右是合適的。

1.5 受電弓氣動(dòng)特性計(jì)算

采用Star CCM+軟件進(jìn)行計(jì)算模型的網(wǎng)格劃分和流場(chǎng)計(jì)算。劃分整個(gè)流體計(jì)算域的網(wǎng)格時(shí),將體網(wǎng)格最大網(wǎng)格尺寸設(shè)為1 024 mm,面網(wǎng)格最小網(wǎng)格尺寸設(shè)為8 mm。為提高數(shù)值模擬的真實(shí)性,需要在列車(chē)表面劃分邊界層,邊界層層數(shù)為7層,邊界層第1層網(wǎng)格距列車(chē)壁面的法向距離為0.5 mm,y+值落在30~100之間,網(wǎng)格總數(shù)為825萬(wàn)。網(wǎng)格劃分完成后,列車(chē)縱剖面附近體網(wǎng)格如圖5所示,受電弓和導(dǎo)流艙面網(wǎng)格如圖6所示。

圖5 列車(chē)縱剖面體網(wǎng)格分布

圖6 受電弓附近面網(wǎng)格

網(wǎng)格劃分、邊界條件設(shè)定完成后,分別進(jìn)行250 km/h和350 km/h速度下的氣動(dòng)模擬仿真計(jì)算。得到受電弓附近表面壓力云圖和其整體氣動(dòng)阻力與氣動(dòng)升力。其中,350 km/h受電弓附近表面壓力云圖如圖7所示。

圖7 受電弓艙表面壓力云圖(v=350 km/h)

由圖7可知,在高速列車(chē)運(yùn)行過(guò)程中,受電弓區(qū)域表面壓力分布由于其導(dǎo)流裝置形狀的不同而不同,在該區(qū)域,最大正壓及負(fù)壓出現(xiàn)位置類(lèi)似,但其分布范圍明顯不同。可以看出矩形導(dǎo)流艙使得受電弓表面壓力最小,在導(dǎo)流艙內(nèi)部壓力分布較其他更為均勻。

圖8為4種內(nèi)置式受電弓艙在兩種速度級(jí)下的氣動(dòng)阻力和升力值,可以明顯看出在兩種速度級(jí)下,矩形導(dǎo)流艙可使受電弓受到的氣動(dòng)阻力與升力最小,因此可以認(rèn)為4種方案中矩形內(nèi)置式受電弓導(dǎo)流艙能夠最大程度地改善受電弓在高速運(yùn)行過(guò)程的氣動(dòng)阻力與升力性能。

圖8 受電弓部分氣動(dòng)力計(jì)算值

由于矩形內(nèi)置式受電弓艙對(duì)于受電弓的氣動(dòng)性能改善程度最大,因此對(duì)其進(jìn)行多目標(biāo)優(yōu)化,進(jìn)一步挖掘其潛在的減阻和升力降低效果。

2 多目標(biāo)氣動(dòng)特性?xún)?yōu)化

2.1 多目標(biāo)優(yōu)化設(shè)計(jì)理論

提取矩形內(nèi)置式導(dǎo)流艙主要結(jié)構(gòu)尺寸參數(shù)為優(yōu)化變量,以350 km/h運(yùn)行時(shí)受電弓整體氣動(dòng)阻力和氣動(dòng)升力為優(yōu)化目標(biāo),需采用多目標(biāo)優(yōu)化理論,其實(shí)質(zhì)為在非單一目標(biāo)的情況下,綜合考慮各優(yōu)化變量變化對(duì)每個(gè)目標(biāo)的影響程度,使各目標(biāo)達(dá)到綜合最優(yōu)。多目標(biāo)優(yōu)化問(wèn)題的數(shù)學(xué)模型為[18]

(2)

本文進(jìn)行多目標(biāo)優(yōu)化時(shí)采用NSGA-Ⅱ算法進(jìn)行計(jì)算,該算法在優(yōu)化計(jì)算過(guò)程中能協(xié)調(diào)各個(gè)目標(biāo)函數(shù)之間的關(guān)系,盡可能找出使各個(gè)目標(biāo)函數(shù)都滿足要求的優(yōu)化變量解集,簡(jiǎn)單有效且在計(jì)算過(guò)程中比較明顯。一般情況下,在多目標(biāo)優(yōu)化過(guò)程中,很難得到一組同時(shí)滿足所有要求的最優(yōu)模型,往往一個(gè)目標(biāo)在變優(yōu)的同時(shí)會(huì)使得另一個(gè)目標(biāo)性能變差。進(jìn)行多目標(biāo)優(yōu)化研究,最終目的是在各個(gè)目標(biāo)的最優(yōu)性能之間進(jìn)行協(xié)調(diào)與權(quán)衡處理,使各優(yōu)化目標(biāo)最大程度達(dá)到最優(yōu)值[19-20]。

2.2 優(yōu)化設(shè)計(jì)變量選取

在三維建模軟件中對(duì)矩形內(nèi)置式受電弓導(dǎo)流艙進(jìn)行參數(shù)化建模,提取主要設(shè)計(jì)參數(shù)為多目標(biāo)優(yōu)化的設(shè)計(jì)變量,矩形內(nèi)置式受電弓導(dǎo)流艙主要設(shè)計(jì)參數(shù)示意如圖9和圖10所示,主要設(shè)計(jì)參數(shù)取值范圍見(jiàn)表2。

表2 主要設(shè)計(jì)參數(shù)值

圖9 矩形導(dǎo)流艙軸測(cè)圖

圖10 矩形導(dǎo)流艙橫截面示意

2.3 近似模型構(gòu)建

利用參數(shù)優(yōu)化軟件Isight,在其中集成軟件CATIA、Star CCM+和軟件本身自帶組件來(lái)搭建優(yōu)化設(shè)計(jì)流程,Isight軟件中搭建的優(yōu)化設(shè)計(jì)流程如圖11所示。

圖11 優(yōu)化設(shè)計(jì)流程

采用最優(yōu)拉丁超立方設(shè)計(jì)的試驗(yàn)設(shè)計(jì)方法(DOE)進(jìn)行樣本空間的樣本點(diǎn)選取,選取60個(gè)隨機(jī)樣本點(diǎn)構(gòu)建響應(yīng)面近似模型(RSM),用來(lái)代替CFD計(jì)算程序進(jìn)行計(jì)算。構(gòu)建的部分響應(yīng)面近似模型如圖12所示,從圖中可以看出,設(shè)計(jì)變量L和H對(duì)優(yōu)化目標(biāo)呈線性的影響。響應(yīng)面模型建立完成之后,應(yīng)對(duì)其進(jìn)行精度驗(yàn)證。文獻(xiàn)[18]中指出,近似模型的精度一般由全相關(guān)系數(shù)R2決定,R2值越接近1即可認(rèn)為所構(gòu)建的代理模型符合要求。現(xiàn)選取其他10組樣本值,進(jìn)行精度誤差分析,經(jīng)過(guò)計(jì)算得到代理模型R2=0.986,接近于要求值,故認(rèn)為該近似模型精度符合要求,可用其替代CFD計(jì)算程序進(jìn)行計(jì)算。

圖12 部分RSM響應(yīng)面三維圖

3 優(yōu)化結(jié)果分析

利用已構(gòu)建的響應(yīng)面近似模型代替CFD進(jìn)行計(jì)算,參考張亮等[21]針對(duì)超高速列車(chē)頭型進(jìn)行多目標(biāo)優(yōu)化的做法,在Isight軟件中選擇已經(jīng)集成的NSGA-Ⅱ算法對(duì)矩形內(nèi)置式導(dǎo)流艙進(jìn)行多目標(biāo)優(yōu)化設(shè)計(jì),將NSGA-II算法的初始采樣點(diǎn)設(shè)定為16個(gè),遺傳更新代數(shù)設(shè)置為20代,共計(jì)完成320次的優(yōu)化迭代計(jì)算,完成后可認(rèn)為矩形內(nèi)置式受電弓導(dǎo)流艙優(yōu)化計(jì)算完畢。

圖13為優(yōu)化過(guò)程中各優(yōu)化設(shè)計(jì)變量的迭代歷史曲線,可以看出變量L、W、H收斂到了不同的位置,這是由于變量L、W、H不能同時(shí)使兩個(gè)優(yōu)化目標(biāo)達(dá)到最優(yōu),即使得一個(gè)目標(biāo)優(yōu)化的同時(shí)另一個(gè)目標(biāo)變差,而Rad有了良好的收斂值,說(shuō)明隨著Rad的變化,兩個(gè)優(yōu)化目標(biāo)能夠同時(shí)變優(yōu)。

圖13 優(yōu)化設(shè)計(jì)變量迭代歷史曲線

圖14為設(shè)計(jì)變量L-Rad與優(yōu)化目標(biāo)之間的相關(guān)性圖,其相關(guān)系數(shù)范圍在[-1,1]之間,圖中柱形越高表示設(shè)計(jì)變量與優(yōu)化目標(biāo)相互影響程度越強(qiáng)。各設(shè)計(jì)變量對(duì)受電弓氣動(dòng)阻力影響效應(yīng)強(qiáng)弱依次為H、W、L、Rad,對(duì)氣動(dòng)升力影響效應(yīng)強(qiáng)弱依次為H、L、W、Rad。其中變量L和H與氣動(dòng)阻力之間的關(guān)系為正相關(guān),說(shuō)明L、H變大會(huì)導(dǎo)致受電弓部位氣動(dòng)阻力呈現(xiàn)增大的趨勢(shì),而設(shè)計(jì)變量W和Rad有著相反的趨勢(shì)。W和氣動(dòng)升力之間的相互關(guān)系為正相關(guān),說(shuō)明設(shè)計(jì)變量W增加會(huì)使得受電弓部位氣動(dòng)升力增大,而設(shè)計(jì)變量L、H和Rad與W對(duì)氣動(dòng)升力的影響相反。

圖14 設(shè)計(jì)變量與優(yōu)化目標(biāo)之間相關(guān)性

圖15為優(yōu)化過(guò)程中各優(yōu)化目標(biāo)的優(yōu)化歷史曲線,可以看出在經(jīng)過(guò)320次迭代之后,兩個(gè)優(yōu)化目標(biāo)值都有了收斂的趨勢(shì),優(yōu)化過(guò)程中出現(xiàn)的受電弓氣動(dòng)阻力最小可達(dá)到1 882 N,表面氣動(dòng)升力最小可達(dá)到74.5 N,較初始值氣動(dòng)阻力最大可減小6.5%,氣動(dòng)升力最大可減小2.7%,但需要注意的是最小的受電弓氣動(dòng)阻力和升力值并不是同一組優(yōu)化變量值所對(duì)應(yīng)的,即兩優(yōu)化目標(biāo)不能同時(shí)達(dá)到最優(yōu)值。

圖15 優(yōu)化目標(biāo)的優(yōu)化歷史曲線

圖16為兩個(gè)優(yōu)化目標(biāo)計(jì)算完成后的Pareto解集圖,可以看出,Pareto解的空間分布合理,三角形“▲”表示原始模型所對(duì)應(yīng)的氣動(dòng)阻力與升力值,“★”表示優(yōu)化計(jì)算完成后得到的Pareto前沿的最優(yōu)解。經(jīng)過(guò)對(duì)比,可以發(fā)現(xiàn)在優(yōu)化完成之后受電弓的氣動(dòng)阻力和氣動(dòng)升力都有了一定的降低。

圖16 Pareto前沿解集

由于優(yōu)化完成后的Pareto前沿解所對(duì)應(yīng)的一系列優(yōu)化目標(biāo)值相差不大,故為使優(yōu)化目標(biāo)同時(shí)最大程度達(dá)到最優(yōu),從其中選出了優(yōu)化模型1~優(yōu)化模型6。原始模型和優(yōu)化模型1~6所對(duì)應(yīng)的優(yōu)化變量和優(yōu)化目標(biāo)的值見(jiàn)表3。

表3 Pareto最優(yōu)解值

從原始模型和優(yōu)化模型1~6所對(duì)應(yīng)的設(shè)計(jì)變量及優(yōu)化目標(biāo)的值可以得出,受電弓艙長(zhǎng)度(L)和深度值(H)變化范圍不大,較原始模型變小,受電弓艙寬度(W)較原始模型變大,受電弓艙傾角(Rad)變化范圍較大,較原始模型增大,說(shuō)明受電弓艙往外傾斜有利于改善兩個(gè)優(yōu)化目標(biāo)。優(yōu)化后受電弓氣動(dòng)阻力最小的是優(yōu)化模型4,受電弓氣動(dòng)阻力大小為1 894.6 N,較原始模型降低5.9%,受電弓氣動(dòng)升力最小的是優(yōu)化模型2,氣動(dòng)升力為74.7 N,較原始模型降低2.5%。

4 結(jié)論

(1)CFD計(jì)算后得出在250 km/h和350 km/h速度時(shí),4種內(nèi)置式受電弓導(dǎo)流艙中矩形導(dǎo)流艙能夠更好地降低受電弓的氣動(dòng)阻力和升力。

(2)構(gòu)建代替CFD計(jì)算程序進(jìn)行計(jì)算的響應(yīng)面近似模型,模型精度為0.986,能很好地?cái)M合輸入與輸出之間的響應(yīng)關(guān)系。

(3)采用Isight中集成的NSGA-Ⅱ算法,對(duì)矩形內(nèi)置式導(dǎo)流艙進(jìn)行多目標(biāo)優(yōu)化設(shè)計(jì),得到的一系列優(yōu)化模型能使受電弓氣動(dòng)阻力最大降低5.9%,氣動(dòng)升力最大降低2.5%。

(4)相較于原始模型,得到的優(yōu)化模型受電弓艙傾角變化較大,較原始模型增大,其余變量變化不大。說(shuō)明受電弓艙往外傾斜有利于改善受電弓氣動(dòng)阻力和升力性能。

猜你喜歡
優(yōu)化模型
一半模型
超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
由“形”啟“數(shù)”優(yōu)化運(yùn)算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 亚洲欧美一级一级a| 亚洲激情99| 久热这里只有精品6| 亚洲人成网站色7777| 99精品国产自在现线观看| 精品成人一区二区| 91九色视频网| 久久久国产精品无码专区| 国产尤物jk自慰制服喷水| 色婷婷亚洲综合五月| 99在线免费播放| 成人年鲁鲁在线观看视频| 丰满的熟女一区二区三区l| 亚洲日韩精品欧美中文字幕| 欧美日韩国产在线观看一区二区三区 | 一本大道无码高清| 中国国产A一级毛片| 国产区精品高清在线观看| 日本高清免费不卡视频| 国产xx在线观看| 国产人人乐人人爱| 九九线精品视频在线观看| 久视频免费精品6| 国产情侣一区二区三区| 欧美福利在线| 香蕉eeww99国产精选播放| 亚洲精品另类| 国产极品美女在线| 国产高清毛片| 亚洲综合激情另类专区| 国产在线拍偷自揄观看视频网站| 亚洲AⅤ波多系列中文字幕 | 中文字幕在线看| 亚洲欧美日韩中文字幕在线| 国产国语一级毛片| 亚洲视频四区| 538国产视频| 欧美三級片黃色三級片黃色1| 亚洲精品爱草草视频在线| 国产欧美日韩va| 日韩欧美中文字幕在线精品| 一级毛片免费播放视频| 久久 午夜福利 张柏芝| 人与鲁专区| 女人毛片a级大学毛片免费| 91青青在线视频| 无码精品一区二区久久久| 日韩av无码精品专区| 亚洲欧美色中文字幕| 青青草原国产精品啪啪视频| 午夜天堂视频| 国产主播喷水| 中文纯内无码H| 97色婷婷成人综合在线观看| 99热最新网址| 91色在线观看| 国产高清自拍视频| 超薄丝袜足j国产在线视频| 在线视频亚洲色图| h视频在线播放| 精品久久香蕉国产线看观看gif | 国产白丝av| 99在线视频精品| 毛片基地美国正在播放亚洲| 亚洲无码熟妇人妻AV在线| 精品国产一二三区| 免费又黄又爽又猛大片午夜| 久久综合伊人 六十路| 露脸真实国语乱在线观看| 天堂成人在线| 国产在线自乱拍播放| 2021亚洲精品不卡a| 国产在线小视频| 免费A级毛片无码免费视频| 免费看一级毛片波多结衣| 乱人伦视频中文字幕在线| 久久国产乱子伦视频无卡顿| 国产全黄a一级毛片| 国产精品无码影视久久久久久久 | 国产三区二区| 日韩欧美网址| 国产在线拍偷自揄拍精品|