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

暴雨環(huán)境下的彈丸氣動(dòng)特性數(shù)值模擬

2022-07-01 10:08:30常思江張競(jìng)文
彈道學(xué)報(bào) 2022年2期

楊 鵬,常思江,張競(jìng)文

(南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 南京 210094)

氣象是影響火炮武器系統(tǒng)作戰(zhàn)效果的重要因素之一。隨著現(xiàn)代及未來(lái)戰(zhàn)場(chǎng)環(huán)境變得日益復(fù)雜,對(duì)彈藥武器裝備在復(fù)雜氣象條件下的作戰(zhàn)適應(yīng)性提出了更高的要求。暴雨、強(qiáng)氣旋等復(fù)雜、惡劣氣象對(duì)彈藥外彈道性能的影響,直接關(guān)系到武器裝備的實(shí)際使用效果。科研過(guò)程中,為了研究武器裝備本身的性能,往往對(duì)復(fù)雜、惡劣氣象予以規(guī)避,盡可能剝離其影響。目前在我國(guó)相關(guān)國(guó)軍標(biāo)中,凡涉及外彈道試驗(yàn)時(shí),均對(duì)氣象采取了條件控制的策略。然而,為了進(jìn)一步提高武器裝備對(duì)戰(zhàn)場(chǎng)環(huán)境的適應(yīng)性,有必要開(kāi)展特殊氣象條件下的彈藥外彈道性能研究。

目前國(guó)內(nèi)外研究人員已就暴雨條件對(duì)一些飛行器的氣動(dòng)性能的影響開(kāi)展了相關(guān)研究。WU等研究了降雨條件對(duì)NACA 0015翼型的影響。結(jié)果發(fā)現(xiàn),降雨對(duì)翼型的氣動(dòng)性能有不利的影響。DUNHAM等根據(jù)液態(tài)水含量來(lái)模擬不同降雨條件。風(fēng)洞實(shí)驗(yàn)表明液態(tài)水含量在13~20 g/m時(shí)機(jī)翼的最大升力下降了15%。TOMPSON等對(duì)NASA 4412翼型進(jìn)行了模擬降雨的風(fēng)洞實(shí)驗(yàn)。發(fā)現(xiàn)翼型氣動(dòng)性能受其表面水膜形成的位置和厚度所決定。此外,該團(tuán)隊(duì)通過(guò)風(fēng)洞實(shí)驗(yàn)研究了可濕表面涂層對(duì)翼型氣動(dòng)性能的影響。結(jié)果表明,在降雨條件下,具有可濕表面涂層的翼型可形成更薄且更均勻的液膜層。降雨環(huán)境中造成的氣動(dòng)損失相比于具有非濕涂層的翼型更小,這對(duì)飛行器表面材料使用具有一定的啟示。

近年來(lái),隨著計(jì)算機(jī)技術(shù)的快速發(fā)展,越來(lái)越多的學(xué)者通過(guò)計(jì)算機(jī)數(shù)值模擬手段對(duì)這一現(xiàn)象進(jìn)行研究。VALENTINE等為了研究降雨對(duì)機(jī)翼的影響,提出了單向和雙向的拉格朗日粒子追蹤算法,模擬了雨滴粒子對(duì)翼型的撞擊作用以及雨滴的濺射效應(yīng),發(fā)現(xiàn)雨滴濺射是導(dǎo)致翼型氣動(dòng)性能下降的原因之一,但并未研究翼型表面水膜的影響。張瑞民等通過(guò)計(jì)算流體力學(xué)(CFD)方法研究了降雨對(duì)翼型氣動(dòng)性能的影響,發(fā)現(xiàn)在小迎角下,翼型表面水膜導(dǎo)致其升力系數(shù)下降和阻力系數(shù)略微變大。同時(shí)該團(tuán)隊(duì)通過(guò)對(duì)CFD軟件中自帶的轉(zhuǎn)捩模型進(jìn)行修正,進(jìn)一步研究降雨對(duì)翼型氣動(dòng)性能的影響機(jī)理,發(fā)現(xiàn)翼型表面存在降雨產(chǎn)生的不均勻水膜,進(jìn)而導(dǎo)致其表面粗糙度增加,使邊界層提前轉(zhuǎn)捩,翼型的氣動(dòng)性能隨即降低。WU等開(kāi)發(fā)了一種雙向動(dòng)量耦合的Eulerian-Lagrangian數(shù)值模擬方法來(lái)研究暴雨中翼型的空氣動(dòng)力學(xué)性能,發(fā)現(xiàn)了雨滴碰撞導(dǎo)致邊界的動(dòng)量損失和不均勻水膜,共同導(dǎo)致了翼型在暴雨環(huán)境中空氣動(dòng)力學(xué)效率降低,認(rèn)為降雨條件對(duì)飛機(jī)氣動(dòng)參數(shù)有不利影響,降低了飛機(jī)穩(wěn)定性和控制性。FATAHIAN等同樣通過(guò)雙向動(dòng)量耦合的Eulerian-Lagrangian數(shù)值模擬方法,研究了降雨對(duì)單段型翼型和板條型翼型的的氣動(dòng)性能的影響。結(jié)果表明,降雨對(duì)板條型翼型的氣動(dòng)影響相比于單段型翼型更大。

