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

船—纜拖曳系統操縱性能分析

2015-04-26 08:07:43孫洪波施朝健林文錦
船舶力學 2015年11期
關鍵詞:船舶系統

孫洪波,施朝健,林文錦

(1上海海事大學商船學院,上海200135;2集美大學航海學院,福建廈門361021;3中海油服,天津300450)

船—纜拖曳系統操縱性能分析

孫洪波1,2,施朝健1,林文錦3

(1上海海事大學商船學院,上海200135;2集美大學航海學院,福建廈門361021;3中海油服,天津300450)

為獲取拖船在拖曳時的操縱性的變化規律。文章采用MMG船舶運動數學模型的建模思想,建立了六自由度拖船運動數學模型,采用有限差分法,建立了拖纜模型。然后,在此基礎上建立將船-纜耦合起來以形成整個系統的運動數學模型,并分別采用龍格庫塔方法對船舶運動積分求解,采用后向差分法對拖纜運動進行求解。通過對比仿真計算分析了水面拖船在拖帶過程中的加速性能、旋回性能及偏轉抑制性能。仿真結果表明在拖船與拖纜的相互影響下,拖船的加速性能和旋回性能有所下降而偏轉抑制性能有所增強。

船舶運動數學模型;有限差分;加速性能;旋回性能;偏轉抑制性能

0 引言

在各種海洋資源勘探及海底電纜鋪設活動中,拖船拖曳是不可避免的環節。在海上,相比水下拖纜部分,水面拖船和拖纜上端靠近水面的部分受到風、波浪和海流的影響很大;而且水面拖船、拖纜之間存在著極為復雜的相互作用,構成了一個復雜的動態響應系統。對這個動態系統的動態特性進行研究,不僅可以更好地對拖纜進行控制,保證拖曳系統設備的可靠性和安全性,而且對保證拖船的安全操縱與航行具有理論上的指導作用和重要的工程應用價值。

國內外已有很多學者分別針對拖船和拖纜的運動建模與仿真預報進行了大量的理論研究。對于水面船體的研究,從上世紀50年代野本[1]建立了描述系統輸入輸出關系的響應模型、上世紀60年代以Abkowitz[2]為代表的整體建模理論的提出,以及上世紀70年代以小川陽弘等[3]為代表的MMG模型的提出,發展到現階段,船舶運動模型雖然在某些方面的模型精度還有提高的空間,但已能滿足工程及船員培訓的需要,并得到廣泛的應用。國內的賈欣樂、楊鹽生[4]對船舶運動建模進行了詳細而系統的總結與研究。水下拖纜自身運動的基礎研究,早在上世紀60年代就已廣泛地展開。拖纜水動力模型建模方法主要包括:集中質量法、有限元法、有限差分法和直接積分法等。例如:Walton和Polchaton[5]首先采用集中質量法研究了海軍武器試驗中水下錨鏈的二維運動響應,并給出詳細的求解算法和公式。中島俊夫等[6]建立了船舶系留錨鏈的三維集中質量法模型,并給出了詳細的算例。Ablow等[7]在建立拖纜運動控制微分方程的基礎上,采用有限差分法求解其的三維動態運動,與集中質量法相比空間與時間步長的選擇具有更大的靈活性。有限元法和直接積分法更適用于拖纜靜力分析,而不適合描述非定常狀態下拖纜水動力性能[8]。國內的王飛、朱軍等[9-10]在拖纜的建模與仿真方面,分別對集中質量法和有限差分法進行了詳細的總結和研究;文獻[11]采用了考慮拖纜彎矩和扭矩的有限差分法,可以對拖纜在低應力情況下進行計算。

通常,拖曳系統由水面拖曳母船、水下拖纜及一些纜載設備等組成,必需考慮船纜之間的相互影響才能準確預報其在各種情況下的運動響應特性。但現有的眾多相關研究,考慮到拖纜相對于拖船質量比較小而忽略了船纜之間的相互作用,這必定會對整個系統的預報精度產生影響。文獻[12]采用數值模擬手段,詳細分析了拖纜對四自由度拖曳船舶操縱性能的影響;文獻[13]采用集中質量法,對三自由度船舶與拖纜完全耦合的運動響應做了詳細分析。本文在此基礎上,首先采用MMG建模思想建立六自由度船舶運動模型,采用有限差分法建立拖纜模型,然后建立船纜相互作用,并考慮風浪影響的船纜拖帶系統模型,最后,對該系統的動態響應進行仿真分析。

