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

底凹結構減阻效應數值分析

2022-04-29 06:37:20李彪王良明楊志偉
北京航空航天大學學報 2022年4期
關鍵詞:結構

李彪,王良明,楊志偉

(南京理工大學能源與動力工程學院,南京 210016)

彈箭阻力小不僅意味著彈箭的射程遠、打擊范圍大,更能使彈箭在落點處保持較大的存速,增強毀傷能力,因此,減阻一直是彈箭設計者的首要任務。彈箭的零升阻力由摩阻、底阻和波阻3部分組成,且在中等超聲速范圍內,底阻占總阻的40% ~50%[1],因此減小底阻對減阻增程意義重大。常見的減小底阻方法有底部排氣法和底凹減阻法。底凹是指在彈丸底部平面上向彈頭方向開一個空腔,具有底凹結構的彈丸被稱作底凹彈。雖然這種減阻方法已經在工程上實施并取得了成功,但是其減阻機理仍沒有被完全理解,尤其是在亞跨聲速范圍內,數值計算和實驗方法所得出的一些結論甚至相反。

Tanner[2]使用實驗方法測得不同攻角下底凹彈和普通彈丸的壓力分布,通過直接分析彈丸底部壓力變化來研究攻角對底凹減阻效應的影響。Sahu等[3-5]分別對超聲速和亞跨聲速范圍內的底凹彈進行了數值仿真,討論了底凹對彈丸飛行性能的影響。文獻[2-5]都是通過比較阻力系數及彈底壓力分布來分析底凹結構的減阻效應,并沒有從彈丸底部流場結構方面來分析底凹減阻效應的產生機理。隨后更高水平的技術,如紋影照相、鐳射多普勒測速被引入到實驗中用于測量超聲速圓柱體[6]和跨聲速錐形板底凹[7-8]的底部流場結構。對于超聲速底部流場,文獻[6]指出,在底部拐角處湍流邊界層折轉并膨脹,與圓柱體發生分離,并自此形成底部的自由剪切層。自由剪切層將外部無黏流和底部回流區隔離開。當自由剪切層到達對稱軸時會再度發生壓縮,之后附著在對稱軸上。自由剪切層下邊緣與對稱軸的交點為駐點,最小斷面處為喉部。對于亞跨聲速底部流場,文獻[7-8]指出,底凹減阻的主要原因是用屈從的流體邊界替代了原有的固體壁面,使得尾渦形成的位置稍稍遠離了底面,而且增大了渦脫離的頻率。這2點可能都是底凹結構減阻的原因,而這2點又與數值計算的結果相反。數值計算的結果指出,底凹的存在會增強尾部渦之間的相互作用,減小渦脫離的頻率,另外尾渦還會部分進入到底凹中。這是因為計算模型是一個完美的二維流場,而實驗永遠不可能完全不受三維效應的影響。文獻[7-8]中使用的風洞規模相對較小,一些三維效應是不可避免的。需要指出的是,文獻[6-8]研究的對象均是錐形平板和圓柱體,而數值計算的對象是二維簡化模型。文獻[9]將粒子圖像測速技術引入到實驗中對亞聲速錐形板底凹流場進行了實驗測量,所得結論與文獻[7-8]基本一致,且指出底凹結構不會對尾部渦街結構造成影響或改變。Fournier等[10-11]對超高聲速下底凹結構對CAN-4彈丸氣動力的影響進行了數值研究,指出超高聲速下底凹對氣動力影響很小。Simon等[12-13]研究了不同湍流模型對底部流場數值計算的影響和效果。

國內也已經對底凹彈丸的減阻效應進行了研究。谷嘉錦[14]通過風洞試驗分析了不同底凹深度對減阻效應的影響,并指出在超聲速條件下,底凹需從側面開孔才能起到減阻的作用,否則減阻效果并不明顯。文獻[15-17]通過對超聲速底凹彈側壁開孔建立數學力學模型,獲得了側壁開孔底凹效應的工程計算方法。西北工業大學的Pan和Cai[18]研究了底凹數目和形狀對底凹減阻效應的影響。近年來,有關底凹效應的研究[19-21]逐漸增多,但對底凹效應的研究往往偏重于實驗而忽視理論分析,并且國內只是對其減阻效果進行了數值模擬驗證,關于底凹減阻的機理探討尚不全面。另外,減小彈丸空氣阻力對改善彈丸的彈道性能(如增大射程和存速、縮短飛行時間等)十分有利,而底阻在彈丸空氣阻力中占比大,因此彈丸的底凹結構研究是十分有意義的。

