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

海洋資源四維勘探拖纜陣列動(dòng)力學(xué)特性仿真研究

2012-01-08 04:59:46吳喆瑩張維競(jìng)張廣磊史斌杰
海洋工程 2012年4期

吳喆瑩,張維競(jìng),劉 濤,張廣磊,史斌杰

(上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海 200240)

海洋資源四維勘探拖纜陣列動(dòng)力學(xué)特性仿真研究

吳喆瑩,張維競(jìng),劉 濤,張廣磊,史斌杰

(上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海 200240)

面對(duì)日益增長(zhǎng)的海洋資源勘探需求,海洋資源四維勘探拖纜陣列的研究越來(lái)越被人們所重視。作為拖纜動(dòng)態(tài)控制策略研究的前提與依據(jù),拖纜陣列的動(dòng)力學(xué)特性始終是該領(lǐng)域在技術(shù)突破進(jìn)程中的重要課題。以Ablow的經(jīng)典模型為基礎(chǔ),采用微元法對(duì)拖纜陣列建立三維模型,并在時(shí)域上采用廣義α算法對(duì)非線性方程進(jìn)行離散求解,通過(guò)程序編譯,對(duì)拖纜在四種工況(即拖船勻速直線航行、回轉(zhuǎn)航行、變速直線航行以及拖點(diǎn)處存在升沉運(yùn)動(dòng)等)下的響應(yīng)進(jìn)行計(jì)算分析,著重考量其深度及張力的變化。仿真計(jì)算結(jié)果表明,廣義α算法對(duì)拖纜陣列姿態(tài)在空間和時(shí)間上離散求解,算法穩(wěn)定性良好,彌補(bǔ)了傳統(tǒng)有限差分法對(duì)于時(shí)域不穩(wěn)定的缺陷。

海洋勘探;拖纜陣列;廣義α算法;動(dòng)力學(xué)特性;仿真

海洋中蘊(yùn)藏著豐富、巨量的油氣資源[1],海洋資源的勘探與開發(fā)已成為能源發(fā)展的重點(diǎn)和熱點(diǎn)。海洋資源四維勘探是目前海底探查成效最高的地球物理技術(shù),它建立在海洋聲學(xué)的基礎(chǔ)上,依靠水下拖曳系統(tǒng)來(lái)實(shí)現(xiàn)勘探作業(yè)[2]。水下拖曳系統(tǒng)通常由水面拖曳母船、拖纜、尾繩及各拖曳設(shè)備等組成一個(gè)復(fù)雜系統(tǒng),掌握其在各種情況下動(dòng)力學(xué)特性,特別是水下拖纜的動(dòng)力學(xué)特性,對(duì)系統(tǒng)工作效率提升、系統(tǒng)穩(wěn)定性保障等有著重要作用。

目前國(guó)內(nèi)外拖纜系統(tǒng)的動(dòng)力學(xué)模型較為成熟,主流的有 Ablow 和 Schechter[3],Sanders[4],Delmer[5],Mili-nazzo[6]等提出的三維模型。其中,Ablow和Schechter以及Milinazzo采用有限差分法來(lái)求解模型;Delmer采用有限元法求解問(wèn)題。近年來(lái)也有不少研究者用集中質(zhì)量法、有限差分法、直接積分法等求解拖纜系統(tǒng)的動(dòng)、穩(wěn)態(tài)運(yùn)動(dòng)[7-12]。基于前人的研究與探索,采用Ablow和Schechter的經(jīng)典模型,采用廣義α算法求解三維拖纜方程,并對(duì)典型工況(包括勻速直線航行、回轉(zhuǎn)航行、變速直線航行、船體升沉運(yùn)動(dòng)等工況)下的拖纜動(dòng)力學(xué)特性進(jìn)行仿真研究,為后期控制策略的深入研究做好鋪墊。

1 拖纜動(dòng)力學(xué)模型

依據(jù)Ablow與Schechter提出并經(jīng)實(shí)驗(yàn)有效驗(yàn)證的模型,假設(shè)拖纜為圓柱型柔性纜,忽略拖纜彎曲應(yīng)力[13]和海流的影響,有如下方程:

圖1 歐拉角(θ,ψ,φ)示意(ψ =90°)Fig.1 Illustration of(θ,ψ,φ)(ψ =90°)

