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

液體重定位推力時(shí)序的優(yōu)化研究

2018-03-21 08:02:08王麗霞周運(yùn)強(qiáng)
關(guān)鍵詞:優(yōu)化

劉 楨,王麗霞,林 宏,林 若,周運(yùn)強(qiáng)

?

液體重定位推力時(shí)序的優(yōu)化研究

劉 楨1,王麗霞2,林 宏3,林 若1,周運(yùn)強(qiáng)1

(1. 北京航天長(zhǎng)征飛行器研究所,北京,100076;2. 北京強(qiáng)度環(huán)境研究所,北京,100076;3. 北京宇航系統(tǒng)工程研究所,北京,100076)

研究了液體重定位的數(shù)值仿真方法,并利用半人馬座在美國(guó)路易斯落塔試驗(yàn)室的落塔試驗(yàn)結(jié)果驗(yàn)證該計(jì)算方法的準(zhǔn)確性,使用該數(shù)值仿真的方法對(duì)某貯箱液體重定位推力時(shí)序進(jìn)行優(yōu)化研究,研究結(jié)果表明通過對(duì)推力時(shí)序的優(yōu)化能夠加速液體重定位并能有效減少該過程需要消耗的能量。

微重力液體管理;液體重定位;推力時(shí)序優(yōu)化

0 引 言

飛行器在軌無(wú)動(dòng)力滑行時(shí),貯箱內(nèi)的液體推進(jìn)劑處于微重力環(huán)境而懸浮于貯箱中的任意位置。若飛行器是在低軌滑行,稀薄空氣的阻力作用甚至?xí)官A箱內(nèi)的液體推進(jìn)劑移動(dòng)到貯箱的上底處。此時(shí),飛行器要完成變軌、軌道調(diào)整或位置保持等動(dòng)作,必須在主發(fā)動(dòng)機(jī)啟動(dòng)前使用姿控發(fā)動(dòng)機(jī)推力將懸浮的液體推進(jìn)劑重新定位到貯箱下底的出口處。液體的這個(gè)流動(dòng)過程稱之為重定位過程,而利用這種方法實(shí)現(xiàn)在軌滑行飛行器多次啟動(dòng)的液體管理方法就是正推重定位管理方式[1]。

正推重定位管理方式中由于貯箱里沒有設(shè)置額外的推進(jìn)劑管理裝置(Propellant Management Devices,PMD),因此不存在推進(jìn)劑與結(jié)構(gòu)材料之間的兼容性問題,經(jīng)常應(yīng)用于長(zhǎng)時(shí)間滑行的低溫上面級(jí)液體管理,例如,美國(guó)通用上面級(jí)——半人馬座系列[2]等。

在正推重定位推力時(shí)序設(shè)計(jì)時(shí),總是希望盡量通過合理控制液體的慣性力,來(lái)抑制重定位過程中的液體動(dòng)能并加速其衰減,達(dá)到貯箱內(nèi)液體快速穩(wěn)定的目的。20世紀(jì)80年代末,美國(guó)提出脈沖式重定位的方 法[3,4],通過合理地對(duì)推力大小、推力作用時(shí)間和周期進(jìn)行優(yōu)化,減少了姿控系統(tǒng)推進(jìn)劑消耗量、提高有效運(yùn)載能力。

在地基上開展液體重定位試驗(yàn)研究需要建造很高的落塔設(shè)備來(lái)提供足夠的失重時(shí)間,研究周期很長(zhǎng)且成本高。而基于天基(如空間站等)的試驗(yàn)研究在試驗(yàn)準(zhǔn)備和試驗(yàn)操作上尚存在很多難以解決的問題。本文旨在對(duì)貯箱內(nèi)液體重定位過程采用有限體積數(shù)值仿真的方法來(lái)研究液體重定位過程推力時(shí)序設(shè)計(jì)的有效優(yōu)化。

1 重定位過程和數(shù)值仿真

