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

彈道槍不同水深下全淹沒式發射膛口流場的數值分析*

2020-10-23 07:37:50張京輝余永剛
爆炸與沖擊 2020年10期

張京輝,余永剛

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

膛口流場一般指槍炮發射時,膛內氣體在膛口外膨脹而形成的隨時間變化的氣流區域[1],是非定常、帶有強激波間斷的復雜流場。眾多學者針對槍炮在空氣中發射的膛口流場進行了大量實驗和數值研究。姜孝海等[2]運用有限體積法,對空氣中膛口流場形成的動力學過程進行了數值模擬。吳偉等[3]基于無網格方法,對包含大位移運動邊界和非平衡化學反應的膛口流場進行了數值模擬。陳川琳等[4]結合實驗與數值模擬,研究了彈丸在膛口流場中的受力情況和運動規律。

而在水下槍炮發射過程中,由于水的物性特點與空氣大不相同,且不同水深下發射的環境壓力也會不同,這不僅使得水下航行體的運動阻力極大,而且會對水下槍炮發射所形成的膛口流場的演化特性產生復雜的影響,這對槍炮等身管武器在水下進行發射提出了極大的挑戰。在水下膛口流場方面,張欣尉等[5-6]運用Fluent 軟件,針對機槍在空氣中和水下密封式發射所形成的膛口流場進行了數值對比分析。水下膛口流場涉及氣液兩相流相互作用問題,圍繞燃氣射流與水的作用,許多學者開展了實驗與數值模擬。Hu 等[7]建立了燃氣射流在受限空間內擴展的三維非穩態數理模型,對貼壁燃氣射流在充液圓管內的擴展特性進行了數值模擬。Zhao 等[8]建立了多股燃氣射流數值求解模型,對錐形分布燃氣射流的干涉、匯聚特性進行了分析。在前人的基礎上,Zhou 等[9]對靜態和準動態條件下氣幕式發射的內彈道及氣液相互作用過程進行了數值模擬。而在水下超聲速射流方面,郝宗睿等[10]采用VOF (volume of fluid)兩相流模型結合動網格方法對系留狀態和運動狀態下的燃氣噴射流場進行數值模擬,分析了兩種狀態下燃氣噴射形成的復雜流場結構。唐云龍等[11]設計了一種解決水下超音速燃氣射流復雜相變過程的計算模型,并對典型工況的水下高溫、高速燃氣射流問題中相變過程進行了數值分析。張煥好等[12]建立了二維軸對稱水下超聲速氣體射流的數值計算模型,并進行了相關的數值模擬,得到了水下超聲速氣體射流的初始流動結構。

已有研究主要針對空氣中膛口流場和水下燃氣射流問題,但對水下全淹沒發射(即身管內充滿水),尤其是不同發射水深條件下的膛口流場演化特性的研究較少。本文中,建立二維軸對稱非穩態膛口流場模型,對水下全淹沒發射膛口流場演變全過程進行數值模擬,并搭建水下可視化射擊實驗平臺,利用高速錄像觀測12.7 mm 口徑彈道槍在水中全淹沒式發射時膛口流場演變的過程。在此基礎上,針對12.7 mm 口徑彈道槍在3 種水深下全淹沒發射膛口流場演化過程進行數值模擬,以期研究結果對深入了解水下膛口流場特性以及降低膛口氣流危害有所幫助,同時對新型水下槍炮設計有一定的參考價值。

1 數理模型與計算方法

1.1 物理模型

根據彈道槍水下全淹沒式發射特點,采用下列物理模型:

(1)火藥燃燒遵循幾何燃燒定律。藥粒均在平均壓力下燃燒,且遵循指數燃速定律。單位質量火藥燃燒所放出的熱量及生成的燃氣溫度均為定值,在以后的膨脹做功過程中,不考慮燃氣組分變化,火藥力fp、余容αp及比熱比kp等均視為常數。

(2)槍管內膛表面的熱散失用減小火藥力fp或增大比熱比kp的方法間接修正,用系數φ來考慮其他的次要功。