本文采用標準k-ε湍流模型,對有無底凹結構的M910彈丸進行數值模擬,從底部流場特性,結合彈丸零升阻力系數對彈丸底凹結構在亞、跨、超聲速范圍下的減阻效應產生機理進行了分析。

1 數值計算方法

1.1 計算模型和網格

本文采用M910彈丸為對象進行數值計算,計算后的結果與實驗數據[22]進行比較以驗證計算的正確性。為M910彈丸引入圓柱底凹結構并將改造后的彈丸命名為M910BC。計算模型的外形和底凹形狀及位置尺寸如圖1所示。

圖1 計算模型Fig.1 Computational model

實驗和數值計算結果都表明,一旦底凹結構深度達到某一臨界值后,繼續增加底凹結構的深度將不會引起底凹結構減阻效應的任何改變,所有實驗中測得的觀測量也沒有發生實質性的變化。文獻[8]還指出,這一臨界深度可能只有彈徑的1/4或1/3,故本文取值為彈徑的1/2以確保超過底凹結構的臨界深度,從而不影響對計算結果的分析。

為了避免文獻[7-8]中提到的計算模型采用二維結構網絡出現的數值計算結果與實驗結果存在差異,本文采用三維結構網格。為了對比分析彈丸有無底凹結構對底部流場的影響,應盡量保證不引入計算網格方面的差異,故本文均采用三維結構網格,并保證除在底凹內部,2種彈丸的網格數目和形狀都完全一致。圖2為M910BC彈丸的底凹計算網格示意圖。

圖2 底凹計算網格Fig.2 Computational mesh of base cavity

為了驗證網格無關性,選取來流速度為1 190.7 m/s。對M910分別取網格數為120萬、200萬和310萬3組網格,對M910BC分別取網格數為130萬、200萬和320萬3組網格進行計算,比較阻力系數如表1所示。可以看出,2組彈丸網格數選取中等粗細和細網格阻力系數差異很小,大大節省了計算時間。因此,通過對網格無關性的研究,獲得了一套既不影響計算精度又能減少計算時間的最優網格。表2為2種彈丸的網格配置參數和最優網格數目,其中徑向邊界距離和前向邊界距離均指相應邊界到彈丸頂點的距離,后向邊界距離指后邊界到彈丸底面的距離,表中距離的單位均為1倍彈徑。

表1 網格無關性驗證Table 1 Grid independence verification

表2 計算網格特性Table 2 Computational mesh characteristics

1.2 Navier-Stokes CFD

1.2.1 控制方程

完整的積分型三維可壓縮Navier-Stokes方程可表示為

式中:W =[ρ ρu ρv ρw ρE]T為守恒變量,ρ為密度,u、v、w分別為3個方向上的速度,E為能量;FC為對流通量;FV為黏性通量;Q為源項;Ω為控制體的體積;S為控制體的各個表面。采用有限體積法解該方程。對流通量采用Roe格式離散,黏性通量采用中心差分格式離散。

1.2.2 湍流模型

標準k-ε湍流模型有較高的穩定性、經濟性和計算精度,應用廣泛[23],故本文選取該湍流模型進行計算。該模型可表示為

該模型的參數設置參照文獻[24]:Cμ=0.09,Cε1=1.44,Cε2=1.92,σK=1.0,σε=1.3,PrT=0.9。

因為壁面附近流場變量的梯度較大,所以壁面對湍流計算的影響很大。因此,在壁面附近要進行特殊處理,常用的方法有壁面函數法和近壁模型法。本文采用壁面函數法,即用半經驗公式將自由流中的湍流與壁面附近的流動連接起來。文獻[25]指出該方法的優點是可以保持雷諾應力與真實湍流一致,對旋流和分離流的模擬結果更符合真實情況。因此,該湍流模型適合于計算彈丸底部流場的分離和尾渦的形成。

近壁面阻尼函數為

1.3 邊界條件

流場的外邊界設置為壓力遠場。來流的壓力為101 325 Pa,溫度為288.5 K。彈體表面為非滑移絕熱壁面,非滑移即在固體邊界上的流體速度等于固體表面速度(本文指固體表面靜止),絕熱即壁面處溫度梯度為0,故非滑移絕熱的邊界條件為

