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

蛟龍號載人潛水器無動力潛浮運動分析系統開發

2012-06-07 10:23:46潘彬彬崔維成劉正元
船舶力學 2012年1期
關鍵詞:深度

潘彬彬,崔維成,葉 聰,劉正元

(中國船舶科學研究中心,江蘇 無錫 214082)

1 引 言

基于氣球原理設計的第一代和第二代載人潛水器,均是采用無動力下潛和上浮的方式[1-3]。這樣做獲得的最大好處就是可以大大地節省潛水器上的能源,但需要解決無動力潛浮過程中的水動力學問

題[4]。

在“蛟龍”號設計階段,劉正元等就開始研究載人潛水器的無動力下潛上浮運動[5-6],給出了一套簡化分析方法并應用于“蛟龍”號載人潛水器的無動力潛浮設計。為了更精確和方便地進行海試和今后應用下潛的配載,我們在“蛟龍”號載人潛水器3 000米級海試數據的基礎上,曾經建立了一個理論上更為精細的載人潛器無動力上浮下潛運動模型,并根據該模型開發了一個運動分析系統[7],用于預報“蛟龍”號5 000米級海試的配載和速度。但是由于3 000米海試所有潛次的海試數據中都存在開動推進器和調節壓載水艙的操作,沒有真正意義上的無動力上浮下潛數據用于標定模型參數,因此,所建立的無動力潛浮運動分析系統并不十分可靠。在2011年7-8月間進行的5 000米級海試過程中發現該模型預報值比實際值偏大了約80kg。在這次的5 000米級海試中,D44潛次首次實現了全程無動力上浮下潛運動,為改進載人潛器無動力上浮下潛運動模型提供了寶貴的數據。本文將對D44潛次的數據進行分析,在此基礎上改進潛器無動力上浮下潛運動模型,期望開發出真正可靠的“蛟龍”號載人潛水器的無動力潛浮運動分析系統,為2012年的7 000米級海試以及今后的應用提供可靠的手段,同時,本文所建立的精細分析方法也可以應用于我國正在進行的4 500米載人潛器無動力潛浮運動的設計。

2 載人潛水器無動力潛浮運動模型

在文獻[7]中詳細給出了載人潛水器的無動力潛浮運動分析精細模型,在該模型中,認為無動力潛浮運動過程中潛器所受的外力只有重力W、浮力B及流體阻力Fd。此外考慮到穩態運動時其角速度為零等因素,模型建立在固定坐標系內,如下圖1所示。

潛器在固定坐標系中oz0方向的力平衡方程為:

其中M=潛器的質量與z0方向附加質量之和;0=潛器垂向(z0方向)加速度;W=重力;B=浮力;Fd=潛器在鉛垂方向運動時受到的 (z0方向)的阻力。潛器無動力下潛上浮的大部分時間中潛器的運動可以近似為平穩運動,即近似認為平穩階段潛器的加速度為零:

圖1 坐標系Fig.1 Coordinate system

但(2)式中各部分力與文獻[7]中的模型有所變化,所以本文將對每個部分力重新介紹。

2.1 阻力項Fd

如圖1所示,潛器運動時存在一定的縱傾角θ,文獻[6]推導了潛器速度方向(攻角)α和縱傾角θ的關系為:

其中:Cx(α)為潛器縱向(x軸方向)運動阻力系數,Cz(α)為潛器垂向(z軸方向)運動阻力系數,兩個方向的阻力試驗結果見文獻[7]。

3 000 米和5 000米級海試的數據顯示縱傾角θ接近零度,而由上式(3)可知,當縱傾角θ接近零度時攻角α可表示為θ的函數,采用數值算法可近似求出該函數,與文獻[7]類似,阻力項寫為:

0

2.2 浮力項B