(3)考慮到在膛口流場中氣相密度遠小于液相密度,重力場影響很小,以二維效應為主,因此本文中采用二維軸對稱模型對膛口流場進行研究,湍流模型采用標準k-ε 模型。

(4)將火藥燃氣看作理想氣體,不考慮燃氣在膛口的二次燃燒,忽略體積力的影響。

(5)假設水為不可壓縮相,密度取998.2 kg/m3。因膛口燃氣射流與水的作用時間短暫,忽略燃氣射流對水的汽化作用。

(6)對于膛口以及彈丸表面附近水的空化現象,采用Schnerr-Sauer 空化模型進行計算。

1.2 數學模型

針對上述物理模型建立下列數學模型:

(1)連續性方程:

式中: Γe和Γc分別為氣泡生長和潰滅時質量傳遞源項; αv為空泡的體積分數,c為單位體積內的氣泡數,ρv和 ρl分別為汽相和液相的密度, ρ 為混合相的密度,pv為蒸汽分壓力,Rb為氣泡的半徑。氣泡半徑的表達式為:

膛口流場的計算中需要耦合內彈道方程組,內彈道方程如下:

(1)形狀函數:

式中:ψ 為火藥燃燒質量百分比,χ、λ 和μc均為火藥形狀函數,Z為火藥燃燒相對厚度,Zk為當顆粒燃盡時燃燒相對厚度, χs和 λs為在火藥表面縮小階段的形狀函數。

(2)燃速方程:

式中:u1為火藥燃速系數,n為火藥燃速指數,e1為火藥半弧厚,p為膛內火藥燃氣平均壓力,t為時間。(3)彈丸運動方程:

式中:pd和pf分別為彈底和彈前壓力,A為彈丸橫截面積,φ為次要功系數,m為彈丸質量,v為彈丸運動速度。

(4)內彈道基本方程:

式中:lψ為藥室自由容積縮徑長,絕熱指數 θ =kp?1 ,w為裝藥質量,x為彈丸行程,Δ為裝填密度,αp為單位質量氣體分子體積有關的修正量,稱為余容。

(5)彈丸速度與行程關系式:

方程(12)~(17)組成內彈道方程組。

1.3 計算方法

本文中借助FLUENT 軟件開展數值模擬,結合用戶自定義函數(user defined function,UDF)技術,耦合內彈道方程組求解得到彈丸速度、位移和膛壓;采用流體體積函數(volume of fluid,VOF)多相流模型來描述氣液相互作用;利用PRESTO!方法對壓力項離散;動量和能量的離散采用一階迎風格式;采用PISO 算法對壓力與速度進行耦合;計算采用的時間步長控制在0.1 μs 內。

2 網格劃分與無關性驗證

2.1 網格劃分

如圖1 所示,計算域分為身管區域和膛口流場區域兩個部分,其中身管區域為長1 m、半徑為6.35 mm的圓柱形區域,流場區域為長1 m、半徑0.3 m 的圓柱形區域。以結構化網格為主,對膛口附近的網格進行局部加密,最小網格尺寸為0.3 mm×0.3 mm。采用了動網格技術中的層鋪法來實現彈丸的運動,隨著彈丸向前運動,當靠近彈底的網格長度被拉長至0.6 mm 時,網格會分裂成兩個網格,同理,當靠近彈頭的網格長度被壓縮至0.2 mm 時,該網格會與相鄰的網格合并。

本文計算的邊界條件為:燃燒室為壓力入口,其初始參數由內彈道方程組求解得出;身管壁面及彈丸為固壁邊界;計算域外邊界為壓力出口,壓力為環境壓力。將彈前身管區域和流場區域初始化為充滿水,溫度和壓力初始化為環境條件,即溫度為300 K,壓力根據不同水深分別設置為111.5、608.0 和1 114.6 kPa。

2.2 網格無關性驗證

