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

基于隨機共振與蝙蝠算法的高速動車組滾動軸承故障診斷

2021-04-20 08:14:20趙鵬振
中國測試 2021年3期
關鍵詞:故障信號

趙鵬振,劉 繼

(同濟大學 鐵道與城市軌道交通研究院 ,上海 201804)

0 引 言

滾動軸承是高速動車組中至關重要的零件,對列車的安全穩(wěn)定運行有著重要的作用。目前,對于滾動軸承的單體研究與實驗數(shù)據(jù)處理已經(jīng)較為成熟,但在高速動車組軸承早期故障的在途診斷方面研究還不完善。

國內(nèi)對于高速動車組的滾動軸承故障研究方面,大多沿用了通用軸承的研究方法,即建立動力學模型求解數(shù)值解。曹青松等[1]建立了CRH1型動車組滾動軸承-車軸耦合系統(tǒng)的非線性動力學模型,采用數(shù)值方法分析了不同工況下動車組軸承-車軸耦合系統(tǒng)的動力學響應與非線性特性。關貞珍等[2]建立了轉(zhuǎn)子-故障滾動軸承-軸承座的非線性振動模型以及軸承外圈、內(nèi)圈、滾動體局部故障非線性動力學模型,運用數(shù)值積分進行仿真分析。東亞斌等[3]對單一軸承外圈故障和滾動體故障建立了動力學模型,并用實際軸承進行了驗證。李佳睿[4]對列車行駛時的滾動軸承故障信號進行了模擬,在模擬故障信號上添加白噪聲后提出相應算法提取原信號。梁瑜[5]對比了目前在噪聲中提取故障軸承加速度信號的理論方法,并提出了在EMD方法上改進的AFD算法。

針對高速動車組滾動軸承故障的研究,學者們的研究大多是對滾動軸承單體,或滾動軸承與少數(shù)部件如車軸、軸承座等的結(jié)合體進行分析;對于從整車角度分析,將故障軸承放在運行的高速動車組中,通過列車運行參數(shù)的變化來診斷滾動軸承故障的研究較少。且在滾動軸承故障信號在復雜背景噪聲中的提取方面,目前研究方法大多是人為施加白噪聲,并未將系統(tǒng)噪聲、軌道不平順等作為噪聲添加到模擬信號中,從而造成了與實際情況中的信號背景噪聲產(chǎn)生差異。

針對這些不足,本文通過對滾動軸承的運動學與動力學分析,研究了滾動軸承外圈故障造成沖擊激勵的頻率與幅值,并據(jù)此模擬軸承故障對列車車體造成的沖擊信號。依據(jù)CRH3型高速動車組動車實際參數(shù),在SIMPACK平臺建立整車模型,模型包含1個車體、2個轉(zhuǎn)向架、4副輪對、8個轉(zhuǎn)臂、4個電機、4個齒輪箱,如圖1所示。

圖1 CRH3型高速動車組動車模型

以中國高速鐵路無砟軌道不平順譜作為軌道激勵[6],驗證模型動力學性能。發(fā)現(xiàn)模型非線性臨界速度為485 km/h;在300 km/h的運行速度下,垂向Sperling指標為1.804 9和1.826 5均符合標準。將沖擊信號輸入到整車模型的軸承位置,使列車在具有軌道不平順的軌道上運行,提取軸箱處的垂向加速度信號。提出基于隨機共振的信號處理方法,增強復雜背景噪聲下有效故障成分占比,并用蝙蝠算法實現(xiàn)對隨機共振的自適應調(diào)參,比較算法處理前后的信噪比,證明算法有效性。并提出一種滾動軸承故障的判斷指標L以及基于L的高速動車組軸承早期故障的在途檢測方法。

1 軸承故障沖擊

工程中,滾動軸承的缺陷往往是由多方面的因素共同造成的,例如摩擦損傷、交變載荷、電化學腐蝕、接觸疲勞等。故障的類型有很多種,包括腐蝕、磨損、疲勞剝落、膠合、塑性變形、斷裂等。這些故障根據(jù)缺陷尺寸與分布規(guī)律可以分為局部式、擴展式和分布式3種。其中,早期以裂紋、凹坑等局部式缺陷為主,中后期滾動體通過缺陷部位時的振動沖擊使缺陷尺寸沿滾道方向逐漸增加[7],進而形成擴展缺陷。分布式缺陷通常由研磨磨損、制造誤差、操作不當?shù)纫穑植加谡麄€滾道,對軸承整體振動響應都有影響。