1 拖船和拖纜系統運動數學模型

研究船舶和拖纜在水中的運動時,一般采用兩種坐標系統:慣性坐標系統和附體坐標系統,如圖1所示。圖中,o0x0y0z0為固定于地球表面的慣性坐標系統,規定x0指向正北,y0指向正東,z0垂直于地表指向上;ogxyz是原點位于船舶重心的附體坐標系。規定x指向船首,y指向右舷,z垂直于龍骨指向上;octnb是原點位于拖纜上的附體坐標系,規定t為拖纜切向,n和b為拖纜的兩個法向;Θ、Φ、ψ為船舶在ogxyz坐標系下的歐拉角,θ、φ為拖纜在octnb坐標系下的歐拉角;vw、ψw、J和ψc分別為風速、風向、流速和流向。

圖1 拖船和拖纜坐標系統Fig.1 Ship and cable coordinates system

1.1 船舶運動數學模型

船舶的實際運動是一個具有6個自由度的復雜的運動。對于大多數情況下的船舶運動及控制而言,可以忽略船舶垂蕩、橫搖和縱搖運動,為了盡可能準確地反應船舶實際的運動姿態,本文采用考慮風、流、浪影響的六自由度船舶運動數學模型,其動力學方程如公式(1)。式中:X、Y、Z、K、M、N分別為作用于船體的力和力矩;m、mx、my、mz分別為船舶質量和各船舶附體坐標軸方向上的附加質量;I、Ia結合下角標分別為繞各船舶附體坐標軸的轉動慣量和附加轉動慣量。下角標H、P、R、cable、wind、current和wave分別表示船體、螺旋槳、舵、纜、風、流和波浪。公式(1)中各參數及力的詳細計算可參考文獻[4,14]。

其中:u、v、w、p、q、r分別為船舶在ogxyz坐標系下的縱向、橫向和垂向速度和繞各船舶附體坐標軸的角速度,與慣性坐標系內船舶運動學物理量之間的關系分別為公式(2)和(3)。

1.2 拖攬運動數學模型

設拖纜是細長、柔性和圓柱形狀。拖纜上受重力、水動力和慣性力作用,由拖船拖動前進,不考慮拖纜所受的扭矩。拖纜微元Δs在t時刻所受張力、慣性力、水中重力和流體力分別為和以s表示未拉伸的拖纜長度,S是拉伸后的拖纜長度。考慮一般拖纜材料的拉伸特性有S′=墜S/墜s=1+eT,其中e=1/EA,E為拖纜彈性模量;A為拖纜橫截面積。由于拖纜幾乎全部位于水面以下,因此忽略了拖纜所受波浪力的影響。根據牛頓第二定律,拖纜微元的動力平衡方程如下:

同時聯立方程(4)和(6),并將其在附體坐標系下展開,可分離得到6個獨立的分量方程。這樣,可以由6個變量T、Vt、Vn、Vb、θ和φ描述拖纜的三維運動。最終可以將拖纜系統的6個分量方程寫成如下矩陣形式,其具體推導過程祥見文獻[15]。

其中:ρ為海水密度;g為重力加速度;A為拖纜橫截面積;m為每米拖纜質量;ρA為每米拖纜的附加質量;w=(m-ρA)g為水中每米纜質量;d=(1+eT)1/2d0為拉伸后纜的直徑,d0為未拉伸纜的直徑;Ct、Cn和Cb為拖纜的水動力系數;考慮均勻流的影響可知Ut=Vt+Jt,Un=Vn+Jn,Ub=Vb+Jb。

2 船—攬拖曳系統運動數學模型

為求解拖纜的三維動態運動,采用有限差分方法離散方程組(7)。設纜總長為s,把纜分成n段,節點數為n+1,然后將方程組按空間和時間兩個方向上作向后差分,得到由6n個含有6n+6個未知變量方程組成的差分方程組。差分方程組只有6n個,故需補充拖曳系統端點邊界條件使方程組得到封閉。

