亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        垂直發(fā)射水下航行體的通氣空化數(shù)值模擬研究

        2017-12-12 19:49:44劉濤濤王國玉張孟杰
        宇航總體技術(shù) 2017年4期
        關(guān)鍵詞:通氣孔旋渦空泡

        劉濤濤,黃 彪,王國玉,張孟杰

        (北京理工大學(xué)機(jī)械與車輛學(xué)院,北京100081)

        垂直發(fā)射水下航行體的通氣空化數(shù)值模擬研究

        劉濤濤,黃 彪,王國玉,張孟杰

        (北京理工大學(xué)機(jī)械與車輛學(xué)院,北京100081)

        基于均相流模型對垂直發(fā)射水下航行體周圍的通氣空化流動進(jìn)行了三維數(shù)值模擬,通過與試驗得到的壓力數(shù)據(jù)比較驗證了數(shù)值方法的可靠性,并討論了通氣時序?qū)ν饪栈鲌龅挠绊憽=Y(jié)果表明:計算結(jié)果與試驗吻合較好,通氣能提高空泡內(nèi)壓力且不同位置壓力一致,空泡末端存在回射高壓;隨著通氣時刻提前,泡內(nèi)壓力和回射高壓增大;空泡的發(fā)展過程伴隨著旋渦的形成與演變,通入空泡的氣體有3種流向,氣體的運(yùn)動導(dǎo)致空泡內(nèi)產(chǎn)生多個旋渦,旋渦的種類可以分為3類。

        數(shù)值模擬;通氣空化;通氣時序;旋渦結(jié)構(gòu)

        0 引言

        水下航行體垂直發(fā)射過程分為3個階段:出筒段、水中段和出水段[1]。其中,在水中段,隨著水下航行體運(yùn)行速度的增大,航行體表面的局部壓力降低,當(dāng)壓力降低至水的飽和蒸汽壓時,水發(fā)生汽化,形成覆蓋航行體的空穴,即自然空化??栈F(xiàn)象的存在對于高速運(yùn)動物體的水動力特性具有決定性的影響,由于自然空化形成的空泡不穩(wěn)定,在航行體運(yùn)行過程中容易受到外界干擾,發(fā)生斷裂脫落等現(xiàn)象,因此會給水下航行體帶來很大的流體動力載荷;同時,由于自然空泡泡內(nèi)壓力很低,航行體出水時,空泡發(fā)生潰滅產(chǎn)生很高的潰滅壓力,給航行體帶來很大的沖擊。

        采用通氣空泡來消除自然空化帶來的不利影響是水下航行體控制的一項新技術(shù)。通過向自然空化形成的空穴中通入不可凝結(jié)的氣體,來增強(qiáng)空泡的穩(wěn)定性,提高泡內(nèi)壓力,改善航行體受到的流體動力載荷和潰滅高壓。通氣空泡是一種非常復(fù)雜的高速流動現(xiàn)象,涉及多相流、湍流、質(zhì)量輸運(yùn)、可壓縮性和非定常性等復(fù)雜的流動機(jī)制[2]。目前國內(nèi)外學(xué)者研究通氣空泡主要通過實驗研究和數(shù)值模擬兩個方面。Semenenko[3]根據(jù)實驗結(jié)果提出了回射流、雙渦以及空泡振蕩3種泄氣方式,并且確定了3種泄氣方式的范圍,提出了通氣空泡的理論計算公式。Kawakami等[4]和Wang等[5]在通氣超空泡的形成過程中觀測到了兩種不同形式的回射流:第一種形式的回射流出現(xiàn)在回射流泄氣向雙渦管泄氣之前,且通入氣體量較小時,空泡整體長度較小,呈現(xiàn)出云霧狀,氣液界面模糊;第二種形式的回射流出現(xiàn)在第一種形式向雙渦管泄氣轉(zhuǎn)變的過程中,空泡整體呈現(xiàn)出透明狀,只有尾部出現(xiàn)小范圍的云霧狀。Karn等[6]綜合利用高速攝像機(jī)和動態(tài)壓力傳感器對通氣超空泡的氣體泄氣方式進(jìn)行了更為詳細(xì)的分析,指出通氣空泡泄氣除了4種基本的穩(wěn)定泄氣方式外,還存在另外5種泄氣方式:泡沫狀空穴、四渦-回射流泄氣、雙渦-四渦泄氣、雙渦-回射流泄氣、空泡界面波動-雙渦泄氣和耦合渦泄氣。賈力平等[7]利用中速可持續(xù)通氣空泡水洞開展了空化器參數(shù)對通氣超空泡形態(tài)影響的研究。張敏弟等[8]采用高速全流場顯示技術(shù)和動態(tài)應(yīng)變式測力系統(tǒng)相結(jié)合研究了繞圓盤空化器通氣空化流動發(fā)展過程和動力特性。在數(shù)值模擬方面,郭建紅等[9]基于輸運(yùn)類空化模型模擬了繞圓盤的軸對稱通氣超空泡流。時素果等[10]采用濾波器湍流模型分析了繞空化器的空泡形態(tài)、流動結(jié)構(gòu)和阻力特性。

        本文采用CFX軟件,對真實水域中垂直發(fā)射的高速水下航行體進(jìn)行了三維數(shù)值模擬,通過在控制方程中耦合重力源項考慮重力效應(yīng)。由于在實際發(fā)射過程中,往往在通氣前航行體周圍已經(jīng)被發(fā)射筒內(nèi)帶出的部分氣體覆蓋,并伴隨著自然空化,因此本文的計算也基于帶空泡狀態(tài)下的航行體進(jìn)行。通過與試驗壓力數(shù)據(jù)的對比分析驗證了計算結(jié)果的正確性,討論了通氣時序?qū)ν饪栈鲌龅挠绊?并對空泡發(fā)展及流場進(jìn)行了分析。

        1 數(shù)值計算方法

        1.1 基本控制方程

        計算中,采用均相流模型,則質(zhì)量方程、Favre平均的N-S方程為

        其中,下標(biāo)i和j分別代表坐標(biāo)方向,ρm、u和p分別為混合介質(zhì)的密度、速度和壓強(qiáng),μ和μt分別為混合介質(zhì)的層流和湍流動力黏性系數(shù)。

        1.2 湍流模型

        本文采用工程計算中廣泛使用的標(biāo)準(zhǔn)k-ε湍流模型,該模型由Launder等于1972年提出[11],對于均相平衡流動的數(shù)值計算,標(biāo)準(zhǔn)k-ε模型的控制方程為:

        其中,P t為湍動能生成項,k、ε分別為湍動能和湍動能耗散率,μT為湍流黏性系數(shù)。模型常數(shù)分別為:C1=1.44,C2=1.92,σε=1.3,σk=1.0。

        1.3 空化模型

        Kubota空化模型結(jié)合泡間兩相流動理論,基于單個空泡生成和發(fā)展時空泡體積的變化,基于Rayleigh-Plesset空泡生長方程推導(dǎo)出了如下蒸發(fā)和凝結(jié)相的表達(dá)式:

        其中,αnuc為汽核體積分?jǐn)?shù),RB為汽泡半徑,pv為汽化壓強(qiáng)。Fe和Fc分別是蒸發(fā)和凝結(jié)經(jīng)驗系數(shù)。計算中,相關(guān)的經(jīng)驗系數(shù)設(shè)定參考Kubota的論文[12]。

        1.4 計算域和參數(shù)設(shè)置

        通氣裝置示意圖如圖1所示,通氣裝置位于航行體肩部,頭部為錐形,計算域的設(shè)定如圖2所示,入口設(shè)定來流速度,按照航行體運(yùn)行速度給定;出口設(shè)定壓力,隨航行體運(yùn)行深度而變化;計算域側(cè)面設(shè)定為開放邊界。網(wǎng)格劃分采用全結(jié)構(gòu)化網(wǎng)格,由于航行體通氣裝置結(jié)構(gòu)復(fù)雜,為了控制網(wǎng)格數(shù)量,分別對通氣裝置和航行體外部流域劃分結(jié)構(gòu)化網(wǎng)格,在計算中通過interface將兩套網(wǎng)格對接進(jìn)行計算,計算總網(wǎng)格為200萬左右。

        圖1 通氣裝置示意圖Fig.1 Sketch of the gas ventilated device

        圖2 計算域及邊界條件Fig.2 Outline of the computational domain with boundary condition

        圖3 網(wǎng)格示意圖 (左邊為外部流域,右邊為通氣裝置)Fig.3 Grids in the computation(Left is out-domain, right is ventilated device)

        計算中選取的通氣率為Q v=34.8,通氣率的定義為為通氣體積流量,v為航行體運(yùn)行的平均速度,S為通氣面積。

        2 結(jié)果與分析

        2.1 壓力分析

        圖4給出了通氣時刻,初始空泡的示意圖,A點為被初始空泡覆蓋的測點,B點、C點為初始空泡外的測點。圖5給出了數(shù)值模擬和試驗中得到的航行體表面壓力隨時間發(fā)展規(guī)律的對比圖。其中縱坐標(biāo)為壓力系數(shù),橫坐標(biāo)為無量綱發(fā)射時間t/T,t為通氣后航行體運(yùn)行的時間,T為通氣時刻至航行體到達(dá)水面的總時間。從圖5可以看出,數(shù)值模擬得到的壓力整體變化趨勢與實驗值吻合較好。

        圖4 測點位置示意圖Fig.4 Positions of monitoring points

        對于A點,由于不可凝結(jié)氣體的通入,氣體的動能瞬間轉(zhuǎn)化為壓能,空泡內(nèi)壓力得到顯著提高,通氣空泡發(fā)展到穩(wěn)定階段后,壓力停止上升,隨著航行體的運(yùn)行,環(huán)境壓力不斷降低,而泡內(nèi)壓力受環(huán)境壓力的影響也逐漸降低。對于B點,空泡未發(fā)展到測點位置時,壓力值隨水深減小而下降,當(dāng)空泡發(fā)展到測點位置時,壓力值上升,出現(xiàn)一個高壓,這是由于空泡閉合造成的回射高壓,圖6(a)中給出了該時刻的壓力場和空泡形態(tài)圖,可以很清楚地看到B點處于空泡閉合處,且該處有一個高壓區(qū)。隨著空泡繼續(xù)向后發(fā)展,測點被空泡覆蓋,壓力迅速下降至泡內(nèi)壓力,并隨水深減小而下降。對于C點,可以觀察到試驗值中出現(xiàn)了2個高壓,第一個高壓是由于航行體完全離開發(fā)射筒時水氣作用劇烈,造成空泡外的測點壓力集體升高,本文數(shù)值模擬中未考慮出筒,因此數(shù)值模擬值沒有觀察到這個高壓,除此之外,C點的變化規(guī)律與B點一致。值得注意的是,B測點和C測點試驗中得到的壓力峰值出現(xiàn)時刻總是滯后于數(shù)值結(jié)果,這可能是試驗過程的不確定外界因素影響了空泡的發(fā)展速度,導(dǎo)致空泡發(fā)展至測點的時間延長,加之壓力傳感器存在一定的響應(yīng)時間造成的。

        圖5 航行體表面測點壓力對比Fig.5 Comparisons of the pressure surrounding the vehicle surface in experimental and numerical results

        圖6給出了計算終了時刻,航行體表面壓力隨軸向位置的變化,其中x為距頭部距離,L為航行體軸向長度。從圖中可以看出壓力變化趨勢大致可以分為4個階段:

        1)壓力下降段,此區(qū)域?qū)?yīng)于航行體頭部,其中壓力最低點位于其肩部位置,這是流動分離造成的;

        2)壓力平直段,此區(qū)域?qū)?yīng)于航行體表面空泡覆蓋位置,其壓力值不隨水深變化且遠(yuǎn)小于無空泡的區(qū)域;

        3)壓力峰值段,此區(qū)域?qū)?yīng)于空泡尾端,這是由于空泡閉合形成的滯止高壓;

        4)壓力二次下降段,此區(qū)域?qū)?yīng)于航行體無空泡區(qū)域,其壓力值保持為環(huán)境壓力。

        對比數(shù)值模擬和試驗結(jié)果發(fā)現(xiàn),兩者在空泡覆蓋段吻合較好,在壓力峰值段中高壓值和高壓值出現(xiàn)位置存在一定偏差,即試驗測量高壓值小于數(shù)值結(jié)果且出現(xiàn)位置相對滯后,這可能是實驗中的測量誤差所致。

        圖6 觸水時刻航行體表面壓力系數(shù)對比Fig.6 The pressure surrounding the vehicle surface when touching on the water surface

        對于垂直發(fā)射的水下航行體,通氣時序?qū)张莅l(fā)展以及泡內(nèi)壓力有顯著影響。本文選取了4個不同時刻t A、t B、t C、t D開始通氣,分別進(jìn)行了計算,其中t A<t B<t C<t D。圖7給出了計算終了時刻航行體的表面壓力分布曲線,壓力變化趨勢與圖6所示基本一致。從圖6中可以發(fā)現(xiàn),隨著通氣時刻提前,回射高壓出現(xiàn)的位置向后移動,說明空泡長度變大,這是由通氣總質(zhì)量決定的;同時,空泡內(nèi)壓力值和回射高壓均增大。由于靠近航行體尾部,水深增大,環(huán)境壓力增大,而空泡內(nèi)壓力和回射高壓均與環(huán)境壓力成正比,故泡內(nèi)壓力和空泡尾部回射高壓均增大。

        圖7 不同通氣時序下航行體的表面壓力Fig.7 The pressure surrounding the vehicle surface of different timing variation of ventilation rate

        2.2 通氣對空泡發(fā)展的影響

        2.2.1 空泡形態(tài)分析

        為了進(jìn)一步分析通氣對空泡形態(tài)的影響,圖8給出了通氣之后空泡形態(tài)的變化過程,t0為通氣時刻,t9為航行體觸水時刻。通氣前,航行體前半部有較小尺度的空泡覆蓋,通氣之后,空泡的中前段首先發(fā)生膨脹,直徑增大,產(chǎn)生一個環(huán)狀的凸起,隨著航行體運(yùn)行,環(huán)狀凸起逐步向空泡后半段移動,同時空泡的長度和直徑增大,t5時刻,環(huán)狀凸起移動到空泡末端,之后空泡形狀穩(wěn)定,空泡長度和直徑繼續(xù)增大。

        圖8 空泡形態(tài)隨時間的演變過程Fig.8 Time-evolution process of cavity shapes

        圖9給出了空泡長度的變化曲線,其中x1為空泡長度。從圖9中可以看出,t0~t3時刻,由于通氣引起的環(huán)狀凸起離空泡末端較遠(yuǎn),空泡增長的速度基本不變,空泡長度呈線性增大;t3~t4時刻,環(huán)狀凸起移動到了空泡后端,在這部分氣體的擠壓作用下,空泡增長的速度有一個瞬時的提升;t4~t6時刻,環(huán)狀凸起繼續(xù)向末端移動,空泡的增速變慢;t6時刻之后,環(huán)狀凸起移動到末端,空泡增長速度基本不變。

        圖9 空泡長度隨時間的變化Fig.9 Time-evolution process of cavity length

        2.2.2 旋渦結(jié)構(gòu)分析

        圖10給出了整個空泡內(nèi)的流線圖,并用黑色實線表示空泡的輪廓。圖11給出了通氣裝置附近的流線圖,以便觀察通入氣體的流向。結(jié)合圖10、圖11可以看出,空泡的發(fā)展過程伴隨著旋渦的形成與演化,通氣空泡內(nèi)大體存在3類旋渦:Ⅰ)靠近通氣孔緊貼壁面,其尺度較小 (圖10藍(lán)色虛線);Ⅱ)靠近空泡閉合端,其尺度較大,延伸至空泡中部(圖10紅色虛線);Ⅲ)對于大尺度的空泡,空泡中部還存在若干個小旋渦,數(shù)量隨空泡尺度增大而增加(圖10黑色虛線)。

        圖10 空泡內(nèi)流線圖Fig.10 Streamline patterns inside the cavity

        通入的氣體在空泡內(nèi)的運(yùn)動分為3種流向: a)從通氣孔流出后反向朝航行體肩部流去,這是因為肩部存在流動分離造成的低壓區(qū),氣體到達(dá)空泡邊緣后再次反向,變?yōu)橘N壁向后流動 (圖11藍(lán)色實線);b)從通氣孔流出后隨主流向后運(yùn)動,到達(dá)空泡尾部后閉合(圖11紅色實線);c)從通氣孔流出后遭遇主流阻礙,貼壁形成反向射流,在通氣孔后部形成旋渦 (圖11綠色實線)。

        圖11 通氣裝置附近流線圖Fig.11 Streamline patterns near the ventilated hole

        旋渦是由主流和貼壁的反向射流的相互作用引起的,當(dāng)主流的強(qiáng)度增加時,旋渦收縮,當(dāng)主流的強(qiáng)度減弱時,旋渦擴(kuò)大。對于通氣孔附近的旋渦,即第1類旋渦,由圖11可以看出,t0~t3時刻,旋渦尺度明顯減小,這是因為氣體通入后,反向朝航行體肩部流去,對旋渦起到壓迫作用。第1類旋渦的收縮使得整個空泡有收縮的趨勢,造成初始階段空泡長度增速較慢。t4時刻之后,第1類旋渦尺度恢復(fù)到初始大小,這是由于通氣孔流出的一部分氣體遭遇主流阻礙,貼壁形成反向射流,在通氣孔后部形成旋渦。當(dāng)主流強(qiáng)度與反向射流強(qiáng)度達(dá)到平衡后,旋渦恢復(fù)初始狀態(tài),不隨空泡尺度增大而改變。對于空泡閉合端的旋渦,即第2類旋渦,此旋渦是由于氣體從通氣孔流出后隨主流向后運(yùn)動,到達(dá)空泡尾部后閉合所形成。由圖10可以看出,其隨著空泡尺度的增大而增大,這是因為空泡尺度增大,空泡閉合端的主流強(qiáng)度減弱。當(dāng)?shù)?類旋渦增大到一定程度后會發(fā)生分裂,形成2個新的旋渦,同時,空泡增大到一定程度后,到達(dá)空泡尾部的氣體隨著反向射流進(jìn)入空泡內(nèi)部,與通氣孔后部出流氣體產(chǎn)生的旋渦相互作用,導(dǎo)致第1、2類旋渦之間的區(qū)域會出現(xiàn)新的小旋渦,兩種因素共同導(dǎo)致了第3類旋渦的生成,如圖10中t7時刻所示。

        3 結(jié)論

        本文采用CFX軟件對垂直發(fā)射的水下航行體進(jìn)行了三維數(shù)值模擬,并對改變通氣時序的影響進(jìn)行了分析,結(jié)果表明:

        1)通氣能提高泡內(nèi)壓力,通氣空泡內(nèi)不同位置的壓力值基本一致;通氣空泡末端存在回射高壓,航行體肩部存在低壓區(qū);隨著通氣時刻提前,終了時刻泡壓和回射壓力增大。

        2)通氣后空泡產(chǎn)生環(huán)狀凸起,環(huán)狀凸起向空泡尾部運(yùn)動,導(dǎo)致空泡長度增加。

        3)通入氣體在空泡內(nèi)有3種流向:向肩部運(yùn)動、隨主流向空泡閉合端運(yùn)動、形成反向射流;空泡內(nèi)旋渦分為3類:通氣孔附近、空泡閉合端、空泡中部。通氣孔附近旋渦不隨空泡增長而變化,空泡閉合端旋渦隨空泡增長而增大,進(jìn)而發(fā)生分裂,空泡中部的旋渦數(shù)目隨空泡增長而增加。

        [1] 黃壽康.流體動力-彈道-載荷-環(huán)境[M].北京:中國宇航出版社,1991.

        [2] 劉筠喬,魯傳敬,李杰,等.導(dǎo)彈垂直發(fā)射出筒過程中通氣空泡流研究 [J].水動力學(xué)研究與進(jìn)展, 2007,22(5):549-554.

        [3] Semenenko V N.Artificial supercavitation:physics and calculation[D].Ukraimian Academy of Sciences,Kiev Inst of Hydromechanics,2001.

        [4] Kawakami E,Arndt R E A.Investigation of the behavior of ventilated supercavities[J].Journal of Fluids Engineering,2011,133(9):091305.

        [5] Wang Z Y,Huang B,Wang G Y,et al.Experimental and numerical investigation of ventilated cavitating flow with special emphasis on gas leakage behavior and re-entrant jet dynamics[J].Ocean Engineering, 2015,108:191-201.

        [6] Karn A,Arndt R E A,Hong J.An experimental investigation into the physics of supercavity closure[J]. Journal of Fluid Mechanics,2015,789:259-284.

        [7] 賈力平,王聰,于開平,等.空化器參數(shù)對通氣超空泡形態(tài)影響的實驗研究[J].工程力學(xué),2007, 24(3):159-164.

        [8] 張敏弟,邵峰,郭善剛,等.繞空化器通氣空化流場的實驗研究[J].工程熱物理學(xué)報,2011,32(10): 1673-1676.

        [9] 郭建紅,魯傳敬,陳瑛,等.基于輸運(yùn)方程類空化模型的通氣空泡流數(shù)值模擬[J].力學(xué)季刊,2009, 30(3):378-384.

        [10] 時素果,王國玉,余志毅,等.FBM湍流模型在非定常通氣超空化流動計算中的評價與應(yīng)用 [J].船舶力學(xué),2012,16(10):1099-1106.

        [11] Launder B E,Spalding D B.The numerical computation of turbulent flows[J].Computational Methods in Applied Mechanics and Engineering,1974,3(2): 269-289.

        [12] Kubota A,Kato H,Yamaguchi H,et al.Unsteady structure measurement of cloud cavitation on a foil section using conditional sampling technique[J]. Journal of Fluids Engineering,1989,111(2): 204-210.

        Numerical Investigation of Ventilated Cavitating Flows around a Vertical Underwater-launched Projectile

        LIU Tao-tao,HUANG Biao,WANG Guo-yu,ZHANG Meng-jie
        (School of Mechanical Engineering,Beijing Institute of Technology,Beijing 100081,China)

        The numerical simulation of ventilated cavitating flow around an underwater-launched projectile is performed based on a homogeneous model.The simulated results are compared with experimental results in terms of the pressure.The influence of the timing variation of ventilation rate on ventilated cavity development and flow characteristics is analyzed.The research results show that the numerical results are consistent with the experimental results.Ventilation can improve the pressure inside the cavity and the pressure of different locations tends to be same.With the ventilation time earlier,the pressure inside the cavity and the high pressure of re-entrant jet increase gradually.In addition,the formation and evolution of vortex are captured in the process of cavity development.For gas flow characteristics,the gas injection moves towards three directions.The motion of gas leads to the formation of vortex and the vortex can be catalogued into three types.

        Numerical simulation;Ventilated cavitation;Timing variation of ventilated rate;Vorticity

        TV131.3

        A

        2096-4080(2017)04-0022-07

        2017-09-19;

        2017-11-01

        國家自然科學(xué)基金(51679005)

        劉濤濤(1989-),男,博士,主要研究方向為多相介質(zhì)復(fù)雜流動特性。Email:liutaotao_0708@126.com

        猜你喜歡
        通氣孔旋渦空泡
        一種耐高溫能力強(qiáng)的橡膠輸送帶
        汽車覆蓋件拉延模排氣截面積研究
        通氣孔對動葉冷卻結(jié)構(gòu)流動和換熱特性的影響
        水下航行體雙空泡相互作用數(shù)值模擬研究
        小心,旋渦來啦
        大班科學(xué)活動:神秘的旋渦
        旋渦笑臉
        山間湖
        筆帽安全通氣孔檢測研究
        基于LPV的超空泡航行體H∞抗飽和控制
        亚洲av精品一区二区| 亚洲av无码久久寂寞少妇| 99视频在线国产| 久久人妻av不卡中文字幕| 中文字幕人妻互换av | 国产精品中文久久久久久久 | 亚洲AV无码乱码1区久久| 日本激情久久精品人妻热 | 久久不见久久见免费影院国语| 中文亚洲av片在线观看不卡| 欧美亚洲尤物久久综合精品| 精品麻豆一区二区三区乱码| 国产激情无码视频在线播放性色| 国产成人av免费观看| 亚洲av乱码专区国产乱码| 国产三级不卡在线观看视频| 成品人视频ww入口| 国产色综合天天综合网| 中国女人a毛片免费全部播放| 日本一区二区三区爱爱视频| 中文字幕亚洲综合久久菠萝蜜| 亚洲人成网7777777国产 | 日韩精品综合在线视频| 中文无码av一区二区三区| 无码人妻精品一区二区三区不卡 | 亚洲男人第一无码av网站| 免费看欧美日韩一区二区三区| 少妇性l交大片免费1一少| 久久精品国产亚洲av果冻传媒| 国精产品一品二品国在线| 亚洲AV无码成人精品区H| 日本精品久久不卡一区二区 | 亚洲av日韩av不卡在线观看| 亚洲色偷偷偷综合网另类小说 | 日韩久久无码免费毛片软件| 天天爽天天爽天天爽| 一区二区三区婷婷中文字幕| 国产日产一区二区三区四区五区| www插插插无码视频网站| 精品亚洲国产探花在线播放| 极品少妇一区二区三区四区视频|