式中:下標w代表壁面處;T為溫度;n為垂直于壁面的法單位向量。

表3顯示了所要計算的工況,還定義了馬赫數所屬的范圍(參照文獻[22]進行馬赫數范圍劃分),為接下來的分析工作帶來了便利。

表3 來流條件與馬赫數的關系Table 3 Relationship between incoming flow conditions and Mach number

2 計算結果與分析

2.1 計算結果與實驗數據對比

圖3為計算所得的M910彈丸阻力系數CD與M910彈丸實驗數據[22]的對比結果。為了便于直觀顯示底凹結構的減阻效應,將M910BC彈丸的計算結果也放于圖中。從圖3可以看出,M910彈丸的計算結果與實驗值在跨聲速段吻合很好,但在亞聲速和超聲速段,計算結果都略大于實驗值。雖然標準k-ε湍流模型對旋流和分離流的模擬具有優勢,但對強逆壓梯度不敏感[24]。隨著馬赫數增大,這種強逆壓梯度越明顯,故計算誤差越大。此外,標準k-ε湍流模型假設流動為完全湍流,忽略了分子黏性的影響。然而在亞聲速段,彈體表面的流動是層流和湍流混合的,此時的計算誤差也會相應較大。

圖3 計算阻力系數與實驗數據的對比Fig.3 Comparison between calculated drag coefficient and experimental data

對于M910和M910BC二種彈丸,強逆壓梯度發生和存在于船尾開始處,有無底凹結構對其影響不大。另外,由于2種彈丸的外形完全一致,相同馬赫數下彈體表面的流動特性也完全相同。綜上所述,2種彈丸在相同馬赫數下由于湍流模型所帶來的計算誤差基本相等,圖3中所顯示的M910BC彈丸阻力系數減少的原因可基本確信為底凹結構的減阻效應,而并不是計算的誤差。

此外,還可根據彈丸底凹結構的減阻效率隨馬赫數變化規律來判斷計算的正確與否。底凹結構減阻效率為

圖4給出了不同馬赫數下底凹結構的減阻效率。可以看出,在亞聲速段底凹結構減阻效率隨馬赫數增大而減小,這與文獻[7-9]的實驗結果規律一致。當馬赫數增大到跨聲速段,底凹結構減阻效應完全消失,這一現象在所有參考文獻中都未曾給出。從圖3可看出,在超聲速段,底凹減阻的差值是先增大,后保持不變,這與文獻[10]中實驗結果一致。然而圖4中減阻效率一直增大的原因是超聲速下阻力系數隨馬赫數的增大而減小。結合底凹結構在亞聲速和超聲速段的減阻效率規律,可以確定在跨聲速段底凹結構的減阻效率會逐漸到某一最小值然后緩慢增大,因此可以判斷本文的計算結果是合理的。綜上所述,本文所采用的數值計算方法基本捕捉到了各馬赫數范圍的流場物理規律,計算結果與實驗結果吻合很好,故本文的計算結果可以用于底凹結構減阻效應的機理性分析。

圖4 底凹結構的減阻效率Fig.4 Drag reduction efficiency of base cavity structure

2.2 亞聲速下的底凹效應分析

亞聲速下的彈丸底部流場結構復雜,不管是采用實驗還是數值計算手段,獲得真實的彈丸底部物理流場都是極具挑戰的任務。大量研究指出,當亞聲速氣流從船尾流向彈底時,船尾處的橫斷面不斷減小,在彈底面處突然消失為0。由于物體橫斷面的減小,由一圈流線所圍成的流管的橫斷面積S必然增大。根據連續方程,流體速度減小。再由伯努利方程可知,當流線上速度減小,壓強p會增大。因此,從船尾起始點處起,后面的流體將被阻滯。當壓強增大到一定程度,阻滯作用會使流體流動停止。隨著壓強繼續增大,就形成了反壓,在反壓的作用下靠近彈體表面的部分流體會形成逆流。此時,附面層將與彈體表面分離,形成渦街。

由于本文采用的是定常CFD計算方法,計算所得的尾渦將是一個等效渦街。該等效渦街只能表征實際渦街的大小和形成位置,關于渦街脫落信息則由彈底橫截面平均壓力系數來推斷。