為盡量避免網格粗細程度對計算結果產生較大影響,分別采用20 萬、15 萬和10 萬的網格數(N)進行試算,選擇彈丸出膛后1 ms 時刻膛口中心軸向速度分布作為網格無關性驗證的特征參數,無關性驗證結果如圖2 所示。相對于20 萬網格數的計算結果,15 萬網格數的最大相對誤差僅為6%,而10 萬網格數的最大相對誤差達到30.6%。綜合考慮計算效率和計算精度,現選擇網格數為15 萬的網格進行計算。

圖1 計算域及邊界條件Fig. 1 Calculation domain and boundary conditions

圖2 不同網格數下膛口中心速度沿軸向的變化(t=1 ms)Fig. 2 Velocity of the muzzle center varying along the axial direction for different grid quantities at t=1 ms

3 模型驗證

為觀察水下全淹沒式膛口流場中氣、汽、液三相的摻混過程以及驗證數值模型的合理性,以12.7 mm口徑的彈道槍為平臺,搭建了全淹沒式發射的可視化實驗測試系統,如圖3 所示。實驗系統由水箱、12.7 mm 口徑彈道槍、水泥平臺以及脈沖電點火器組成,其中水箱上開有射擊孔,供身管伸入,水箱側面開有觀察窗。采用FASTCAM-Ultima APX 高速攝像機和計算機進行流場演化過程的觀察和記錄,其中高速攝像機最大分辨率為1 024 ×512 像素,頻率為4 000 s?1。

圖3 實驗裝置示意圖Fig. 3 Schematic diagram of the experimental device

在實驗準備階段,將彈丸填充在槍膛底部,通過彈帶進行密封以阻止水進入彈后燃燒室;為保證全淹沒發射環境,將身管以及水箱內注滿水;調整高速攝像拍攝的焦距與角度,使其能透過觀察窗清楚的拍攝到膛口流場區域;當儀器調整好后,將稱量好的火藥裝入燃燒室,裝好電底火,并將其連接到脈沖電點火器上。所有裝置就緒后,用脈沖電點火器放電進行點火。實驗后,用計算機將高速攝像拍攝的結果進行截取和后處理。

本次模型驗證所選取的實驗工況水深為0.3 m,對應的數值模擬工況水深也為0.3 m。采用的實心彈丸質量為45 g,裝藥采用6 g 的4/7 單基火藥,實驗測得彈丸初速為264.2 m/s。實驗陰影照片與數值模擬相圖的對比情況如圖4 所示,將彈丸出膛口時刻定義為t=0 時刻。

圖4 實驗照片陰影圖與模擬相圖對比情況Fig. 4 Comparison of experimental shadow photos and simulated phase diagrams

由圖4 可見,在全淹沒發射條件下,身管內與外界環境都充滿了水。圖4(a)~(b)對應點火后,彈丸在膛內如同活塞般不斷推動彈前水柱的過程,其中圖4(a)對應彈丸出膛前0.8 ms 時刻。水柱在彈丸推動下不斷加速,使得膛口附近的局部壓力低于水的飽和蒸汽壓力,膛口附近的水不斷汽化為水蒸汽,在膛口逐漸堆積成一個蒸汽團。圖4(b)對應彈丸即將出膛口瞬間,彈丸頭部被側面的蒸汽團完全包裹。圖4(c)~(e)對應彈丸逐漸飛離膛口的過程。在圖4(c)中,高壓火藥燃氣由膛內噴出,與蒸汽團迅速摻混,初步形成了泰勒空腔,同時提高了膛口附近的壓力,使得蒸汽團迅速潰滅。隨著彈丸運動到射流頭部,對射流在軸向上的約束減弱,使射流在軸向上迅速擴展,而在徑向上擴展仍受周圍水介質的限制,因此頭部呈錐形,同時氣液之間的強烈湍流摻混引起的Kelvin-Helmholtz 不穩定效應使得射流邊界變得不規則。彈丸穿過泰勒空腔后相對水介質做高速運動,使得彈丸頭部表面的水的壓強降低,初步形成自然空泡。在圖4(d)和(e)中彈丸已完全超越射流,彈丸出膛1.2 ms 后,隨著彈丸繼續向前飛行,彈丸表面的空泡長度逐漸變長,形成超空泡,包裹著彈丸前進,同時在彈丸尾部會留下細長的氣柱。