綜上所述,有關(guān)飛行器在降雨環(huán)境中的氣動(dòng)特性的研究處于起步階段,過(guò)去的大部分研究都集中于翼型,對(duì)炮彈的研究則更少。本文采用計(jì)算流體力學(xué)方法,選取雙向動(dòng)量耦合的Eulerian-Lagrangian模型對(duì)暴雨環(huán)境進(jìn)行模擬,以美國(guó)M910彈丸和某120 mm迫擊炮彈為研究對(duì)象,分別研究其氣動(dòng)性能在暴雨環(huán)境中的變化,并開(kāi)展機(jī)理分析。

1 數(shù)值模擬方法

1.1 多相流模型

暴雨環(huán)境實(shí)質(zhì)為氣液兩相流流動(dòng),對(duì)其進(jìn)行數(shù)值模擬可采用的模型主要有兩種:Euler-Euler模型和Eulerian-Lagrangian模型。Euler-Euler模型是將不同的相處理成互相貫穿的連續(xù)介質(zhì);Eulerian-Lagrangian模型是將多相流流動(dòng)分為連續(xù)相和離散相,適用于一般離散相的體積分?jǐn)?shù)較小的情況。其中連續(xù)相通過(guò)直接求解N-S方程得到,而離散相則通過(guò)計(jì)算流場(chǎng)中的粒子運(yùn)動(dòng)得到。

具體而言,連續(xù)相是通過(guò)采用標(biāo)準(zhǔn)的-湍流模型求解質(zhì)量和動(dòng)量的定常雷諾平均N-S方程,求解彈丸的繞流場(chǎng)。而雨滴粒子被認(rèn)為是離散相,采用Eulerian-Lagrangian模型進(jìn)行模擬,該模型通過(guò)向連續(xù)相中注入雨滴顆粒形成。雨滴的拉格朗日運(yùn)動(dòng)方程可表示為

(1)

(2)

(3)

(4)

式中:p分別為雨滴相和空氣在方向上的速度;(-p)為單位質(zhì)量粒子所受的阻力;是雨滴粒子的阻力系數(shù);是雨滴粒子的雷諾數(shù);和為雨滴相和空氣相的密度,為雨滴粒徑的大小。

1.2 暴雨條件模擬

模擬暴雨環(huán)境需要考慮的因素通常有降雨率、雨滴粒徑和雨滴末速度。而對(duì)于大氣中雨量的測(cè)量,通常以液態(tài)含水量衡量,表示單位體積中包含的液體質(zhì)量,計(jì)算公式如下式:

=0054084

(5)

式中:為地面上的降雨率,單位mm/h。

降雨條件下,雨滴的粒徑分布非常復(fù)雜。為了簡(jiǎn)化計(jì)算,可選取在一定降雨強(qiáng)度下雨滴數(shù)量最多的粒徑值作為代表粒徑:

(6)

式中:和的經(jīng)驗(yàn)系數(shù)分別取13和0232;為雨型常數(shù),取225。

雨滴假設(shè)為不蒸發(fā)、不變形以及不相互作用的球體,在空中只受阻力和重力的作用。當(dāng)受力達(dá)到平衡時(shí),雨滴勻速下降,此時(shí)雨滴的速度稱(chēng)為雨滴的末速度,可采用下式估算:

(7)

式中:為雨滴的末速度。

