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

基于近似CMS法及模態(tài)測(cè)試的失諧葉盤結(jié)構(gòu)動(dòng)力學(xué)特性研究

2016-10-18 03:20:50亮,李
中國(guó)測(cè)試 2016年6期
關(guān)鍵詞:模態(tài)振動(dòng)結(jié)構(gòu)

張 亮,李 欣

(遼寧工業(yè)大學(xué)機(jī)械工程與自動(dòng)化學(xué)院,遼寧 錦州 121001)

基于近似CMS法及模態(tài)測(cè)試的失諧葉盤結(jié)構(gòu)動(dòng)力學(xué)特性研究

張亮,李欣

(遼寧工業(yè)大學(xué)機(jī)械工程與自動(dòng)化學(xué)院,遼寧 錦州 121001)

為解決工作轉(zhuǎn)速下失諧葉盤結(jié)構(gòu)動(dòng)力學(xué)特性分析計(jì)算量大的問(wèn)題,提出一種近似CMS超單元法。該方法采用先預(yù)應(yīng)力子結(jié)構(gòu)后整體的分析方式,分析葉盤結(jié)構(gòu)在工作轉(zhuǎn)速下的動(dòng)力學(xué)特性。通過(guò)循環(huán)對(duì)稱分析法驗(yàn)證該方法的分析精度。采用模態(tài)測(cè)試及有限元法識(shí)別葉片失諧參數(shù),基于近似CMS超單元方法分析某失諧葉盤結(jié)構(gòu)在工作轉(zhuǎn)速下的動(dòng)力學(xué)特性。結(jié)果表明:諧調(diào)葉盤結(jié)構(gòu)無(wú)量綱動(dòng)頻相對(duì)誤差最大值為3.07%,滿足分析精度要求。與循環(huán)對(duì)稱分析法相比,近似CMS超單元法可用于諧調(diào)或周期失諧葉盤結(jié)構(gòu)在工作轉(zhuǎn)速下的動(dòng)力學(xué)特性分析,并且適用于任意失諧葉盤結(jié)構(gòu)。

失諧葉盤結(jié)構(gòu);CMS超單元法;子結(jié)構(gòu);模態(tài)測(cè)試;動(dòng)力學(xué)特性

0 引 言

葉盤結(jié)構(gòu)作為航空發(fā)動(dòng)機(jī)的重要零部件,其動(dòng)力學(xué)特性直接影響發(fā)動(dòng)機(jī)的工作可靠性及結(jié)構(gòu)完整性。目前主要采用試驗(yàn)測(cè)量、有限元數(shù)值仿真研究葉盤結(jié)構(gòu)的動(dòng)力學(xué)特性。由于葉盤結(jié)構(gòu)較大且葉片形狀復(fù)雜,整體試驗(yàn)測(cè)試時(shí),工作量很大,且有些測(cè)點(diǎn)振動(dòng)很難測(cè)量,因此通常采用有限元法研究葉盤結(jié)構(gòu)等大型復(fù)雜結(jié)構(gòu)的動(dòng)力學(xué)特性,由于分析精度的要求,葉盤模型劃分網(wǎng)格及節(jié)點(diǎn)數(shù)較多,導(dǎo)致其節(jié)點(diǎn)自由度較大,對(duì)于普通計(jì)算機(jī),直接求解葉盤模型動(dòng)力學(xué)方程組是幾乎不可能的。采用大型數(shù)值仿真工作站即使能夠求解,也比較耗時(shí)。近年來(lái)國(guó)內(nèi)外很多學(xué)者已經(jīng)對(duì)復(fù)雜大型結(jié)構(gòu)的動(dòng)力學(xué)問(wèn)題進(jìn)行了廣泛研究,其主要目的是建立高效的降階方法,既有效減縮系統(tǒng)模型的自由度,又能滿足實(shí)際結(jié)構(gòu)分析精度要求。在國(guó)外,Guyan[1]和Irons[2]提出了自由度靜態(tài)縮減法;Paz[3]提出了自由度動(dòng)態(tài)縮減法。由于靜態(tài)縮減法簡(jiǎn)單實(shí)用,一些學(xué)者在此基礎(chǔ)上提出了幾種模態(tài)綜合法,如 Hurty[4]提出了固定界面法、Hou[5]與Goldman[6]分別提出了自由界面法。Wei[7]在模態(tài)綜合法的基礎(chǔ)上提出了模態(tài)綜合超單元法。在國(guó)內(nèi),王建軍等[8]利用循環(huán)對(duì)稱分析法對(duì)錯(cuò)頻周期失諧葉盤系統(tǒng)的振動(dòng)局部化特性進(jìn)行了分析。周傳月等[9-10]利用循環(huán)對(duì)稱分析法對(duì)某燃?xì)廨啓C(jī)諧調(diào)帶冠葉片進(jìn)行了耦合振動(dòng)分析;采用模態(tài)綜合超單元法對(duì)燃?xì)廨啓C(jī)失諧葉盤系統(tǒng)進(jìn)行了靜頻分析。秦飛等[11]采用直接分析法建立失諧葉盤系統(tǒng)整體有限元模型,分析其在工作轉(zhuǎn)速下的動(dòng)頻。王培屹等[12]和臧朝平等[13]各提出了一種失諧葉盤動(dòng)力響應(yīng)預(yù)測(cè)的減縮計(jì)算方法,并通過(guò)實(shí)例驗(yàn)證了方法的分析精度,但上述兩種縮減方法均未考慮轉(zhuǎn)速(預(yù)應(yīng)力)的影響。