從圖4 中可以看出實驗陰影照片與數值模擬相圖中的膛口流場的外部輪廓基本吻合,為進一步驗證數值模型的合理性,做出了射流軸向、徑向最大位移對比圖,如圖5 所示。高速錄像在進行拍攝時,由于光路傳播原因,照片上水箱前后透明觀察窗固定螺釘的間距為等比例縮小。在實驗中,槍管位于水箱對稱面位置,根據梯形中位線特征,在進行膛口數據處理時,前后觀察窗的螺釘間距可以直接從照片中量取,兩者的平均值與實際螺釘間距的比值即為量取長度與實際長度的比值,再結合照片中量取的射流軸向、徑向最大位移從而得到實際射流的最大位移。為了減小誤差,對每種工況實驗數據進行了多次測量求平均值。由圖5 可知,模擬結果中不同時刻燃氣射流的軸向、徑向最大位移與實驗測量值吻合較好,最大相對誤差為6.4%,說明本文中采用的數值模型基本合理。

圖5 射流軸向、徑向最大位移對比Fig. 5 Comparison of the maximum axial and radial displacements of the jet flow between experimental and simulated results

4 結果與分析

為研究發射水深對全淹沒式發射膛口流場的影響,針對12.7 mm 口徑彈道槍在3 種水深h下全淹沒發射膛口流場演化過程進行了數值模擬和對比分析。3 種工況的初始參數如表1所示,表中l為身管長度,v0為彈丸初速,pk0為膛口壓力。圖6 給出了不同水深條件下彈丸的膛外行程對比曲線。

表1 膛口初始參數Table 1 Initial parameters for the muzzle

圖6 不同水深下彈丸膛外行程Fig. 6 Displacement-time curves of the projectile at different water depths

由表1 可見,隨著水深增加,彈丸初速降低明顯,膛口壓力也隨之升高。由圖6 可見,水深越深,彈丸速度越慢。對膛口流場影響范圍內,彈丸行程隨時間變化特性進行擬合,發現其在不同水深條件下均滿足指數函數規律,即:

式中:xt為彈丸膛外行程,m;t為時間,ms;x0、x1和t1為彈丸膛外行程隨時間變化的擬合參數,如表2 所示。由表2 可見,x0、x1和t1都會隨著水深增大而減小。

表2 擬合參數Table 2 Fitting parameters

全淹沒發射條件下,由于膛內的火藥燃氣處于高溫高壓狀態,當燃氣從膛口噴出后,會迅速膨脹與高密度水介質以及彈丸相互作用,從而形成復雜的波系結構。為了便于對不同發射環境條件下的膛口流場流譜結構進行對比分析,圖7 給出了空氣中槍炮發射時的膛口流場流譜示意圖[1], 圖8 給出了全淹沒發射時膛口流場流譜示意圖。為研究水深對膛口流場演化特性的影響,圖9~11 分別給出了3 種水深條件下的膛口中心剖面的馬赫數云圖和等值線圖。

圖7 空氣中發射時典型膛口流場流譜示意圖Fig. 7 Schematic diagram of the typical muzzle flow field during air launch

圖8 全淹沒發射時膛口流場流譜示意圖Fig. 8 Schematic diagram of the muzzle flow field during submerged launch

由圖7 可見,空氣中膛口流場將射流理論邊界包圍的整個流場區域劃分為了5 個區域,分別為:(1)區,即瓶狀激波內的自由膨脹區,氣流主要在此瓶狀區內膨脹,壓力劇降、速度激增;(2)區,即相交激波與射流邊界之間的超聲速區;(3)區,即馬赫盤下游的亞聲速區,氣流經過正激波后速度劇降、壓力和溫度陡增;(4)區,經過兩次斜激波后,流動情況變得復雜,壓力與(3)區相同,但為超聲速區,兩區之間形成類似于拉瓦爾噴管形狀的切向間斷面;(5)區,即湍流混合過渡區。