定義 Y(s,t)=(T,Vt,Vn,Vb,φ,θ)T,則上述方程(1)~ (6)有矩陣形式:

其中,

2 邊界條件

2.1 拖纜首端邊界條件

首端纜上拖點(diǎn)的速度與拖船速度相同,即

式中:Vx,Vy,Vz為拖船航速在固定坐標(biāo)系 x,y,z方向上的分量。

2.2 尾端邊界條件

文獻(xiàn)[12]、[14]等中對(duì)于拖纜系統(tǒng)尾端邊界條件設(shè)為自由端并定義張力T=0。實(shí)際工況中,拖纜系統(tǒng)尾端有尾浮,其一般在拖纜末端產(chǎn)生一個(gè)較大的預(yù)張力以增加系統(tǒng)的穩(wěn)定性及可控性。由于尾浮與拖纜間的耦合作用,拖纜的運(yùn)動(dòng)將變得非常復(fù)雜,而尾浮在側(cè)向的偏移很小,其側(cè)向與縱向的耦合可以忽略,僅在速度方向上取近似阻力[15]。

參考Egil Pedersen[16]的實(shí)船試驗(yàn),在穩(wěn)態(tài)情況下簡(jiǎn)化假設(shè)尾浮對(duì)拖纜系統(tǒng)提供一個(gè)定常的張力T=T0。

3 拖纜方程的離散

3.1 空間離散

拖纜方程在空間上的離散參考Ablow和Schechter采用的基于二階精度的有限差分方法,假設(shè)將拖纜劃分為n個(gè)節(jié)點(diǎn),記j為第j個(gè)空間點(diǎn)(j=1,2,…,n),則拖纜的動(dòng)態(tài)方程可離散為n-1個(gè)矩陣方程:

3.2 時(shí)間離散

然而,有限差分方法在時(shí)間離散求解上并不穩(wěn)定[17]。對(duì)此,Gobat曾將廣義α算法用于二維拖纜方程的求解中,并已獲得良好的精度和穩(wěn)定性[17-18]。經(jīng)驗(yàn)算,將廣義α算法用到了三維拖纜方程的求解中,對(duì)方程在時(shí)間上進(jìn)行離散。

廣義α算法在時(shí)間上的離散形式:

式中:Δt表示時(shí)間步長(zhǎng),αm、αk和γ則構(gòu)成了廣義α算法,i為第i個(gè)時(shí)間點(diǎn)。

由于廣義α算法在求解線性問(wèn)題時(shí),系數(shù)矩陣M、K恒定不變,因此可以用式(12)進(jìn)行計(jì)算。但是對(duì)于非線性問(wèn)題系數(shù)矩陣會(huì)隨時(shí)間的變化而改變,算法的穩(wěn)定性將隨時(shí)間有條件收斂。為避免該問(wèn)題以提高算法的穩(wěn)定性,采用系數(shù)矩陣與速度矢量、位移矢量具有相同權(quán)重的方法。此時(shí),算法的表達(dá)形式變?yōu)?/p>

4 典型工況仿真計(jì)算結(jié)果

在前期的研究中,通過(guò)與Ablow在1983年所做實(shí)驗(yàn)的數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)算,并結(jié)合近期在拖曳水池所做的模型試驗(yàn)的趨勢(shì)分析結(jié)果,表明利用廣義α算法求解拖纜動(dòng)力學(xué)模型是可靠的[19]。在此基礎(chǔ)上,進(jìn)行進(jìn)一步動(dòng)力學(xué)特性仿真分析。

4.1 拖纜系統(tǒng)穩(wěn)態(tài)仿真

4.1.1 考慮尾拖作用的拖纜系統(tǒng)穩(wěn)態(tài)仿真

1)仿真參數(shù)

在當(dāng)前的實(shí)際勘探作業(yè)中,陣列纜的長(zhǎng)度一般為3 000~8 000 m長(zhǎng),勘探船一般以4~5節(jié)的船速航行。參考文獻(xiàn)[11]、[15]、[16]等中的拖纜參數(shù),這里以250 m長(zhǎng)的引導(dǎo)拖纜、3 000 m長(zhǎng)陣列纜和50 m長(zhǎng)的尾纜作為研究對(duì)象。