綜上,諧調(diào)或周期失諧葉盤結(jié)構(gòu)可采用循環(huán)對(duì)稱法分析其動(dòng)力學(xué)特性,而任意失諧葉盤結(jié)構(gòu)可采用模態(tài)綜合法或直接法分析其動(dòng)力學(xué)特性。由于直接法比較耗時(shí),通常采用模態(tài)綜合法。模態(tài)綜合法要把葉盤結(jié)構(gòu)劃分成若干個(gè)子結(jié)構(gòu),而失諧葉盤結(jié)構(gòu)的動(dòng)頻及動(dòng)態(tài)響應(yīng)分析是有預(yù)應(yīng)力的子結(jié)構(gòu)模態(tài)綜合分析,經(jīng)常采用先整體后子結(jié)構(gòu)的分析方式,但這種分析方式僅適用中小模型,分析失諧葉盤結(jié)構(gòu)等大模型非常耗時(shí)。基于以上不足,提出一種近似的模態(tài)綜合分析方法,采用先預(yù)應(yīng)力子結(jié)構(gòu)后整體的分析方式,基于模態(tài)測(cè)試及有限元法識(shí)別葉片失諧參數(shù),運(yùn)用該近似模態(tài)綜合法分析了某失諧葉盤結(jié)構(gòu)的動(dòng)力學(xué)特性。

1 近似CMS超單元法

葉盤結(jié)構(gòu)被劃分成N個(gè)子結(jié)構(gòu),每個(gè)葉盤基本扇區(qū)為一個(gè)子結(jié)構(gòu),葉盤結(jié)構(gòu)中第i個(gè)子結(jié)構(gòu)在工作轉(zhuǎn)速下的靜力學(xué)(預(yù)應(yīng)力)方程為

式中:[Ki]——第i個(gè)子結(jié)構(gòu)的剛度矩陣;

{Fi}——由工作轉(zhuǎn)速引起的載荷向量;

{ui}——第i個(gè)子結(jié)構(gòu)的位移向量。

將位移向量分成主自由度(界面自由度,用下標(biāo)m表示)和從自由度(非界面自由度,用下標(biāo)s表示),即:

將式(2)代入式(1)得到:

對(duì)于無(wú)阻尼葉盤結(jié)構(gòu),其第i個(gè)子結(jié)構(gòu)的動(dòng)力學(xué)方程為

式中:[Mi]——子結(jié)構(gòu)質(zhì)量矩陣;

將式(2)代入式(6),并考慮結(jié)構(gòu)自由振動(dòng)時(shí),界面以外節(jié)點(diǎn)所受作用力式(6)變?yōu)?/p>

采用自由界面法子結(jié)構(gòu)模態(tài)綜合超單元法分析,即釋放界面主自由度約束界面自由的無(wú)阻尼子結(jié)構(gòu)自由振動(dòng)方程為

基于自由界面法時(shí),第一變換矩陣為

式中:[I]——單位矩陣;

把N個(gè)子結(jié)構(gòu)的動(dòng)力學(xué)方程全部轉(zhuǎn)換成相應(yīng)模態(tài)坐標(biāo)下的動(dòng)力學(xué)方程,完成了N個(gè)超單元的生成。