由圖9~11 可見,水下全淹沒式發射時,由于燃氣射噴出膛口后迅速膨脹,膨脹波系在氣液邊界以及彈底發生反射形成激波,約在t=0.4 ms 時刻,演化成了如圖8 所示的波系結構。類似于空氣中膛口流場典型結構,入射激波、馬赫盤以及反射激波交匯形成了三波點,通過正激波后的亞聲速氣流與經過斜激波后速度仍為超聲速的氣流之間有速度差,因此在三波點沿流向形成切向間斷。瓶狀激波雖然已經形成,但受高密度水介質的影響,流場波系結構的分布特點存在一些不同:(1)區(自由膨脹區)受氣液邊界的影響,形狀尚不飽滿,入射斜激波在徑向上向內收縮,馬赫盤在軸向上向外突出且直徑較小;受氣液邊界和(1)區影響,(2)區被擠壓呈狹長狀被未能充分發展至將(1)區包裹起來;(3)區也受高密度水介質的壓縮影響,在徑向上大大被壓縮。在后續的發展中(t>0.4 ms),隨著彈丸完全脫離泰勒空腔,對膛口近場的波系結構不再有影響,同時膛內后續氣體不斷噴出,使得泰勒空腔在徑向上擴展明顯,對瓶狀激波在徑向上的限制減小。整個瓶狀激波系逐漸飽滿,最終形成如圖7 所示的膛口流場典型波系結構。

從圖9~11 的對比中可以看出,在3 種水深下,彈丸出膛后0.4 ms,包括入射激波、反射激波、馬赫盤和三波點等結構的瓶狀激波都已經形成。當水深為1 m 時,在t=0.4~1.4 ms 期間,瓶狀激波在徑向上趨于飽滿,馬赫盤直徑也相應增大,最終在t=1.4 ms 形成與如圖所示的膛口流場典型波系結構。而在水深為50 m 時,在t=0.4~1.4 ms 期間,瓶狀激波在徑向上有一定擴展,馬赫盤沿軸向逐漸向下游移動,使得瓶狀激波的“瓶口”仍向外突出。在t=1.4 ms 時,馬赫盤已經到達距膛口最大位移并逐漸向膛口方向移動,在t=1.8 ms 時,瓶狀激波才趨于飽滿,形成膛口流場典型波系結構。而在水深為100 m 的條件下,在t=2.4 ms 時,流場內才形成膛口流場典型波系結構。

圖9 h=1 m 時膛口中心剖面馬赫數云圖與等值線圖Fig. 9 Mach number cloud map and contour map at h=1 m

圖10 h=50 m 時膛口中心剖面馬赫數云圖和等值線圖Fig. 10 Mach number cloud map and contour map at h=50 m

圖11 h=100 m 時膛口中心剖面馬赫數云圖和等值線圖Fig. 11 Mach number cloud map and contour map at h=100 m

這說明相同裝藥條件下,水深越大,對應的彈丸初速越低,脫離膛口流場區域所需的時間越長,對流場內的波系結構演化的限制時間更長,同時也由于水深處的環境壓力較大,使得水深越大,瓶狀激波在徑向上需要更長時間趨于飽滿,形成膛口流場典型波系結構。

為研究流場內壓力變化特性,圖12 給出了水深為50 m 的條件下膛口中心剖面的壓力云圖,圖13 則給出了不同水深下,壓力沿膛口中心軸線的分布曲線和瓶狀激波內(x=1.02 m 處)壓力沿徑向的分布曲線。

圖12 h=50 m 是膛口中心剖面壓力云圖Fig. 12 Pressure cloud map of the muzzle center section at h=50 m

