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

基于EFFD參數(shù)化的風(fēng)扇/壓氣機(jī)葉片-端壁一體化伴隨優(yōu)化設(shè)計(jì)

2023-07-07 09:48:02李偉偉季路成
航空發(fā)動(dòng)機(jī) 2023年3期
關(guān)鍵詞:變形優(yōu)化設(shè)計(jì)

李 鑫 ,張 韜 ,李偉偉 ,周 玲 ,季路成

(清華大學(xué)航空發(fā)動(dòng)機(jī)研究院1,車輛與運(yùn)載學(xué)院2:北京 100084;3.北京動(dòng)力機(jī)械研究所,北京 100024;4.北京理工大學(xué)宇航學(xué)院,北京 100081)

0 引言

現(xiàn)代先進(jìn)航空發(fā)動(dòng)機(jī)的內(nèi)部流動(dòng)非線性程度高[1]、設(shè)計(jì)參數(shù)眾多,完全依靠設(shè)計(jì)師憑經(jīng)驗(yàn)進(jìn)行手動(dòng)設(shè)計(jì)已難以滿足日益提高的發(fā)動(dòng)機(jī)性能發(fā)展需求[2]。因此,依托計(jì)算機(jī)和計(jì)算技術(shù)對以眾多參數(shù)描述的葉片通道進(jìn)行精細(xì)設(shè)計(jì)與優(yōu)化[3]是發(fā)展現(xiàn)實(shí)與趨勢。在計(jì)算機(jī)優(yōu)化算法中,全局隨機(jī)類優(yōu)化算法與常規(guī)梯度類優(yōu)化算法均無法應(yīng)對設(shè)計(jì)參數(shù)眾多的精細(xì)化設(shè)計(jì)需求。

伴隨方法具有計(jì)算成本與設(shè)計(jì)變量數(shù)目無關(guān)的特性,是最為契合精細(xì)化設(shè)計(jì)概念的合理選擇。該方法由Lions 于1971 年首次提出[4],2000 年由Liu[5]、Dreyer 等[6]用于葉輪機(jī)氣動(dòng)外形優(yōu)化。歷經(jīng)近20 年的發(fā)展,伴隨優(yōu)化已成為葉輪機(jī)精細(xì)化設(shè)計(jì)與優(yōu)化的重要工具,在該方法基礎(chǔ)上,多點(diǎn)優(yōu)化[7]、多排葉輪機(jī)優(yōu)化[8-9]、多學(xué)科一體化優(yōu)化[10]、魯棒性設(shè)計(jì)[11]及氣動(dòng)不確定性分析[11]等問題均得到了充分研究。在伴隨優(yōu)化系統(tǒng)流程中,參數(shù)化方法描述了離散的設(shè)計(jì)變量取值與連續(xù)的外形幾何間一一對應(yīng)的映射關(guān)系,因而成為連接流動(dòng)數(shù)值仿真與目標(biāo)函數(shù)梯度分析、梯度尋優(yōu)的重要橋梁。2017 年Xu 等[12]指出,參數(shù)化方法應(yīng)具有全自動(dòng)、空間大、無震蕩、便于施加幾何約束等系統(tǒng)特征。歷經(jīng)近20 年的發(fā)展,以擾動(dòng)控制為主的擾動(dòng)參數(shù)化方法具有更高的設(shè)計(jì)空間與靈活性,更符合葉輪機(jī)精細(xì)化設(shè)計(jì)的需求,也因此被廣泛應(yīng)用,如Hicks-henne 函數(shù)、Cst 函數(shù)等[13]。作為3 維空間的隨機(jī)擾動(dòng)施加方法,自由變形技術(shù)(Free-form Deforma?tion,F(xiàn)fd)[14]自2018 年由John 等[15]首次用于葉輪機(jī)伴隨優(yōu)化后,便因其設(shè)計(jì)空間廣、3 維連續(xù)性強(qiáng)等優(yōu)勢獲得研究人員的關(guān)注。然而,F(xiàn)fd 應(yīng)用過程卻無法與葉輪機(jī)幾何分布特征完全契合。從原理上講,F(xiàn)fd 變形過程受限于長方體控制體,難以與圓柱形流道完全匹配,且盡管其能夠保證控制體內(nèi)部實(shí)體的曲線連續(xù),于邊界位置也難以施加約束,難以保證回轉(zhuǎn)幾何的周期性。