本研究聚焦于高速動車組軸箱軸承早期故障的特征信號提取與故障診斷。我國高速動車組軸箱軸承均采用雙列圓錐滾子軸承。因此,本文以雙列圓錐滾子軸承外圈滾道的早期局部式缺陷為對象展開研究。提出的信號提取算法屬于從復雜背景噪聲中提取小幅規(guī)律故障信號的方法,可同樣應用于擴展式缺陷等中后期軸承故障。

如圖2所示,對滾動軸承結(jié)構(gòu)進行分析,假設滾子均勻地分布在保持架中且不相對于滾道移動。設軸承外圈的旋轉(zhuǎn)頻率為fo,內(nèi)圈的旋轉(zhuǎn)頻率為fi,滾子公轉(zhuǎn)速度為vC,軸承節(jié)圓直徑為D,滾子直徑為d,接觸角為 α,滾子個數(shù)為Z。

圖2 滾動軸承結(jié)構(gòu)圖

設某一滾動體中心點為C,與內(nèi)圈滾道的接觸點為A,與外圈滾道的接觸點為B。則A、B點繞軸承中心轉(zhuǎn)動直徑分別為D?dcosα,D+dcosα。

軸承運轉(zhuǎn)時,滾子公轉(zhuǎn)速度為:

則單個滾子的轉(zhuǎn)動頻率為:

滾子在轉(zhuǎn)過外圈滾道時的頻率為:

由于軸承有Z個滾子,且在高速動車組滾動軸承運行過程中,軸承外圈是相對固定的,即fo=0。

則圓錐滾子軸承外圈故障特征頻率為:

已知高速動車組直線運行速度為350 km/h,輪對滾動圓半徑為0.46 m,CRH3型高速動車組采用FAG軸承,主要型號為TAR01.130-240-B-TVP,軸承參數(shù)如表1所示。不考慮蠕滑影響,可計算得到輪對滾動的頻率為33.64 Hz,由式(4)得,外圈故障頻率為246.5 Hz,角速度為1 548.8 rad/s。

表1 圓錐滾子軸承TAR01.130-240-B-TVP參數(shù)

在李治強[7]研究中建立的軸承模型基礎上,本文以雙列圓錐滾子軸承為研究對象,建立了如圖3所示的軸承坐標系統(tǒng)。模型將軸承內(nèi)外圈考慮為剛體,將滾動體與滾道的接觸考慮為彈簧-阻尼系統(tǒng)。考慮滾動體與滾道之間的Hertz接觸變形。本文聚焦于軸承故障振動響應的研究,為了簡化分析,忽略對軸承振動響應特征影響較小的次要因素,模型簡化條件如下:

圖3 滾動軸承動力學模型

1)軸承運行時環(huán)境條件不變;

2)不考慮滾動體的質(zhì)量;

3)軸承內(nèi)圈與軸為剛性連接,外圈與軸承座通過彈簧-阻尼系統(tǒng)連接;

4)假設在保持架作用下,滾動體之間的角間距保持不變,忽略滾動體的慣性及離心力影響;

5)忽略潤滑油膜及位移點振動。

首先,對整個雙列圓錐滾子軸承分析,定義滾動軸承承受的力和轉(zhuǎn)矩為矢量F:

式中:Fx、Fy、Fz——水平徑向x方向、軸向y方向、垂直徑向z方向的力,N;

Mx、My、Mz——繞x、y、z方向的轉(zhuǎn)矩,N·m。

將滾動軸承的內(nèi)外圈考慮為剛體,以內(nèi)外圈之間的相對運動定義軸承系統(tǒng)的接觸變形,內(nèi)外圈之間的相對轉(zhuǎn)動和位移矢量q為:

式中: δx、δy、δz——內(nèi)外圈在水平徑向x、軸向y、垂直徑向z方向的相對位移,μm;

θx、θy、θz——內(nèi)外圈繞x、y、z方向的相對轉(zhuǎn)動,(°)。