將葉盤結(jié)構(gòu)中N個(gè)超單元模態(tài)坐標(biāo)下的動(dòng)力學(xué)方程組合起來(lái)為

假設(shè)各超單元界面為剛性連接,即界面節(jié)點(diǎn)自由度耦合,則第i個(gè)超單元和與其相連接的第j個(gè)超單元界面之間滿足位移協(xié)調(diào)條件和力平衡方程,位移協(xié)調(diào)條件為即:

力平衡方程為

式中[α]為第二變換矩陣,其值由式(14)確定。將式(16)代入到式(13)中獲得廣義坐標(biāo){p}下的綜合方程為

式(18)為考慮旋轉(zhuǎn)預(yù)應(yīng)力影響的整個(gè)葉盤結(jié)構(gòu)在廣義坐標(biāo){p}下的自由振動(dòng)方程。通過(guò)求解式(18)獲得整個(gè)葉盤結(jié)構(gòu)的動(dòng)頻及廣義坐標(biāo){p}下的振型,完成超單元的使用。

通過(guò)代入坐標(biāo)變換式(16)和式(10)可以獲得葉盤結(jié)構(gòu)在物理坐標(biāo){u}下的振型,完成超單元的擴(kuò)展。

2 分析精度驗(yàn)證

采用循環(huán)對(duì)稱分析法驗(yàn)證近似CMS超單元法分析精度。由于循環(huán)對(duì)稱分析法只適用于循環(huán)對(duì)稱結(jié)構(gòu)或周期失諧結(jié)構(gòu),因此假設(shè)葉盤結(jié)構(gòu)是諧調(diào)的,即每個(gè)葉片具有相同的材料參數(shù)。分別采用上述兩種方法對(duì)某壓氣機(jī)第一級(jí)葉盤結(jié)構(gòu)進(jìn)行動(dòng)頻分析,葉盤基本扇區(qū)及葉盤結(jié)構(gòu)有限元模型如圖1所示。葉盤工作轉(zhuǎn)速為11 000 r/min,葉片數(shù)為38。葉片材料參數(shù)為:彈性模量 E0=1.135×1011Pa、泊松比 μ0=0.3、密度 ρ0=4 380 kg/m3;輪盤材料參數(shù)為:彈性模量E1=1.15×1011Pa、泊松比μ1=0.3、密度ρ1=4640kg/m3。

圖1 葉盤基本扇區(qū)及葉盤結(jié)構(gòu)有限元模型

定義無(wú)量綱動(dòng)頻為葉盤結(jié)構(gòu)各節(jié)徑(或各階)動(dòng)頻與諧調(diào)葉片靜頻的比值。無(wú)量綱動(dòng)頻相對(duì)誤差為

式中:j——節(jié)徑數(shù);

采用以上兩種方法獲得的諧調(diào)葉盤結(jié)構(gòu)在工作轉(zhuǎn)速下的無(wú)量綱動(dòng)頻及相對(duì)誤差如表1所示,隨節(jié)徑數(shù)增加,無(wú)量綱動(dòng)頻的相對(duì)誤差逐漸減小,其最大值是3.07%,符合精度要求。

表1 近似CMS超單元法分析精度驗(yàn)證

3 失諧葉盤結(jié)構(gòu)動(dòng)力學(xué)特性分析

3.1葉片失諧參數(shù)識(shí)別

葉盤結(jié)構(gòu)通常被認(rèn)為是循環(huán)對(duì)稱或諧調(diào)的,但由于材料不均、加工誤差、使用磨損或抑制顫振,實(shí)際中各葉片存在小量差異,稱為葉片失諧[14]。假設(shè)輪盤是諧調(diào)結(jié)構(gòu),對(duì)葉片引入不同的彈性模量擾動(dòng)系數(shù)Pj模擬葉片失諧,如下式:

式中:E0——諧調(diào)葉片彈性模量;

Ej——第j個(gè)失諧葉片的彈性模量;

Pj——第j個(gè)失諧葉片彈性模量的擾動(dòng)系數(shù),

Pj可與動(dòng)力學(xué)方程矩陣關(guān)聯(lián)[15]。