由圖12 可見,彈丸出膛后,膛內的火藥燃氣噴出膛口,迅速膨脹,在t=0.4 ms 已經形成了包含瓶狀激波在內的復雜激波結構,從而影響了流場內的壓力分布。瓶狀激波內是燃氣射流的主要膨脹區,射流在瓶區內迅速膨脹,使得瓶區內的壓力劇降。而在馬赫盤的下游,由于射流受到周圍高密度水介質與彈底的壓縮作用,使得壓力驟升,形成高壓區域。隨后由于燃氣自身被壓縮后的膨脹作用,以及受到的壓縮作用,使得馬赫盤后高壓區域與低壓區域交替出現。在t=1.4 ms 后,由于彈丸遠離膛口區域,以及后續燃氣噴出,使得馬赫盤下游的壓力趨于平均,不再有突變,而自由膨脹區內的壓力仍然較低,瓶狀激波結構仍清晰可見。

圖13 膛口壓力分布曲線Fig. 13 Pressure distribution curves of muzzle

由圖13(a)~(b)可見,在水深為1 m 時,在彈丸出膛0.4 ms 時刻,此時彈丸仍處在膛口近場區域內,燃氣從膛口噴出后迅速膨脹,在軸向上壓力降至約0.2 MPa。隨后由于水介質與彈底的壓縮作用,壓力在馬赫盤處分別驟升至1.6 MPa。由于射流自身被壓縮后的膨脹作用,以及水介質與彈底的壓縮作用,流場內的壓力在馬赫盤至彈底區間內多次震蕩,最終在彈底驟升至約60 MPa。而在徑向上,燃氣壓力在小幅度下降后,于瓶狀激波的斜激波處壓力驟升至1.3 MPa,隨后在斜激波至氣液邊界區間內震蕩,最終緩慢下降。

隨著彈丸逐漸遠離膛口(t=1.0 ms),膨脹區內的最低溫度與最低壓力有所降低,馬赫盤處的溫度與壓力峰值也相應減小。在馬赫盤下游,壓力在軸向上的震蕩幅度有所減弱,但仍十分明顯,而在徑向上,壓力在斜激波處到達壓力峰值后呈緩慢下降趨勢,不再震蕩。彈丸出膛1.4 ms 后,馬赫盤后壓力的震蕩幅度進一步減弱,幾乎消失,壓力呈緩慢下降趨勢。

對比不同發射水深下的壓力曲線,可以發現,在軸向上,在彈丸出膛后相同時刻,相對于1 m 水深時,50 和100 m 水深發射時,在軸向上,射流出膛口后在馬赫盤處的壓力峰值偏小,且馬赫盤后壓力震蕩程度大大減弱,震蕩次數也相對減少,在彈丸出膛1.4 ms 后,馬赫盤后的壓力幾乎不再震蕩,沿程緩慢下降。在徑向上,由于環境壓力的增大,與1 m 水深時不同,50 和100 m 水深發射時,彈丸出膛1.0 ms 后,燃氣壓力于瓶狀激波的斜激波處壓力驟升后,在斜激波至氣液邊界區間內仍有較大幅度的震蕩,直到t=1.4 ms,壓力才趨于平穩,幾乎不再震蕩。

整體而言,在不同發射水深條件下,燃氣射流出膛后迅速膨脹,壓力在軸向與徑向上迅速降低,隨后分別在馬赫盤與斜激波處驟升。隨后軸向與徑向壓力都多次震蕩,且震蕩幅度隨時間逐漸減小,最終呈緩慢下降趨勢。當發射水深越深時,在軸向上,燃氣在馬赫盤處的溫度與壓力峰值相應越低,壓力震蕩程度也越小,且更快趨于平穩,而在徑向上,壓力的震蕩持續的時間更長,需要更長時間趨于平穩。

5 結 論

針對彈道槍水下全淹沒發射膛口流場問題建立了二維軸對稱非穩態模型,借助Fluent 軟件,結合動網格和用戶自定義函數技術,對3 種水深下的膛口流場演變過程進行數值模擬,得到如下結論:

