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

        ?

        基于扭擺的中低頻響應(yīng)譜修正方法研究

        2020-07-21 10:06:44惠安民閆明馮麟涵楊寧
        兵工學(xué)報(bào) 2020年6期
        關(guān)鍵詞:振子乘法修正

        惠安民,閆明,馮麟涵,楊寧

        (1.沈陽工業(yè)大學(xué) 機(jī)械工程學(xué)院,遼寧 沈陽 110870;2.海軍研究院,北京 100161)

        0 引言

        實(shí)船或浮動(dòng)平臺(tái)水下非接觸爆炸試驗(yàn)是檢測(cè)艦艇及其艦載設(shè)備抗沖擊能力的重要試驗(yàn)方法。試驗(yàn)中需要準(zhǔn)確測(cè)量沖擊響應(yīng)譜,以標(biāo)定被測(cè)試設(shè)備的沖擊環(huán)境。據(jù)國(guó)內(nèi)外公開文獻(xiàn)報(bào)道,測(cè)量沖擊響應(yīng)譜的方法,一般為通過實(shí)船水下非接觸爆炸測(cè)量沖擊時(shí)刻的加速度時(shí)間歷程,再應(yīng)用Duhamel積分獲得沖擊響應(yīng)譜[1-4],但該方法很難獲取高質(zhì)量的沖擊信號(hào)。究其原因,主要是加速度信號(hào)中含有的趨勢(shì)項(xiàng)誤差會(huì)使沖擊響應(yīng)譜在中低頻段出現(xiàn)較明顯的失真現(xiàn)象,對(duì)艦艇抗沖擊設(shè)計(jì)產(chǎn)生較大誤導(dǎo)[5]。為從試驗(yàn)數(shù)據(jù)中獲取較可靠的結(jié)果,諸多學(xué)者提出了各種數(shù)學(xué)修正方法。Grillo[6]提出基于最小二乘法和樣條擬合曲線法修正趨勢(shì)項(xiàng)誤差,其原理主要是對(duì)測(cè)量得到的加速度信號(hào)進(jìn)行多次積分?jǐn)M合再微分。Goebel等[7]、李海廣等[8]、于大鵬等[9]、李慧浩等[10]、梁兵等[11]分別依據(jù)經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)法、低通濾波法、小波變換法等不同修正方法,對(duì)加速度的趨勢(shì)項(xiàng)誤差進(jìn)行修正研究。但這些數(shù)學(xué)方法會(huì)引入主觀誤差,導(dǎo)致修正后的結(jié)果還需要其他手段進(jìn)行驗(yàn)證;除上述數(shù)學(xué)方法外,還可以通過給加速度傳感器加裝機(jī)械濾波器來修正沖擊響應(yīng)譜測(cè)量時(shí)中低頻段出現(xiàn)的失真現(xiàn)象[3],但該方法對(duì)沖擊響應(yīng)譜的修正效果有限;或根據(jù)沖擊響應(yīng)譜定義,應(yīng)用中低頻彈簧振子和簧片儀對(duì)沖擊響應(yīng)譜進(jìn)行直接測(cè)量。但兩種設(shè)備都無法單獨(dú)完成測(cè)量工作。由于簧片儀自身結(jié)構(gòu)特點(diǎn),在測(cè)量低頻譜時(shí)位移響應(yīng)較大,極易因根部應(yīng)力過大產(chǎn)生斷裂現(xiàn)象,且測(cè)量過程中產(chǎn)生的響應(yīng)值為多階模態(tài)疊加的總和,因此僅適用于5 Hz以上的測(cè)量工作[12]。而低頻彈簧振子在測(cè)量5 Hz以下低頻譜時(shí),由于自身固有頻率低,彈簧拉伸、壓縮位移形變量過大、質(zhì)量塊重等缺點(diǎn),均不利于沖擊譜的修正工作。

        本文針對(duì)趨勢(shì)項(xiàng)誤差的作用機(jī)理進(jìn)行分析,提出一種扭擺模型用以修正沖擊響應(yīng)譜中低頻段趨勢(shì)項(xiàng)誤差。對(duì)該扭擺進(jìn)行縱向沖擊機(jī)沖擊試驗(yàn),并根據(jù)趨勢(shì)項(xiàng)誤差產(chǎn)生機(jī)理設(shè)計(jì)一種基于EMD與最小二乘法聯(lián)合修正趨勢(shì)項(xiàng)誤差的數(shù)學(xué)方法。分別對(duì)比彈簧振子、扭擺以及EMD與最小二乘法聯(lián)合修正方法的修正數(shù)據(jù),驗(yàn)證所提出的扭擺修正沖擊響應(yīng)譜中低頻段的可行性。

        1 沖擊響應(yīng)譜簡(jiǎn)介

        圖1 沖擊響應(yīng)譜模型Fig.1 Shock response spectrum model

        Gaberson[15]為更加全面地反映沖擊信息,發(fā)明了一種對(duì)數(shù)四坐標(biāo)沖擊譜繪制方法,將上述3種沖擊響應(yīng)在同一個(gè)圖中表達(dá)出來,如圖2所示。圖2中:與縱軸呈+45°、-45°的線分別代表位移譜和加速度譜;對(duì)數(shù)四坐標(biāo)系中的沖擊譜形似三折線,左側(cè)折線表示彈簧振子的極限譜位移,中間直線表示極限譜速度,右側(cè)折線表示極限譜加速度。依據(jù)三折線譜,只要確定設(shè)備的安裝頻率,即可以確定該安裝頻率設(shè)備在沖擊載荷作用下產(chǎn)生的最大相對(duì)位移D、相對(duì)速度v和絕對(duì)加速度a[16].

        圖2 沖擊譜與三折線譜Fig.2 Shock spectrum and trilinear spectrum

        2 趨勢(shì)項(xiàng)誤差對(duì)沖擊譜的作用機(jī)理與修正

        2.1 趨勢(shì)項(xiàng)誤差的作用機(jī)理

        在實(shí)際測(cè)量沖擊譜時(shí),含有趨勢(shì)項(xiàng)誤差的加速度信號(hào)所積分得到的速度以及位移曲線在沖擊結(jié)束后,其數(shù)值并不歸0,甚至趨于無限大。這與艦艇在承受水下非接觸爆炸沖擊后,最終穩(wěn)定靜止于水面,其速度值為0的事實(shí)并不相符。此時(shí)含有趨勢(shì)項(xiàng)的加速度信號(hào)繪制成的沖擊譜在低頻區(qū)的譜位移與譜速度高于正常值幾個(gè)數(shù)量級(jí),呈現(xiàn)明顯的失真現(xiàn)象,對(duì)艦船的抗沖擊設(shè)計(jì)產(chǎn)生極其嚴(yán)重的誤導(dǎo)。沖擊譜在低頻區(qū)域的數(shù)值,一般取決于殘余沖擊響應(yīng)譜[17]。殘余譜是指沖擊輸入結(jié)束后,系統(tǒng)自由振動(dòng)階段的響應(yīng)。由文獻(xiàn)[17-19]可知,單自由度彈簧- 質(zhì)量- 阻尼系統(tǒng)的相對(duì)位移沖擊響應(yīng)可由(1)式表示:

        (1)

        (2)

        (3)

        再假設(shè)單自由度的彈簧- 質(zhì)量系統(tǒng)的阻尼為0,即有ξ=0,則(3)式可化簡(jiǎn)為(4)式:

        (4)

        由三角函數(shù)性質(zhì),該方程可進(jìn)一步化簡(jiǎn)為(5)式:

        (5)

        式中:φ為相位角。

        由此可求出此時(shí)的最大位移響應(yīng)為

        (6)

        圖3 最大殘余位移響應(yīng)ze、沖擊結(jié)束時(shí)刻速度響應(yīng)對(duì)系統(tǒng)最大位移響應(yīng)zmax的影響Fig.3 Effects of maximum residual displacement response ze and velocity response at the end of shock on system maximum displacement response zmax

        2.2 趨勢(shì)項(xiàng)誤差的EMD與最小二乘法聯(lián)合修正方法

        通過2.1節(jié)趨勢(shì)項(xiàng)誤差對(duì)沖擊譜作用機(jī)理的分析,可以為修正趨勢(shì)項(xiàng)誤差方法提供有效的理論根據(jù),即消除趨勢(shì)項(xiàng)的同時(shí),還要保證采樣結(jié)束時(shí)刻速度歸零。因此根據(jù)以上理論,本文設(shè)計(jì)一種基于EMD與最小二乘法的聯(lián)合修正方法。該修正方法經(jīng)過多次擬合與數(shù)學(xué)變換,仍無法解決人為引入的主觀誤差,但在一定程度上可以驗(yàn)證試驗(yàn)結(jié)果的正確性。

        EMD可把原始信號(hào)中的各種頻率成分以固有模態(tài)的函數(shù)形式,由高頻到低頻顯示出來[20-23],如(7)式所示:

        (7)

        式中:x(t)表示原信號(hào);cj(t)表示EMD得到的固有模態(tài)函數(shù)分量;r(t)表示殘余信號(hào),通常代表原信號(hào)的趨勢(shì)或均值;N表示頻率成分個(gè)數(shù)。因此,對(duì)含有趨勢(shì)項(xiàng)誤差的信號(hào)進(jìn)行EMD法分解,得到的殘余信號(hào)r(t)可以作為原信號(hào)中的趨勢(shì)項(xiàng)誤差。雖然該方法在消除趨勢(shì)項(xiàng)誤差方面很有優(yōu)勢(shì),但卻無法保證修正后的速度信號(hào)在采樣結(jié)束時(shí)刻的值為0 m/s,使得修正效果并不完美。因此,需要進(jìn)一步對(duì)數(shù)據(jù)進(jìn)行修正。最小二乘法是一種通過最小化誤差的平方和來尋找數(shù)據(jù)最佳的函數(shù)匹配方法,可以保證修正后的時(shí)域信號(hào)在采樣結(jié)束時(shí)刻數(shù)值為零,在趨勢(shì)向誤差修正中得到一定應(yīng)用。因此,通過以EMD方法為主、最小二乘法進(jìn)行二次修正的聯(lián)合修正方法,可以對(duì)趨勢(shì)項(xiàng)誤差進(jìn)行較良好的修正。

        3 扭擺模型

        由于數(shù)學(xué)修正方法中的不確定性,工程上常采用彈簧振子與簧片儀直接修正由趨勢(shì)項(xiàng)誤差引起的中低頻沖擊響應(yīng)譜的失真現(xiàn)象,并將修正結(jié)果與數(shù)學(xué)修正方法進(jìn)行相互印證,以保證在中低頻范圍內(nèi)得到較準(zhǔn)確且完整的沖擊譜。但彈簧振子與簧片儀本身也有著局限性,如低頻彈簧振子體積、質(zhì)量巨大,極不利于工程技術(shù)人員攜帶;而簧片儀雖然體積較小,但由于本身為懸臂梁結(jié)構(gòu),測(cè)量得到的數(shù)據(jù)為多階模態(tài)相互疊加的綜合響應(yīng),且在強(qiáng)沖擊環(huán)境下,固定端易發(fā)生破損、折斷等現(xiàn)象。

        根據(jù)傳統(tǒng)中低頻沖擊測(cè)量裝置中存在的不足,本文提出一種用以測(cè)量沖擊譜的扭擺裝置(以下簡(jiǎn)稱扭擺)。該扭擺在結(jié)構(gòu)上與簧片儀相似,均為集中質(zhì)量的懸臂梁結(jié)構(gòu),但扭擺為獲得更大的量程,改變簧片儀固定端完全約束的連接方式,改為鉸鏈連接方式,并在鉸鏈段加裝扭簧作為彈性元件。扭擺的工作原理為:扭擺鉸鏈與基礎(chǔ)剛性連接,當(dāng)受到?jīng)_擊載荷的沖擊時(shí),其集中質(zhì)量的懸臂端以固定頻率圍繞鉸鏈開始往復(fù)振動(dòng)。其中彈簧振子、簧片儀以及扭擺的原理圖如圖4所示。

        圖4 彈簧振子、簧片儀以及扭擺的原理圖Fig.4 Schematic diagram of spring oscillater, reed gage and torsional pendulum

        3.1 扭擺運(yùn)動(dòng)方程及其近似解析解

        設(shè)旋轉(zhuǎn)中心到質(zhì)量塊重心的距離為l,擺桿質(zhì)量為MRod,扭簧剛度為kTP,其擺角為θ,扭擺原理圖如圖5所示。

        圖5 扭擺的原理圖Fig.5 Schematic diagram of torsional pendulum

        建立扭擺的運(yùn)動(dòng)微分方程為

        (8)

        式中:g為重力加速度。

        引入sinθ的泰勒級(jí)數(shù)展開式,并取其中的前兩項(xiàng)代入原方程,可得(9)式:

        (9)

        對(duì)(9)式進(jìn)行無因次化后,可得方程

        (10)

        對(duì)(10)式應(yīng)用Ritz-Galerkin平均法求解其近似解析解,其解的形式為

        (11)

        式中:κ1,κ2,…,κn為待定的權(quán)系數(shù);φ1(t),φ2(t),…,φn(t)為假定時(shí)間函數(shù),n為正整數(shù)。假設(shè)上述微分方程的1階近似解為(12)式:

        (12)

        式中:A0為1階近似解的振幅;ω為1階近似解的角頻率。則(10)式可表述為

        (13)

        (14)

        (13)式代入(14)式中并整理,可到非線性方程為

        (15)

        即扭擺非線性方程的解為

        (16)

        (16)式中,對(duì)于非平凡解而言A0≠0,故該方程的解可表示為

        (17)

        由(17)式可知,扭擺振幅A0是影響扭擺系統(tǒng)非線性的關(guān)鍵因素。在較小擺幅下,扭擺固有圓頻率穩(wěn)定,且呈現(xiàn)弱非線性性質(zhì),線性性質(zhì)較明顯。該段區(qū)域可以很好地應(yīng)用于低頻段沖擊譜的修正工作中。但隨著沖擊輸入的增大,振幅A0超過一定限度后,扭擺系統(tǒng)進(jìn)入強(qiáng)非線性區(qū)域,導(dǎo)致固有頻率隨著振幅而改變,此時(shí)修正低頻譜值不可信。為分析扭擺響應(yīng)的近似線性范圍,本文進(jìn)行如下沖擊試驗(yàn)。

        3.2 扭擺幾何參數(shù)與樣機(jī)

        根據(jù)(17)式,設(shè)計(jì)固有頻率分別為6 Hz、10 Hz和20 Hz的扭擺,各扭擺的幾何參數(shù)如表1所示。

        表1 扭擺樣機(jī)的幾何參數(shù)Tab.1 Geometric parameters of torsional pendulum

        如圖6所示為20 Hz扭擺樣機(jī)圖。

        圖6 20 Hz扭擺樣機(jī)圖Fig.6 20 Hz torsional pendulum prototype

        4 彈簧振子與扭擺的沖擊試驗(yàn)

        由于扭擺受到?jīng)_擊后,其響應(yīng)值為擺動(dòng)角度,對(duì)沖擊譜而言,擺動(dòng)角度并無實(shí)際物理意義。為分析扭擺對(duì)中低頻段沖擊響應(yīng)譜的趨勢(shì)項(xiàng)誤差修正數(shù)據(jù)的可行性與可靠性,需先行分析相同沖擊條件下中低頻彈簧振子的最大位移響應(yīng)與同頻率下扭擺的最大擺角響應(yīng)間的關(guān)系,同時(shí)還可以驗(yàn)證擺桿的近似線性擺角范圍。

        4.1 沖擊試驗(yàn)原理

        沖擊試驗(yàn)中,首先通過加速度傳感器測(cè)量彈簧振子與扭擺安裝基座的加速度信號(hào),并通過位移傳感器直接獲取彈簧振子在沖擊過程中的最大位移響應(yīng)。其次通過角位移傳感器測(cè)量沖擊過程中扭擺的最大角度響應(yīng),并建立彈簧振子最大位移響應(yīng)與扭擺最大角度響應(yīng)的函數(shù)關(guān)系。最后通過該函數(shù)關(guān)系,應(yīng)用扭擺修正趨勢(shì)項(xiàng)誤差,與彈簧振子修正趨勢(shì)項(xiàng)誤差和EMD聯(lián)合最小二乘法修正趨勢(shì)項(xiàng)誤差進(jìn)行對(duì)比。

        4.2 沖擊試驗(yàn)

        沖擊試驗(yàn)采用500 kg級(jí)垂向沖擊試驗(yàn)機(jī)提供沖擊載荷,采用中國(guó)ECON公司生產(chǎn)的MI7016型16通道數(shù)據(jù)采集儀進(jìn)行數(shù)據(jù)采集,采用中國(guó)東華公司生產(chǎn)的DH131E型壓電式加速度傳感器記錄彈簧振子的瞬時(shí)加速度變化,采用美國(guó)MTS公司生產(chǎn)的CS系列磁致非接觸型位移傳感器對(duì)彈簧振子的位移響應(yīng)進(jìn)行測(cè)量,采用中國(guó)Robo Brain公司生產(chǎn)的RB100LA型角位移傳感器對(duì)扭擺的擺角進(jìn)行測(cè)量。如圖7所示為沖擊試驗(yàn)圖。圖7中在沖擊臺(tái)面上還布置有加速度傳感器,用以測(cè)量沖擊臺(tái)的瞬時(shí)加速度變化量。

        圖7 沖擊試驗(yàn)圖Fig.7 Impact test

        圖8 某次沖擊載荷加速度時(shí)域曲線波形Fig.8 Time domain wave forms of acceleration of a certain impact load

        沖擊載荷是根據(jù)德國(guó)BV043/85沖擊標(biāo)準(zhǔn)推薦使用的組合半正弦波形沖擊載荷,該組合半正弦波為沖擊載荷的加速度時(shí)域曲線,分為正波和負(fù)波,正波波峰較高、脈寬較短,負(fù)波波峰較低、脈寬較長(zhǎng),且正波對(duì)時(shí)間積分的面積等于負(fù)波對(duì)時(shí)間積分的面積。圖8(a)所示為某次沖擊載荷與經(jīng)帶通濾波后加速度時(shí)域曲線,圖8(b)所示為該沖擊載荷理論加載波形和經(jīng)濾波后實(shí)際加載波形對(duì)比。圖8(b)中,a+為沖擊載荷加速度正弦正波峰值,a-為沖擊載荷加速度正弦負(fù)波峰值,t+為正波脈寬,t-為負(fù)波脈寬,v+為正波對(duì)時(shí)間積分的面積,v-為負(fù)波對(duì)時(shí)間積分的面積,該積分面積代表沖擊速度。

        在該沖擊載荷下,彈簧振子與扭擺的響應(yīng)如圖9所示,其中,圖9(a)為彈簧振子位移響應(yīng)時(shí)域曲線,圖9(b)為扭擺擺角響應(yīng)時(shí)域曲線。

        圖9 沖擊載荷下彈簧振子與扭擺時(shí)域響應(yīng)曲線Fig.9 Time domain response curves of spring oscillater and torsion pendulum under impact load

        沖擊試驗(yàn)臺(tái)液壓壓強(qiáng)由小逐漸增大,彈簧振子的最大沖擊位移響應(yīng)與對(duì)應(yīng)頻率下扭擺的最大擺角響應(yīng)如圖10所示。通過觀察圖10中的數(shù)據(jù)點(diǎn)發(fā)現(xiàn),二者在一定范圍內(nèi)存在較良好的線性關(guān)系,因此基于最小二乘法,應(yīng)用MATLAB軟件,計(jì)算得到擬合曲線的線性段關(guān)系式如表2所示。

        圖10 彈簧振子最大位移響應(yīng)與扭擺最大擺角響應(yīng)的關(guān)系圖Fig.10 Relationship between the maximum displacement response of spring oscillater and the maximum swing angle response of torsional pendulum

        通過圖10中擬合曲線可以看出,在扭擺擺角小于20.06°(0.35 rad)的范圍內(nèi),可以把扭擺的最大擺角響應(yīng)轉(zhuǎn)化為沖擊譜中具有物理意義的彈簧振子最大位移響應(yīng),并對(duì)中低頻段的沖擊譜進(jìn)行測(cè)量與修正工作,且固有頻率較高的扭擺保持線性區(qū)域的能力高于低頻扭擺。

        4.3 扭擺修正沖擊譜趨勢(shì)項(xiàng)誤差

        通過擬合出的線性段方程,可得到扭擺最大擺角響應(yīng)與彈簧振子最大位移響應(yīng)的關(guān)系。為驗(yàn)證該關(guān)系的準(zhǔn)確性以及扭擺修正中低頻沖擊譜方法的可行性,取上述試驗(yàn)線性區(qū)間內(nèi)彈簧振子、扭擺的測(cè)量值與沖擊臺(tái)臺(tái)面加速度,通過Duhamel積分獲得沖擊響應(yīng)譜,將EMD和最小二乘法聯(lián)合修趨勢(shì)項(xiàng)誤差后的沖擊譜進(jìn)行對(duì)比,如圖11所示。

        表2 彈簧振子與扭擺最大響應(yīng)關(guān)系擬合線性曲線Tab.2 Fitting curve of maximum response relationship between spring oscillator and torsional pendulum at sample point frequency

        圖11 彈簧振子、扭擺與EMD和最小二乘法聯(lián)合修正中低頻沖擊譜圖Fig.11 Middle and low frequency impact spectra corrected by spring oscillater, torsional pendulum, EMD and least square method

        通過圖11可以看出:EMD與最小二乘法聯(lián)合修正、扭擺對(duì)沖擊譜修正的結(jié)果,與彈簧振子對(duì)沖擊譜結(jié)果較相似,三者相互印證修正了結(jié)果的準(zhǔn)確性;在小擺幅下扭擺對(duì)沖擊譜有較良好的修正與驗(yàn)證作用。由于該沖擊臺(tái)臺(tái)面的橫向固有頻率約在2~3 Hz之間,修正后的數(shù)據(jù)在2~3 Hz間存在大于譜位移的現(xiàn)象。

        5 結(jié)論

        2)通過縱向沖擊臺(tái)沖擊試驗(yàn)發(fā)現(xiàn),本文提出的扭擺在采樣點(diǎn)分別為6 Hz、10 Hz和20 Hz條件下且保證扭擺擺桿振動(dòng)幅度在0~20.06°范圍內(nèi),對(duì)中低頻沖擊譜的修正較合理,該修正數(shù)據(jù)通過與彈簧振子修正數(shù)據(jù)、EMD和最小二乘法聯(lián)合修正數(shù)據(jù),可互相印證修正結(jié)果,同時(shí)也驗(yàn)證了扭擺對(duì)修正中、低頻沖擊譜的可行性與準(zhǔn)確性。

        猜你喜歡
        振子乘法修正
        算乘法
        Some new thoughts of definitions of terms of sedimentary facies: Based on Miall's paper(1985)
        修正這一天
        快樂語文(2021年35期)2022-01-18 06:05:30
        我們一起來學(xué)習(xí)“乘法的初步認(rèn)識(shí)”
        彈簧振子問題的分析與求解
        《整式的乘法與因式分解》鞏固練習(xí)
        合同解釋、合同補(bǔ)充與合同修正
        法律方法(2019年4期)2019-11-16 01:07:28
        把加法變成乘法
        非線性Duffing擾動(dòng)振子共振機(jī)制的研究
        軟件修正
        少妇高潮惨叫久久久久电影| 亚洲愉拍99热成人精品热久久| 成人麻豆视频免费观看| 无码gogo大胆啪啪艺术| 久久精品人人做人人爽| 欧美日韩性高爱潮视频| 我也色自拍俺也色自拍| 亚洲综合网国产精品一区| 一区二区三区国产亚洲网站| 亚洲AV无码精品色欲av| 蜜臀久久久精品国产亚洲av| 国产精品成人亚洲一区| 久久精品麻豆日日躁夜夜躁| 欧美a级在线现免费观看| av男人的天堂手机免费网站| 最新露脸自拍视频在线观看| 性色av无码久久一区二区三区| 久久人人97超碰超国产| 亚洲国产免费公开在线视频 | 亚洲av成人永久网站一区| 免费大片黄国产在线观看| 中文字幕熟妇人妻在线视频| 国产成人一区二区三区高清| 男女视频一区二区三区在线观看| 女同同志熟女人妻二区| 国产精品久久久久aaaa| 国产啪精品视频网给免丝袜| 在线视频一区二区三区中文字幕| 性感美女脱内裤无遮挡| 亚洲人成无码区在线观看| 香蕉久久夜色精品国产2020| 白色月光在线观看免费高清 | 国产在线一区二区三区av| 99热国产在线| 在线看亚洲一区二区三区| 亚洲av网一区二区三区| 亚洲另类激情综合偷自拍图 | 精品蜜臀国产av一区二区| 国产精品久人妻精品老妇| 色一情一区二| 在线播放中文字幕一区二区三区|