采用上述近似CMS超單元法時(shí),式(20)表示把每個(gè)葉片的整體剛度作為子矩陣,并引入一個(gè)失諧彈性模量擾動(dòng)參數(shù)Pj[16]。本文采用基于葉片模態(tài)測(cè)試及有限元法相結(jié)合的葉片失諧參數(shù)識(shí)別方法,對(duì)某壓氣機(jī)一級(jí)葉盤結(jié)構(gòu)各葉片失諧參數(shù)進(jìn)行識(shí)別,具體識(shí)別流程及模態(tài)測(cè)試系統(tǒng)參數(shù)詳見(jiàn)文獻(xiàn)[15]。表2給出了采用上述方法識(shí)別的各葉片失諧彈性模量。

表2 失諧彈性模量

3.2葉盤結(jié)構(gòu)動(dòng)力學(xué)特性分析

采用近似CMS超單元法對(duì)表2中給出的某失諧葉盤結(jié)構(gòu)進(jìn)行動(dòng)頻分析。圖2為諧調(diào)葉盤結(jié)構(gòu)的振型圖。由于葉盤結(jié)構(gòu)是諧調(diào)的,因此其振動(dòng)形式為節(jié)徑振動(dòng)。圖2(a)為0節(jié)徑振動(dòng)即節(jié)圓振動(dòng),即存在一個(gè)圓心為系統(tǒng)中心的圓線,其上的振幅為0。圖2(b)為3節(jié)徑振動(dòng),即存在3條通過(guò)系統(tǒng)中心的徑向線,其上的振幅為0。

圖2 諧調(diào)葉盤結(jié)構(gòu)節(jié)徑振型圖

圖3為失諧葉盤結(jié)構(gòu)振型圖。失諧葉盤結(jié)構(gòu)第6、7階振型是與諧調(diào)葉盤結(jié)構(gòu)的3節(jié)徑振型相對(duì)應(yīng)的,由于葉片存在失諧,第6、7階振型不再呈現(xiàn)3節(jié)徑振動(dòng),振動(dòng)能量集中在少數(shù)葉片上,出現(xiàn)嚴(yán)重的振動(dòng)局部化現(xiàn)象。

圖4為諧調(diào)與該失諧葉盤結(jié)構(gòu)1彎族各階無(wú)量綱動(dòng)頻。從圖中可看出,諧調(diào)葉盤系統(tǒng)出現(xiàn)重頻現(xiàn)象,即除第1階(0節(jié)徑)及第38階(19節(jié)徑)外,頻率成對(duì)出現(xiàn);而失諧葉盤結(jié)構(gòu)由于失諧的存在,出現(xiàn)了頻率分離的現(xiàn)象,即原來(lái)成對(duì)出現(xiàn)的重頻不再相等;由于各葉片失諧彈性模量均小于諧調(diào)葉片的彈性模量,因此該失諧葉盤結(jié)構(gòu)各階無(wú)量綱動(dòng)頻均小于諧調(diào)葉盤結(jié)構(gòu)。

圖4 諧調(diào)與失諧葉盤結(jié)構(gòu)1彎族各階無(wú)量綱動(dòng)頻

對(duì)比近似CMS超單元法與循環(huán)對(duì)稱分析法,循環(huán)對(duì)稱分析法只適用于循環(huán)對(duì)稱(周期對(duì)稱)結(jié)構(gòu),即諧調(diào)或周期失諧葉盤結(jié)構(gòu)。而本文方法可對(duì)每個(gè)葉片引入相同或不同的材料參數(shù),如葉片彈性模量的失諧,既可分析諧調(diào)、周期失諧葉盤結(jié)構(gòu)的動(dòng)力學(xué)特性,又可分析任意失諧葉盤結(jié)構(gòu)的動(dòng)力學(xué)特性。

4 結(jié)束語(yǔ)

本文針對(duì)現(xiàn)有方法分析有預(yù)應(yīng)力失諧葉盤結(jié)構(gòu)等復(fù)雜大型模型時(shí)存在的不足,提出了一種近似CMS超單元法。通過(guò)循環(huán)對(duì)稱分析法驗(yàn)證了該方法的分析精度;采用模態(tài)測(cè)試及有限元法識(shí)別了葉片失諧參數(shù);基于該近似CMS超單元方法分析了某失諧葉盤結(jié)構(gòu)在工作轉(zhuǎn)速下的動(dòng)力學(xué)特性。結(jié)果表明:諧調(diào)葉盤結(jié)構(gòu)無(wú)量綱動(dòng)頻相對(duì)誤差最大值為3.07%,符合精度要求。與循環(huán)對(duì)稱法相比,本文方法不僅可分析諧調(diào)、周期失諧葉盤結(jié)構(gòu)在工作轉(zhuǎn)速下的動(dòng)力學(xué)特性,還可分析任意失諧葉盤結(jié)構(gòu)在工作轉(zhuǎn)速下的動(dòng)力學(xué)特性。