當滾動軸承內(nèi)圈和軸以角速度ws旋轉(zhuǎn)時,保持架的旋轉(zhuǎn)角速度為:

式中:rp——滾動軸承分度圓直徑;

rd——滾子半徑;

α——滾動軸承接觸角。

則滾動軸承中,第q列第j個滾子角位置 θjq為:

其中N是滾子個數(shù)。

第q列第j個滾子在徑向的位移 δrjq為:

其中c0是滾動軸承徑向間隙。

當出現(xiàn)軸承外圈故障時,根據(jù)劉永強等提出的故障模型[8],如圖4所示。

圖4 外圈故障示意圖

設缺陷處于外圈滾道底部區(qū)域,缺陷對應的圓心角為 θcr,損傷寬度為L0,圓錐滾子質(zhì)心處的截面半徑為rg,滾子通過缺陷時,滾子在缺陷處的接觸變化量為:

則變化后的接觸變形量為:

由赫茲接觸理論,單個滾動體與滾道間接觸力為:

其中K為滾動體與滾道間的有效接觸剛度系數(shù)。當滾動軸承系統(tǒng)承受水平徑向載荷Fx和垂直徑向載荷Fz時,雙列圓錐滾子軸承動力學方程如下:

式中:c——接觸阻尼;

Fix、Fiz、Fox、Foz——內(nèi)外圈在x方向和z方向的負載。

定義公式(13)中的狀態(tài)矢量w(t)為:

式中:xi(t)、zi(t)、xo(t)、zo(t)——內(nèi)外圈在x和z方向的位移,是時間t的函數(shù)。

δx(t)、 δz(t)——軸承內(nèi)外圈相對位移在x、z方向的分量,定義為:

公式(13)中4容量矩陣Rjq(t)為正交坐標系向徑向坐標系轉(zhuǎn)換的系數(shù),定義式為:

公式(13)中質(zhì)量矩陣M,剛度矩陣K,阻尼矩陣C分別為:

式中:mi——軸承內(nèi)圈與軸的質(zhì)量,kg;

mo——軸承外圈與軸承座的質(zhì)量,kg;

kox、koz——軸承外圈剛度系數(shù),MN/m;

cox、coz——軸承外圈阻尼系數(shù),N·s/m。

以TAR01.130-240-B-TVP型雙列圓錐滾子軸承為研究對象,進行數(shù)值計算。與實際情況結(jié)合,令滾動軸承系統(tǒng)承受z向靜載荷73.5 kN,即二分之一軸重。滾子與內(nèi)外滾道的接觸剛度系數(shù)K=14.5 GN/m1.5,阻尼系數(shù)c=800 N·s/m,mi=1 627 kg,mo=84 kg,其余參數(shù)如表1所示。

滾動軸承中滾子與滾道間接觸力的突然變化是產(chǎn)生沖擊的原因。經(jīng)計算可知,當雙列圓錐滾子軸承單列外圈滾道發(fā)生故障,故障深度為50 μm,缺陷圓周角為15°時,由外圈振動產(chǎn)生的沖擊信號如下如圖5所示。

圖5 軸承故障沖擊信號模擬

2 信號增強算法

2.1 隨機共振算法

隨機共振是基于雙穩(wěn)態(tài)模型通過噪聲增強微弱周期沖擊成分的一種算法,可以用非線性朗之萬方程表示為:

式中:x(t)——系統(tǒng)的輸出信號;

s(t)——系統(tǒng)的輸入信號;

ξ(t)——均值為0方差為1的高斯白噪聲;

D——噪聲強度;

U(x)是雙穩(wěn)態(tài)系統(tǒng)的勢函數(shù),表達式為:

U(x)有兩個穩(wěn)態(tài)點和一個非穩(wěn)態(tài)點x=0,勢壘高度為U=a2/4b。 系統(tǒng)輸出x(t)是周期性輸入信號s(t)和噪聲n(t)共同作用下布朗粒子的運動軌跡。當s(t)和n(t)同時存在時,粒子可以積累一定的能量完成在兩個勢阱間的躍遷,從而產(chǎn)生于輸入信號頻率相同的振動,即隨機共振。此時部分噪聲能量轉(zhuǎn)換成為信號能量,提高系統(tǒng)輸出信噪比,進而實現(xiàn)在噪聲中提取微弱信號的功能。