1.3 離散隨機(jī)游走模型

在計(jì)算中本文還考慮了由于連續(xù)相的湍流引起的顆粒湍流擴(kuò)散。對(duì)于顆粒的湍流擴(kuò)散計(jì)算,使用了積分時(shí)間尺度的概念,其表示為顆粒沿著其運(yùn)動(dòng)軌跡d所經(jīng)歷的湍流運(yùn)動(dòng)狀態(tài)的時(shí)間。可采用下式表達(dá):

(8)

1.4 壁膜模型

采用Fluent中的Wall Film模型模擬雨滴和彈丸表面的碰撞作用。雨滴撞擊到彈丸表面形成壁膜模型可以分成四類(lèi):顆粒撞擊與壁面的相互作用、顆粒撞擊壁膜后的跟蹤,壁膜變量的計(jì)算以及相間耦合。

除上述外,還需考慮顆粒相(雨滴)和連續(xù)相(空氣)之間的相互作用。這種雙向耦合作用交替求解,直到兩相均達(dá)到穩(wěn)定狀態(tài),具體計(jì)算過(guò)程如圖1所示。

圖1 計(jì)算過(guò)程流程圖

2 計(jì)算和結(jié)果討論

2.1 計(jì)算驗(yàn)證

為了驗(yàn)證以上計(jì)算方法,本文通過(guò)對(duì)文獻(xiàn)[17]中的NACA 64-210翼型開(kāi)展數(shù)值模擬,計(jì)算液態(tài)含水量=25 g/m時(shí)翼型升力系數(shù)的變化,并與文獻(xiàn)進(jìn)行對(duì)比,結(jié)果如圖2所示。

圖2 NACA 64-210升力系數(shù)對(duì)比

從圖中可以看出,在無(wú)雨條件下,模擬值稍大于參考值,誤差在可接受的范圍之內(nèi);在降雨條件下,模擬值和參考值符合情況良好,計(jì)算得出暴雨對(duì)翼型升力系數(shù)的影響和文獻(xiàn)中基本一致,誤差均小于5%,驗(yàn)證了本文計(jì)算方法的正確性。

2.2 彈丸模型

2.2.1 M910彈丸模型

本節(jié)采用M910彈丸,其外形結(jié)構(gòu)與具體參數(shù)如圖3所示。對(duì)其外流場(chǎng)進(jìn)行網(wǎng)格結(jié)構(gòu)化劃分,如圖4所示。

圖3 M910彈丸外形圖

圖4 M910彈丸的彈體及其外流場(chǎng)網(wǎng)格

為了驗(yàn)證該模型求解的精確性,在不同馬赫數(shù)條件下對(duì)該彈丸進(jìn)行數(shù)值模擬。數(shù)值模擬的阻力系數(shù)與文獻(xiàn)[20]的對(duì)比結(jié)果如表1所示,可以看出兩者基本一致,誤差小于5%。

表1 M910彈丸阻力系數(shù)計(jì)算值與文獻(xiàn)[20]的對(duì)比情況

2.2.2 某120 mm迫彈模型

為了研究暴雨環(huán)境對(duì)不同彈型是否會(huì)產(chǎn)生不同的影響,本文還計(jì)算了某120mm迫彈(如圖5所示)在暴雨環(huán)境下阻力系數(shù)和升力系數(shù)的變化,該迫彈帶有八片尾翼。本文計(jì)算的馬赫數(shù)為0.7,為滿足亞聲速壓力遠(yuǎn)場(chǎng)邊界足夠遠(yuǎn)的條件,遠(yuǎn)場(chǎng)區(qū)域前方和徑向距離彈體20(為彈體直徑),遠(yuǎn)場(chǎng)后方距離彈體35,并對(duì)彈體附近網(wǎng)格進(jìn)行加密,整個(gè)區(qū)域網(wǎng)格數(shù)約為445萬(wàn)。其網(wǎng)格劃分如圖6所示。

圖5 迫彈模型圖

圖6 120 mm迫彈彈體及外流場(chǎng)網(wǎng)格

2.3 邊界條件

①空氣相邊界條件設(shè)置為壓力遠(yuǎn)場(chǎng)邊界條件,保證流場(chǎng)內(nèi)所研究的流體在流過(guò)模型時(shí)產(chǎn)生的擾動(dòng)不被邊界反射,取標(biāo)準(zhǔn)大氣海平面壓力和溫度作為入口壓力和溫度,為0.7,攻角取0°、5°、10°、15°和20°。