[1]GUYAN R J.Reduction of stiffness and mass matrices[J]. AIAA J,1965,3(2):380.

[2]IRONS B M.Structural eigenvalue problems:elimination of unwanted variables[J].AIAA J,1965,3(5):961-962.

[3]PAZ M.Dynamic condensation[J].AIAA J,1984,22(5):724-727.

[4]HURTY W C.Vibration of structural system by component mode synthesis[J].Journal of the Engineering Mechanics Division ASCE,1960(86):51-69.

[5]HOUS.Reviewof mode synthesis techniques and a newapproach[J].Shock and Vibration Bulletin,1969,40(4):25-30.

[6]GOLDMANRL.Vibrationsanalysisbydynamic partitioning[J].AIAA J,1969,7(6):1152-1154.

[7]WEI S T,PIERRE C.Statistical analysis of the forced response of mistuned cyclic assemblies[J].AIAA J,1990,28(5):861-868.

[8]王建軍,于長(zhǎng)波,李其漢.錯(cuò)頻葉盤結(jié)構(gòu)振動(dòng)模態(tài)局部化特性分析[J].航空動(dòng)力學(xué)報(bào),2009,24(4):788-792.

[9]周傳月,鄒經(jīng)湘,聞雪友,等.燃?xì)廨啓C(jī)帶冠葉片耦合振動(dòng)分析[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào),2001,33(1):129-133.

[10]周傳月,鄒經(jīng)湘,聞雪友,等.燃?xì)廨啓C(jī)失調(diào)葉盤系統(tǒng)的振動(dòng)特性分析[J].燃?xì)廨啓C(jī)技術(shù),2000,13(3):42-46.

[11]秦飛,陳立明.失調(diào)葉片-輪盤系統(tǒng)耦合振動(dòng)分析[J].北京工業(yè)大學(xué)學(xué)報(bào),2007,33(2):126-128.

[12]王培屹,李琳.用于失諧葉盤動(dòng)力學(xué)特性分析的減縮計(jì)算方法[J].航空動(dòng)力學(xué)報(bào),2014,29(6):1395-1402.

[13]臧朝平,段勇亮,PETROV E P.失諧葉片輪盤的減縮建模及動(dòng)力響應(yīng)預(yù)測(cè)方法[J].航空學(xué)報(bào),2015,36(10):3305-3315.

[14]張亮,袁惠群,韓清凱,等.基于微動(dòng)滑移摩擦模型的失諧葉盤系統(tǒng)振動(dòng)分析[J].振動(dòng)工程學(xué)報(bào),2012,25(3):289-293.

[15]張亮,李欣,袁惠群.基于模態(tài)測(cè)試及有限元法的葉片失諧參數(shù)識(shí)別[J].中國(guó)測(cè)試,2015,41(11):16-19.

[16]張亮.航空發(fā)動(dòng)機(jī)葉片輪盤系統(tǒng)振動(dòng)特性及多長(zhǎng)耦合力學(xué)特性研究[D].沈陽(yáng):東北大學(xué),2013.

(編輯:劉楊)

Study on structural dynamics of mistuned bladed disk system based on approximate CMS method and modal testing

ZHANG Liang,LI Xin
(Faculty of Mechanical Engineering and Automation,Liaoning University of Technology,Jinzhou 121001,China)

In order to save the great computation effort for the structural dynamic characteristics analysis of mistuned bladed disk at working speed,an approximate CMS super-element method is proposed.The analysis mode from the prestressed substructure to the whole model is applied to analyze the dynamic characteristics of bladed disk system at working speed.Analysis accuracy of this method is verified by the cyclic symmetry analysis method.The mistuning parameter of blades is identified by modal testing and finite element method,and dynamic characteristics of mistuned bladed disk system at working speed are analyzed based on the approximate CMS super-element method.The results showed that the maximum relative error of dimensionless dynamic frequency of the tuned system is 3.07%,which is in conformity with the requirements of accuracy.Compared with the cyclic symmetry analysis method,this method not only can be used for the analysis of dynamic characteristics of tuned and cyclic mistuned bladed disk system at working speed,but also can be applied to the arbitrary mistuned bladed disk system.