隨機共振系統(tǒng)通常會受到小參數(shù)的限制,即系統(tǒng)輸入信號s(t)的頻率、幅值等遠小于1。在高速動車組軸承故障檢測中,故障特征頻率較高,無法利用小參數(shù)隨機共振來實現(xiàn)。結(jié)合文獻[9-11]提出的理論與處理方法,本文采用在大參數(shù)信號處理中常用的變尺度信號處理方法,壓縮軸箱垂向加速度信號,使其滿足隨機共振小參數(shù)的要求。尺度變換的方法如下:

設定頻率壓縮比K,根據(jù)K重新定義一個壓縮后的采樣頻率fks=fs/K,而后設置步長h=1/fks。

本文中,使用龍格庫塔法來實現(xiàn)隨機共振的計算,同時采樣頻率采用5 000 Hz,取K=50,則步長為h= 1/100。

2.2 蝙蝠算法

蝙蝠算法(bat algorithm,BA)[12]是 Yang在2010年提出的一種優(yōu)化算法。Ali,Rodrigues D等的研究結(jié)果[13-14]表明:蝙蝠算法具有參數(shù)少,效率高的特點,對于高速動車組軸承故障診斷這種大計算量的工程項目具有較高的性價比與實用價值。與其他最優(yōu)算法相比,在蝙蝠算法中,當某一解超出了約束條件,并非簡單的將其舍棄,而是在其基礎上進行下一步處理,從而加速全局搜索到局部搜索的進程,極大地提升了算法的收斂效率。

蝙蝠i發(fā)出的聲波頻率fi,自身的速度和位置的更新公式如下式:

式中:β——[0,1]之間的一個隨機數(shù);

x——使目標函數(shù)最優(yōu)時蝙蝠所處的位置;

fmin和fmax——聲波頻率的最小值和最大值。

位置x、脈沖響度Ai和脈沖發(fā)射率Ri更新公式為:

式中:ε——[-1,1]中的隨機數(shù);

At——所有蝙蝠在t時刻的平均響度;

β——脈沖響度增加系數(shù);

γ——脈沖頻度衰減系數(shù)。

本文中,使用BA算法對隨機共振過程的系數(shù)a、b進行自適應調(diào)節(jié),從而最大限度地提取故障信號中的有效成分。

2.3 故障信號采集

將軸承外圈故障沖擊信號輸入SIMPACK模型中第一輪對右側(cè)軸承位置,模擬列車在具有軌道不平順的直線軌道上(采用中國高速鐵路無砟軌道不平順譜)以350 km/h的速度直線運行,進行仿真計算。以5 000 Hz的采樣頻率提取第一輪對上的右側(cè)軸箱處垂向加速度信號,采樣長度為8 s,作為故障信號,進行算法處理及故障診斷。

本文對正常工況與3種不同程度的軸承外圈故障工況進行仿真計算,如表2所示。并基于故障信號驗證信號處理方法的有效性。

表2 4種工況

2.4 信號處理算法

為了在復雜噪聲中判斷列車軸承是否發(fā)生故障,需要對原始的軸箱垂向加速度信號進行處理,放大其中的有效成分。本文采取的信號增強算法分為三步:1)對原始信號進行高通濾波;2)進行希爾伯特變換并取其解析信號;3)對解析信號進行自適應隨機共振處理。其中,第二步對信號進行希爾伯特變換并取其解析信號是為了使信號正頻段幅值增大,平穩(wěn)性增強,便于第三步的隨機共振處理。第三步中運用蝙蝠算法對隨機共振處理進行自適應調(diào)參,目標函數(shù)設置為信噪比,使蝙蝠參數(shù)的更新向使目標函數(shù)最小的方向發(fā)展。流程圖如圖6所示。

圖6 信號處理流程圖

2.5 信號增強效果

為驗證信號增強算法的有效性,在故障工況下,比較原始信號與信號增強算法處理后的輸出信號中,有效成分的信噪比。

在故障工況2下(外圈故障缺陷圓周角為15°),仿真采集的原始故障信號如圖7所示。

圖7 故障工況2下原始故障信號

此信號為復雜背景噪聲下的故障信號。仿真采集相同故障狀態(tài)下,列車在無軌道不平順的光滑軌道上運行時的軸箱處垂向加速度信號,作為無背景噪聲的有效故障信號。經(jīng)計算可知,原始故障信號的信噪比為-19.56 dB。