隨著潛器所處深度不同,海水密度和壓力都會不同。一方面密度的變化會使潛器浮力發生相應變化;另一方面潛器的各耐壓部件、密封部位和壓力補償部件的排水體積也會隨著海水壓力變化而變化。所以為了較準確地計算潛水器在各個深度的浮力,需要采用對應深度處的密度代入計算,而且需要較準確地計算潛器排水體積的變化量,潛器的浮力項可寫為:

B=ρg (V-ΔV)

(5)其中:ρ=ρ (depth)為深度depth處的海水密度;ΔV為深度depth處潛器排水體積減少量,ΔV的變化分為兩個大的階段,第一階段是潛器壓載水箱注水或者排水階段,此階段ΔV變化速度較快;第二階段是潛器壓載水箱注滿水后的階段,此階段潛器的排水體積變化主要是由于海水壓力導致的,本文主要分析第二階段的ΔV;所以V為潛器的壓載水箱注滿水后的排水體積。

根據水池均衡實驗結果[7],潛器在水池平衡時總重20 807kg,而此時除了潛器本身的排水體積外,還有27kg鉛塊(密度11 700kg/m3)。此時潛水器在水池中平衡,有:

其中:Vm=潛器本體的排水體積。由(6)式可求得:

除了潛器本體,潛器攜帶的固定壓載、可拋壓載和作業工具等也產生一定的排水體積,這部分排水體積稱為可變排水體積Vc。則潛器剛潛入水面以下時的排水體積(未考慮體積變化量)為:

對于排水體積變化量ΔV,如前所述,在海水壓力下潛器排水體積發生變化的主要是耐壓結構、密封部位和壓力補償部件,下面分兩部分進行分析。

2.2.1 耐壓結構排水體積變化

“蛟龍”號載人潛器上的主要耐壓結構按外形可以分為兩類:

(1)球殼:載人艙、可調壓載水艙、高壓氣罐;

(2)圓柱殼:計算機罐、水聲通信機罐、測深側掃聲納罐、配電罐。

根據設計規范,可知潛器的耐壓結構主體的材料在工作壓力下還屬于彈性范圍內,所以可以采用理論公式計算排水體積變化量。

2.2.1.1 球殼耐壓結構

表1 “蛟龍”號主要球殼耐壓結構Tab.1 Large spherical pressure hull of “JiaoLong”

球殼在下潛開始時內壓已經存在,此時球殼的排水體積為:

其中Δr0為球殼徑向變形量:

(10)式中:P1=球罐的內壓,Ri=球殼內半徑,t=殼厚,E=材料楊氏模量,μ=材料泊松比,各參數值可見表1。

當潛器下潛到某深度時,設海水壓力為P2,此時潛器上的各個球殼同時受外壓P2和內壓P1,球殼的排水體積為:

其中Δr1為球殼徑向變形量:

則各個球殼的排水體積減小量為:

2.2.1.2 圓柱殼耐壓結構

表2 “蛟龍”號主要圓柱殼耐壓結構Tab.2 Large cylindrical pressure hull of “JiaoLong”

圓柱殼在下潛開始時內壓已經存在,此時球殼的排水體積為:

其中Δr0和Δl0分別是圓柱殼的徑向和軸向變形量:

(15)式和(16)式中:r2=圓柱殼外半徑;r1=球殼內半徑;l為圓柱殼軸線長度;其他參數的定義同(10)式。各參數的取值見表2。

當潛器下潛到某深度時,設海水壓力為P2,此時潛器上的各個圓柱殼同時受外壓P2和內壓P1,圓柱殼的排水體積為:

其中Δr1和Δl1分別是圓柱殼的徑向和軸向變形量:

則各個圓柱殼的排水體積減小量為:

則整個潛水器的耐壓結構排水體積減小量為:

需要注意的是,圓柱筒形耐壓結構的封蓋可能是平封頭,也可能是半球封頭,且封頭上一般有穿艙件和密封結構,所以封蓋的變形量不能用理論公式方便地求出,這部分的變形量將并入下一類。

2.2.2 其他排水體積變化