液體重定位問題是一種無(wú)法用數(shù)學(xué)公式描述的復(fù)雜氣、液兩相流問題,要準(zhǔn)確地對(duì)貯箱內(nèi)液體重定位過程進(jìn)行數(shù)值仿真,必須先要了解重定位過程各階段的物理特征。

重定位過程中的姿控推力一般比較小,屬于低重力的過載環(huán)境范疇,此時(shí),表面張力對(duì)流動(dòng)作用還不能忽略。因此,液體重定位流動(dòng)的控制力主要包括姿控推力產(chǎn)生的慣性力、液體的黏性力和流動(dòng)相的表面張力等。

從流動(dòng)特征來(lái)看重定位過程可分為4個(gè)階段[5](見圖1):a)液體呈薄層狀從箱頂沿箱壁向箱底加速流動(dòng)過程;b)薄層狀流動(dòng)液體前沿在箱底匯合以“噴泉”狀涌起;c)“噴泉”消失,所有液體在箱底聚集晃動(dòng);d)箱內(nèi)液體基本平穩(wěn),液體中混夾的氣泡逐漸上升析出。

圖1 重定位過程

數(shù)值計(jì)算的初始狀態(tài),需先構(gòu)造盡量與實(shí)際狀態(tài)一致的流體區(qū)域,然后以某一特定條件(如“邦德數(shù)=15[6]”等)維持貯箱內(nèi)液面的低重形狀達(dá)到平衡,以保證重定位開始時(shí)所有液體沿箱壁向下流動(dòng)。由重定位過程特征可知,該問題是一個(gè)典型的氣、液兩相流問題,重定位流動(dòng)開始之后,在整個(gè)重定位流動(dòng)過程中箱內(nèi)氣枕不僅影響著前期薄層的流動(dòng)形態(tài),在中期氣枕分裂成氣泡、氣泡破碎成小氣泡以及后期的小氣泡合并成大氣泡析出都是液體重定位過程重要的組成部分。

因此,在數(shù)值計(jì)算過程中不能簡(jiǎn)單地只考慮表面張力和貯箱整體的氣液界面,可以采用絕熱氣泡模型模擬流動(dòng)過程中氣體內(nèi)壓力、位置的變化導(dǎo)致氣泡的體積變化以及分裂、合并析出等現(xiàn)象。重定位過程中液體流速比較低,此時(shí),可以忽略液體的湍流效應(yīng)而采用層流模型來(lái)考慮液體黏性力對(duì)運(yùn)動(dòng)的影響。

2 人馬座落塔試驗(yàn)和數(shù)值仿真的驗(yàn)證

人馬座重定位落塔試驗(yàn)研究[6]是在美國(guó)路易斯(Lewis)零重力設(shè)備上完成的,該落塔設(shè)備的自由落體高度約為142 m,可以提供大約5 s的失重時(shí)間。此零重力設(shè)備的自由落體高度處于世界領(lǐng)先水平,即便如此,其所能提供的失重時(shí)間與實(shí)際的人馬座貯箱重定位的時(shí)間還是不在一個(gè)量級(jí),而必須采用縮比模型進(jìn)行模擬試驗(yàn),通過犧牲模型的幾何尺度來(lái)延長(zhǎng)模擬試驗(yàn)的時(shí)間尺度。試驗(yàn)中用于模擬氧化劑箱和燃燒劑箱的模型貯箱半徑為5.5 cm和7.0 cm,長(zhǎng)度分別為 21.8 cm和27.8 cm,分別以Freon-TF和FC-78為試驗(yàn)液體,并通過反推噴氣系統(tǒng)調(diào)節(jié)下落時(shí)的推力來(lái)達(dá)到重定位初期以及重定位流動(dòng)開始之后的邦德數(shù)相似。