經(jīng)信號增強算法處理后,輸出結(jié)果如圖8所示。

圖8 處理后輸出結(jié)果

基于相同的無背景噪聲的有效故障信號,計算處理后信號的信噪比為-0.06 dB。比較可看出,本文提出的信號處理方法可以顯著增強復雜背景噪聲下的有效成分占比,使之更易被觀測。

3 故障診斷指標L

3.1 指標L計算方法

若列車軸承存在外圈故障,則在軸箱垂向加速度信號中,外圈故障特征頻率fout的整數(shù)倍成分應比周圍成分功率高。若要判斷列車滾動軸承是否存在外圈故障,需要判斷輸出信號中,在外圈故障特征頻率fout整數(shù)倍的成分的功率是否突出。通常,判斷信號特定頻率成分功率是否突出的指標是信噪比,信噪比是指有用信號的功率與噪聲信號的功率的比。用于本文情況,即將信號特定頻率成分視為有效信號,其余成分視為噪聲信號,進行信噪比判斷。

在本文中,由于經(jīng)自適應隨機共振算法處理后,輸出信號中各成分的功率并不確定,所以若采用比較信噪比的方式,需加入濾波器提取出有效信號成分,再用輸出信號減去有效信號得到噪聲信號,再計算功率比。這種方法會增加誤差并增加計算量。所以,本文將信噪比概念與實際情況結(jié)合,提出了一種基于輸出信號頻率密度的故障指標L。計算流程如圖9所示。

圖9 故障診斷指標計算流程

3.2 指標L有效性分析

在正常工況下,仿真采集的原始故障信號在200~300 Hz區(qū)域的頻譜圖像如圖10所示,此時故障指標L為1.471 8。

圖10 正常工況200~300 Hz區(qū)間頻域原始信號

信號增強算法處理后,輸出結(jié)果如圖11所示。可以看到,在外圈故障特征頻率(246.5 Hz)附近也無明顯峰值。此時故障指標L為1.837 8。

圖11 正常工況200~300 Hz區(qū)間處理后結(jié)果

在故障工況2下(外圈故障缺陷圓周角為15°),仿真采集的原始故障信號在200~300 Hz區(qū)域的頻譜圖像如圖12所示,此時故障指標L為2.549 8。

圖12 故障工況2下200~300 Hz區(qū)間頻域原始信號

通過信號增強算法對原始故障信號進行處理,輸出信號在200~300 Hz區(qū)域的頻譜圖像如圖13所示,此時故障指標L為5.995 5。

比較可看出,故障工況下,采集的故障信號在外圈故障特征頻率(246.5 Hz)附近峰值明顯。經(jīng)信號增強算法處理后,峰值更加突出,故障工況的指標L值由2.549 8增高至5.995 5,遠大于正常工況下的1.837 8。

圖13 故障工況2下200~300 Hz區(qū)間處理后結(jié)果

分別對4種工況下原始故障信號、經(jīng)信號增強算法處理后的輸出信號的故障診斷指標L進行計算,計算結(jié)果如表3所示。

表3 4種工況下處理前后的故障診斷指標L

可以看出,3種故障工況下,原始故障信號的故障指標L值均大于正常工況,且經(jīng)信號增強算法處理后,故障指標L值均有較大提升,與正常工況間的差別明顯。正常工況下故障指標L的值較小,算法處理前后均不超過2,而在故障工況下,經(jīng)信號增強算法處理后信號的故障指標L值均超過2.5。因此,故障指標L能顯著區(qū)分正常工況與故障工況。與信噪比等傳統(tǒng)指標相比,故障指標L既降低了運算量,又可以直觀地顯示軸承運行狀態(tài),具有較強的實用價值。

結(jié)合信號處理算法與故障診斷指標L,可以提出一種判別動車組滾動軸承外圈故障的新方法。即根據(jù)指標L設定故障報警閾值Lo,如對本文研究的350 km/h運行的CHR3動車組,可令Lo=2.3。當列車運行過程中,提取軸箱處垂向加速度信號,經(jīng)信號增強算法處理后,計算故障指標L值,若L值持續(xù)大于閾值Lo,且逐漸加大時,則認為滾動軸承發(fā)生外圈故障,應立即報警并采取安全措施。