除了主要的耐壓結構,潛器上存在大量的密封部件和壓力補償部件,例如載人艙觀察窗有機玻璃和觀察窗基座之間存在摩擦,需要通過有限元接觸分析才能計算觀察窗有機玻璃的變形量。又例如潛器的耐壓結構上的封頭和耐壓結構主體的裝備處往往采用O型圈等密封措施。在這些密封部件處,裝配體之間往往存在間隙,當海水壓力小于某個值時,裝配部位通過密封圈接觸而裝備體互相之間不接觸,由于密封圈材料比較軟(通常密封圈材料為橡膠),此時變形量隨海水壓力增加而有較大的增加;當海水壓力大于某個值時,密封圈變形量已經足夠大,使得裝配體之間開始互相直接接觸,此后變形量隨海水壓力增加的速度明顯減小。這樣的現象在壓力筒試驗中已經觀察到了數次。例如模型球壓力筒試驗中密封端蓋的變形量在約25MPa處存在明顯的拐點。此外潛器上存在的許多壓力補償的部件和設備,例如充油電池箱、充油管路和電纜線等,這些設備和部件的內部承壓部件往往外形復雜且包含多種材料。例如充油電池箱中承壓的是一塊塊的電池,而這些電池的外形比較復雜,電池的等效楊氏模量也未知,所以采用結構分析來計算變形量變得不切實際。本文將所有不能通過球殼和圓柱殼理論公式計算變形量的部件的體積變化歸為“其他排水體積變化”,且假設這一變化量為潛器所受的海水壓力的函數,即:

則潛器總的體積變化量為:

2.3 重力項W

根據水池均衡試驗和海試的經驗,潛器的重量可分為潛器本體重量和可變重量兩部分組成。

根據水池均衡試驗,潛器在水池中均衡時重量為20 807kg[7],此時除了潛器本體重量外,還包含有潛航員重量240kg、艙內耗材(包括二氧化碳吸收劑)13kg和可調壓載水艙中140kg的淡水,則潛器本體重量為:

需要注意的是正常情況下潛器上安裝的鈦制采樣籃是不會變動的,采樣籃將算作潛器本體的一部分,故M0中包括了20kg的鈦合金制采樣籃,而D44潛次中改為重量約25kg的鋼制采樣籃,所以將D44潛次中鋼制采樣籃比原裝鈦制采樣籃多出的5kg歸入作業工具。因此,在計算作業工具的排水體積時,需要扣除更換采樣籃帶來的排水體積變化。潛器的重量可以劃分為如表3所示的幾個部分。將除了M0外的潛器可變重量之和記為MC,則潛器下潛時的重量M為:

表3 重量組成Tab.3 Weight components

2.4 潛器無動力潛浮運動平衡方程

潛器無動力潛浮運動遵守牛頓第二定律,運動方程如(1)式所示,在文獻[7]中,潛器的下潛和上浮運動劃分為以下幾個階段:

(1)壓載水箱注水階段;

(2)平穩下潛階段;

(3)拋掉P1后減速階段;

(4)潛器靠底階段;

(5)拋掉P2后加速階段;

(6)平穩上浮階段;

(7)壓載水箱排水階段。

在文獻[9]中,對本文模型的7個階段均進行了詳細的分析,出于篇幅考慮,此處只給出平穩下潛階段和平穩上浮階段的分析結果。

2.4.1 平穩下潛階段

壓載水箱注水完成后潛器排水體積不再有大的變化 (除了由于海水壓力導致的排水體積減小),此階段潛器的運動也是動態平衡的,即隨著重力與浮力阻力之差程減速—加速—減速—加速……的過程,但是通過分析3 000米海試和5 000米海試的數據可見每個減速—加速循環的周期不長,且潛器的速度隨著下潛深度的增加會逐漸降低,D44潛次從潛航員關閉壓載水箱上閥門到拋載P1的下潛過程的速度—深度曲線如圖2所示。

圖2 D44平穩下潛階段速度—深度曲線Fig.2 The velocity curve of the steady descent stage of dive 44