圖5和圖6分別為彈丸底部中心線上的壓力系數曲線及彈丸底部橫截面上的平均壓力系數曲線。壓力系數Cp的定義如式(7)所示。彈丸底部橫截面是指圓心在底部中心線上、以彈徑為直徑且平行于彈底面的一系列圓面。

圖5 彈底中心線壓力系數曲線Fig.5 Centerline pressure coefficient curves of projectile base

圖6 彈底橫截面平均壓力系數曲線Fig.6 Cross-sectional average pressure coefficient curves of projectile base

式中:∞表示“遠場”處。

彈底中心線可從壓力系數最大點處分為2段,如圖5所示。壓力系數最大點的左側流場速度方向與來流相反,屬于逆流,而右側則屬于順流。對于M910彈丸,逆流是從壓力高的地方向壓力低的地方加速流動,在遇到固體底面后流體運動受到阻滯,流速降低,壓力系數增大。這種流體運動遇到固體壁面受到阻滯的作用是在距壁面很小一段距離內發生的,從圖5可以看出,在彈底尾部中心線上,壓力系數在很短距離內快速減小,之后又迅速增大。而對于M910BC彈丸,由于底凹結構的存在,以流體邊界取代原有的固體底面,逆流所受到的阻滯作用就沒那么強,速度減小的少,在底面上,M910BC的壓力系數比M910彈丸要小。

需要指出的是,圖5顯示的只是彈底中心線上的壓力變化規律,而圖6則顯示了整個彈底橫截面上的平均壓力變化規律。由圖6可知,相同馬赫數下,底凹底面上的平均壓力系數明顯大于固體底面上的平均壓力系數,表明底凹結構的存在可以增大彈底壓力,減小彈丸阻力。

圖7為亞聲速下彈底壓力分布云圖。由于計算結果在豎直平面內是軸對稱的,為了便于比較分析,將相同馬赫數下的2種彈丸計算結果各取一半放在圖中。

圖7 亞聲速下彈底壓力對比Fig.7 Comparison of projectile base pressure at subsonic speed

結合圖7可知,對于M910彈丸,由于固體底面的阻滯作用,在底面中心處形成了一塊面積不大的高壓區。該高壓區使得渦街內部的逆流方向發生改變,而并不是像固體底面那樣產生阻滯作用,因此方向改變的逆流速度并沒有減小。當逆流的方向改變90°后,流體又加速從彈底中心流向彈底四周,這個過程中流體的壓力在減小,因此在彈底四周形成了一片低壓區,如圖7所示。

底凹結構的引入給彈底流場帶來了2個變化:①在底凹內部形成了高壓“死水區”;②以“屈從”的流體邊界代替了原來的固體底面邊界。“屈從”的流體邊界對逆流的阻滯作用不如固體底面那么強,因此在底面中心所形成的高壓區不明顯,面積也較小。與M910彈丸一樣,渦街內部的逆流也會在底凹的流體邊界上發生轉向,流體由底面中心加速流向四周,從而導致沿彈徑方向的壓力減小。這樣,底凹內外壓差在拐點A(見圖7(a))處達到最大。當內外壓差達到一定值后,流體邊界會被破壞,導致渦街部分進入底凹結構中,而邊界破壞就發生在A點附近。由于在底面中心處流體邊界內外的壓差并不大,A點處邊界不會被破壞,渦街也不能由此進入到底凹結構內,如圖7所示。因此,渦街只是部分進入到底凹結構中,這一結論與文獻[7-9]中的實驗結果一致,反駁了底部渦街會進入到底凹內部且渦街運動會受到底凹限制從而提高底部壓力的假設。

圖5和圖6還顯示出,M910BC彈丸在底凹內部中心線上的壓力系數變化不大,從約0.076 m處開始減小,這也證明了文獻[8]中的結論,即底凹結構深度達到某一臨界值后,底凹結構減阻效果不再隨底凹結構深度的增加而變化。另外,“屈從”的流體邊界在遇到渦街逆流時,在邊界兩側會發生壓力波動,如圖5和圖6所示,這種壓力波動正是由流體邊界“屈從”的特性所引起的。