以10 m為一段,將引導(dǎo)拖纜分為25段,陣列纜分為300段,尾纜分為5段,共330個(gè)分段;以尾端為初始點(diǎn)(第1點(diǎn))標(biāo)記,直至拖纜與拖船結(jié)合點(diǎn)(即拖點(diǎn))為末位點(diǎn)第331點(diǎn),總共331個(gè)節(jié)點(diǎn)。其中,陣列纜在水中的凈浮力為零(拖纜參數(shù)見表1)。

表1 拖纜系統(tǒng)的物理特性Tab.1 Physical properties of the towed cable

拖曳船以一定速度拖帶拖纜系統(tǒng)沿直線穩(wěn)定航行,拖曳系統(tǒng)尾部考慮拖船作用時(shí),將其近似為尾部張力T=T0。參考并引用Egil Pedersen[16]實(shí)驗(yàn)結(jié)果,在拖船速度為 2.0 m/s 和2.6 m/s時(shí)分別取尾端水阻力(Fd)為406 N和743 N,則拖纜尾端邊界條件如表2所示。

2)仿真計(jì)算結(jié)果

分別計(jì)算該條件下拖船以2.0 m/s和2.6 m/s速度沿直線穩(wěn)定航行工況下,仿真結(jié)果如圖2和圖3。

表2 拖纜尾端邊界條件Tab.2 Boundary conditions at the end of towed cable

圖2 考慮尾拖船時(shí)拖纜在勻速直線拖曳工況下深度響應(yīng)Fig.2 Position of towed cable with tail boat in steady straight course

圖3 考慮尾拖船時(shí)拖纜勻速直線拖曳工況下張力響應(yīng)Fig.3 Tension of towed cable with tail boat in steady straight course

結(jié)果顯示:2.0 m/s拖曳速度下,由于引導(dǎo)纜與尾繩凈浮力不為零且拖船速度較慢,陣列纜深度在-8.8~-13 m之間,當(dāng)拖船速度增大時(shí),拖纜陣列有所上浮且整條纜姿態(tài)更趨于平穩(wěn)(見圖2),同時(shí)全纜各處的張力均勻變化,在即定速度拖帶下,拖纜拖點(diǎn)處所承受的張力最大,尾端張力最小(見圖3)。

4.1.2 不考慮尾拖作用的拖纜系統(tǒng)穩(wěn)態(tài)仿真

1)仿真參數(shù)

[12]、[14]、[15]中,當(dāng)不考慮尾浮時(shí),尾端可近似看做自由端,且有T=0。由于此時(shí)拖纜尾端為自由端,為了使尾繩在水中的重力影響拖纜原有的姿態(tài),將尾繩調(diào)整水中凈浮力為零的纜,參數(shù)詳見表3。

表3 拖纜系統(tǒng)的物理特性Tab.3 Physical properties of the towed cable

2)仿真結(jié)果

結(jié)果顯示:2.0 m/s拖曳速度下,陣列纜深度維持在-8.8 m左右;2.6 m/s拖曳速度下,陣列纜深度維持在-5.4 m左右(圖4)。由于陣列纜為零浮力纜,其在水中姿態(tài)平穩(wěn)。全纜各處的張力變化趨勢(shì)不變(對(duì)比圖5與圖3)。

圖4 不考慮尾拖船時(shí)勻速直線拖曳工況下拖纜深度圖Fig.4 Position of towed cable without tail boat in steady straight course

圖5 不考慮尾拖船時(shí)拖纜勻速直線拖曳工況下張力響應(yīng)Fig.5 Tension of towed cable without tail boat in steady straight course

4.2 拖纜系統(tǒng)動(dòng)態(tài)仿真

拖纜分段方式及節(jié)點(diǎn)選取方式同穩(wěn)態(tài)情況下的仿真計(jì)算一樣,拖纜參數(shù)詳見表3,各工況下的仿真不計(jì)尾拖船的影響。

4.2.1 拖船回轉(zhuǎn)航行

船舶以一定速度拖帶系統(tǒng)回轉(zhuǎn)。計(jì)算工況:初始時(shí)刻,拖船系統(tǒng)以vx=2 m/s的速度勻速穩(wěn)定直航;此后,進(jìn)入半徑為1 000 m的回轉(zhuǎn)運(yùn)動(dòng)并完成180°回轉(zhuǎn);最后,仍以vx=2 m/s的速度開出,并持續(xù)航行至穩(wěn)定狀態(tài)。仿真6 000 s后,結(jié)果如圖6和圖7所示。