2013 年研究團(tuán)隊(duì)研發(fā)了基于簡單H 網(wǎng)格、凍結(jié)粘性法、薄層削尖模型的連續(xù)型伴隨優(yōu)化系統(tǒng)Opti?turbo[16],經(jīng)后續(xù)不斷改進(jìn),已具備工程可行的伴隨優(yōu)化應(yīng)用能力[17];為解決由網(wǎng)格及湍流模型簡化引發(fā)的數(shù)值精度問題,2017 年進(jìn)一步研發(fā)了葉輪機(jī)氣動(dòng)外形離散型伴隨優(yōu)化系統(tǒng)Turbosim_Un[18],實(shí)現(xiàn)了流動(dòng)數(shù)值仿真及伴隨敏感性分析精度的顯著改進(jìn)。然而,其仍使用基于Hicks-henne 函數(shù)[19-20]的擾動(dòng)參數(shù)化方法,不僅操作復(fù)雜,也嚴(yán)重限制了設(shè)計(jì)空間。為此,從實(shí)用的需求出發(fā),順應(yīng)行業(yè)領(lǐng)域發(fā)展趨勢,本文以拓展自由變形技術(shù)(Extended Free-form Deformation,EFFD)為基礎(chǔ)開發(fā)了具有葉片-端壁協(xié)調(diào)控制能力的全新參數(shù)化方法,并對跨聲速壓氣機(jī)轉(zhuǎn)子Rotor 67 進(jìn)行優(yōu)化,在驗(yàn)證TurboSim_Un 優(yōu)化能力的同時(shí)保證參數(shù)化方法的有效性。

1 基于拓展自由變形技術(shù)的參數(shù)化方法

拓展自由變形技術(shù)是Coquillart[21]于1900 年對FFD 方法的改進(jìn)拓展,其主要改進(jìn)在于對長方體控制體進(jìn)行變形以構(gòu)造任意形狀控制體,并結(jié)合局部變形概念改善控制體邊界附近曲線的連續(xù)性。文獻(xiàn)[20]已對EFFD 理論進(jìn)行了詳細(xì)描述,本文僅以Rotor 67轉(zhuǎn)子通道為例展示其實(shí)施過程:

(1)以通道形狀為控制體,繪制均勻分布的控制點(diǎn)。與流域空間不對等的長方體控制體是FFD 技術(shù)無法實(shí)現(xiàn)葉片-端壁一體化設(shè)計(jì)的主要限制,而EFFD 方法中通過控制體變形可構(gòu)造任意形狀控制體。為此,直接以通道流域?yàn)榭刂企w進(jìn)行變形處理,并均勻分布控制體內(nèi)部控制點(diǎn),如圖1(A)所示。

圖1 拓展自由變形過程

(2)對于任意實(shí)體位置,應(yīng)用局部變形理念確定其所在的局部控制體位置。為解決控制體邊界附近的曲線變形連續(xù)性等問題,EFFD 方法需使用局部變形,即僅實(shí)體位置附近的控制點(diǎn)干擾方能改變實(shí)體位置。局部控制體位置計(jì)算需利用Bezier 超曲面性質(zhì),以判斷包含該實(shí)體位置的具體子控制體微元。

(3)對控制點(diǎn)施加擾動(dòng),實(shí)現(xiàn)實(shí)體變形。施加某隨機(jī)擾動(dòng)后的邊界控制點(diǎn)位置如圖1(B)所示。

(4)自動(dòng)計(jì)算局部控制體內(nèi)虛擬控制點(diǎn)位置。為保證不同子控制體間曲線變形的斜率連續(xù),應(yīng)在子控制體內(nèi)施加虛擬控制點(diǎn),其不參與參數(shù)化變形,而根據(jù)真實(shí)控制點(diǎn)的變化自動(dòng)調(diào)整位置,如圖1(C)所示。

(5)使用牛頓迭代,計(jì)算實(shí)體位置在局部控制體內(nèi)的相對坐標(biāo)。