圖2為人馬座落塔試驗(yàn)和數(shù)值仿真對(duì)比。圖2a是半徑為7.0 cm,長(zhǎng)度為27.8 cm的人馬座縮比貯箱重定位試驗(yàn)結(jié)果。人馬座落塔試驗(yàn)初期使用邦德數(shù)=15維持貯箱內(nèi)液面的低重形狀。試驗(yàn)中典型重定位過程的前2個(gè)階段的特征非常明顯,由于邦德數(shù)較大,流體的慣性作用影響顯著,薄層狀流動(dòng)液體前沿在箱底匯合后涌起強(qiáng)健有力的柱狀“噴泉”,噴泉沖上貯箱的頂部后再沿箱壁往下流動(dòng)形成流動(dòng)循環(huán)。試驗(yàn)后期流體動(dòng)能減少逐漸使噴泉變細(xì),最后消失。由于落塔時(shí)間的限制,重定位的過程中噴泉消失后的晃動(dòng)以及液面基本平穩(wěn)之后的氣泡析出過程無(wú)法通過試驗(yàn)?zāi)M。圖2b是采用有限體積的數(shù)值計(jì)算方法對(duì)人馬座縮比貯箱重定位試驗(yàn)進(jìn)行數(shù)值仿真得到的在相同特征時(shí)刻的流型圖。

圖2 人馬座落塔試驗(yàn)和數(shù)值仿真對(duì)比

從圖2中可以看出,數(shù)值仿真的液體重定位過程與試驗(yàn)研究過程一致性很好,在同一特征時(shí)刻,數(shù)值仿真的液體的運(yùn)動(dòng)形態(tài)與試驗(yàn)結(jié)果非常一致。從而說明數(shù)值計(jì)算中采用的物理模型基本上符合液體重定位過程的物理特征。由于數(shù)值計(jì)算相比于試驗(yàn)研究方法來(lái)說具有不受自由落體時(shí)間限制的優(yōu)勢(shì),非常適用于液體重定位推力時(shí)序設(shè)計(jì)的優(yōu)化。

3 重定位推力時(shí)序優(yōu)化研究

在現(xiàn)役飛行器的實(shí)際應(yīng)用中,通常采用小推力發(fā)動(dòng)機(jī)推力讓液體沉底來(lái)減少液體在重定位過程對(duì)貯箱的反向沖擊,同時(shí)達(dá)到減少姿控系統(tǒng)燃料消耗的目的,等液體基本聚集在箱底后再采用大推力發(fā)動(dòng)機(jī)推力抑制貯箱內(nèi)液體晃動(dòng),加速氣泡析出,確保能給主發(fā)動(dòng)機(jī)啟動(dòng)輸送不夾氣的推進(jìn)劑。典型液體重定位推力時(shí)序如圖3所示[1]。

圖3 液體重定位典型推力時(shí)序

沉底段使用小推力發(fā)動(dòng)機(jī)需要的重定位時(shí)間長(zhǎng),一方面影響飛行器的機(jī)動(dòng)特性;另一方面也不一定能最大限度地減少姿控系統(tǒng)燃料的消耗。因此在沉底段考慮使用大推力發(fā)動(dòng)機(jī)加速重定位過程。但是大推力發(fā)動(dòng)機(jī)的過早介入,有時(shí)不僅達(dá)不到抑制貯箱內(nèi)液體晃動(dòng)的目的,還會(huì)使貯箱內(nèi)的液體晃動(dòng)加劇。因此,選擇變推力時(shí)間點(diǎn)讓大推力發(fā)動(dòng)機(jī)作用適時(shí)地介入到重定位過程中是時(shí)序優(yōu)化的關(guān)鍵。

在本次重定位推力時(shí)序優(yōu)化研究中,假設(shè)姿控系統(tǒng)僅具有300 N和600 N 2種正向推力的發(fā)動(dòng)機(jī),且發(fā)動(dòng)機(jī)不具有脈沖的工作能力。

試驗(yàn)參數(shù)包括:模型貯箱直徑為0.88 m;柱段長(zhǎng)度為0.57 m;四氧化二氮密度為1458 kg/m3;表面張力系數(shù)為25.61×10-3N/m;動(dòng)力粘度系數(shù)為0.4189×10-3kg/m?s。該模型參數(shù)在300 N推力時(shí)邦德數(shù)為150左右。