為了簡化分析,本文將這一階段潛器的運動近似為平穩運動,即認為加速度項為0,加速度為零則不考慮附加質量,此階段潛器的各項力為:

2.4.2 潛器平穩上浮階段

潛器的重力、浮力和阻力達到動態平衡后,潛器的上浮運動可看成平穩運動,此階段前期的運動稱為平穩上浮階段,此階段潛器的各項受力為:

則潛器這一階段的運動可以近似描述為平衡方程:

則潛器這一階段的運動可以近似描述為平衡方程:

3 D44潛次數據分析

D44潛次實現了完全的無動力下潛上浮運動,其速度曲線相對沒有實現無動力上浮下潛運動的潛次來說規律性比較明顯,為了確定已經建立的無動力運動模型中的未知參數,本文對D44潛次的配載情況、航行記錄和潛浮運動數據進行了詳細的分析。

3.1 D44配載情況

D44潛次實現了無動力上浮下潛運動,其配載較合理,根據水池均衡試驗結果和D44潛次的記錄,D44潛次潛器的配載情況如下表4所示。

表4 D44配載情況Tab.4 Weight of dive 44

3.2 D44無動力潛浮運動數據

與無動力潛浮運動相關的D44潛次的海水密度—深度曲線、深度—時間曲線、速度—深度曲線如下圖3-5所示。

圖3 海水密度—深度曲線Fig.3 Density curve

圖4 D44下潛和上浮的深度—時間曲線Fig.4 The depth curve of submerge and ascent stage of dive 44

需要說明以下幾點:(a)從圖3中可看到深度大于1 000米后,海水密度隨深度的變化呈線性,所以本文經將D44潛次的海水密度曲線外推至7 500米,且將該曲線和3 000米海試的海水密度曲線比較后發現,1000米以深的海水密度曲線差別較小。(b)圖4和圖5中的下潛曲線包括了從壓載水艙開始注水到拋掉P1的整個下潛過程;而上浮曲線包括了拋掉P2到壓載水艙開始排水的整個上浮過程。(c)與上浮至3 600米左右控制潛器上浮的操作記錄相對應,可看到圖5中平穩上浮階段的速度曲線在3 600米附近有一個“尖刺”。

而本文主要分析潛器下潛和上浮的平穩運動階段,故將加速度值較大的曲線段去除,僅對平穩階段的曲線進行分析,由于潛器是動態平衡的,所以速度曲線存在較多小幅的加速—減速振蕩,本文對D44平穩下潛階段和平穩上浮階段的速度曲線進行了光順處理,以抓住平穩運動的規律,D44平穩下潛階段的速度曲線和光順后的曲線見圖6所示,同理D44平穩上浮階段的速度曲線和光順后的曲線見圖7所示。

圖5 D44下潛和上浮的速度—深度曲線Fig.5 Velocity curve of submerge and ascent stage of dive dive 44

圖6 D44平穩下潛階段速度—深度曲線及光順后曲線Fig.6 The smoothed velocity curve of the steady descent stage of dive 44

圖7 D44平穩上浮階段速度—深度曲線及光順后曲線Fig.7 The smoothed velocity curve of the steady asscent stage of dive 44

4 模型參數辨識

4.1 平穩下潛階段

由平穩下潛階段的平衡方程式(26)可知,無動力運動模型中的“其他排水體積變化量—ΔVother”是未知的,需要由5 000米海試D44潛次的平穩下潛階段的曲線求得。D44下潛的ΔVother對海水壓力P的曲線見圖8所示,由光順后曲線可看到在約20-26MPa之間潛器的排水體積變化量存在明顯的拐點,可認為潛器上大部分設備的密封部位和補償部位在20-26MPa時達到緊實,此后ΔVother隨P增長的速度降低(表現為圖上的曲線斜率降低)。

圖8 D44平穩下潛階段ΔVother隨海水壓力P變化曲線Fig.8 The ΔVothercurve of the steady descent stage of dive 44

4.2 平穩上浮階段