②離散相邊界條件:對(duì)離散相粒子采用非穩(wěn)態(tài)追蹤的方法,并定義一個(gè)粒子射入面。為了保證雨滴粒子能夠全部落到彈體上,本次研究取彈體上方5倍彈徑的平面為雨滴粒子射入面。彈丸表面設(shè)置為Wall-film邊界條件,表示雨滴落到彈丸表面會(huì)形成一個(gè)液膜,雨滴與彈丸表面的碰撞采用Stanton-Rutland模型;整個(gè)外流場(chǎng)表面設(shè)置為逃逸邊界。

③降雨條件:描述暴雨環(huán)境的具體參數(shù)如表2所示。

表2 本文研究采用的降雨參數(shù)

2.4 計(jì)算結(jié)果

2.4.1 M910彈丸計(jì)算結(jié)果

M910彈丸在=07、=30 g/m,不同攻角下的氣動(dòng)系數(shù)的變化情況如表3及表4所示。

表3 M910彈丸阻力系數(shù)計(jì)算結(jié)果(Ma=0.7)

表4 M910彈丸升力系數(shù)計(jì)算結(jié)果(Ma=0.7)

由表3和表4可知,在暴雨條件下,阻力系數(shù)整體變化不大,均有略微的降低。相較于阻力系數(shù),可以明顯觀察到升力系數(shù)受降雨的影響。升力系數(shù)的最大變化出現(xiàn)在5°攻角時(shí),較無(wú)雨條件下降了約14.5%。并且隨著攻角的增大,升力系數(shù)下降的幅度逐漸降低,在10°攻角時(shí)下降6.7%,而在15°以及20°攻角時(shí),升力系數(shù)的下降約為5%。

正如文獻(xiàn)[7]所提及,暴雨對(duì)翼型氣動(dòng)參數(shù)的影響是由翼型表面不均勻的水膜引起的。因此,本文通過(guò)雙向動(dòng)量耦合的方法對(duì)彈丸表面這種現(xiàn)象進(jìn)行研究。圖7為液態(tài)含水量=30 g/m時(shí),M910彈丸在不同攻角下的水膜高度情況。圖8為彈體表面場(chǎng)應(yīng)力分布。

如圖7所示,彈頭前端的水膜層高度較高,彈身部分的水膜層厚度較低。這是由于在亞聲速環(huán)境下(=0.7),彈頭表面所受氣流的切應(yīng)力較小,大部分的水膜集中在彈頭,而彈身表面的切應(yīng)力較大,超過(guò)了水膜附在壁面上的黏附力,導(dǎo)致水膜產(chǎn)生分離。此外,彈體表面的切應(yīng)力隨著攻角增大而增大。圖7中可以觀察到彈丸表面由于雨滴撞擊產(chǎn)生的厚度不同的水膜層。攻角為0°時(shí),較高水膜層幾乎都聚集在彈頭前表面,這是因?yàn)榍氨砻孀矒舢a(chǎn)生的可塑性雨滴最多,水膜層也較厚。隨著攻角增大,空氣側(cè)速度的向上分量對(duì)雨滴粒子軌跡產(chǎn)生影響,水膜的最大厚度逐漸向下表面移動(dòng),彈頭上表面附近水膜開(kāi)始減少,下表面附近的水膜開(kāi)始增加。攻角為20°時(shí),彈體下表面的水膜層厚度高于上表面,并且側(cè)面的水膜像一層薄片覆蓋在彈體表面。這是由于彈頭前表面的水坑不再穩(wěn)定,在氣流的作用下,不同雨滴粒子向下緩慢移動(dòng),并和其他雨滴粒子合并,直到合并了大部分的雨滴粒子,從而形成一層層均勻的水膜。

圖7 Ma=0.7、不同攻角下的M910彈體水膜高度分布(單位:m)

圖8 Ma=0.7、不同攻角下M910彈體表面切應(yīng)力分布(單位:Pa)