4 結(jié)束語

本文首先研究了滾軸承外圈故障引起的沖擊激勵,通過分析滾動軸承運動特點計算沖擊頻率,并建立滾動軸承動力學模型,計算外圈故障引起的沖擊信號。在SIMPACK平臺建立CRH3高速動車組動車整車模型,將滾動軸承外圈故障沖擊激勵添加到故障位置,從而模擬出故障列車在具有軌道不平順的軌道上直線行駛時的運行工況。

提出了基于蝙蝠算法的自適應隨機共振算法,對軸箱處垂向加速度信號進行處理,增強故障信號中的有效成分占比。通過對比故障工況下,信號處理前后故障信號的信噪比,驗證了信號增強算法對于在復雜背景噪聲中提取滾動軸承外圈故障信號有明顯效果。

提出動車組滾動軸承外圈故障診斷指標L。基于處理后的加速度信號與故障診斷指標L,提出了一種新的高速動車組滾動軸承外圈故障判別方法。通過對比多工況下的故障診斷指標L值,說明其能準確反映高速動車組軸承的外圈故障狀況,明顯區(qū)分正常工況與故障工況。從而驗證了故障診斷指標L在高速動車組早期軸承外圈故障的檢測上是十分有效的,有較高的應用價值。

猜你喜歡
故障信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
故障一點通
孩子停止長個的信號
奔馳R320車ABS、ESP故障燈異常點亮
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
故障一點通
故障一點通
故障一點通
主站蜘蛛池模板: 婷婷激情五月网| 人人91人人澡人人妻人人爽| 久久精品视频一| www.av男人.com| 中国国产A一级毛片| 日本www在线视频| 亚洲av无码成人专区| 国产拍揄自揄精品视频网站| 中文字幕啪啪| 日本91视频| 高清亚洲欧美在线看| 国产凹凸视频在线观看| 久久青草精品一区二区三区| 亚洲精品无码不卡在线播放| 九九视频免费在线观看| 欧美精品啪啪| 久久久久久国产精品mv| 在线国产你懂的| 久久香蕉国产线| 中文成人在线视频| 亚洲精品男人天堂| 99这里只有精品在线| 国产a v无码专区亚洲av| 亚洲丝袜中文字幕| 久久99这里精品8国产| 久久这里只有精品8| av在线5g无码天天| 毛片免费高清免费| 国产一区二区视频在线| 亚洲an第二区国产精品| 一本色道久久88| 亚洲一区二区三区香蕉| 欧美日韩国产在线播放| 日韩一区二区三免费高清| 欧美自拍另类欧美综合图区| 欧美va亚洲va香蕉在线| 欧美成在线视频| 国产白浆视频| 91免费精品国偷自产在线在线| 永久在线精品免费视频观看| 男女性午夜福利网站| 欧美色综合网站| 国产成人精品高清不卡在线| 成人综合网址| 97综合久久| 伊人查蕉在线观看国产精品| 激情成人综合网| 亚洲精品男人天堂| 在线无码九区| 色偷偷男人的天堂亚洲av| 亚洲黄网视频| 熟妇无码人妻| 一本一道波多野结衣一区二区 | 国产高清在线精品一区二区三区 | 亚洲αv毛片| 亚洲成人网在线播放| 永久免费AⅤ无码网站在线观看| 真实国产乱子伦视频| 一级做a爰片久久免费| 国产区福利小视频在线观看尤物| 欧美无遮挡国产欧美另类| 精品国产电影久久九九| 亚洲无限乱码| 国产综合另类小说色区色噜噜| 日韩免费毛片视频| 国产鲁鲁视频在线观看| 国产亚洲精品精品精品| 99热国产这里只有精品无卡顿" | 亚洲男人的天堂在线| 国产成人精彩在线视频50| 日韩天堂视频| 99热这里只有精品2| 亚洲天堂网视频| 999在线免费视频| 91九色国产porny| 亚洲精品视频免费看| 97精品国产高清久久久久蜜芽| 日韩国产综合精选| 99爱视频精品免视看| 亚洲色图欧美激情| 国产亚洲精品在天天在线麻豆| 国产av一码二码三码无码|