平穩上浮階段的平衡方程見(27)式,ΔVother同樣是未知的,需要由5 000米海試D44潛次的平穩上浮階段的曲線求得,需要注意的是由于壓縮量恢復有一定的滯后,上浮階段的ΔVother和下潛階段的ΔVother可能不同。ΔVother隨P變化的曲線見圖9,可見當海水壓力大于6MPa時潛器ΔVother隨P的變化可近似為直線,而P在4.5MPa至6MPa之間曲線存在明顯的拐點,當P小于4.5MPa時ΔVother減小的速度明顯加快,可認為在4.5~6MPa之間潛器密封部件的密封面脫離閉合狀態,當然,僅由D44潛次一個潛次的數據還不能下充分的結論。

圖9 D44平穩上浮階段ΔVother隨海水壓力P變化曲線Fig.9 The ΔVothercurve of the steady ascent stage of dive 44

4.3 參數取值

將圖8和圖9中的D44下潛和上浮求出的ΔVother隨海水壓力P變化的光順后的曲線取出放在同一張圖中比較,如圖10所示。可見:隨著下潛深度增大,下潛和上浮的ΔVother值逐漸接近,這個現象和實際情況較符合,因為D44潛次最大潛深為5180米,在較深的深度處,例如5 000米處,潛器下潛和上浮至該深度處的壓縮量應該差別不大,因為潛器下潛超過5 000米深度只有180米,在最深的5 180米處作業接近2小時,這2個小時內潛器相當于處在保壓狀態,ΔVother將隨著作業時間增加而略有增加,但是從5 180米上浮至5 000米時的壓縮量和下潛至5 000米時的壓縮量不應有大的差別。

圖10 D44航次下潛和上浮ΔVother隨P變化曲線(光順后)比較Fig.10 Comparison of smoothed ΔVothercurve of the steady descent stage and ascent stage of dive 44

為了將5 000米海試的結果外推至7 000米,需要對圖10中的曲線進行分析和外插。由圖10,下潛的ΔVother-P曲線可近似為三段:

(1)P<13MPa時,此階段潛器的ΔVother隨P的變化可近似為直線;

(2)13MPa<P<33MPa時,此階段潛器的ΔVother隨P的變化出現拐點,直接取由D44下潛數據反求出的光順后的曲線段;

(3)P>33MPa時,ΔVother隨P的變化也可近似為直線,將該直線段外插至76.7MPa(約7 500米),以便于估算7 000米處的ΔVother。

圖11 下潛ΔVother-P曲線的處理和外插Fig.11 Extrapolation of ΔVothercurve of the steady descent stage of dive 44

近似結果如上圖11所示。

類似的,上浮曲線也可分為三段:

(1)P<3MPa時,此階段潛器的ΔVother隨P的變化可近似為直線;

(2)3MPa<P<7.6MPa時,此階段潛器的ΔVother隨P的變化出現拐點,直接取由D44上浮數據反求出的光順后的曲線段;

(3)P>7.6MPa時,ΔVother隨P的變化也可近似為直線,將該直線段外插至76.7MPa(約7 500米)。

上浮近似結果如下圖12所示。

將處理并外插后的下潛和上浮ΔVother-P曲線置于同一張圖中,如圖13所示,可看到下潛和上浮曲線的交點在約5 736米處,與最大下潛深度5 180米差556米,一方面可能是由于潛器在高壓下作業的過程中ΔVother會存在少量增加,另一方面也可能是由于上浮和下潛的數據存在誤差等導致。

圖12 上浮ΔVother-P曲線的處理和外插Fig.12 Extrapolation of ΔVothercurve of the steady ascent stage of dive 44

圖13 上浮和下潛的ΔVother-P曲線Fig.13 Used ΔVothercurve