為了選擇合理的變推力點(diǎn),先采用300 N推力完成整個(gè)重定位過程,并提取計(jì)算過程中能夠表征液體運(yùn)動(dòng)幅度的軸向重心位置和液體動(dòng)能隨時(shí)間的變化曲線,如圖4所示。分析中統(tǒng)一采用動(dòng)能變化曲線穩(wěn)定低于0.0005時(shí)作為重定位運(yùn)動(dòng)結(jié)束的判據(jù)。

圖4 沉底段采用300N推力液體動(dòng)能和液體重心軸向變化

從圖4中可以看到,重定位過程大約耗時(shí)75 s。液體重心軸向變化反映液體的勢(shì)能變化,當(dāng)動(dòng)能具有峰值時(shí),重心位置基本上處于極小值位置,符合能量守恒規(guī)律。

優(yōu)化的推力時(shí)序,可以從以下4種方式考慮變推力時(shí)間點(diǎn):a)液體基本聚集在箱底后再采用大推力發(fā)動(dòng)機(jī)推力,這也是飛行器實(shí)際應(yīng)用中通常采用的方式。從圖4中可以看出,在44 s左右是重心位置的極小值點(diǎn),若在此刻變大推力有利于減小液體整體勢(shì)能,圖5為300 N正推重定位過程中變推力時(shí)刻的液體流型圖。此時(shí)刻液體流型圖如圖5d所示,液體已經(jīng)基本上都聚集在箱底;b)大推力發(fā)動(dòng)機(jī)提前介入,從減小整體液體勢(shì)能的角度考慮在重心變化曲線第1個(gè)波谷位置 (20 s)施加變大推力,此時(shí)的貯箱內(nèi)流型圖如圖5c所示,噴泉已經(jīng)涌起;c)在薄層狀流動(dòng)液體前沿在箱底匯合時(shí)刻(14 s)施加大推力,貯箱內(nèi)流型圖如圖5b所示,此時(shí)液體鋪滿在貯箱內(nèi)表面,重心位置位于貯箱的中心位置左右;d)作為對(duì)比分析,從一開始就采用600 N大推力沉底,重定位初始狀態(tài)設(shè)定為與 300 N狀態(tài)一致,起始狀態(tài)貯箱內(nèi)流型圖如圖5a所示。

圖5 300N正推重定位過程中變推力時(shí)刻的液體流型

2種恒定推力方式和3種變推力方式重定位過程中的動(dòng)能變化如圖6所示,時(shí)間和能量消耗的對(duì)比結(jié)果如表1所示。

圖6 各種重定位推力方式下的液體動(dòng)能變化

表1 各種重定位推力方式下的時(shí)間和能量消耗對(duì)比

Tab.1 Comparision of Time and Energy in Various Thrust Modes

參數(shù)推力方式重定位時(shí)間/s各工況與恒定300N時(shí)的重定位時(shí)間比值重定位消耗的能量/J各工況與恒定300N時(shí)消耗能量的比值 恒定300N751225001 恒定600N600.8360001.6 44s變推力720.96300001.33 20s變推力570.76282001.25 14s變推力500.67258001.15

從重定位的時(shí)間上考慮,恒定300 N小推力方式重定位需要的時(shí)間最長(zhǎng),但是,恒定600 N大推力重定位的時(shí)間卻不是最短的,這也驗(yàn)證了上文所述的大推力過早介入重定位帶來(lái)負(fù)面效應(yīng)的觀點(diǎn)。在重心變化曲線波谷位置,變大推力來(lái)減少液體整體勢(shì)能加速重定位過程的方式具有一定的效果,尤其是20 s變推力方式的重定位時(shí)間比恒定600 N大推力方式更少。而14 s變推力方式的重定位過程中,液體動(dòng)能在經(jīng)歷2次波動(dòng)后迅速下降并達(dá)到穩(wěn)定,其50 s的重定位時(shí)間是所有方式中最短的。