隨著馬赫數的增大,彈底整體壓力減小,形成的渦街會更強。渦街形成的位置也會更遠離彈底面,如表4所示。無論是固體底面還是流體邊界底面,渦街的逆流對其影響都會減弱。對于固體底面,渦街的遠離就相當于在固體底面附近自然形成一個高壓“死水區”,如圖7中紅色框圖區域所示。這樣因底凹結構所形成的高壓“死水區”的作用就會減弱,在亞聲速范圍內,隨著馬赫數的增大,底凹結構減阻效率會降低,如圖4所示。

表4 亞聲速下渦街中心信息Table 4 Information of vortex street center at subsonic speed

2.3 跨聲速下的底凹效應分析

跨聲速是亞聲速和超聲速之間的過渡階段,不僅同時具備了二者的特點,還具有自己的特點。從圖4看,底凹結構的減阻效率在跨聲速范圍內先減弱至零,后慢慢增強。在該范圍內低馬赫數下,流場規律還是與亞聲速流場規律一致,即隨馬赫數增大,底凹結構的減阻效果減弱。在該過程中,形成的尾部渦街越來越大,渦街形成位置也越來越遠離彈底面,如表3和表4所示。

圖8為跨聲速下彈底壓力對比。可以看出,當Ma=0.90時船尾處就已經形成了膨脹波,此時的馬赫角幾乎接近90°,附面層無法完成重新附著。隨著馬赫數的增大,馬赫角逐漸減小,附面層向彈軸方向偏轉的角度就越大。除此之外,經過膨脹波的流體會加速,附面層會快速流向彈底中心線完成重新附著,形成封閉的回流區。尾部渦街被“困”于回流區內,將不再脫離。此時定常計算的結果將與非定常計算的結果一樣[22]。

圖8 跨聲速下彈底壓力對比Fig.8 Comparison of projectile base pressure at transonic speed

由表5可知,對于2種彈型,當馬赫數大于0.98后,渦街的中心位置不再向彈底后方移動,而是轉而向彈頭方向移動。這是由于回流區形成以后,渦街無法脫離便被其后更高的壓力壓向彈頭方向,在渦街和固體壁面的相互作用下,最終渦街中心會反向移動。在該過程中,因渦街位置越來越遠而形成的彈底中心的高壓“死水區”被回流區所破壞并成為回流區的一部分,最終形成新的穩定回流。

表5 跨聲速下渦街中心信息Table 5 Information of vortex street center at transonic speed

對于M910彈丸,彈底中心高壓區與渦街中心的壓差在亞聲速段越來越大,而進入跨聲速后壓差又會減小,圖8中的結果證明了該結論。對于M910BC彈丸,彈丸底部的高壓“死水區”也隨著渦街的位置變化而變化,然而這些變化并未影響底凹內的“死水區”,因此底凹結構并不能起到減阻作用,如圖4所示。

2.4 超聲速下的底凹效應分析

超聲速下的彈丸固體底面底部流場結構已經清楚,下面討論底凹結構如何影響超聲速底部流場,從而實現減阻。

對于底凹彈丸,當回流區形成以后,回流區里的流體將與底凹結構中的流體混合,該混合隨著馬赫數的增大而加深。如圖9(a)、(b)所示,雖然自由剪切層和底凹壁面已經形成一個封閉的區域,但該區域內的流體并不全是回流。圖9(c)、(d)又顯示,在高馬赫數下流體充分混合。根據底排彈的工作原理可知,只要向低壓回流區中“添質加能”,就能減小底阻。因此,回流區里的流體與底凹中的流體混合可以看成一種“添質”行為,底凹結構在超聲速下也可以實現減阻。

圖9 超聲速下彈底壓力對比Fig.9 Comparison of projectile base pressure at supersonic speed

一方面,底凹減阻與2種流體的混合程度有關,當完全混合后,減阻效果不會繼續增大,如圖4所示;另一方面,在確定馬赫數的情況下,2種流體混合所造成的“添質”行為有一個上限,即無法將底凹結構中所有的流體都視作“添質”。這又證明了增加底凹結構深度到某一極限值后,底凹結構的減阻效果就不再改變。

圖9為超聲速下彈底壓力對比云圖。可以發現,底凹結構確實將駐點后移,自由剪切層的折轉角度減小,自由剪切層外緣抬高,喉部的高度增大,這些現象都說明了底凹結構起到了減小底阻的作用[26]。

3 結 論