圖6 回轉(zhuǎn)航行纜上各點(diǎn)深度隨時(shí)間變化曲線Fig.6 Response of depth on different points of towed cable under circular manoeuver

結(jié)果顯示:拖纜上的點(diǎn)在進(jìn)入回轉(zhuǎn)狀態(tài)時(shí)先有所下沉,在回轉(zhuǎn)運(yùn)動(dòng)的尾聲接近最低點(diǎn),之后再緩慢上浮并趨近初始平衡位置(如圖6所示);對(duì)比圖6(a)與圖6(b),纜上后端位置的點(diǎn)在回轉(zhuǎn)過(guò)程中深度的波動(dòng)幅度小于纜上前端位置的點(diǎn);拖點(diǎn)處張力在此過(guò)程中,先緩慢變小并于回轉(zhuǎn)運(yùn)動(dòng)的尾聲到達(dá)最低點(diǎn),之后張力緩慢增加并趨近初始平衡位置水平(如圖7)。這與文獻(xiàn)[15]所得的結(jié)果趨勢(shì)相符。

4.2.2 拖船沿直線變速航行

船舶開始以vx=2.6 m/s的的速度沿直線行駛,假定某一時(shí)刻沿航行方向(即vx方向)船速做幅值為0.1 m/s,周期為1 000 s的正弦變化來(lái)模擬實(shí)船直線變速航行工況。仿真3 000 s后得結(jié)果如圖8和圖9所示。

圖7 回轉(zhuǎn)航行拖點(diǎn)處張力隨時(shí)間變化曲線Fig.7 Response of tension at towed point of the towed cable under circular manoeuver

圖8 直線變速航行拖纜各位置深度沿時(shí)間變化Fig.8 Response of depth at certain points of towed cable during speed change of straight run

圖9 直線變速航行拖點(diǎn)張力沿時(shí)間變化Fig.9 Response of tension at towed point of towed cable during speed change of straight run

結(jié)果顯示:1)在船舶變速直線行駛過(guò)程中,纜上各點(diǎn)的深度波動(dòng)變化不大,較之纜長(zhǎng)可以忽略不計(jì);2)對(duì)于低頻擾動(dòng),拖纜深度變化的影響難以傳到尾端,故尾端在深度上基本沒有變化;3)拖纜上各點(diǎn)張力變化卻非常大,在±0.1 m/s速度變化的情況下拖點(diǎn)處張力變化范圍超過(guò)2 kN。

4.2.3 拖點(diǎn)處有升沉運(yùn)動(dòng)

船舶以vx=2.6 m/s的速度沿直線航行,在拖點(diǎn)處沿深度方向(即vz方向)增加正弦變化速度以模擬拖點(diǎn)處升沉運(yùn)動(dòng),其幅值為1.59 m,周期為10 s。仿真結(jié)果如圖10和圖11。

圖10 升沉運(yùn)動(dòng)時(shí)陣列纜首端深度變化Fig.10 Response of depth at certain points of towed cablein heaving motion

圖11 升沉運(yùn)動(dòng)時(shí)拖點(diǎn)張力變化Fig.11 Response of tension at certain pointss of towed cable in heaving motion

結(jié)果顯示:拖纜在深沉運(yùn)動(dòng)激勵(lì)下深度變化較為明顯;然而張力變化相對(duì)直線變速工況而言變化較小。

通過(guò)上述計(jì)算,可以看出在時(shí)間步長(zhǎng)從0.625 s到10 s,仿真時(shí)間從100 s到6 000 s過(guò)程,工況從直線穩(wěn)定到回轉(zhuǎn)運(yùn)動(dòng)、變速直線運(yùn)動(dòng)以及拖點(diǎn)升沉運(yùn)動(dòng)的狀況下都未出現(xiàn)不穩(wěn)定情況。再次表明,這里所用廣義α算法穩(wěn)定性良好,進(jìn)一步驗(yàn)證了此算法在求解拖纜動(dòng)力學(xué)模型是正確可行的。仿真結(jié)果可為下一步制定拖纜在上述工況下的動(dòng)態(tài)控制策略給出了良好的分析依據(jù)。