從能量消耗的角度分析,恒定300 N小推力重定位過程雖然重定位時(shí)間最長(zhǎng),但是消耗的能量卻是最少的。而恒定600 N大推力方式最大的能量消耗沒有換來(lái)最短的重定位時(shí)間。20 s變推力方式和44 s變推力方式在能量消耗方面表現(xiàn)優(yōu)于600 N大推力方式。14 s變推力方式由于重定位時(shí)間上的較大優(yōu)勢(shì),消耗的能量?jī)H僅為恒定300 N方式消耗能量的115%,優(yōu)于其它推力方式。

由于受網(wǎng)格大小的影響,計(jì)算過程中只能描述較大氣泡(大于3 cm)的析出情況,而更小氣泡的析出情況無(wú)法模擬。考慮到大推力時(shí)氣泡析出快于小推力時(shí)氣泡析出的特點(diǎn),在計(jì)算的重定位結(jié)束時(shí)刻,14 s變推力方式液體中的氣泡析出情況應(yīng)該優(yōu)于300 N推力方式,也就是說恒定300 N的推力方式要達(dá)到相同的氣泡的析出程度還需更長(zhǎng)的時(shí)間,說明14 s變推力的方式消耗的能量低于恒定300 N重定位消耗的能量。

綜合以上分析,在300 N的推力作用14 s之后,變600 N大推力的重定位推力時(shí)序不僅能讓主發(fā)動(dòng)機(jī)最快啟動(dòng),同時(shí)減少姿控系統(tǒng)的燃料消耗,是比較合理的重定位推力時(shí)序。通過對(duì)14 s變推力重定位時(shí)序的重定位過程中的速度云圖進(jìn)行分析,可以得到優(yōu)化變推力介入時(shí)間點(diǎn)的規(guī)律,這對(duì)重定位時(shí)序設(shè)計(jì)有重要的意義。

圖7是14 s變推力方式重定位過程中形成的噴泉時(shí)的軸向速度云圖。速度云圖以液體推力方向?yàn)檫\(yùn)動(dòng)正方向。

圖7 14s變推力方式的噴泉速度云圖

從圖7中可以看到,此時(shí)噴泉的上部已經(jīng)開始下落,而噴泉的下部則保持向上運(yùn)動(dòng),兩股液體的逆向運(yùn)動(dòng)在箱體中部形成對(duì)沖。從運(yùn)動(dòng)區(qū)域分布來(lái)看,對(duì)沖的兩股流動(dòng)動(dòng)量基本相當(dāng),方向相反。兩股流動(dòng)對(duì)沖混合之后,箱內(nèi)液體的總動(dòng)量將迅速降低,此刻大推力的介入將有效地抑制液體勢(shì)能和動(dòng)能之間的轉(zhuǎn)化。通過以上分析可知,當(dāng)重定位過程重心位于貯箱中心位置時(shí),施加大推力可以在第1次噴泉涌起之后的下墜過程中的一個(gè)時(shí)間點(diǎn)使箱內(nèi)的液體具有最小總動(dòng)量。

如本例計(jì)算結(jié)果所示,在重定位推力時(shí)序設(shè)計(jì)時(shí)將液體重心位于貯箱中心時(shí)刻設(shè)計(jì)為大推力介入時(shí)刻可以快速完成重定位并減小能量消耗。

4 時(shí)序優(yōu)化的應(yīng)用分析

以上計(jì)算分析工作是基于某一種貯箱特定的工況,分析在說明通過合理的優(yōu)化推力時(shí)序,可實(shí)現(xiàn)快速節(jié)能完成液體重定位的一種途徑。當(dāng)液體重定位過程重心位于貯箱中心位置時(shí),施加大推力可以平衡噴泉上涌和下落液體的動(dòng)量,推力時(shí)序優(yōu)化方法對(duì)于貯箱快速節(jié)能地完成液體重定位具有一般性。由于推力介入時(shí)刻的確定對(duì)推進(jìn)劑本身的物理屬性并不敏感[5],因此認(rèn)為該推力時(shí)序優(yōu)化方法也可以用于燃燒劑重定位過程。