(6)對子控制體進(jìn)行FFD 變形,計(jì)算公式見文獻(xiàn)[14],最終變形結(jié)果如圖1(d)所示。

2 優(yōu)化案例

2.1 流動(dòng)數(shù)值仿真介紹

Rotor 67 轉(zhuǎn)子葉片是20 世紀(jì)70 年代由美國Nasa設(shè)計(jì)的低展弦比跨聲速風(fēng)扇/壓氣機(jī)轉(zhuǎn)子葉片,1980年Anthony 等使用激光測速等方法獲得了其詳細(xì)的內(nèi)部流場數(shù)據(jù)[22],包括特性線、截面馬赫數(shù)云圖及流向速度分布等,該轉(zhuǎn)子葉片已成為葉輪機(jī)領(lǐng)域廣為應(yīng)用的Cfd 數(shù)值校驗(yàn)、流動(dòng)機(jī)理分析和優(yōu)化設(shè)計(jì)的典型案例。Rotor 67轉(zhuǎn)子部分設(shè)計(jì)及試驗(yàn)參數(shù)見表1。

表1 Rotor 67轉(zhuǎn)子部分設(shè)計(jì)及試驗(yàn)參數(shù)

優(yōu)化過程Rotor 67 轉(zhuǎn)子計(jì)算流域進(jìn)口與葉根前緣的距離為0.10 M,流域出口與葉根尾緣的距離為0.12 M,子午面如圖2 所示。采用Autogrid 5 繪制網(wǎng)格,基于Sa(Spalart And Allmaras)模型的網(wǎng)格無關(guān)性校驗(yàn)結(jié)果如圖3所示。

圖2 Rotor 67子午面

圖3 網(wǎng)格無關(guān)性校驗(yàn)

經(jīng)過選擇,后續(xù)Rotor 67轉(zhuǎn)子優(yōu)化計(jì)算使用58萬網(wǎng)格方案,其壁面y+值如圖4 所示,基于此網(wǎng)格的Ro?tor 67 轉(zhuǎn)子特性線如圖5 所示,其中下標(biāo)choke 表示堵塞流量,本文后續(xù)堵塞流量均以圖5為基準(zhǔn)。

圖5 Rotor 67轉(zhuǎn)子特性線

2.2 優(yōu)化設(shè)置

優(yōu)化過程以效率為目標(biāo)函數(shù),并采用罰函數(shù)法對流量及壓比進(jìn)行約束

式中:η為絕熱效率;ω1為流量約束權(quán)值,取值約為1000;ω2為壓比約束權(quán)值,取值為200;M為流量;Pr為總壓比;下標(biāo)0表示原型的參數(shù)。

在采用EFFD 參數(shù)化方法的優(yōu)化過程中,軸向控制點(diǎn)數(shù)為10,周向控制點(diǎn)數(shù)為7,徑向控制點(diǎn)數(shù)為8,除保證曲線連續(xù)型的若干位置控制點(diǎn)、排除葉尖復(fù)雜流動(dòng)干擾的尖部控制點(diǎn)不變外,總變形控制點(diǎn)數(shù)為462,總設(shè)計(jì)參數(shù)的數(shù)量為1386。

3 結(jié)果與分析

3.1 優(yōu)化進(jìn)程

在上述設(shè)置下,本文Rotor 67 轉(zhuǎn)子伴隨優(yōu)化過程中相關(guān)參數(shù)變化如圖6所示。

圖6 優(yōu)化過程中性能參數(shù)變化

圖中,橫軸ITER 表示優(yōu)化后所獲得的幾何模型次序編號(hào),左縱軸為流量相對偏差δm及總壓比相對偏差δpr,右縱軸為效率絕對偏差ζη

由此可知,經(jīng)16 次伴隨優(yōu)化,Rotor 67 轉(zhuǎn)子效率提升了0.74%,而流量及總壓比偏差較小,分別為0.42%及0.11%。

3.2 優(yōu)化結(jié)果性能變化