2.1 拖纜自由端邊界條件

在拖纜自由端,即s=0處,拖纜張力與歐拉角的變化率均為零。因此,在自由端的邊界條件可以直接表示為:

2.2 隨船拖帶點邊界條件

拖帶點處的邊界條件由在慣性坐標系中,拖纜拖帶端的三個速度分量VCx0、VCy0和VCz0與拖船上系固點的三個速度分量VTx0、VTy0和VTz0相等得到。前者具體形式可由公式(6)得到。而后者可根據拖船拖帶點P(xP,yP,zP)距船舶重心的位置,通過解算船舶運動模型得到,其具體表達如下:

2.3 船—纜耦合條件

在拖曳系統運動過程中,拖船與拖纜是相互作用、相互影響的。拖船提供張力拖帶著拖纜向前運動,反過來拖纜也影響了拖船的運動。因此,為進一步提高船纜拖曳系統在不同情況下的預報精度,必須將拖船和拖纜視為一個相互影響的整體進行考慮。船—纜耦合的主要方法有兩種,一種是文獻[16]采用的方法,即將拖船的運動方程離散成拖纜的方程形式,然后聯立成方程組一并求解。該方法雖然能夠保證船與纜之間的相互影響在時間步長上的同步性,但使整個系統的運動方程復雜化,并降低了解算效率。

本文的做法是將拖船與拖纜兩個部分分別求解。拖船部分采用龍格庫塔法,而拖纜部分采用有限差分法。這種處理方法思路清楚,物理意義明確,求解算法簡單明了。雖然在解算過程中,是將拖纜前一時刻的張力作為作用在拖船上的影響力,但對于船纜拖曳這種低頻運動系統,在較小的時間步長下仍能保證足夠的精度。

具體的方法是根據(7)、(8)式和(9)式離散后的拖纜非線性運動方程組,在給定了拖點速度條件下采用Newton迭代法解算[15],求得拖纜的形態及張力。在拖帶點P(xP,yP,zP)已知的情況下,將求得的拖點張力經慣性坐標和船體坐標變換后,得到各船體坐標軸上的力和力矩,如下式:

將(10)式代入到(1)式中,用龍格庫塔法解算(1)式的非線性微分方程組,從而得到拖船在拖纜影響下的操縱性能。纜索的動力學方程式(7)是一個一階的非線性偏微分方程組。對其數值求解采用后向差分格式,如以i和j分別表示時間節點和空間節點,先在時域內對該方程利用后向差分格式進行離散,再在空間域內對方程仍采用中心差分格式離散。最后整理可得后向差分格式為:

對此非線性代數方程組(11),利用Newton-Raphson迭代法在時域內求解,即可求出任意時刻水中拖纜的運動姿態。

3 仿真計算

拖纜運動的計算涉及解算大型差分方程組,利用MATLAB軟件結合并行算法對所建船—纜拖曳系統進行求解仿真。仿真所用拖船及拖纜的主尺度如表1和表2。

表1 拖船主尺度Tab.1 Towing ship main dimensions

表2 拖纜主尺度Tab.2 Towed cable main dimensions

為獲取拖船在拖曳時的操縱性的變化規律,以及進一步提高水下拖曳系統在不同情況下的預報精度,分別對船—纜拖曳系統進行了加速、旋回和改向仿真試驗。仿真結果如圖2~4所示。

圖2 系統加速過程仿真Fig.2 Acceleration simulation of towing system

變速性能是度量船舶慣性的技術指標,其包括加速、減速、停船和倒航性能,本文僅對拖曳系統的加速性能進行仿真。圖2為系統仿真初始速度為零,螺旋槳RPM為50的仿真結果。從圖2(a)中可以看出在拖纜的影響下,與無拖纜拖船相比,在相同的400 s加速過程中,隨拖纜長度的增加,拖船的最終速度有所減小。圖2(b)為1 000 m拖纜分為5段,6個節點的加速過程中的姿態變化。可以看出,隨拖船速度的增加,各節點的浸水深度逐漸減小。圖2(b)也給出了拖纜末端的沉深變化,同時可以看出船纜相互作用對拖纜末端沉深的影響隨船速的增加差異越明顯,加速末期,拖纜末端沉深相差近30 m。因此,為提高水下拖曳系統的預報精度,在加速階段必須考慮船纜相互影響。