(1)全淹沒發射條件下,在彈丸出膛前,身管內水柱在彈丸推動下不斷加速,在膛口空化形成蒸汽團;高壓火藥燃氣由膛內噴出后,與蒸汽團迅速摻混,初步形成了泰勒空腔;穿過泰勒空腔后,彈丸頭部初步形成自然空泡;隨著彈丸繼續向前飛行,彈丸表面形成超空泡,包裹著彈丸前進;實驗結果與數值模擬結果吻合較好,驗證了數值模型的合理性。

(2)同空氣中膛口流場典型波系結構相比,全淹沒發射條件下,彈丸出膛后0.4 ms 瓶狀激波雖已形成,但在徑向上尚不飽滿;在后續發展中,瓶狀激波在徑向上逐漸擴展,最終形成膛口流場典型波系結構,且水深越深,形成膛口流場典型波系結構所需時間越長。

(3)在膛口流場影響范圍內,彈丸膛外行程隨時間在不同水深條件下均滿足指數函數規律,即:x(t)=x0+x1e?t/t1,且水深越深燃氣在膛口軸向馬赫盤處的溫度與壓力峰值越低,壓力振蕩幅度也越小,更快趨于平穩,但在徑向上,水深越深,壓力振蕩的持續時間越長。

主站蜘蛛池模板: 2022国产无码在线| 99er精品视频| 国产成人综合在线观看| 久久这里只有精品23| 一级做a爰片久久免费| 亚洲熟女中文字幕男人总站| 欧美黄网站免费观看| 日韩欧美中文在线| 国产在线拍偷自揄观看视频网站| 国产xxxxx免费视频| 91蜜芽尤物福利在线观看| www.99在线观看| AV在线麻免费观看网站 | 亚洲AⅤ无码国产精品| 深夜福利视频一区二区| 国产专区综合另类日韩一区| 国产自视频| 色综合久久无码网| 国产精品永久不卡免费视频| 国产三级国产精品国产普男人| 精品国产自在现线看久久| 在线毛片免费| 人人91人人澡人人妻人人爽| 日韩美毛片| 久久久91人妻无码精品蜜桃HD| 亚洲swag精品自拍一区| 国产美女在线观看| 国产精品色婷婷在线观看| 日韩午夜片| 国产精品午夜福利麻豆| 91蝌蚪视频在线观看| 91精品国产自产在线观看| 国产香蕉在线视频| 野花国产精品入口| 亚洲人成高清| av一区二区三区高清久久| 欧美日韩第二页| 青草娱乐极品免费视频| 国产精品无码翘臀在线看纯欲| 九色在线视频导航91| 91福利片| 老司机精品一区在线视频| 国产男女XX00免费观看| 在线看片国产| 亚洲精品国产精品乱码不卞| 国产91久久久久久| 精品综合久久久久久97超人| 亚洲国产欧美中日韩成人综合视频| 欧美精品在线免费| 啪啪永久免费av| 久久久久久久久亚洲精品| 黄色免费在线网址| 亚洲av综合网| 免费看a级毛片| 精品国产免费第一区二区三区日韩| 91无码视频在线观看| 亚洲天堂日韩在线| 久久香蕉国产线看精品| 天天躁狠狠躁| 国产精品美女自慰喷水| 免费午夜无码18禁无码影院| 国产一区二区三区在线精品专区| 亚洲成肉网| 五月激激激综合网色播免费| 四虎成人在线视频| 亚洲伊人天堂| 五月婷婷伊人网| 亚洲 成人国产| 日本欧美在线观看| 99中文字幕亚洲一区二区| 国产日韩欧美在线播放| 欧美、日韩、国产综合一区| 伊人成色综合网| 日韩精品成人网页视频在线| 国产精品免费久久久久影院无码| 日韩成人在线网站| 亚洲精品国偷自产在线91正片| 国产精品无码作爱| 丁香综合在线| 亚洲第一区欧美国产综合| 日本成人不卡视频| 中文字幕亚洲第一|