由于液體重定位是復(fù)雜的非線性過程,液體重定位流動(dòng)形態(tài)與推進(jìn)劑的初始容積以及貯箱的尺寸形狀相關(guān),而各飛行器的推力系統(tǒng)組成也不盡相同。因此,在實(shí)際應(yīng)用中重定位推力時(shí)序優(yōu)化設(shè)計(jì)時(shí),本文中的推力時(shí)序優(yōu)化方法可以作為參考。

對(duì)于雙組元推進(jìn)劑系統(tǒng),燃燒劑和氧化劑的貯箱尺寸和充液比都存在差別,各貯箱液體重定位的特征時(shí)刻并不完全同步,需要多方面綜合權(quán)衡燃燒劑和氧化劑貯箱的具體情況制定重定位推力時(shí)序。若燃燒劑和氧化劑質(zhì)量相差較大時(shí),單從能量的角度來(lái)說以質(zhì)量較大的貯箱作為重定位推力時(shí)序優(yōu)化主要研究對(duì)象是比較合適的。

5 結(jié)束語(yǔ)

液體重定位由于其運(yùn)動(dòng)的復(fù)雜性,要通過試驗(yàn)手段比較完整地模擬整個(gè)過程,無(wú)論是通過地基還是天基試驗(yàn)設(shè)備都幾乎是不可能完成的。而不受過載和時(shí)間約束的數(shù)值分析方法逐漸成為液體重定位運(yùn)動(dòng)研究的重要手段。

本文通過數(shù)值分析方法,對(duì)多種推力方式的重定位效果進(jìn)行對(duì)比,證明了恒定的使用大推力或是小推力來(lái)完成重定位的方式都不是最理想的,通過合理的優(yōu)化推力時(shí)序,使用變推力的沉底方式,可以達(dá)到快速完成重定位同時(shí)又達(dá)到最小能量消耗的目的。

[1] 劉楨, 褚桂敏, 李紅, 王麗霞. 運(yùn)載火箭上面級(jí)微重力環(huán)境下的推進(jìn)劑管理[J]. 導(dǎo)彈與航天運(yùn)載技術(shù), 2012(4): 20-26.

Liu Zhen, Chu Guimin, Li Hong, Wang Lixia. Propellant management of rocket upper stage in microgravity environment[J]. Missiles and Space Vehicles, 2012(4): 20-26.

[2] 王永志, 王丹陽(yáng). 同步通信衛(wèi)星的發(fā)射[M]. 北京: 國(guó)防工業(yè)出版社, 2005.

Wang Yongzhi, Wang Danyang. The launch of synchronous communication satellite[M]. Beijing: National Defense Industry Press, 2005.

[3] Jaekle D E, Jr. Propellant management device conceptual design and analysis: traps and troughs[R]. AIAA 95-2531, 1995.

[4] Jaekle D E, Jr. Propellant management device conceptual design and analysis: Vanes[R]. AIAA 91-2172, 1991.

[5] 黃懷德, 褚桂敏, 張秀清. 振動(dòng)工程(下)[M]. 北京: 宇航出版社, 1995.

Huang Huaide, Chu Guimin, Zhang Xiuqing. Vibration engineering(Ⅱ) (M). Beijing: Aerospace Press, 1995.

[6] Jack A S, William J M, Raymond F L. Low-gravity reorientation in a scale-model centaur liquid-hydrogen tank[R]. NASA TN D-7166, 1973.

Thrust Time-series Optimization for Reorientation of Propellant

Liu Zhen1, Wang Li-xia2, Lin Hong3, Lin Ruo1, Zhou Yun-qiang1

(1. Beijing Institute of Space Long March Vehicle, Beijing,100076; 2. Beijing Institute of Structure and Environment Engineering, Beijing, 100076; 3. Beijing Institute of Astronautical Systems Engineering, Beijing, 100076)