船舶的旋回性是船舶最基本的重要操縱性能之一,通常采用旋回初徑DT與船長L之比DT/L,即相對旋回初徑來衡量。圖3為系統仿真初始速度為12 kns,螺旋槳RPM為124的仿真結果。從圖3(a)中可以看出在拖纜的影響下,隨拖纜長度的增加,拖船的旋回圈增大,且旋回舵角越小差異越明顯,橫距比進距變化明顯。相對旋回初徑如表3所示。

圖3 系統旋回過程仿真Fig.3 Turning simulation of towing system

表3 拖帶系統相對旋回初徑Tab.3 Relative tactical diameter of towing system

圖3(b)為500 m拖纜20°舵角旋回過程中的姿態變化。可以看出,由于拖船在旋回過程中存在明顯的降速,因此,隨拖帶點速度的下降,各節點的浸水深度逐漸增加。圖3(b)給出了拖纜末端的沉深變化,變化幅度近150 m,如直接認為旋回過程中速度不變,這顯然是不符合實際情況的,將會引起不小的誤差。同時可以看出船纜相互作用對拖纜末端沉深的影響,由于旋回過程中,拖船在有拖纜的情況下,由于船尾在拖纜的拉力及力矩的作用下降速明顯。因此,旋回末期受拖纜影響減小,但相對于加速過程,拖纜末端的沉深有近5 m的差異。

為分析拖纜對拖船偏轉抑制性即保向性的影響,本文采用10°/10°Z形仿真試驗。系統仿真初始速度為12 kns,螺旋槳RPM為124,仿真結果如圖4。從圖4(a)中可以看出,無拖纜情況下的航向超越角與拖帶500 m拖纜的情況差別不大,而與拖帶1 000 m拖纜的情況差別明顯。根據Z形仿真試驗結果,可以方便地計算出上述三種仿真試驗工況下拖船的旋回性指數K和追隨性指數T,具體計算方法見文獻[17]。計算結果如表4所示,其中:K′=KL/Vs;T′=KVs/L為K和T的無量綱化,L為船長,Vs為船速。計算結果表明拖船在拖纜的影響下,旋回性指數K減小說明拖船的旋回性能變差,且拖纜越長旋回性越差。這一點同旋回仿真試驗的結果是一致的;追隨性指數T減小說明船舶應舵快,易于保向。總體來說,拖船的旋回性變差,追隨性變好,但操舵后,雖然應舵較快,但定常旋回角速度小,旋回圈大。

從圖4(b)中可以看出,由于拖船在改向時存在降速而引起的拖纜末端沉深的變化。同時可以看出船纜相互作用對拖纜末端沉深的影響比較明顯,Z形仿真試驗末期,拖纜末端沉深相差近2.5 m。

圖4 系統改向過程仿真Fig.4 Altering course simulation of towing system

表4 拖帶系統操縱性指數Tab.4 Maneuverability indices of towing system

4 結論

本文為獲取拖船在拖曳時的操縱性的變化規律,以及進一步提高水下拖曳系統在不同情況下的預報精度,將水面拖船和水下拖纜視為一個相互影響的整體進行考慮。其中,水面拖船的運動采用六自由度的MMG模型,而水下拖纜則采用有限差分法進行建模,在此基礎上建立了考慮拖船與拖纜的相互影響船-纜拖曳系統運動數學模型,并采用四階龍格庫塔方法和后向差分法對拖船與拖纜兩個部分分別求解。最后,分析了包括水面拖船在內的整個系統的操縱運動響應。仿真試驗結果表明,受拖纜影響,拖船加速性能下降,拖纜長度越長速度下降越明顯;拖船的旋回直徑增大,旋回性能隨纜長的增加而下降,且旋回舵角越小差異越明顯;拖船的偏轉抑制性有所增強。此外,由于拖船在旋回和改向過程中存在降速的影響,拖纜的姿態變化更為符合實際的運動特征。