5 結(jié)語(yǔ)

采用微元法對(duì)海洋資源四維勘探拖纜陣列建立三維動(dòng)力學(xué)模型進(jìn)行動(dòng)力學(xué)特性研究,并用廣義α算法對(duì)拖纜陣列姿態(tài)在空間和時(shí)間上離散求解,彌補(bǔ)了傳統(tǒng)有限差分法對(duì)于時(shí)域不穩(wěn)定的缺陷。通過(guò)計(jì)算機(jī)仿真程序編譯,算出拖纜在勻速直線航行、回轉(zhuǎn)航行、變速直線航行以及拖點(diǎn)升沉運(yùn)動(dòng)四種工況下拖纜的動(dòng)力學(xué)特性,為后續(xù)進(jìn)一步探究拖纜動(dòng)態(tài)控制策略及纜上控制器設(shè)計(jì)等提供了依據(jù)。

參考文獻(xiàn):

[1]趙政璋,趙賢正,李景明,等.國(guó)外海洋深水油氣勘探發(fā)展趨勢(shì)及啟示[J].中國(guó)石油勘探,2005(6):71-76.

[2]陳小宏,牟永光.四維地震油藏監(jiān)測(cè)技術(shù)及其應(yīng)用[J].石油地球物理勘探,1998(6):707-715.

[3]Ablow C M,Schechter.Numerical simulation of undersea cable dynamics[J].Ocean Engineering,1983(10):443-457.

[4]Sanders J V.A three-dimensional dynamic analysis of a towed system[J].Ocean Engineering,1982,9(5):483-499.

[5]Delmer T N,Stephens T C,Coe J M.Numerical simulation of towed cables[J].Ocean Engineering,1983,10(2):119-132.

[6]Milinazzo F,Wilkie M,Latchman S A.An efficient algorithm for simulating the dynamics of towed cable systems[J].Ocean Engineering,1987,14(6):513-526.

[7]Wang F,Huang G L,Deng D H.Steady state analysis of towed marine cable[J].Journal of Shanghai Jiao Tong University:Science,2008,13(2):239-244.

[8]Sun Y,Leonard J W,Chiou R B.Simulation of unsteady oceanic cable deployment by direct integration with suppression[J].O-cean Engineering,1994,21(3):243-256.

[9]Huang S.Dynamic analysis of 3-dimensional marine cables[J].Ocean Engineering,1994,21(6):587-605.

[10]張維競(jìng),張小卿,陳 峻.基于嵌入式水鳥的海洋地震拖纜運(yùn)動(dòng)狀態(tài)仿真研究[J].海洋工程,2009,27(4):81-86

[11]王春杰,張維競(jìng),劉 濤,等.橫向海流作用下海洋地震拖纜姿態(tài)控制[J].大連海事大學(xué)學(xué)報(bào),2011,37(1):9-17.

[12]鄧德衡,黃國(guó)樑,樓連根.拖曳線列陣陣形與姿態(tài)數(shù)值計(jì)算[J].海洋工程,1999,17(1):17-26.

[13]王 飛,陳錦標(biāo),涂興華.低速低應(yīng)力時(shí)拖纜運(yùn)動(dòng)仿真[J].上海交通大學(xué)學(xué)報(bào),2010,44(6):828-832.

[14]羅 薇,張 攀.拖纜系統(tǒng)運(yùn)動(dòng)仿真[J].武漢理工大學(xué)學(xué)報(bào),2005,29(5):724-726

[15]王 飛.海洋勘探拖曳系統(tǒng)運(yùn)動(dòng)仿真與控制技術(shù)研究[D].上海:上海交通大學(xué),2006.

[16]Egil Pedersen.A Nautical Study of Towed Marine Seismic Streamer Cable Configurations[D].Doctoral Thesis in Nautical Engineering of Norweigian University of Science and Technology,1996.

[17]Jason I Gobat.The Dynamics of Geometrically Compliant Mooring Systems[D].Massachusetts Institute of Technology,Woods Hole Oceanographic Institution,2000.

[18]Gobat J I,Grosenbaugh M A.Time-domain numerical simulation of ocean cable structures[J].Ocean Engineering,2005(33):1373-1400.