相比于無(wú)雨條件下的流動(dòng),降雨使得彈丸表面形成貼體的水膜層,當(dāng)后續(xù)雨滴粒子下落時(shí),會(huì)對(duì)彈丸表面的水膜進(jìn)行撞擊,形成凹坑,也就是產(chǎn)生了不均勻的水膜層。這一定程度上改變了彈丸表面的光滑度,增大了彈丸表面對(duì)氣體邊界層的擾動(dòng),從而使得氣體邊界層過(guò)早的從層流轉(zhuǎn)變?yōu)橥牧鬟吔鐚印2⑶液罄m(xù)的雨滴撞擊與湍流邊界層相互作用,降低了彈丸的氣動(dòng)性能。

此外,=0.7,=10°和20°下,M910彈丸的速度流場(chǎng)圖如圖9和圖10所示。

圖9 Ma=0.7、10°攻角下M910彈丸速度流場(chǎng)

圖10 Ma=0.7、20°攻角下M910彈丸速度流場(chǎng)

從圖9和圖10中可以看出,在降雨條件下,彈丸邊界附近的速度明顯低于無(wú)雨條件,這是由于水膜層的存在以及雨滴粒子與彈丸表面碰撞、回濺后引起彈丸邊界動(dòng)量的損失。并且從彈丸速度流場(chǎng)中可以看出,無(wú)雨條件下的彈丸上下表面的壓差要大于降雨條件下的壓差,因此降低了彈丸的升力系數(shù)。

綜上所述,彈丸表面不均勻的水膜層和雨滴顆粒的碰撞引起的邊界動(dòng)量損失是彈丸氣動(dòng)性能降低的主要原因。

2.4.2 120 mm迫彈計(jì)算結(jié)果

迫彈在=07、=30 g/m,不同攻角下的氣動(dòng)系數(shù)如表5和表6所示。和M910彈丸類(lèi)似的是,暴雨對(duì)迫彈阻力系數(shù)的影響程度不大,但對(duì)升力系數(shù)的影響顯著。從表中數(shù)據(jù)可以看出,該迫彈升力系數(shù)的最大下降出現(xiàn)在5°攻角條件下,下降了21.9%,并且隨著攻角的增大,暴雨對(duì)升力系數(shù)的影響程度逐漸降低,在20°攻角時(shí),升力系數(shù)較無(wú)雨條件下降了約12.9%。

表5 120 mm迫彈阻力系數(shù)計(jì)算結(jié)果(Ma=0.7)

表6 120 mm迫彈升力系數(shù)計(jì)算結(jié)果(Ma=0.7)

不同攻角下迫彈表面的水膜層分布如圖11所示。由圖可見(jiàn),水膜層基本都分布在迫彈上表面,與M910彈丸不同的是,水膜層較高的部分不僅僅分布于彈頭部,尾翼部分也存在較高的水膜層。并且隨著攻角的增大,彈體前表面的水膜變得不穩(wěn)定,在氣流作用下,不同雨滴粒子向下緩慢移動(dòng),導(dǎo)致水膜層逐漸向尾翼聚集,當(dāng)彈丸攻角達(dá)到20°時(shí),較厚的水膜層基本聚集在尾翼。

圖11 Ma=0.7、不同攻角下迫彈表面水膜層的分布

=0.7,=10°和20°下,迫彈的速度流場(chǎng)圖如圖12和圖13所示。

圖12 Ma=0.7、10°攻角下迫彈速度流場(chǎng)圖

圖13 Ma=0.7、20°攻角下迫彈速度流場(chǎng)圖

降雨條件下,迫彈邊界層附近的速度要明顯低于無(wú)雨條件,這是因?yàn)槠葟棻砻嫠拥拇嬖谝约坝甑瘟W优c彈表面碰撞、回濺后引起彈體邊界動(dòng)量的損失。與M910彈不同的是,迫彈上表面水膜層比M910彈丸分布更高,而且由于尾翼的存在,較高的水膜層不僅在迫彈頭部處聚集,在迫彈尾翼處也分布較高的水膜層。因此,后續(xù)雨滴粒子的撞擊使得迫彈表面形成更多的凹坑,對(duì)氣體邊界層產(chǎn)生了更大的擾動(dòng)。并且從圖12和圖13可以看出,暴雨環(huán)境中,迫彈上表面的速度要明顯低于M910彈丸,迫彈尾翼部分也存在較大的邊界動(dòng)量損失,這兩個(gè)方面導(dǎo)致暴雨對(duì)迫彈(尾翼彈)氣動(dòng)性能的影響程度要高于M910彈丸(旋轉(zhuǎn)彈)。

