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

        ?

        直升機槳葉鳥撞試驗方法的數(shù)值模擬?

        2012-12-26 09:08:26林長亮王浩文陳仁良
        振動工程學報 2012年6期
        關鍵詞:模型

        林長亮,王浩文,陳仁良

        (1.南京航空航天大學直升機旋翼動力學國家級重點實驗室,江蘇 南京 210016;2.清華大學航天航空學院,北京 100084)

        引 言

        鳥撞事故一般發(fā)生在 600 m以下的空域,而低空飛行是直升機的顯著使用特點,因此直升機發(fā)生鳥撞的可能性很大,并且隨著飛行速度的提高,鳥撞事故的危害性也在逐步加大。Dolbeer在第8屆鳥撞會議上指出[1],從 1999至 2005年,直升機鳥撞事故共 370起,其中 186(大于 50%)起引起了事故,67(18%)起引起了嚴重事故。從這些數(shù)據(jù)可以看出鳥撞事故對直升機的飛行安全有很大影響。發(fā)達國家越來越重視民用直升機使用過程中的鳥撞危害性問題,要求研制單位在民用直升機研制過程中必須滿足抗鳥撞技術要求,并在適航條例中加入了相關條款。2002年 7月 2日中國民航總局頒布了新版的CCAR29,在 CCAR29.631條款中明確規(guī)定民用直升機抗鳥撞設計要求。新條例的頒布一方面使中國民航適航條例與國際接軌,有利于國產(chǎn)民用直升機盡快進入國際市場,另一方面也對民用直升機的研制提出了新的技術要求。

        對于鳥撞問題的研究,在飛機鳥撞方面,國內(nèi)外學者作了大量的研究工作。 Wilbeck和 Barber采用10%孔隙率的明膠代替真鳥進行試驗研究[2,3],發(fā)現(xiàn)鳥撞過程可被描述成一個非恒定的流體動力學過程。Meguid等通過計算研究了鳥體形狀對計算結(jié)果的影響[4],發(fā)現(xiàn)鳥體和目標的最初接觸面積對接觸力的峰值有巨大影響。張志林、姚衛(wèi)星等結(jié)合鳥撞風擋試驗[5],討論了鳥彈質(zhì)量及風擋厚度變化對鳥撞風擋動響應的影響規(guī)律。朱書華等對鳥體形狀對風擋動態(tài)響應的影響進行了研究[6]。萬小朋等對飛機機翼前緣的抗鳥撞性能進行了分析[7]。陳偉、關玉璞等對發(fā)動機葉片鳥撞的瞬態(tài)響應計算以及損傷試驗進行了研究[8,9]。

        而在直升機鳥撞方面,國外僅有西科斯基公司在 S-92直升機抗鳥撞設計中對垂尾前緣、尾槳鳥撞進行了動力學仿真,同時開展了鳥撞試驗驗證工作[10,11]。目前,中國對直升機結(jié)構(gòu)的鳥撞問題的研究還很少。王益鋒等采用有限轉(zhuǎn)動梁理論處理槳葉的彈性變形[12],引入位移協(xié)調(diào)方程識別彈性碰撞載荷,運用Hamilton原理建立了槳葉正碰撞的動力學方程,溫海濤等對直升機主槳葉鳥撞的過程進行了數(shù)值模擬[13,14],而在直升機結(jié)構(gòu)鳥撞試驗方面還沒有見到相關報道。

        本文通過對旋翼旋轉(zhuǎn)狀態(tài)下的鳥撞分析,提出采用集中載荷作用下的槳葉靜態(tài)鳥撞試驗替代旋轉(zhuǎn)狀態(tài)下槳葉鳥撞試驗。采用MSC.Dytran軟件,基于 ALE流固耦合方法對兩種狀態(tài)下的直升機槳葉鳥撞進行了數(shù)值模擬。通過數(shù)值模擬發(fā)現(xiàn)兩種狀態(tài)下的葉鳥撞動態(tài)響應一致性很好,說明了替代方法是有效的,從而為直升機旋翼鳥撞試驗提供了一種理論依據(jù)。

        1 鳥撞試驗

        1.1 槳葉鳥撞試驗方法

        由于直升機旋翼直徑較大,并且在旋轉(zhuǎn)過程中槳葉具有揮舞、擺振、扭轉(zhuǎn)等運動,邊界條件復雜,受試驗條件和試驗場地的限制,主槳葉全尺寸鳥撞試驗很難實現(xiàn),本文提出采用集中載荷作用下的槳葉靜態(tài)鳥撞試驗替代旋轉(zhuǎn)狀態(tài)下槳葉鳥撞試驗。

        槳葉在旋轉(zhuǎn)狀態(tài)與靜止狀態(tài)下發(fā)生鳥撞的主要區(qū)別在于以下 3點:1.槳葉以恒定的角速度旋轉(zhuǎn),由于離心載荷的作用,會使槳葉的剛度有所增加;2.槳葉上會產(chǎn)生氣動升力;3.在旋轉(zhuǎn)過程中槳葉具有揮舞、擺振、扭轉(zhuǎn)等運動。

        在靜止條件下模擬旋轉(zhuǎn)槳葉的鳥撞過程,可以通過:1.采用集中載荷替代離心載荷;2.由于撞擊是瞬態(tài)發(fā)生的,其產(chǎn)生的撞擊力要遠遠大于撞擊部位上槳葉的氣動力,因此忽略氣動力的影響以及槳盤入流的變化;3.由于鳥撞持續(xù)時間很短(大約3 ms),鳥體沖量有限,認為在鳥撞過程中槳葉的運動姿態(tài)沒有明顯改變,即保持鳥撞之前的狀態(tài)。在試件安裝過程中,可以根據(jù)槳葉撞擊部位在運動中的姿態(tài)進行調(diào)整。

        1.2 試驗裝置

        飛機結(jié)構(gòu)部件鳥撞試驗技術已經(jīng)很成熟,并制定了相關標準[15]。鳥撞試驗裝置主要由空氣炮和動態(tài)數(shù)據(jù)測量系統(tǒng)組成。圖 1為鳥撞試驗裝置示意圖[16]。試件固定在試驗臺上,然后通過空氣炮將鳥彈發(fā)射出去撞擊試件。測量系統(tǒng)一般包括位移傳感器、超動態(tài)應變儀以及高速攝影機等。多臺高速攝影機可以從不同方向拍攝鳥與試件的撞擊過程。

        圖1 鳥撞試驗裝置示意圖Fig.1 Schematic diagram of bird impact facility

        在飛機鳥撞試驗裝置上修改試驗臺試件的夾持機構(gòu),使之可以施加集中載荷,便可開展槳葉鳥撞試驗。

        2 槳葉鳥撞數(shù)值模擬

        CCAR29.631條款要求,旋翼航空器設計必須滿足:在 2 440 m高度,速度等于Vne或Vh(取較小者)時,受到 1.0 kg(2.2 lb)的鳥擊后能繼續(xù)安全飛行和著陸(對 A類)或安全著陸(對 B類)。 本文以 1 kg的鳥體以298 m/s(其中選取鳥飛行速度為15 m/s,直升機前飛速度為85 m/s,撞擊槳葉0.9R處的切向速度為198m/s)的速度沿水平方向撞擊某型直升機槳尖段前緣為例,對兩種狀態(tài)下的鳥撞過程進行數(shù)值模擬,通過槳葉動態(tài)響應結(jié)果的對比驗證本文提出的試驗方法有效性。

        2.1 ALE流固耦合方法

        最初用于求解沖擊問題的計算機程序主要分為兩類:拉格朗日法和歐拉法[17]。在拉格朗日法中,計算網(wǎng)格固定在物質(zhì)上,網(wǎng)格點與物質(zhì)點在物體的變形過程中始終保持重合;在歐拉法中,計算網(wǎng)格固定在空間,在物體變形過程中保持不變,而材料相對于網(wǎng)格運動。近 20年發(fā)展起來的任意拉格朗日-歐拉(ALE)耦合方法,將拉格朗日和歐拉法結(jié)合,材料在歐拉網(wǎng)格中移動的同時歐拉網(wǎng)格節(jié)點本身也在運動,使得歐拉網(wǎng)格的位置和形狀在不斷調(diào)整。使得ALE方法容易處理網(wǎng)格畸變,能夠精確地描述接觸滑移面,適合于處理接觸碰撞的大變形與材料破壞等復雜力學問題。

        2.2 鳥體模型

        2.2.1 鳥體幾何模型

        國家軍用標準GJB2464-95規(guī)定:鳥撞試驗中的鳥彈外形長徑比為2∶1的圓柱體。故本文采用形狀長徑比為 2∶1的圓柱體來模擬鳥體,圓柱體底面直徑為 88 mm,長度為 176 mm。

        2.2.2 鳥體材料模型

        鳥體的材料特性是鳥撞仿真分析的重點和難點。真實鳥體的本構(gòu)方程很難描述,目前尚未見到“準真實”鳥體的模擬報道,在實際操作中往往根據(jù)撞擊速度的大小采用彈性體、彈塑性體、理想流體等模型來模擬鳥體本構(gòu)模型。在鳥體與槳葉的碰撞過程中,相對速度很大,鳥體迅速減速而產(chǎn)生的應力遠遠大于鳥體材料的屈服強度,鳥體發(fā)生流變,呈現(xiàn)了流體特性,因此在數(shù)值模擬過程中,鳥體材料采用理想流體模型,鳥體密度d=930 kg/m3,體積模量K=2 200 MPa。

        2.2.3 鳥體模型驗證

        為了驗證本文采用鳥體模型的準確性及可靠性,針對文獻 [18]中鳥撞鋁板進行了數(shù)值計算。鋁板尺寸為 410 mm×500 mm×10 mm,4邊固支,鳥體以 91.1 m/s的速度與水平方向成 30°角撞擊鋁板中央。圖2為鳥撞鋁板有限元模型。圖中網(wǎng)格部分為鳥體,底部為鋁板。

        圖2 鳥撞鋁板有限元模型Fig.2 Finite element model of bird impact on Aluminum plate

        鋁板材料采用彈塑性模型,具體參數(shù)如表 1所示。

        表 1 鋁板材料Tab.1 Material properties of Aluminum plate

        鋁板中心處位移時間歷程的數(shù)值結(jié)果與試驗結(jié)果的對比關系如圖 3所示[18],從對比圖可以看出數(shù)值計算結(jié)果與測試值吻合較好。因此,驗證了所采用鳥體模型的有效性。

        2.3 槳葉模型

        2.3.1 槳葉有限元模型

        直升機旋翼槳葉結(jié)構(gòu)比較復雜,圖 4為“海豚”直升機復合材料槳葉結(jié)構(gòu)示意圖[19]。槳葉結(jié)構(gòu)主要由大梁、蒙皮、后段件和接頭等組成,此外還有槳尖罩、后緣調(diào)整片、平衡配重等結(jié)構(gòu)。槳葉外形也比較復雜,為了提高槳葉的氣動性能,槳葉沿展向進行了負扭轉(zhuǎn)處理,槳尖處為了延緩激波出現(xiàn),通常采取后掠或斜削。

        圖3 鋁板中心點位移-時間歷程Fig.3 Displacement-time curve of the center point of the Aluminum plate

        圖4 “海豚”直升機槳葉結(jié)構(gòu)示意圖Fig.4 Schematic diagram of"Dolphin"helicopter blades

        根據(jù)數(shù)值模擬研究的重點,對槳葉模型進行了簡化處理,采用一個等截面矩形槳葉,主要由大梁、蒙皮、填充泡沫3部分組成。通過簡化可以極大減小建模工作量、單元數(shù)量,從而提高計算效率。槳葉翼型為OA209,蒙皮采用四節(jié)點殼單元劃分,大梁和泡沫采用八節(jié)點六面體單元劃分。

        2.3.2 槳葉材料模型

        蒙皮材料為G827/3234碳纖維織物與759/3233阻燃玻璃纖維織物復合材料,其室溫下的材料力學性能參數(shù)來自文獻[20],其中碳纖維布4層,單層厚度為 0.46 mm,鋪層角為+45°/-45°,玻璃布 1層 ,厚度為0.14 mm,鋪層角為 0°。復合材料的失效及破壞模式采用Chang-Chang準則預測。

        大梁材料為鈦合金,考慮應變率對材料性能的影響,采用 Johnson-Cook屈服模型描述材料屈服特性,失效采用最大應變失效準則預測,本構(gòu)參數(shù)如表2所示。

        泡沫材料采用各向同性材料模型,密度 50 kg/m3,體積模量為 8 MPa,二次體積粘性系數(shù)1.6,線形體積粘性系數(shù) 0.1。

        表 2 鈦合金 J-C本構(gòu)參數(shù)Tab.2 J-C constitutiverelation parameters of Titanium alloy

        2.3.3 槳葉邊界條件

        在旋轉(zhuǎn)狀態(tài)中,升力采用壓力形式作用于槳葉下表面。由于槳葉在旋轉(zhuǎn)狀態(tài)過程中,受離心力的影響,槳葉上具有了一定的預應力,因此在數(shù)值模擬中,首先采用預應力分析將應力場施加到槳葉上,然后采用 ALE流固耦合算法進行槳葉鳥撞動響應計算。槳葉轉(zhuǎn)速為 358 r/min。

        在靜態(tài)數(shù)值計算中,槳葉根部固支,在槳葉尖端施加集中載荷(等于槳葉0.9R處的離心力)。槳葉鳥撞有限元計算模型如圖5所示。其中圓柱體網(wǎng)格部分代表鳥體。

        圖5 槳葉鳥撞有限元計算模型Fig.5 FEM modeling of bird impact on rotor blade

        4 數(shù)值計算結(jié)果與分析

        4.1 槳葉鳥撞過程

        在圖 6中,顯示了鳥體在撞擊槳葉被撕碎的變化過程。圖中的長方體為歐拉區(qū)域,圓柱體為鳥體。可以看出,由于相對撞擊速度比較高,鳥體被切分成兩塊,分別沖擊槳葉的上下表面。在撞擊過程中,鳥體發(fā)生了流變,成碎片狀向外飛濺,體現(xiàn)出一定的流體特性。圖7為美國西科斯基公司在 S-92直升機尾槳鳥撞試驗拍攝的照片[10],可以清晰地看出撞擊過程中鳥體的變化,通過對比,計算結(jié)果準確反映了這一過程。

        圖6 鳥撞過程后處理顯示Fig.6 The picture of changes in the bird body during theimpact course

        圖7 明膠鳥撞試驗Fig.7 The picture of gelatin bird impact on tail

        4.2 接觸壓力

        圖8為鳥體撞擊區(qū)域單元的壓力時間曲線,在撞擊初始階段,鳥體局部相對于物體表面滯流,受到壓縮,在該區(qū)域產(chǎn)生了極高的壓力,最大值近 200 MPa,使得后繼的鳥體材料承受了非常高的壓力梯度,于是鳥體材料就會向外徑方向擴張,形成弱的膨脹波,導致了鳥體內(nèi)部壓力逐漸降低,膨脹波在鳥體中經(jīng)過幾次反射,進入了恒定流動階段,這時,鳥體具有相對恒定的壓力場和速度場。隨著鳥體的消耗,壓力場也發(fā)生變化,當鳥體消耗完時,恒定流動也結(jié)束,撞擊表面壓力快速下降至零,整個鳥撞過程結(jié)束。

        圖8 接觸壓力時間歷程曲線Fig.8 Contact pressure time history curve

        Barber和Wilbeck認為鳥撞過程可被描述成一個非恒定的流體動力學過程[2,3]。S.T.Jenq等將鳥撞過程分為初始撞擊、壓力衰減、恒定流動及流動終止4個階段[21]。并給出了初始沖擊壓力、恒定流動壓力、臨界長徑比以及加載持續(xù)時間等計算公式,圖 9為理想鳥體撞擊接觸壓力示意圖。通過對比,數(shù)值計算結(jié)果的壓力曲線趨勢與理論曲線相一致,也是由這 4個階段組成。

        圖9 理想鳥體撞擊接觸壓力示意圖Fig.9 A schematic representation of ideal bird strike induced contact pressure curve composed of four phases

        4.3 兩種狀態(tài)下的數(shù)值結(jié)果對比分析

        撞擊節(jié)點處的位移隨時間變化的曲線如圖 10所示,在 2.5 ms后兩種狀態(tài)的節(jié)點位移趨于一致,最大位移為 17.5 mm。

        圖11給出了撞擊節(jié)點處的撞擊力的時間曲線,在撞擊初始階段,由于鳥體前端相對于槳葉表面滯流,速度突然降為零,從而在該區(qū)域產(chǎn)生了極高的壓力,撞擊力在很短的時間迅速達到峰值,旋轉(zhuǎn)狀態(tài)的最大值為 21.304 k N,靜止狀態(tài)的最大值為 20.715 k N,相對于旋轉(zhuǎn)狀態(tài)下的差異為 2.8%。

        槳葉根部的應變時間歷程如圖 12,13所示,可以看出應變響應也吻合較好。

        圖10 撞擊節(jié)點位移-時間曲線Fig.10 Displacement-timecurves of the impact node

        圖11 撞擊節(jié)點的載荷-時間曲線Fig.11 Load-time curves of theimpact node

        圖12 槳葉根部 X方向主應變-時間曲線Fig.12 X direction principal strain-time curves of blade root

        5 結(jié) 論

        (1)采用鳥撞鋁板試驗驗證了本文采用鳥體本構(gòu)模型的有效性。

        圖13 槳葉根部Y方向主應變-時間曲線Fig.13 Y direction principal strain-timecurves of blade root

        (2)通過后處理顯示的鳥體撞擊槳葉過程與試驗照片的對比,以及接觸壓力曲線的趨勢,說明數(shù)值方法能夠模擬鳥撞過程。

        (3)位移、撞擊力以及應變對比結(jié)果表明旋轉(zhuǎn)與靜態(tài)的撞擊結(jié)果一致性很好,說明采用靜態(tài)當量辦法替代旋轉(zhuǎn)槳葉鳥撞試驗的方法是可行的,從而為直升機旋轉(zhuǎn)部件的碰撞試驗提供了一種理論依據(jù)。

        (4)由于本文采用槳葉撞擊位置處的離心載荷作為集中載荷施加在槳葉的端部來替代整片槳葉實際承受的離心載荷,集中載荷在槳葉內(nèi)是均勻分布的,而離心載荷是隨槳葉半徑變化的。集中載荷與離心載荷在槳葉撞擊部位對槳葉剛度的作用是相同的,但是在其他位置是不同的。因此,靜動結(jié)果具有較好的一致性,但仍然存在一些差異。

        [1] Dolbeer R A,Wright S E,Cleary E C.Bird strikes to civil helicopters in the United States[A].1990-2005.the 8th meeting of Bird Strike Committee-USA/Canada[C].Missouri USA,August 21-24,2006:45—50.

        [2] Barber JP,Taylor H R,Wilbeck J S.Characterization of bird impact on a rigid plate: partⅠ [J].AFFDL-TR-75-5,AD A021142,1975:1—44.

        [3] Barber J P,Wilbeck J S,Taylor H R.Bird impact forces and pressures on rigid and compliant targets[J].AFFDL-TR-77-60,AD A 061313,1978:1—78.

        [4] Meguid S,Mao R,Ng T.FE analysis of geometry effects of an artificial bird striking an aeroengine fan blade[J].International Journal of Impact Engineering,2008,35(6):487—498.

        [5] 張志林,姚衛(wèi)星.飛機風擋鳥撞動響應分析方法研究[J].航空學報,2004,25(6):577— 580.Zhang Zhi-lin,Yao Wei-xing.Research on dynamic analysis of bird impact on aircraf t windshield[J].Acta Aeronautica et Astronautica Sinica,2004,25(6):577—580.

        [6] 朱書華,童明波.鳥體形狀對飛機風擋鳥撞動響應的影響 [J].南京航空航天大學學報,2008,40(4):551—555.Zhu Shuhua,Tong MingBo.Bird shape sensitivity to dynamic response of bird strike on aircraft windshield[J].Journal of Nanjing University of Aeronautics&Astronautics,2008:40(4):551—555.

        [7] 萬小朋,龔倫,趙美英.基于ANSYS/LS-DYNA的飛機機翼前緣抗鳥撞分析[J].西北工業(yè)大學學報,2007,25(2):285—289.Wan Xiao-Peng,Gong Lun,Zhao Mei-Ying.A method for calculating anti-bird impact capability of an aircraft wing using LS-DYNA software[J].Journal of Northwestern Polytechnical University,2007,25(2):285—289.

        [8] 陳偉,關玉璞,高德平.發(fā)動機葉片鳥撞擊瞬態(tài)響應的數(shù)值模擬 [J].航空學報,2003,24(6):531— 533.Chen Wei,Guan Yu-Pu,Gao De-Ping.Numerical simulation of the transient response of blade due to bird impact[J].Acta Aeronautica et Astronautica Sinica,2003,24(6):531—533.

        [9] 關玉璞,張在坤,趙振華,等.粒子分離器渦流葉片鳥撞擊損傷試驗 [J].航空動力學報,2007,22(12):2 094—2 100.Guan Yu-Pu,Zhang Zai-Kun,Zhao Zhen-Hua.Experimental study on bird damage of vortex vane in the particle separator[J].Journal of Aerospace Power,2007,22(12):2 094—2 100.

        [10]Alan Dobyns,Frank federici,Robert Young.Bird strike analysis and test of a spinning S-92 tail rotor[A].American Helicopter Society 57th Annual forum[C].May 9-11,2001.

        [11]Alan Dobyns.Bird strike analysis of S-92 vertical tail cover using DYTRAN[A].American Helicopter Society Affordable Composite Structures Conference[C].October 7,1998.

        [12]王益鋒,王浩文 ,高正.直升機旋翼槳葉的彈性碰撞動力學建模 [J].航空動力學報,2009,24(05):2 046—2 050.Wang Yi-Feng,Wang Hao-Wen,Gao Zheng.Dynamic model of elastic impact of helicopter rotor blades[J]. Journal of Aerospace Power,2009,24(05):2 046—2 050.

        [13]溫海濤 ,關玉璞,高德平.直升機主槳葉的鳥撞有限元數(shù)值模擬 [J].航空動力學報,2009,24(05):1 150—1 158.Wen Hai-Tao,Guan Yu-Pu,Gao De-Ping.Numerical simulations of bird impact on helicopter rotor blades[J]. Journal of Aerospace Power,2009,24(05):1 150—1 158.

        [14]林長亮 ,王浩文,陳仁良,等.采用流固耦合方法的直升機槳葉鳥撞數(shù)值模擬 [J].科學技術與工程,2012,12(1):1—6.Lin Chang-Liang,Wang Hao-Wen,Chen Ren-Liang.Numerical simulation of bird impact on the helicopter blade by fluid-structure coupling method[J].Science Technology and Engineering,2012,12(1):1—6.

        [15]中華人民共和國國家軍用標準.GJB 2462—95飛機玻璃鳥撞試驗方法[S].1996.

        [16]朱書華 ,童明波,王躍全.某型飛機風擋鳥撞試驗與數(shù)值模擬 [J].應用力學學報,2009,26(03):444— 450.Zhu Shu-hua,Tong Ming-Bo,Wang Yue-Quan.Experiment and numerical simulation for bird impact on aircraft windshield[J].Chinese Journal of Applied Mechanics,2009,26(03):444—450.

        [17]張雄,王天舒.計算動力學 [M].北京:清華大學出版社,2007.

        [18]白金澤.基于神經(jīng)網(wǎng)絡方法的鳥撞飛機風擋反問題研究 [D].西安:西北工業(yè)大學,2003.Bai Jin-Ze. Inverse issue study of bird-impact to aiecraft windshield based on neural network method[D].Xi′an:Northwestern Polytechnical University,2003.

        [19]楊乃賓,倪先平.直升機復合材料結(jié)構(gòu)設計 [M].北京:國防工業(yè)出版社,2008.

        [20]陳祥寶.聚合物基復合材料手冊 [M].北京:化學工業(yè)出版社,2004.

        [21]Jenq S T,Hsiao F B,Lin I C.Simulation of a rigid platehit by a cylindrical hemi-spherical tip-ended soft impactor[J].Computational Materials Science,2007,39(1):518—526.

        猜你喜歡
        模型
        一半模型
        一種去中心化的域名服務本地化模型
        適用于BDS-3 PPP的隨機模型
        提煉模型 突破難點
        函數(shù)模型及應用
        p150Glued在帕金森病模型中的表達及分布
        函數(shù)模型及應用
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
        3D打印中的模型分割與打包
        射精情感曰妓女色视频| 狠狠色噜噜狠狠狠8888米奇| 国产农村熟妇videos| 少妇人妻偷人精品一区二区| 一区二区久久不射av| 一区二区三区成人av| 国产午夜在线视频观看| 欧美 国产 综合 欧美 视频| 久久久精品3d动漫一区二区三区 | 一区二区三区激情免费视频| 亚洲精品无码久久久影院相关影片 | 日韩在线精品视频一区| 亚洲国产精品久久人人爱| 欧美成人一级视频| 国产精品一区二区久久毛片| 开心五月天第四色婷婷| 黑人巨大精品欧美一区二区| 欧美在线播放一区二区| 亚洲素人av在线观看| 日本一区二区三区免费精品| 国产麻豆成人精品av| 国产哟交泬泬视频在线播放| 亚洲无av码一区二区三区| 国产成人无码av一区二区在线观看| 秋霞午夜无码鲁丝片午夜精品| 中字亚洲国产精品一区二区| 亚洲国产91高清在线| 男人和女人做爽爽视频| 激情97综合亚洲色婷婷五| 亚洲黄片高清在线观看| 就爱射视频在线视频在线| 日韩精品无码中文字幕电影| 百合av一区二区三区| 中文字幕精品久久一区二区三区| 国产日韩精品suv| 乌克兰少妇xxxx做受6| 99精品国产成人一区二区在线| 久久婷婷综合激情五月| 欧美大成色www永久网站婷| 老汉tv永久视频福利在线观看| 亚洲av一二三四五区在线|