mistunedbladeddisksystem;CMSsuper-elementmethod;substructure;modal testing;dynamic characteristics

A

1674-5124(2016)06-0117-05

10.11857/j.issn.1674-5124.2016.06.025

2016-01-27;

2016-03-03

國(guó)家自然科學(xué)基金項(xiàng)目(51505206);遼寧省教育廳科學(xué)研究一般項(xiàng)目(L2014246);遼寧工業(yè)大學(xué)教師科研啟動(dòng)基金(X201202)

張亮(1983-),男,遼寧葫蘆島市人,副教授,博士,研究方向?yàn)闄C(jī)械振動(dòng)、轉(zhuǎn)子動(dòng)力學(xué)。

猜你喜歡
模態(tài)振動(dòng)結(jié)構(gòu)
振動(dòng)的思考
《形而上學(xué)》△卷的結(jié)構(gòu)和位置
振動(dòng)與頻率
論結(jié)構(gòu)
中立型Emden-Fowler微分方程的振動(dòng)性
論《日出》的結(jié)構(gòu)
國(guó)內(nèi)多模態(tài)教學(xué)研究回顧與展望
創(chuàng)新治理結(jié)構(gòu)促進(jìn)中小企業(yè)持續(xù)成長(zhǎng)
基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
UF6振動(dòng)激發(fā)態(tài)分子的振動(dòng)-振動(dòng)馳豫
主站蜘蛛池模板: 亚洲天堂日本| 青草视频在线观看国产| 暴力调教一区二区三区| 亚洲欧美激情小说另类| 成人综合久久综合| 日韩无码真实干出血视频| 精品91视频| 亚洲日本一本dvd高清| 老司机久久99久久精品播放| 亚洲日韩国产精品无码专区| 国产美女无遮挡免费视频网站 | 欧美日韩精品一区二区在线线| 亚洲一区免费看| 久久免费视频播放| 97久久免费视频| 日韩精品成人网页视频在线| 精品久久久久久久久久久| 97久久超碰极品视觉盛宴| 波多野结衣一区二区三区AV| 欧洲欧美人成免费全部视频| 97人人做人人爽香蕉精品| 亚洲日韩精品无码专区97| 欧美在线天堂| 午夜视频免费试看| 九色视频一区| 亚洲性色永久网址| 91久久精品国产| 亚洲一区二区视频在线观看| 99热这里只有成人精品国产| 亚洲一级毛片| 亚洲精品777| 欧美a在线| 久久6免费视频| www.狠狠| 日韩少妇激情一区二区| 免费观看精品视频999| 日本成人一区| 亚洲AV无码乱码在线观看代蜜桃 | 中文字幕无码av专区久久| 午夜视频在线观看区二区| 亚洲伊人天堂| 久久久精品国产SM调教网站| 久久香蕉欧美精品| 欧美成人看片一区二区三区| 伊人久久青草青青综合| 欧美色伊人| 国产激情无码一区二区免费| 免费女人18毛片a级毛片视频| 久久永久视频| 国产专区综合另类日韩一区| 成人在线观看一区| 色综合国产| 无码AV高清毛片中国一级毛片 | 日韩黄色精品| 亚洲男人天堂久久| 欧美一级片在线| 精品小视频在线观看| 国产精品久久自在自线观看| 国产青榴视频| 免费人成视网站在线不卡| 国产精品七七在线播放| 久久午夜夜伦鲁鲁片无码免费| 99精品国产高清一区二区| 真实国产乱子伦高清| 婷婷综合缴情亚洲五月伊| 思思热在线视频精品| 国产视频久久久久| 玖玖免费视频在线观看| 日日碰狠狠添天天爽| 99草精品视频| 狠狠亚洲婷婷综合色香| 国产青青操| 亚洲无码一区在线观看| 欧美福利在线观看| 国产91丝袜在线播放动漫 | 久久国产V一级毛多内射| 福利姬国产精品一区在线| 日韩精品成人在线| 欧美区在线播放| 国产精品国产三级国产专业不 | 免费国产好深啊好涨好硬视频| 国产一级视频在线观看网站|