[19]張廣磊,張維競(jìng),劉 濤.廣義α算法在拖曳陣列動(dòng)態(tài)仿真中的應(yīng)用研究[J].哈爾濱工程大學(xué)學(xué)報(bào),2012,33(5):1-6.

Simulation study on dynamic characteristics of towed array for four-dimensional exploration of marine resources

WU Zhe-ying,ZHANG Wei-jing,LIU Tao,ZHANG Guang-lei,SHI Bin-jie
(National Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

In the face of the growing demand of marine resources exploration,studies are increasingly concentrated on towed array for four-dimensional exploration of marine resources.As the premise and basis of further study on dynamic control strategy of towed array,dynamic characteristics of towed array are always the crucial task of technology development in this field.Based on the classic Ablow model,a three-dimensional model is established by applying micro-element method to towed array in this paper and a generalized-α algorithm is introduced to temporally discrete and solve the nonlinear equations.By programming,the computation analysis of response of towed cable,especially the variation in depth and tension,under four different motions,such as steady straight motion,circular manoeuver motion,speed changing motion and heaving motion,is made,which provide an important basis for the following research in making dynamic control strategy of towed array.

marine resources exploration;towed cable;generalized-α;dynamic characteristics;simulation

P741

A

1005-9865(2012)04-0118-07

2011-11-18

國(guó)家自然科學(xué)基金資助項(xiàng)目(51079083)

吳喆瑩(1986-),女,上海人,碩士生,主要從事輪機(jī)工程方向的研究。E-mail:tracywzy@sjtu.edu.cn

張維競(jìng)。E-mail:wjzhang@sjtu.edu.cn

主站蜘蛛池模板: 91精品国产麻豆国产自产在线| 亚洲欧美国产五月天综合| 国产精品无码作爱| 人人澡人人爽欧美一区| 亚洲国产一区在线观看| 成人在线不卡视频| 亚洲一级毛片免费看| 91久久精品国产| 免费 国产 无码久久久| 欧美天堂久久| 国产高清在线精品一区二区三区 | 亚洲人妖在线| 午夜无码一区二区三区在线app| 国产人妖视频一区在线观看| 91在线精品免费免费播放| 四虎永久免费地址在线网站| 亚洲精品国产乱码不卡| 免费观看男人免费桶女人视频| 麻豆精品在线视频| 亚洲中久无码永久在线观看软件 | 亚洲毛片一级带毛片基地| 欧美精品亚洲二区| 男人的天堂久久精品激情| 国产 在线视频无码| 福利小视频在线播放| 国产精品男人的天堂| 亚洲an第二区国产精品| 美女无遮挡被啪啪到高潮免费| 亚洲无线一二三四区男男| 国产精品hd在线播放| 中文国产成人精品久久一| 国产成人精品男人的天堂下载| 欧美中文字幕第一页线路一| 国产一区二区三区在线精品专区 | aⅴ免费在线观看| 久久久国产精品无码专区| 99久久国产综合精品2023| 欧美色伊人| 国产男人天堂| a色毛片免费视频| 91成人免费观看| 国产精品精品视频| 熟女成人国产精品视频| 妇女自拍偷自拍亚洲精品| 亚洲欧美国产五月天综合| 亚洲Av激情网五月天| 亚洲精品免费网站| 亚洲欧美日韩色图| 国产成在线观看免费视频| 国产成人成人一区二区| 亚洲第一视频免费在线| 国产二级毛片| 亚洲经典在线中文字幕| 午夜日b视频| 成人福利在线视频| 国产精品太粉嫩高中在线观看| 国产真实乱人视频| 日韩在线观看网站| 国产丝袜精品| 天堂成人av| 久久久久免费看成人影片| 中文精品久久久久国产网址| 国产91高清视频| 日韩欧美国产三级| 久久成人18免费| 国产成人91精品免费网址在线| 色婷婷亚洲综合五月| 成人午夜视频网站| 三上悠亚一区二区| 亚洲品质国产精品无码| 精品一区国产精品| 先锋资源久久| 国产网站一区二区三区| 色综合a怡红院怡红院首页| 日韩最新中文字幕| 欧美日本一区二区三区免费| 无码日韩人妻精品久久蜜桃| 欧美高清日韩| 成人午夜免费视频| 台湾AV国片精品女同性| 亚洲第一天堂无码专区| 久久精品aⅴ无码中文字幕|