亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        水下航行體舵板張開特性研究與分析

        2016-10-12 07:13:05李四超
        海軍航空大學(xué)學(xué)報 2016年4期
        關(guān)鍵詞:張開航行耦合

        李四超

        (海軍駐鄭州軍事代表室,鄭州450015)

        水下航行體舵板張開特性研究與分析

        李四超

        (海軍駐鄭州軍事代表室,鄭州450015)

        在對流固耦合仿真計算方法分析的基礎(chǔ)上,建立起水下航行體舵板張開過程的仿真計算模型,采用該模型對舵板張開過程的角速度、載荷以及響應(yīng)等變化規(guī)律進(jìn)行了研究。在水下航行體發(fā)射試驗中,對舵板張開角度、舵板應(yīng)變等參數(shù)進(jìn)行了測試,對測試獲取的數(shù)據(jù)與仿真結(jié)果進(jìn)行了對比分析,研究表明二者一致性較好,能夠為水下航行體的水彈道分析和舵板的結(jié)構(gòu)設(shè)計、材料選擇等提供指導(dǎo)。

        舵板;流固耦合;數(shù)值仿真;角速度

        舵板是某水下航行體上的重要部件,當(dāng)水下航行體以一定速度從發(fā)射管進(jìn)入水中時,其尾部的舵板在水流沖擊作用下及時、可靠張開,如圖1所示,實(shí)現(xiàn)對水下航行體姿態(tài)的有效操控,同時承受張開瞬間的巨大沖擊載荷。

        圖1 舵板張開過程示意圖Fig.1 Sketch of rudder opening

        在水下航行體從發(fā)射管發(fā)射入水時,管內(nèi)燃?xì)饬鳟a(chǎn)生的后效、航行體頭部產(chǎn)生的空泡效應(yīng)以及舵板周圍流體形成的紊流等因素耦合作用在舵板結(jié)構(gòu)周圍,形成復(fù)雜的多相流場。因此,舵板的張開過程是在一個極其復(fù)雜的力學(xué)環(huán)境中完成的[1]。

        對舵板張開過程進(jìn)行深入研究和分析,掌握其轉(zhuǎn)動規(guī)律和張開到位時的沖擊響應(yīng)特性,對產(chǎn)品的設(shè)計具有重要意義。而要深入分析舵板張開特性,就必須采用仿真計算和試驗相結(jié)合的方法,才能獲得較好的研究效果。

        1 基于流固耦合的仿真計算

        1.1流固耦合計算模型分析

        舵板張開過程受到航行體運(yùn)動速度、舵板表面壓強(qiáng)以及航行體尾流場的影響;反之,舵板的轉(zhuǎn)動又影響周圍流場的變化,這二者相互作用、相互影響,因而這一過程屬于典型的流固耦合問題。流固耦合的概念在1933年由H.M.Westergaard提出,但真正發(fā)展卻是從20世紀(jì)80年代之后。在近30年中,流固耦合數(shù)值計算經(jīng)歷了解耦—弱耦合—強(qiáng)耦合的發(fā)展過程[2-4],相應(yīng)建立流固耦合模型的方法也分為以下3種。

        1)古典分析法。將內(nèi)部耦合的非線性問題分解為2個獨(dú)立的解耦問題分別求解。因此,從嚴(yán)格意義上講,古典分析法還不是一種真正的流固耦合方法。

        2)交錯積分分析法。首先對流體作用下的結(jié)構(gòu)響應(yīng)進(jìn)行積分,推至下一時刻,將結(jié)構(gòu)邊界位移和運(yùn)動傳遞給流體系統(tǒng),更新流體動態(tài)網(wǎng)格,然后進(jìn)行流體積分計算新的壓力場,再將流體壓力轉(zhuǎn)換成結(jié)構(gòu)載荷,傳遞給結(jié)構(gòu)重新進(jìn)行計算,如此往復(fù)循環(huán)[5]。

        3)整體積分法。將流體和結(jié)構(gòu)看作通過耦合界面連接的單一連續(xù)介質(zhì),用單一的算子描述整個計算場的控制方程。

        以上3種方法中,整體積分法在2個場中的時間積分完全同步,不存在時間滯后和能量不守恒現(xiàn)象,是一種完全的強(qiáng)耦合方法,也是解決流固耦合問題最理想的方法。但在實(shí)際應(yīng)用中,將流體和固體的控制方程完全結(jié)合起來尚存在一定難度,且同步求解的收斂性差,占用機(jī)時資源大,因而該方法目前只應(yīng)用于一些簡單問題的求解,尚未得到廣泛應(yīng)用。

        本文采用交錯積分法建立舵板張開過程的仿真計算模型,其計算過程如圖2所示。

        圖2 舵板張開過程的流固耦合計算Fig.2 Fluid-structure calculation of rudder opening

        1.2仿真計算模型的建立

        1)控制方程。由于本文采用了交錯積分方法開展仿真計算,因而舵板張開運(yùn)動及動力學(xué)計算和流場計算分別采用相應(yīng)的控制方程。舵板運(yùn)動學(xué)方程為:

        式(1)、(2)中:M1、M2、M3分別為流體的推動力矩、緩沖管的反作用力矩和銷軸的摩擦力矩;J為舵板轉(zhuǎn)動慣量。

        舵板運(yùn)動到位后的動力學(xué)方程為:

        綜合權(quán)衡計算對象特點(diǎn)以及計算效率,流場計算采用一方程模型對湍流流動現(xiàn)象進(jìn)行模擬。一方程模型是在時均連續(xù)方程和雷諾方程組成的方程組的基礎(chǔ)上,再建立一個湍動能k的輸出方程,而湍動黏度μt表示成k的函數(shù),從而使得方程組封閉。湍動能k的輸運(yùn)方程可寫為[6-7]:

        式(4)中:ui、uj為速度分量;μ為流體動力粘度;σk、CD為經(jīng)驗常數(shù);l為長度比尺。

        2)幾何模型的簡化。水下航行體及舵板結(jié)構(gòu)復(fù)雜,建模時需對舵板的局部區(qū)域做合理簡化:去掉舵板上凹槽及初始分離機(jī)構(gòu),將航行體表面假設(shè)為平整表面,對舵板根部做適當(dāng)切除,但應(yīng)確保簡化后的舵板轉(zhuǎn)動軸心與實(shí)際轉(zhuǎn)軸中心保持相同,且轉(zhuǎn)動慣量發(fā)生的變化可以被忽略。

        3)計算域模型的確定。計算域是包圍水下航行體和舵板的大圓柱體(水域)。根據(jù)經(jīng)驗,當(dāng)水域直徑為水下航行體直徑20倍時,可滿足計算精度要求,不會因邊界設(shè)置對流場計算產(chǎn)生影響。水域上表面為海水表面,即水域高度為水下航行體的發(fā)射水深。

        4)計算網(wǎng)格的劃分。為保證計算準(zhǔn)確性和收斂性,需合理進(jìn)行網(wǎng)格劃分。水下航行體附近網(wǎng)格較密,其中舵板與水下航行體之間的網(wǎng)格最為密集;遠(yuǎn)離水下航行體、流動變化小的部位網(wǎng)格較稀疏。因采用動網(wǎng)格技術(shù),模型中的網(wǎng)格均為四面體網(wǎng)格[8-9]。網(wǎng)格生成后,應(yīng)用Smoother工具進(jìn)行了優(yōu)化,得到質(zhì)量較高的計算網(wǎng)格,如圖3所示,網(wǎng)格數(shù)量共43萬。

        圖3 計算模型網(wǎng)格示意圖Fig.3 Grid sketch of calculation module

        5)初始條件的確定。計算中,將初始時刻確定為舵板初始張開瞬間。由于舵板模型簡化時已去掉初始分離機(jī)構(gòu),使得舵板在0°時不存在驅(qū)動舵板轉(zhuǎn)動的載荷,故采用初始時刻試算的方法,確定一個合理的初始位置以啟動分析計算。以一個初始的角度進(jìn)行流場迭代,當(dāng)計算能夠收斂即流場狀態(tài)穩(wěn)定時,將此初始角確定為計算初始時刻。經(jīng)仿真試驗,計算初始時刻確定為舵板張開夾角5°。

        6)邊界條件的設(shè)定。將位于水下航行體尾部的發(fā)射筒口設(shè)為壓強(qiáng)入口邊界條件,給定發(fā)射筒噴氣的總壓,隨著航行體向上運(yùn)動,噴氣總壓隨之不斷下降。計算域的上表面(自由表面)為壓強(qiáng)出口條件,靜壓值為一個標(biāo)準(zhǔn)大氣壓。由于流體粘性的存在,固體表面設(shè)定無滑移邊界條件,其余邊界,如計算域的側(cè)面和底面均是由假想表面與海水相交組成,設(shè)定為滑移邊界條件[10-11]。

        1.3仿真計算結(jié)果分析

        采用建立起的仿真計算模型,得到不同時刻舵板的流體載荷及其旋轉(zhuǎn)驅(qū)動力矩,據(jù)此計算出舵板的角加速度和角速度,進(jìn)而確定舵板的運(yùn)動規(guī)律,同時計算出舵板上的應(yīng)力分布,如圖4、5所示。

        圖4 舵板正面應(yīng)力云圖Fig.4 Stress cloud graph of rudder front

        圖5 舵板張開過程仿真計算Fig.5 Calculation of rudder opening

        在舵板張開到位瞬間,受到一個較大的沖擊載荷。此時,在舵板正面兩側(cè)邊緣、舵板耳座等部位出現(xiàn)較大的應(yīng)力,見圖4,最大應(yīng)力達(dá)到1 141 MPa,最小應(yīng)力為422 MPa。經(jīng)分析,舵板兩側(cè)邊緣、舵板耳座應(yīng)力較大與實(shí)際工況是相符的,但因這些部位均處于邊緣或夾角處,故計算誤差可能造成計算值偏大。

        由圖5可知,舵板張開到位時間約為43 ms,張開角度約為76°。在舵板張開過程中呈現(xiàn)明顯的加速運(yùn)動趨勢,前20 ms舵板張開角度約為15°,后20 ms舵板張開角度約為56°,與實(shí)際物理規(guī)律一致:即舵板張開初期,作用在舵板上的轉(zhuǎn)動力矩較小,隨著舵板張開角度的增大,作用在舵板上的水動力開始增加,相應(yīng)的轉(zhuǎn)動力矩也逐漸變大,使得舵板作加速運(yùn)動。舵板張開過程實(shí)測數(shù)據(jù)見圖6。

        圖6 舵板張開過程實(shí)測數(shù)據(jù)Fig.6 Experiment data of rudder opening

        2 舵板張開過程的試驗數(shù)據(jù)分析

        2.1測試方法

        舵板轉(zhuǎn)動測試采用角位移傳感器。角位移傳感器的實(shí)質(zhì)是一個旋轉(zhuǎn)變壓器,其轉(zhuǎn)子固定在測試銷軸上,該測試銷軸和舵板運(yùn)動規(guī)律保證一致;在航行體表面安裝固定基座,角位傳感器的定子安裝在此固定基座上,如圖7所示。當(dāng)舵板由于水流沖擊張開時,舵板測試銷軸的轉(zhuǎn)動就會帶動角位移傳感器轉(zhuǎn)子轉(zhuǎn)動,由于角位移傳感器定子部分被固定,因而在定子和轉(zhuǎn)子之間會產(chǎn)生磁力線切割,角位移傳感器產(chǎn)生電壓信號,經(jīng)數(shù)據(jù)轉(zhuǎn)換器變換和處理后,獲得舵板轉(zhuǎn)動角度-時間曲線。

        圖7 角位移傳感器安裝示意Fig.7 Fixing sketch of angular distance sensor

        依據(jù)仿真計算結(jié)果,在舵板正面兩側(cè)邊緣、中間以及舵板耳座布置了應(yīng)變片。其中,舵板耳座位置為三向應(yīng)變片,其他位置為單向應(yīng)變片,單向應(yīng)變片方向與舵板軸向一致[12-13]。應(yīng)變片布置如圖8所示。

        圖8 應(yīng)變片布置示意圖Fig.8 Placement sketch of strain chip

        2.2試驗數(shù)據(jù)分析

        某發(fā)次試驗中,2個角位移傳感器獲取的數(shù)據(jù)如圖6所示。2個傳感器測得的舵板最大轉(zhuǎn)角分別為81°和74°,與仿真計算結(jié)果相符較好。舵板張開時間分別為83 ms和77 ms,這一數(shù)據(jù)大于仿真計算的43 ms。經(jīng)分析,造成這一現(xiàn)象的原因主要有2點(diǎn):一是在仿真計算中預(yù)置了5°的初始角度,造成舵板實(shí)際張開時間大于計算值,根據(jù)測試數(shù)據(jù)或仿真計算結(jié)果,舵板初始5°轉(zhuǎn)動中消耗的時間大約為20 ms左右;二是仿真計算得出的轉(zhuǎn)動力矩大于實(shí)際值,從舵板的實(shí)際角位移曲線也可以看出,計算得到的舵板角速度要大于實(shí)際舵板運(yùn)行角速度。

        試驗中,1#和3#應(yīng)變片、4#三向應(yīng)變片測試獲取的曲線如圖9所示。

        圖9 舵板應(yīng)力-時間歷程Fig.9 Stress-time course of rudder

        1#、3#應(yīng)變片測得的最大應(yīng)力分別為732 MPa和689 MPa,二者在整個時間歷程中吻合性較好,這說明測試數(shù)據(jù)準(zhǔn)確可信。測試值小于計算值充分驗證了1.3節(jié)中的分析:在計算中,舵板建模時未對其邊緣進(jìn)行倒角處理,導(dǎo)致計算中出現(xiàn)明顯的應(yīng)力集中和相應(yīng)的計算誤差;在實(shí)際工況中,舵板邊緣存在一個較大倒角;此外,實(shí)際測試中應(yīng)變片與舵板邊緣有10mm距離,也使得測試中舵板邊緣的應(yīng)力明顯小于計算值。

        從4#三向應(yīng)變片可以看出,在舵板耳座表面位置存在兩向應(yīng)力,但在一個方向主應(yīng)力明顯大于另外一個方向。

        按式(1)、(2)可分別進(jìn)行計算主應(yīng)力大小和方向:

        式(5)、(6)中:ε1、ε2、ε3為0°、45°、90°方向上測得的應(yīng)變;E為材料彈性模量(MPa);μ為波松比。

        計算得耳座最大主應(yīng)力約為1 100 MPa,超過舵板材料的屈服強(qiáng)度,且接近強(qiáng)度極限。經(jīng)進(jìn)一步分析可知,在數(shù)據(jù)處理中按應(yīng)力-應(yīng)變?yōu)榫€性關(guān)系進(jìn)行處理,即將應(yīng)變片測試值直接乘彈性模量得到該點(diǎn)應(yīng)力。這一處理方法在結(jié)構(gòu)彈性變形范圍內(nèi)是適用的,但若結(jié)構(gòu)發(fā)生塑性變形,這樣的數(shù)據(jù)處理方法就會大大放大應(yīng)力值[14-15]。為取得真實(shí)的應(yīng)力值,重新測試舵板材料的應(yīng)力-應(yīng)變曲線,獲取該材料破斷之前完整的應(yīng)力-應(yīng)變關(guān)系曲線,根據(jù)該曲線得出測試應(yīng)變所對應(yīng)的實(shí)際應(yīng)力。經(jīng)再次處理后,4#位置最大應(yīng)力為864 MPa,略高于屈服強(qiáng)度。這說明該舵板耳座與本體連接處會產(chǎn)生一定塑性變形,但不會影響舵板安全工作。舵板張開到位后正常工作角度時,舵板上的應(yīng)力約為100~200 MPa,遠(yuǎn)小于舵板屈服強(qiáng)度。

        3 結(jié)論

        通過仿真計算和試驗結(jié)果,可得出以下結(jié)論。

        1)舵板到位后張開角度約為80°左右,完全能夠滿足對水彈道的姿態(tài)控制要求,舵板張開呈加速運(yùn)動趨勢,張開到位時間大約在70~80 ms之間,和水下航行體在水中的整個運(yùn)動時間相比是一個小量,因而舵板張開過程不會對水下航行體的水彈道控制造成影響。

        2)舵板張開位移-時間歷程的仿真計算和實(shí)驗實(shí)測數(shù)據(jù)吻合較好。

        3)實(shí)驗結(jié)果表明:舵板的張開過程工作在一個復(fù)雜的紊流環(huán)境中。不同發(fā)次中,舵板相同部位的應(yīng)力存在一定離散。

        4)仿真與實(shí)驗測量均表明:舵板張開到位正常工作角度時,應(yīng)力不足200 MPa,能夠確保其安全可靠工作;但當(dāng)舵板和緩沖裝置接觸時,在舵板上產(chǎn)生一個最大值約為860 MPa的應(yīng)力沖擊,會造成舵板產(chǎn)生一定塑性變形,這就要求在舵板材料選擇時,應(yīng)重點(diǎn)考慮材料的沖擊韌性和塑性指標(biāo)。

        [1]羅金玲,何海波.潛射導(dǎo)彈的空化特性研究[J].戰(zhàn)術(shù)導(dǎo)彈技術(shù),2004(3):14-17. LUO JINLING,HE HAIBO.Research on cavitation for under-water launching missile[J].Tactical Missile Technology,2004(3):14-17.(in Chinese)

        [2]張阿漫,戴紹仕.流固耦合動力學(xué)[M].北京:國防工業(yè)出版社,2011:2-3. ZHANG AMAN,DAI SHAOSHI.Fluid solid interaction dynamics[M].Beijing:National Defense Industry Press,2011:2-3.(in Chinese)

        [3]任弘,李范春,杜玲.流固耦合作用對螺旋槳強(qiáng)度影響的數(shù)值計算[J].武漢理工大學(xué)學(xué)報,2015,39(1):43-52. REN HONG,LI FANCHUN,DU LING.Numerical calculation for the effect of fsi on marine propeller strength[J]. Journal of Wuhan University of Technology,2015,39(1):43-52.(in Chinese)

        [4]YOUNG Y L.Fluid-structure interaction analysis of flexible composite propellers[J].Journal of Fluids and Structures,2008,24:799-818.

        [5]宋學(xué)官,蔡林,張華.ANSYS流固耦合分析與工程實(shí)例[M].北京:中國水利水電出版社,2012:5-9. SONG XUEGUAN,CAI LIN,ZHANG HUA.ANSYS fluid solid interaction analysis and engineering examples [M].Beijing:China WaterPower Press,2012:5-9.(in Chinese)

        [6]周俊杰,徐國權(quán),張華俊.FLUENT工程技術(shù)與實(shí)例[M].北京:中國水利水電出版社,2013:8-13. ZHOU JUNJIE,XU GUOQUAN,ZHANG HUAJUN. FLUENT engineering technology and case[M].Beijing:China Water Power Press,2013:8-13.(in Chinese)

        [7]張磊,郎進(jìn)花,王松嶺.流固耦合問題數(shù)值模擬算法研究進(jìn)展[J].熱力發(fā)電,2015,44(1):1-7. ZHANG LEI,LANG JINHUA,WANG SONGLING. Recnt development of numerical simulation methods for fluid-structure interaction[J].Thermal Power Generation,2015,44(1):1-7.(in Chinese)

        [8]都軍民.潛基導(dǎo)彈發(fā)射動力系統(tǒng)推力能控制技術(shù)研究[D].大連:大連理工大學(xué),2008:72-73. DU JUNMIN.Available energy control in submarinebased missile launch power system[D].Dalian:Dalian University of Technology,2008:72-73.(in Chinese)

        [9]劉維偉,張定華,王軍偉.葉片造型萬個扭曲的校正方法研究[J].機(jī)床與液壓,2004(1):65-67. LIU WEIWEI,ZHANG DINGHUA,WANG JUNWEI. Research on rectification technology of parametric grid distortion in blades modeling process[J].Machine Tool and Hydraulics,2004(1):65-67.(in Chinese)

        [10]DU JUNMIN,HU LIZHONG,KANG NING,et al.Numerical simulation of the add mass and drag of an accelerated motion ball[C]//2014 International Conference on Mechanics and Materials Engineering.Xi’an,2014:104-109.

        [11]李秋實(shí),徐飛,李志平.一種包含運(yùn)動邊界的高精度流場數(shù)值計算方法[J].航空學(xué)報,2014,35(7):1815-1824. LI QIUSHI,XU FEI,LI ZHIPING.A numerical method for simulation flow involving moving boundaries with high order accuracy[J].Acta Aeronautica et Astronautica Sinica,2014,35(7):1815-1824.(in Chinese)

        [12]孫曉丹,侯鋼領(lǐng),王月敏,等.基于靈敏度的平板結(jié)構(gòu)多類型傳感情優(yōu)化布置[J].工程力學(xué),2015,32(4):77-84. SUN XIAODAN,HOU GANGLING,WANG YUEMIN,et al.Optimal placement of multi-type sensor based on sensitivity for plate structure[J].Engineering Mechanics,2015,32(4):77-84.(in Chinese)

        [13]TANG Z Z,LIANG J.Three dimensional digital image correlation system for deformation measurement in experimental mechanics[J].Optical Engineering,2010,49:013601-1-9.

        [14]張川,郭楠.超大變形應(yīng)變測量方法的研究[J].中國測試,2014,40(S1):90-93. ZHANG CHUAN,GUO NAN.Research on large deformation strain measurement method[J].China Measurement&Test,2014,40(S1):90-93.(in Chinese)

        [15]郭芳,楊錄,張艷花.應(yīng)用超聲波技術(shù)對材料應(yīng)力特征的提取與分析[J].聲學(xué)技術(shù),2008,27(2):217-220. GUO FANG,YANG LU,ZHANG YANHUA.Measurement and analysis of material stress with ultrasonic technology[J].Technical Acoustics,2008,27(2):217-220.(in Chinese)

        Research and Analysis for Opening Property of the Rudder Plank on Underwater Vehicle

        LI Sichao
        (Military Representatives Office of Navy in Zhengzhou,Zhengzhou 450015,China)

        Based on the analysis of fluid-structure interaction method,the numerical simulation modular for rudder plank opening was built.With this modular,the research for change rule of the angular velocity,the load and the responding on the opening rudder plank of a certain underwater vehicle was carried out.In the underwater vehicle launching test,the angular displacement-time curve and the strain-time curve were measured.Meanwhile,the test data and simulation calculation data were compared and analyzed.The research showed that the calculation data coincided with the test data,so these data could offer advice for the analysis of trajectory analysis,structure design and material choice.

        rudder plank;fluid-structure interaction;numerical simulation;angular velocity

        TB21

        A

        1673-1522(2016)04-0475-05

        10.7682/j.issn.1673-1522.2016.04.012

        2016-02-28;

        2016-06-03

        國家部委科研基金資助項目(010202)

        李四超(1977-),男,工程師,大學(xué)。

        猜你喜歡
        張開航行耦合
        非Lipschitz條件下超前帶跳倒向耦合隨機(jī)微分方程的Wong-Zakai逼近
        到慧骃國的航行
        開花
        詩潮(2019年10期)2019-11-19 13:58:55
        小舟在河上航行
        航行
        青年歌聲(2017年6期)2017-03-13 00:57:56
        基于“殼-固”耦合方法模擬焊接裝配
        大型鑄鍛件(2015年5期)2015-12-16 11:43:20
        Dynamical Properties of a Diluted Dipolar-Interaction Heisenberg Spin Glass?
        求解奇異攝動Volterra積分微分方程的LDG-CFEM耦合方法
        非線性耦合KdV方程組的精確解
        视频一区视频二区亚洲| 国产精品亚洲二区在线看| 精品一区二区三区芒果| 国产高清av首播原创麻豆| 99久久精品免费看国产情侣| 国产欧美日本亚洲精品一5区| 看国产亚洲美女黄色一级片| 中文字幕亚洲综合久久| 中国农村妇女hdxxxx| 国产在线一区观看| 成人影院免费观看在线播放视频| 国产精品人妻熟女男人的天堂 | 狠狠精品久久久无码中文字幕| 少妇被粗大的猛进69视频| 亚洲AV无码一区二区三区少妇av| 精品成人av人一区二区三区| 欧美丰满熟妇bbb久久久| 久久99精品国产99久久| 精品黑人一区二区三区| 国产一区二区长腿丝袜高跟鞋| 久久久久无码精品国产app| 亚洲综合国产精品一区二区99| 国产精品成人有码在线观看| 久久综合噜噜激激的五月天| 白又丰满大屁股bbbbb| 综合91在线精品| 男女打扑克视频在线看| 亚洲精品视频1区2区| 久久久久国色av免费观看性色| 精品国内自产拍在线观看| 亚洲欧美久久婷婷爱综合一区天堂| av天堂网手机在线观看| 亚洲性无码一区二区三区| 亚洲狠狠网站色噜噜| 一区二区三区在线观看视频免费| 大陆老熟女自拍自偷露脸| 在线看片免费人成视频电影| 中文字幕久无码免费久久| 亚洲欧洲综合有码无码| 国产性色av一区二区| 巨大巨粗巨长 黑人长吊|