In this paper, a numerical simulation method was adopted to study the reorientation of propellant in microgravity environment, and verify the model in correct using the result of drop tower test of the centaur in Lewis Zero-Gravity facility .Based on the numerical simulation method, investigate the thrust time-series optimization of the reorientation of propellant .The result show that, the thrust time-series optimization can accelerate the reorientation and reduce the energy consumption.

Propellant management in microgravity environment; Reorientation of propellant; Thrust time-series optimization

1004-7182(2018)01-0106-05

10.7654/j.issn.1004-7182.20180121

V23

A

2016-03-03;

2017-07-03

劉 楨(1982-),男,高級(jí)工程師,主要研究方向?yàn)橐后w動(dòng)力學(xué)

猜你喜歡
優(yōu)化
超限高層建筑結(jié)構(gòu)設(shè)計(jì)與優(yōu)化思考
PEMFC流道的多目標(biāo)優(yōu)化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設(shè)計(jì)優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
由“形”啟“數(shù)”優(yōu)化運(yùn)算——以2021年解析幾何高考題為例
圍繞“地、業(yè)、人”優(yōu)化產(chǎn)業(yè)扶貧
事業(yè)單位中固定資產(chǎn)會(huì)計(jì)處理的優(yōu)化
4K HDR性能大幅度優(yōu)化 JVC DLA-X8 18 BC
幾種常見的負(fù)載均衡算法的優(yōu)化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 久久无码高潮喷水| 日本福利视频网站| 婷婷在线网站| 丁香婷婷久久| 色婷婷在线播放| 亚洲精品无码日韩国产不卡| 国产拍揄自揄精品视频网站| 99久久这里只精品麻豆| 国产欧美日韩综合在线第一| 午夜免费小视频| 免费在线不卡视频| 九九视频免费在线观看| 日韩中文无码av超清| 欧美日韩中文国产| 欧美一区二区啪啪| 日本午夜在线视频| 欧美一级在线看| 亚洲精品国产首次亮相| 98超碰在线观看| 在线观看视频一区二区| 国产成人91精品| 国产理论最新国产精品视频| 粗大猛烈进出高潮视频无码| 国产免费黄| 2020久久国产综合精品swag| 无码日韩精品91超碰| 日本午夜影院| 国内精品一区二区在线观看| 久久天天躁狠狠躁夜夜2020一| 无码人妻热线精品视频| 无码区日韩专区免费系列| 国产亚洲精品97AA片在线播放| 中文字幕日韩久久综合影院| 手机看片1024久久精品你懂的| 欧美a级在线| 午夜免费小视频| 九色视频在线免费观看| 亚洲视频二| 91福利片| 国产十八禁在线观看免费| 伊大人香蕉久久网欧美| 亚洲无线观看| 欧美精品啪啪| 热这里只有精品国产热门精品| 亚洲热线99精品视频| 四虎国产精品永久一区| 色综合日本| 精品伊人久久久久7777人| 亚洲国产精品不卡在线| 日韩AV无码免费一二三区| 国产精品无码AV中文| 91成人在线观看视频| 无套av在线| 国产视频久久久久| 国产成人一区在线播放| 欧美日本在线| 亚洲天堂伊人| 欧美不卡视频在线观看| 特级做a爰片毛片免费69| 免费毛片全部不收费的| 久久国语对白| 欧美日韩免费观看| 四虎在线观看视频高清无码| 99久久这里只精品麻豆| 国产成人精品日本亚洲| 国产精品开放后亚洲| 中国一级特黄视频| 91丨九色丨首页在线播放| 人妻出轨无码中文一区二区| 国产成年女人特黄特色毛片免| 好久久免费视频高清| 伊人久久精品无码麻豆精品| 亚洲中久无码永久在线观看软件 | 欧美一级大片在线观看| 欧美日韩久久综合| 国模极品一区二区三区| 亚洲一级毛片免费观看| 在线免费看片a| 国内精品视频区在线2021 | 免费播放毛片| 国产精品99久久久久久董美香| 美女一级免费毛片|