[1]野本謙作,田口賢士,本田啓之輔,平野進.船の操性に就いて(1)[J].日本造船協會論文集,1956,99:75-82. Nomoto Kensaku,Taguchi Kenshi,Honda Keinosuke,Hirano Susumu.On the steering qualities of ships(1)[J].The Japan Society of Naval Architects and Ocean Engineers,1956,99:75-82.

[2]Abkowitz M A.Lectures on ship hydrodynamics steering and manoeuvability[R].Hydro-and Aerodynamics Laboratory Report No.Hy-5,Lyngby,Denmark.1964.

[3]小川陽弘,小山健夫,貴島勝郎.MMG報告-I操縱運動の數學モデルにつぃて[J].日本造船學會志,1977,575:192-198. Ogawa Youko,Koyama Keno,Takashima Katsuro.MMG Rep-1 mathematical model for motions of ships[J].Journal of The Society of Naval Architects of Japan,1977,575:192-198.

[4]賈欣樂,楊鹽生.船舶運動數學模型—機理建模與辨識建模[M].大連:大連海事大學出版社,1999:234-246. Jia Xinle,Yang Yansheng.Ship motion mathematical model[M].Dalian:Press of Dalian Marine University,1999:234-246.

[5]Walton T S,Polachech H.Calculation of transient motion of submerged cables[J].Mathematics of Computation,1960,14: 27-46.

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

[8]苑志江,金良安,田恒斗,盧祎斌.海洋拖曳系統的水動力理論與控制技術研究綜述[J].科學技術與工程,2013,13 (2):408-420. Yuan Zhijiang,Jin Liang’an,Tian Hengdou,Lu Yibin.Comments on the research of hydrodynamic and controltechnology of underwater towed system[J].Science Technology and Engineering,2013,13(2):408-420.

[9]王飛.海洋勘探拖曳系統運動仿真與控制技術研究[D].上海:上海交通大學,2007. Wang Fei.Simulation and control research of marine towed seismic system[D].Shanghai:Shanghai Jiao Tong University, 2007.

[10]朱軍,熊鷹,王志國.拖纜系統直線定常運動仿真計算[J].海軍工程大學學報,2001,13(2):17-20. Zhu Jun,Xiong Ying,Wang Zhiguo.Towed cable behaviour during straight moving[J].Journal of Naval University of Engineering,2001,13(2):17-20.

[11]Liu Tao,Zhang Weijing,Ma Jie,Zhang Guanglei.Transient dynamic analysis of towed low-tension cable with experimental verification[J].Journal of Ship Mechanics,2013,17(3):197-213.

[12]朱軍,黃若波,胡忠平.拖曳系統對艦船操縱性影響計算[J].船海工程,2002,145(2):5-10. Zhu Jun,Huang Ruobo,Hu Zhongping.Calculation of towing system effects on ship maneuvering[J].Ship&Ocean Engineering,2002,145(2):5-10.

[13]王飛,黃國樑,伍生春.水下拖曳系統纜-船耦合運動模擬[J].上海:上海交通大學學報,2011,45(4):571-575. Wang Fei,Huang Guoliang,Wu Shengchun.Dynamic research on the coupling response of cable-towing ship system[J]. Journal of Shanghai Jiao Tong University,2011,45(4):571-575.

[14]李積德.船舶耐波性[M].哈爾濱:哈爾濱船舶工程學院出版社,1991. Li Jide.Seakeeping[M].Harbin:Press of Harbin Engineering University,1991.

[15]鄧德衡,黃國樑,樓連根.拖曳線列陣陣形與姿態數值計算[J].海洋工程,1999,17(1):18-26. Deng Deheng,Huang Guoliang,Lou Liangen.Numerical calculation of towed array shape and attitude[J].Ocean Engineering,1999,17(1):18-26.

[16]李英輝,連璉.拖曳聲納系統仿真運算研究[C].第十二屆中國海岸工程學術討論會論文集,2005. Li Yinghui,Lian Lian.The towed sonar system simulation computing research[C].The 12th Academic Conference on China Coastal Engineering,2005.

[17]古文賢.船舶操縱[M].大連:大連海事大學出版社,1993:53-55. Gu wenxian.Ship Handling[M].Dalian:Press of Dalian Marine University,1993:53-55.

Analysis of maneuverability of towing cable ship system