通過對M910和M910BC(M910彈丸引入圓柱底凹結構)2種彈丸全馬赫數下的三維定常CFD仿真計算,比較2種彈丸的計算結果,分析了底凹結構減阻效應的產生機理。本文結論可為彈丸底凹結構的設計提供參考,還可為彈丸的彈道優化設計提供依據,具有一定的工程應用價值。本文得出以下結論:

1)馬赫數一定時,底凹結構深度達到極限值后,底凹結構減阻效率將保持不變,且該極限深度隨馬赫數的增大而減小。

2)底凹結構在亞跨聲速和超聲速范圍內的減阻機理不同。在亞跨聲速下,底凹結構以“屈從”的流體邊界代替了固體底面,且為彈丸引入了高壓“死水區”從而實現減阻。

3)在跨聲速下,隨著馬赫數的增大,尾部渦街會遠離彈底面,而在彈底面與渦街之間形成高壓區。該部分高壓區與底凹結構中的“死水區”作用相似,底凹結構減阻作用逐漸減小,最后消失。隨著馬赫數進一步增大,彈底上下附面層會在彈底中心線上形成封閉區域,即底部回流區。底部回流區的形成又使得渦街中心前移從而影響之前生成的高壓區,底凹減阻效果又將出現。

4)在超聲速下,底凹結構通過底凹內部流體與回流區中的流體混合,以向低壓回流區增加質量的方式來減阻,這種混合會隨著馬赫數的增大而加深,直到完全混合后底凹結構的減阻效果將不變。

猜你喜歡
結構
DNA結構的發現
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
循環結構謹防“死循環”
論《日出》的結構
縱向結構
縱向結構
我國社會結構的重建
人間(2015年21期)2015-03-11 15:23:21
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 成人午夜免费观看| 激情六月丁香婷婷四房播| 日韩中文精品亚洲第三区| 无码免费视频| 亚洲国产欧洲精品路线久久| 手机在线国产精品| 色妞www精品视频一级下载| 亚洲国产成人综合精品2020| 精品视频一区在线观看| 精品少妇人妻无码久久| 波多野结衣在线一区二区| 日本一区中文字幕最新在线| 99re这里只有国产中文精品国产精品 | 国产在线97| 久久久四虎成人永久免费网站| 拍国产真实乱人偷精品| 国产亚洲精久久久久久久91| 免费jjzz在在线播放国产| 国产成人亚洲综合a∨婷婷| 欧美视频二区| 亚洲综合色吧| 免费A∨中文乱码专区| 精久久久久无码区中文字幕| 国产爽歪歪免费视频在线观看| 日韩久久精品无码aV| 国产一区二区三区免费观看 | 国产一级小视频| 国产成人精品2021欧美日韩| 青青草原国产精品啪啪视频| 亚洲天堂视频在线观看免费| 成人福利在线观看| 天天视频在线91频| 99在线观看免费视频| 久久综合干| 亚洲资源站av无码网址| 亚洲高清在线播放| 一区二区三区高清视频国产女人| 日韩在线播放中文字幕| 久久久久久久97| 人禽伦免费交视频网页播放| 任我操在线视频| 国产91九色在线播放| 欧美一区二区自偷自拍视频| 久久精品波多野结衣| 久久99久久无码毛片一区二区| 2021最新国产精品网站| 久热中文字幕在线| 亚洲aaa视频| 亚洲欧美一区在线| 国产又大又粗又猛又爽的视频| 在线无码九区| 国产无码性爱一区二区三区| 久久精品人妻中文视频| 久久天天躁狠狠躁夜夜2020一| 国产精品第5页| 色婷婷天天综合在线| 国产亚洲欧美日本一二三本道| 亚洲无码电影| 国产精品女在线观看| 91av成人日本不卡三区| 国产免费久久精品99re不卡| 亚洲欧美成人影院| 国产区免费| 欧美中文一区| 久久精品国产999大香线焦| 亚洲国产成人久久精品软件 | 无码aⅴ精品一区二区三区| 久久不卡国产精品无码| 日韩无码黄色| 亚洲欧美日本国产综合在线| 99精品一区二区免费视频| 婷婷激情亚洲| 国产福利影院在线观看| 国产成人a毛片在线| 色网在线视频| 国产青榴视频在线观看网站| 免费毛片视频| 色香蕉影院| 欧美国产日韩另类| 国产欧美日韩免费| 国产日韩丝袜一二三区| 亚洲中文字幕在线观看|