可見,在前述模型的基礎上,為了更準確地描述潛器上浮階段的運動,需要增加新的變量:作業階段(保壓階段)體積壓縮量ΔVwork—反應潛器在底部作業過程中,即保壓狀態下排水體積減小量。對指定的潛器,ΔVwork是作業壓力和作業時間的函數,由于“蛟龍”號的作業時間基本固定,且還沒有足夠的數據支持,本文在此暫時假設ΔVwork只是作業壓力的線性函數:

則平穩上浮階段平衡方程改為:

4.4 程序修改

根據以上計算模型和參數取值,對原有程序進行了修改,界面也進行了改進,如圖14所示。計算輸入分為兩個部分,第一部分是預計作業深度Z0、上浮拋載、下潛拋載及密度曲線選擇;第二部分是潛器重量和排水情況。計算結果將輸出Z0處的均衡情況、下潛速度曲線和上浮速度曲線,同時也給出下潛/上浮的速度均值以及時間估算。

程序的大部分操作和原程序類似,但是存在以下不同:

(1)程序要求較準確地輸入潛器各部分重量以及排水體積情況。注意“艙內耗材”、“艙內人員”和“壓載水箱注水”三項的排水體積為0;浮拋載和下潛拋載需要根據計算均衡結果人工進行調整;

圖14 程序界面圖14 GUI of the program

(2)程序的“參考密度曲線”為5 000米海試得到的密度—深度曲線拓展至10 000米的密度曲線;“實測密度曲線”為下潛海區的實測密度—深度曲線,數據格式見程序目錄下的 “深度密度曲線.csv”文件,且要求“深度密度曲線.csv”中的第一個深度是0,最后一個深度大于程序輸入的Z0(大深度下密度-深度曲線可近似為直線,所以操作者可將實測密度-深度曲線的后半段 (例如對于Z0=7 000米,而實測曲線只有6 000米的數據,則操作者可以對3 000米以深的密度—深度曲線進行線性擬合,操作者可以根據具體曲線形狀在excel中方便地進行線性擬合和外插操作,本程序不再提供該功能)線性外插拓展至需要的深度 (一般可外插至10 000米),并將拓展后的曲線保存至 “深度密度曲線.csv”中)。

(3)程序的輸出調整為下潛深度Z0處的均衡情況。如果要求Z0處潛器達到均衡,則應調整下潛拋載P2(因為潛器下潛至Z0處時P1已經被拋掉,所以P1不影響均衡位置)使結果接近于0;如果要求Z0處潛器具有足夠的坐底力,則應調整下潛拋載P2使結果接近于設計坐底力(“蛟龍 ” 號 設 計 坐 底 力 ≈40)。

4.5 D44潛次計算結果

采用修改后的程序,將D44潛次的載荷數據填入后采用程序進行計算,計算結果如圖15所示。

可見5 180處潛器重力浮力差(即坐底力)為35kg,與40kg比較接近,且速度曲線和圖6及圖7類似。

5 結 語

本文根據“蛟龍”號5 000米海試的D44潛次無動力潛浮海試數據對“蛟龍”號無動力潛浮運動模型進行了再分析,改進了3 000米海試后提出的分析模型,主要是采用了可變的體積壓縮率的概念,且對程序進行了修改,開發出了更為可靠的蛟龍號載人潛器無動力潛浮運動分析系統,使得計算結果更接近實際情況,程序界面更加符合實際應用,為2012年的7 000米級海試做好了準備,也為4 500米載人潛器無動力潛浮運動的設計打下了基礎。

[1]Busby,Frank.Manned Submersibles[M].Arlington,VA:Office of the Oceanographer of the Navy,1976.

[2]Allmendinger,Eugene E.Submersible Vehicle Systems Design[M].ill.Jersey City,NJ:Society of Naval Architects and Marine Engineers,1990.

[3]朱繼懋主編.潛水器設計[M].上海:上海交通大學出版社,1992.

[4]崔維成,馬 嶺.潛水器設計中所要解決的水動力學問題[C].第九屆全國水動力學學術會議暨第二十二屆全國水動力學研討會文集,吳有生,劉 樺,許唯臨,周連第,楊顯成主編,海洋出版社,北京,2009:9-29.