SUN Hong-bo1,2,SHI Chao-jian1,LIN Wen-jin3
(1.Merchant Marine College,Shanghai Maritime University,Shanghai 200135,China;2.Navigation College, Jimei University,Xiamen 361021,China;3.COSL,Tianjin 300450,China)

For obtaining the regular rules of a cable towing vessel,the following methods are adopted.According to the MMG ship motion mathematical modeling theory,a 6dof ship motion mathematical model was built.And a 6dof cable motion mathematical model was established with finite difference method.Then,a coupling ship-cable system model was established.Finally,the system motion simulations containing acceleration,turning and zigzag maneuvering were carried out,and the results indicate that the abilities of acceleration and turning are decline,but yaw checking ability is enhanced.

ship motion mathematical model;finite difference method;acceleration ability; turning ability;yaw checking ability

U661.33

A

10.3969/j.issn.1007-7294.2015.11.005

1007-7294(2015)11-1325-09

2015-06-13

國家自然科學基金資助項目(51109090);李尚大集美大學學科建設基金資助項目(ZC2010011)

孫洪波(1977-),男,講師,博士研究生,E-mail:sunhongbo1977@126.com;

施朝建(1951-),男,教授,博士生導師。

猜你喜歡
船舶系統
計算流體力學在船舶操縱運動仿真中的應用
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
基于改進譜分析法的船舶疲勞強度直接計算
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
船舶!請加速
基于PowerPC+FPGA顯示系統
BOG壓縮機在小型LNG船舶上的應用
半沸制皂系統(下)
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
主站蜘蛛池模板: 亚洲第一成网站| 中文字幕av一区二区三区欲色| 免费国产在线精品一区| 日本不卡免费高清视频| 国产成人亚洲欧美激情| 国产又粗又爽视频| 欧美一区二区福利视频| 亚洲中文字幕23页在线| 久久精品亚洲中文字幕乱码| 99re在线观看视频| 最新午夜男女福利片视频| 99中文字幕亚洲一区二区| 国产在线观看人成激情视频| 久久精品国产精品一区二区| 国产福利拍拍拍| 91av国产在线| 国产JIZzJIzz视频全部免费| 国产精品冒白浆免费视频| 日韩成人在线一区二区| 成人年鲁鲁在线观看视频| 这里只有精品在线| 日韩福利在线观看| 伦伦影院精品一区| 最新痴汉在线无码AV| 国产网友愉拍精品| 亚洲天堂免费在线视频| 又黄又爽视频好爽视频| A级毛片高清免费视频就| 六月婷婷精品视频在线观看| 国产流白浆视频| 在线一级毛片| 亚洲视频三级| 成人欧美在线观看| 亚洲AV无码一二区三区在线播放| 97青青青国产在线播放| 国产高清不卡视频| 国产成人在线无码免费视频| 福利国产在线| 国产欧美日韩专区发布| 欧美日韩精品一区二区在线线| 亚洲最大福利视频网| 中日韩一区二区三区中文免费视频| 天天综合网站| 波多野结衣无码中文字幕在线观看一区二区| 亚洲第七页| 国产高清在线观看91精品| 国产精品九九视频| 99re66精品视频在线观看| 国产精品自在线拍国产电影| 亚洲精品无码久久毛片波多野吉| 久久99精品久久久大学生| 99这里只有精品免费视频| 国产成人超碰无码| 青青青国产免费线在| 日韩免费毛片| 97se亚洲综合不卡| 无码网站免费观看| 欧美特黄一级大黄录像| 国产丝袜啪啪| 国产网站免费| 久草青青在线视频| 国产精品流白浆在线观看| 欧美爱爱网| 国产精品尤物铁牛tv | 毛片一区二区在线看| 国产一区二区三区在线观看免费| 色AV色 综合网站| 婷五月综合| 99免费视频观看| 91亚瑟视频| 一本综合久久| 亚洲AV无码不卡无码 | 992Tv视频国产精品| 精品综合久久久久久97超人| 一区二区三区精品视频在线观看| 亚洲日韩AV无码一区二区三区人 | 欧美19综合中文字幕| 亚洲综合色婷婷中文字幕| 思思99思思久久最新精品| 国产av剧情无码精品色午夜| 国产日本视频91| 国产一区二区网站|