為檢驗(yàn)優(yōu)化前后Rotor 67 轉(zhuǎn)子的整體性能,分別展示了原型及優(yōu)化后結(jié)果在70%、90%及100%轉(zhuǎn)速下的特性線分布,優(yōu)化前后特性線變化如圖7 所示。從圖中可見,本次伴隨優(yōu)化已實(shí)現(xiàn)轉(zhuǎn)子氣動(dòng)性能的顯著改進(jìn),并未發(fā)生工況漂移問題,且100%轉(zhuǎn)速狀態(tài)由于近峰值效率點(diǎn)的工況約束程度較高,堵塞流量僅增大約0.24%,而總壓比特性線與原型的仍保持較高的吻合度。具體而言,在100%轉(zhuǎn)速下轉(zhuǎn)子等熵效率提升達(dá)到近0.8%;此外,在90%和70%轉(zhuǎn)速下的變工況特性中,峰值效率及近堵塞工況的效率指標(biāo)同樣有所提升,且隨著轉(zhuǎn)速的降低,流量及壓比偏差減小,但效率提升也減小。

圖7 優(yōu)化前后特性線變化

3.3 優(yōu)化結(jié)果幾何分析

為描述伴隨優(yōu)化前后Rotor 67 轉(zhuǎn)子葉片發(fā)生的幾何變動(dòng),展示了若干展向位置的葉型截面,如圖8 所示。從圖中可見,葉片變形集中在前緣及前緣吸力面附近,而葉片后方變形極小,確保在多排一體化應(yīng)用時(shí)優(yōu)化結(jié)果不會(huì)對下游工況產(chǎn)生過大干擾。具體而言,前緣均向壓力面?zhèn)绕D(zhuǎn),且往往伴隨著前緣厚度的顯著減小,越靠近根部則絕對葉片厚度減小越明顯。

圖8 葉型變化

厚度是保證優(yōu)化結(jié)果結(jié)構(gòu)強(qiáng)度的重要指標(biāo),則進(jìn)一步提取了相應(yīng)展向位置的葉片厚度及中弧線彎角分布,如圖9 所示。從圖中可見,最大厚度位置前的葉片厚度減小是各截面葉型的共同變化特征,但最大厚度及其位置仍與原型的保持較高的吻合度,基本保證結(jié)構(gòu)強(qiáng)度不會(huì)大幅降低;另一方面,中弧線在保持彎角分布規(guī)律一致的前提下,在前緣后方局部區(qū)域內(nèi)彎角減小而后整體增大,且隨展高的增大變化更為明顯。

圖9 葉型厚度與彎角分布

優(yōu)化結(jié)果葉表變形量如圖10 所示,更形象地展示了葉片變化位置。從圖中可見,盡管葉根附近變形量絕對值更大,但考慮輪轂徑向變化的參與,其葉型變化并不十分明顯(圖8);除葉尖幾何未參與氣動(dòng)優(yōu)化外,80%展高附近葉片變化最為劇烈,且變化以前緣后方位置為主;尾緣附近變形均整體偏小,與圖8變化對應(yīng)。

圖10 優(yōu)化結(jié)果葉表變形量分布

從圖10 中還可見,葉根附近較大的變形量表明輪轂已發(fā)生顯著改變,作為葉片-端壁一體化優(yōu)化的另一主要調(diào)節(jié)對象,輪轂面的徑向變化如圖11 所示。從圖中可見,經(jīng)伴隨優(yōu)化,輪轂面具有典型“波浪形”分布特征:在前緣前、尾緣前葉表兩側(cè)存在3 組較為明顯的通道下壓區(qū),而在吸力面前緣附近、尾緣前通道中部則明顯上抬;整體而言,通道變形形式較為復(fù)雜,在任意軸向位置處輪轂壁面型線不再保持等徑或單調(diào)分布,S3截面部分葉片及輪轂形狀如圖12所示。

圖11 模型6輪轂面變形量

圖12 S3截面部分葉片及輪轂形狀

圖中,SS表示吸力面,PS表示壓力面。從圖中可見,經(jīng)伴隨優(yōu)化,除前緣附近外,葉片中后部區(qū)域內(nèi)原本近乎90°的葉片及端壁夾角轉(zhuǎn)化為鈍角,且尤以吸力面角區(qū)更為明顯;此外,端壁曲線呈現(xiàn)典型的“下壓-上抬-下壓”變化:吸力面及端壁交匯區(qū)域以二面角變化為主,通道中部輪轂上抬形成尖峰,而在壓力面與輪轂一體化區(qū)域并未顯著改變二面角,整體表現(xiàn)為平移下沉。