[5]劉正元.7000米載人潛器無動力潛浮運動研究[R].無錫:中國船舶科學研究中心科技報告,2004.

[6]劉正元.潛水器大攻角范圍內運動的仿真[J].船舶力學,2005,9(2):54-59.Liu Zhengyuan.Simulation of submersible motion in large attack angle[J].Journal of Ship Mechanics,2005,9(2):54-59.

[7]劉正元,潘彬彬,崔維成.7000m載人潛器海試試驗數據整理和分析[R].無錫:中國船舶科學研究中心科技報告,2010.

[8]范欣珊.壓力容器的應力分析與強度設計[M].北京:原子能出版社,1979.

[9]潘彬彬,崔維成,葉 聰,劉正元.7000m載人潛器無動力潛浮運動分析[R].無錫:中國船舶科學研究中心科技報告,2011.

猜你喜歡
深度
深度理解不等關系
四增四減 深度推進
深度理解一元一次方程
深度觀察
深度觀察
深度觀察
深度觀察
芻議深度報道的深度與“文”度
新聞傳播(2016年10期)2016-09-26 12:14:59
提升深度報道量與質
新聞傳播(2015年10期)2015-07-18 11:05:40
微小提議 深度思考
主站蜘蛛池模板: 一级毛片免费的| 91在线播放国产| 99精品视频播放| 特级毛片免费视频| 国产在线观看高清不卡| 久久精品无码中文字幕| 久久永久精品免费视频| 久久婷婷五月综合色一区二区| 久久久久国色AV免费观看性色| 亚洲性一区| 国产精品性| 麻豆精品在线播放| 中文字幕天无码久久精品视频免费| 欧美性精品| 一级毛片在线免费视频| 福利小视频在线播放| 亚洲欧美激情另类| 婷婷99视频精品全部在线观看| 国产无码性爱一区二区三区| 一级毛片免费观看不卡视频| 日韩精品高清自在线| 波多野结衣一区二区三区88| 米奇精品一区二区三区| 欧美亚洲日韩中文| 亚洲另类国产欧美一区二区| 99re在线免费视频| 久久精品欧美一区二区| 999精品视频在线| 亚洲无码高清免费视频亚洲| 国产精品福利社| AV天堂资源福利在线观看| 免费观看国产小粉嫩喷水 | 国产鲁鲁视频在线观看| 亚洲精品午夜无码电影网| 国产无码精品在线| 欧美狠狠干| 国产正在播放| 97影院午夜在线观看视频| 国产亚洲视频免费播放| 在线视频亚洲色图| 99久久国产自偷自偷免费一区| 国产成人午夜福利免费无码r| 青青草综合网| 伊人激情综合| 午夜精品国产自在| 老汉色老汉首页a亚洲| 久久精品日日躁夜夜躁欧美| 亚洲国产成人麻豆精品| 欧洲亚洲欧美国产日本高清| 99精品伊人久久久大香线蕉| 热伊人99re久久精品最新地| 亚洲综合婷婷激情| 亚洲AV无码一二区三区在线播放| 波多野结衣在线一区二区| 久久亚洲国产一区二区| 尤物特级无码毛片免费| 欧美黄色网站在线看| 免费毛片视频| 九九视频免费在线观看| 日韩精品无码一级毛片免费| 91成人在线免费视频| 久久情精品国产品免费| 天堂久久久久久中文字幕| 中国特黄美女一级视频| 91一级片| 国产成人av一区二区三区| 国产91无码福利在线| 91伊人国产| 国产高清在线精品一区二区三区| 亚洲国产一区在线观看| 狠狠色噜噜狠狠狠狠奇米777| 国产精品自拍合集| 亚洲精品第一页不卡| 亚洲第一成年网| 国产AV毛片| a天堂视频在线| AV色爱天堂网| 国产欧美精品一区aⅴ影院| 精品1区2区3区| 国产午夜在线观看视频| 亚洲精品第五页| 国产91精品久久|