3 結(jié)論

本文以M910彈丸和某120 mm迫彈為對(duì)象,采用計(jì)算流體力學(xué)方法數(shù)值模擬了暴雨對(duì)這兩種彈丸氣動(dòng)特性的影響。利用雙向動(dòng)量耦合的Eulerian-Lagrangian模型對(duì)雨滴顆粒進(jìn)行追蹤,采用了Wall-Film模型來(lái)模擬雨滴粒子與彈丸表面的相互作用,并選用離散隨機(jī)游走模型模擬由于連續(xù)相湍流引起的雨滴粒子隨機(jī)擴(kuò)散。得到如下初步結(jié)論:

①暴雨對(duì)兩種彈丸的升力系數(shù)有較大影響,當(dāng)=0.7時(shí),M910彈丸升力系數(shù)最大下降了14.5%,120 mm迫彈升力系數(shù)最大下降了21.9%,且下降程度隨著攻角增大而逐漸降低;兩種彈丸的阻力系數(shù)受暴雨影響較小;由于尾翼的存在,暴雨對(duì)尾翼彈的影響程度要高于旋轉(zhuǎn)彈。

②根據(jù)數(shù)值模擬結(jié)果,觀察到暴雨導(dǎo)致彈丸表面形成不均勻水膜層;雨滴與彈丸表面碰撞所導(dǎo)致彈丸邊界的動(dòng)量損失。這兩個(gè)因素共同作用引起彈丸氣動(dòng)性能的變化。

主站蜘蛛池模板: 成年看免费观看视频拍拍| 日日噜噜夜夜狠狠视频| 国产一区二区网站| 国产欧美成人不卡视频| 一本一本大道香蕉久在线播放| 91麻豆精品国产91久久久久| 国产毛片高清一级国语| 国产成人高清精品免费软件 | 久久久久88色偷偷| 91精选国产大片| 日韩欧美网址| 91成人在线免费视频| 国产91九色在线播放| 国内毛片视频| 亚洲二区视频| 91香蕉视频下载网站| 日本国产精品一区久久久| 色哟哟精品无码网站在线播放视频| 国产91导航| 国产一区免费在线观看| 国产a在视频线精品视频下载| 国内精品一区二区在线观看| 免费人成又黄又爽的视频网站| 91无码视频在线观看| 波多野结衣爽到高潮漏水大喷| 国产一区二区三区在线观看视频 | 亚洲无卡视频| 国产永久在线视频| 久久窝窝国产精品午夜看片| 色丁丁毛片在线观看| 91福利免费| 欧美一级在线| 成人福利在线视频| 91九色国产在线| 亚洲日韩精品无码专区| 又黄又湿又爽的视频| 国产真实乱人视频| 国产精品白浆在线播放| 成人综合在线观看| 国产日韩欧美中文| 91丨九色丨首页在线播放| 天天综合网色| 国产精品毛片一区视频播| 国产精品一区二区不卡的视频| 波多野结衣二区| 国产免费羞羞视频| 99视频在线看| 中文字幕无码中文字幕有码在线 | 亚洲毛片一级带毛片基地| 久久精品人人做人人爽| 波多野结衣一级毛片| 99久久国产综合精品2020| 国模沟沟一区二区三区| aaa国产一级毛片| 朝桐光一区二区| 一本色道久久88综合日韩精品| 久久久黄色片| 国产精品55夜色66夜色| 九九热精品视频在线| 免费看av在线网站网址| 日韩欧美国产中文| 99热国产这里只有精品9九| 日韩在线影院| 欧美在线黄| 日韩 欧美 国产 精品 综合| 亚洲第一网站男人都懂| 欧美一级高清片欧美国产欧美| 亚洲免费毛片| 精品一区二区三区中文字幕| 国产亚洲第一页| 免费99精品国产自在现线| 亚洲欧美日韩综合二区三区| 在线播放91| 欧美亚洲一二三区| 1769国产精品免费视频| 国产黑丝视频在线观看| 51国产偷自视频区视频手机观看| 一区二区三区精品视频在线观看| 国产永久在线观看| 国产在线专区| 丁香婷婷激情综合激情| 国产成人高清精品免费|