3.4 優(yōu)化結(jié)果流動(dòng)分析

為判斷流動(dòng)性能提升的內(nèi)在機(jī)理,提取了100%轉(zhuǎn)速峰值效率點(diǎn)位置的出口截面內(nèi)若干參數(shù)展向分布,如圖13所示。其中,各相對高度下絕熱效率為

圖13 出口截面若干氣動(dòng)性能參數(shù)展向分布

式中:Pt為總壓;Tt為總溫;γ為比熱比;下標(biāo)in表示進(jìn)口,out表示出口。

經(jīng)伴隨優(yōu)化,近乎在全展范圍內(nèi)效率指標(biāo)均顯著提升,而60%展高以下區(qū)域總壓均顯著提高,在壓比約束下尖部總壓略有降低。

以20%及80%展高為例,從S1截面的視角展示了優(yōu)化前后流動(dòng)的變化,如圖14 所示。從圖中可見,盡管20%展高尚未形成激波結(jié)構(gòu),但吸力面加速、尾跡均是流道內(nèi)主要流動(dòng)現(xiàn)象。在葉根附近,經(jīng)伴隨優(yōu)化,前緣吸力面?zhèn)雀咚賲^(qū)范圍明顯縮減,并結(jié)合葉片前半段的厚度調(diào)整,吸力面加速區(qū)范圍同樣顯著減??;在流速變化影響下,尾跡范圍明顯減小,表征尾跡損失顯著改善。而在80%展高截面內(nèi)盡管尾跡、激波結(jié)構(gòu)變化并不明顯,但借助厚度調(diào)整,激波前高馬赫數(shù)的范圍和數(shù)值均略有減小,同樣表征激波損失的改善。

圖14 若干展向位置相對馬赫數(shù)分布

優(yōu)化前后葉表吸力面靜壓分布如圖15 所示。從圖中可見,優(yōu)化結(jié)果的低壓區(qū)范圍明顯縮減,對應(yīng)葉表高等熵馬赫數(shù)區(qū)域減小,與前述結(jié)論一致。

圖15 優(yōu)化前后葉表靜壓分布

角區(qū)和輪轂形狀的變化及端區(qū)附近流動(dòng)性能的顯著提升,表明輪轂附近流動(dòng)已大幅改善,為此,展示了圖12中截取位置2的周向S3截面馬赫數(shù)分布,如圖16 所示。從圖中可見,相較于原型,優(yōu)化結(jié)果的主流區(qū)流速變化較小,而端區(qū)流速略有提高;以相同的速度等值線與端壁型線間相對面積作為衡量附面層損失的近似判據(jù)判斷其符合二面角原理,增大的二面角分布對改善角區(qū)附面層交匯損失具有正向增益作用。

圖16 截取位置2的S3截面相對馬赫數(shù)分布

優(yōu)化結(jié)果的輪轂面型線具有典型的“下壓-上抬-下壓”分布特征,其中部上抬的輪轂可有效抑制角區(qū)橫向二次流傳播,進(jìn)而改善吸力面分離渦結(jié)構(gòu),吸力面?zhèn)任簿壐浇谋诿鏄O限流線如圖17 所示。從圖中可見,經(jīng)伴隨優(yōu)化,尾緣后方的渦結(jié)構(gòu)甚至基本消失,吸力面后方的渦尺寸同樣略有減小,表明端區(qū)性能有極大改善。優(yōu)化前后通道出口截面的總壓分布如圖18 所示。從圖中可見,根部區(qū)域的低總壓區(qū)范圍明顯減小,再次表明根部區(qū)域流動(dòng)性能顯著改善。

圖17 吸力面?zhèn)任簿壐浇谋诿鏄O限流線

圖18 優(yōu)化前后通道出口截面總壓分布

然而,考慮壁面極限流線(圖17)為基于SA 湍流模型的計(jì)算結(jié)果,其在復(fù)雜湍流狀態(tài)下往往誤差較大。為此,以CFX 商業(yè)軟件為工具,使用剪切應(yīng)力輸運(yùn)(Shear Stress Transport, SST)湍流模型在相同邊界條件下對原型及優(yōu)化結(jié)果進(jìn)行氣動(dòng)仿真,并提取其相應(yīng)位置的流線分布,分離渦結(jié)構(gòu)計(jì)算結(jié)果如圖19所示。從圖中可見,經(jīng)伴隨優(yōu)化,尾緣后方渦結(jié)構(gòu)獲得了極大改善,整體表現(xiàn)為展向尺寸及軸向范圍大幅縮減。

圖19 分離渦結(jié)構(gòu)SST模型CFX計(jì)算結(jié)果

4 結(jié)論

(1)經(jīng)伴隨優(yōu)化,在100%轉(zhuǎn)速下Rotor 67 轉(zhuǎn)子近峰值效率點(diǎn)效率提升0.74%,且工況約束較高,特性線同樣得到全面改善。

(2)葉片與端壁間二面角的調(diào)整,對調(diào)節(jié)端區(qū)附面層交匯進(jìn)而調(diào)節(jié)角區(qū)分離具有重要影響,是改善端區(qū)氣動(dòng)性能時(shí)不可忽視的重要設(shè)計(jì)空間。

(3)拓展自由變形技術(shù)具備葉片通道的全3維調(diào)節(jié)能力,結(jié)合伴隨方法,可有效探索葉片通道的全3維設(shè)計(jì)空間。

猜你喜歡
變形優(yōu)化設(shè)計(jì)
超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
瞞天過海——仿生設(shè)計(jì)萌到家
“我”的變形計(jì)
設(shè)計(jì)秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設(shè)計(jì)叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
例談拼圖與整式變形
主站蜘蛛池模板: 国产综合欧美| 欧美影院久久| 国产精品第一区| 日韩毛片在线视频| 无码AV日韩一二三区| 国产精品一区二区久久精品无码| 成人国产小视频| AV天堂资源福利在线观看| 午夜老司机永久免费看片| 色噜噜狠狠狠综合曰曰曰| 99久久精品免费视频| 中文纯内无码H| 最新亚洲人成网站在线观看| 亚洲成在人线av品善网好看| 免费jizz在线播放| 午夜啪啪网| 无码AV高清毛片中国一级毛片| 欧洲极品无码一区二区三区| 精品国产www| 国产成人精品一区二区秒拍1o| 伊人色综合久久天天| 91啦中文字幕| 91精品专区国产盗摄| 色综合中文字幕| 欧美精品v日韩精品v国产精品| 一个色综合久久| 偷拍久久网| 精品一区二区三区水蜜桃| 国产一区免费在线观看| 尤物视频一区| a级毛片在线免费| 色婷婷视频在线| 99re精彩视频| 黄色网址免费在线| 再看日本中文字幕在线观看| 国产福利大秀91| 国产成人亚洲欧美激情| 精品偷拍一区二区| 天天色天天操综合网| 欧美一区二区福利视频| 丝袜美女被出水视频一区| 国产91蝌蚪窝| 亚洲男女在线| 国产欧美中文字幕| 亚国产欧美在线人成| 国产地址二永久伊甸园| 亚洲成A人V欧美综合天堂| 亚洲妓女综合网995久久| 国产剧情国内精品原创| 久久久久国产一级毛片高清板| a毛片免费观看| 综合网天天| 久久一色本道亚洲| 成人精品午夜福利在线播放 | 一本二本三本不卡无码| 亚洲日本韩在线观看| 久久久四虎成人永久免费网站| 国产精品久久久久婷婷五月| 国产精品一区不卡| 国产福利2021最新在线观看| 精品国产一区91在线| 亚洲天堂日本| 国产成人精品视频一区二区电影| 久久国产亚洲偷自| 亚洲一欧洲中文字幕在线| 免费播放毛片| 18禁黄无遮挡免费动漫网站| 亚洲欧美成人影院| 欧美国产综合视频| 国产SUV精品一区二区6| 亚洲乱亚洲乱妇24p| 亚洲日韩高清无码| 国产高清在线丝袜精品一区| 四虎成人在线视频| 国产超碰一区二区三区| 欧美色视频日本| 99视频在线免费| 综合五月天网| 欧美一级在线播放| 四虎成人在线视频| 高潮毛片无遮挡高清视频播放| 国产主播喷水|