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

        ?

        一種基于系統(tǒng)保形解耦的力源識別方法

        2017-05-04 03:00:48王淑娟王國巧沈繼紅韓偉民
        船舶力學(xué) 2017年1期
        關(guān)鍵詞:力源二階時域

        王淑娟,王國巧,沈繼紅,韓偉民

        (哈爾濱工程大學(xué) 理學(xué)院,哈爾濱150001)

        一種基于系統(tǒng)保形解耦的力源識別方法

        王淑娟,王國巧,沈繼紅,韓偉民

        (哈爾濱工程大學(xué) 理學(xué)院,哈爾濱150001)

        基于模態(tài)分析法的力源時域識別方法中涉及利用模態(tài)矩陣將系統(tǒng)方程轉(zhuǎn)化為非耦合形式,這僅對比例阻尼或經(jīng)典阻尼系統(tǒng)才能實(shí)現(xiàn)。文章基于二階系統(tǒng)保結(jié)構(gòu)解耦方法建立了新的力源識別數(shù)學(xué)模型。首先利用Lancaster結(jié)構(gòu)建立系統(tǒng)的解耦模型,并求得解耦變換;其次,利用解耦變換推導(dǎo)建立系統(tǒng)非耦合響應(yīng)模型,并在微小時段內(nèi)力源為線性變化的假設(shè)下,推導(dǎo)出具體的力源識別公式。最后,采用精細(xì)逐步積分方法,對模型進(jìn)行解算,由結(jié)構(gòu)動態(tài)響應(yīng)反求力源的時間歷程。數(shù)值實(shí)驗不僅驗證了所提方法可以提高識別精度,而且也驗證了所提方法對非比例阻尼系統(tǒng)也是有效的,這也是文中方法的重點(diǎn)。

        Lancaster結(jié)構(gòu);解耦;精細(xì)逐步積分;線性假設(shè)

        0 引 言

        力源識別是根據(jù)已知系統(tǒng)的動態(tài)特性和實(shí)測的動力響應(yīng)反算結(jié)構(gòu)所受的激勵,是一個較難的問題,但又是結(jié)構(gòu)動態(tài)設(shè)計的關(guān)鍵之一。力源識別問題的研究在各個工程領(lǐng)域都具有重要的意義,如隨著IMO組織對船舶噪聲尤其是艙室噪聲的要求越來越嚴(yán),機(jī)電設(shè)備作為船舶主要的激勵源,其振動控制越來越受到人們的關(guān)注。設(shè)備的低噪聲設(shè)計是解決設(shè)備振動問題的最根本方法,而力源識別又是設(shè)備低噪聲設(shè)計的前提條件之一。

        力源識別方法分為頻域法和時域法。頻域法發(fā)展較早,計算方法較為成熟,精度較高。此方法可確定力譜的均值與方差,但對于力的時間歷程還有一定困難,特別是可能會出現(xiàn)奇異值和不穩(wěn)定現(xiàn)象。而時域識別是20世紀(jì)80年代中期發(fā)展起來的一種方法,直接依據(jù)結(jié)構(gòu)的響應(yīng)時程來識別未知力源的時間歷程,便于工程應(yīng)用,受到工程界的歡迎。其中現(xiàn)有的一類力源時域識別方法利用結(jié)構(gòu)的模態(tài)模型將系統(tǒng)方程轉(zhuǎn)化為非耦合形式建立反演模型,具有一定的工程應(yīng)用價值。

        Desanghere等[1]、Ory等[2]于1985年首次將模態(tài)坐標(biāo)變換方法引入到動態(tài)力源時域識別過程中,將運(yùn)動微分方程轉(zhuǎn)化為非耦合方程形式,將實(shí)際載荷的反演轉(zhuǎn)化為對模態(tài)力的反演,在微小時間間隔內(nèi)力源為階躍力的假設(shè)下,推導(dǎo)出已知系統(tǒng)響應(yīng)等條件下的力源反演公式。唐秀近等[3-4]提出基于模態(tài)分解Duhamel積分的一種較為有效的時域載荷識別模式,并較為系統(tǒng)地對該方法的識別精度、模態(tài)截斷的影響作了初步的研究。初良成等[5]利用正分析的手法來求解反問題,建立了階躍力假設(shè)下的時域載荷識別新方法,克服了識別精度受測量信息數(shù)目限制的缺點(diǎn)。徐倩等[6]應(yīng)用精細(xì)逐步積分建立了一種基于階躍力假設(shè)下的力源識別模式,由結(jié)構(gòu)動態(tài)響應(yīng)求出動態(tài)力的時間歷程。張運(yùn)良等[7]在微小時段內(nèi)力源為線性變化的假設(shè)下,推導(dǎo)了僅由結(jié)構(gòu)的一種響應(yīng)(位移、速度或加速度)進(jìn)行力源識別的遞推連鎖計算公式,以便于實(shí)際應(yīng)用。

        這一類力源識別方法無論是階躍力還是線性力假設(shè),在最初都是利用模態(tài)坐標(biāo)變換將系統(tǒng)多自由度微分方程轉(zhuǎn)化為非耦合形式,然后建立具體的反演公式。但利用坐標(biāo)變換將三個系數(shù)矩陣同時對角化,且僅當(dāng)比例阻尼或經(jīng)典阻尼系統(tǒng)才能夠?qū)崿F(xiàn),這就使得算法的應(yīng)用受到了限制。數(shù)值代數(shù)領(lǐng)域認(rèn)為,通過主坐標(biāo)變換來實(shí)現(xiàn)三個矩陣的同時對角化一般是無法實(shí)現(xiàn)的,因此本文將力源識別中涉及的三矩陣同時對角化轉(zhuǎn)化為Lancaster擴(kuò)展系統(tǒng)的塊陣同時對角化,建立系統(tǒng)解耦模型,求得相應(yīng)解耦變換,并通過理論推導(dǎo)建立系統(tǒng)的非耦合方程表示形式。最后,在微小時間段線性力假設(shè)下推導(dǎo)出力源時域識別公式,并利用精細(xì)逐步積分進(jìn)行模型計算。本文數(shù)值試驗將涉及無阻尼系統(tǒng)和比例阻尼系統(tǒng),驗證本文方法較現(xiàn)有基于模態(tài)分析載荷識別模型精度高;并將利用兩組船舶運(yùn)動數(shù)據(jù)驗證本文提出算法對非比例阻尼系統(tǒng)也適用,且精度較滿意。

        1 基于模態(tài)分析法的力源識別模型

        一個n自由度二階系統(tǒng)往往可由下式來描述:

        其中:M、C、K和F分別表示為質(zhì)量矩陣、阻尼矩陣、剛度矩陣和外力向量。由實(shí)測或有限元分析可得系統(tǒng)的模態(tài)參數(shù):固有頻率ωr、阻尼比ξr及模態(tài)矩陣Φ,其中,利用模態(tài)矩陣做變換x(t)=Φq(t),得

        左乘矩陣ΦT

        若上式為無耦合方程,設(shè)Φ= [Φ1,Φ2,…,Φn]則可寫成分解式

        其中:r=1,2,…,n,qr(t)為第r階模態(tài)坐標(biāo)。設(shè)(2)式右端項為fr(t),將(2)式降階為一階微分方程

        其中:[T(τ)]=eHτ,τ=t-tj,{vp(t )}為方程(3)的特解[5]。利用精細(xì)逐步積分法計算指數(shù)矩陣 [T (τ)][5],便可得模態(tài)坐標(biāo)qr(t)與外力源F(t)之間的關(guān)系式。再由x(t)=Φq(t)推導(dǎo)出結(jié)構(gòu)位移響應(yīng)x(t)與外力源F(t)的關(guān)系,進(jìn)而實(shí)現(xiàn)由結(jié)構(gòu)響應(yīng)反向推導(dǎo)力源。

        顯然,從(1)式到(2)式的轉(zhuǎn)變意味著主坐標(biāo)變換x(t)=Φq(t)將三個矩陣M、C、K同時對角化,將多自由度系統(tǒng)簡化為多個單自由度子系統(tǒng),這個過程也被稱為二階系統(tǒng)解耦[8]。眾所周知,當(dāng)M和K都為對稱陣,且其中之一為正定陣時,廣義特征值分解可實(shí)現(xiàn)M和K的同時對角化,相關(guān)的研究方法已經(jīng)很成熟[8]。三個矩陣的同時對角化不僅要求M、C、K同為對稱陣和某個矩陣為正定陣,還需要滿足KM-1C=CM-1K,但這個條件是很難滿足的,因此工程上涉及到二階系統(tǒng)解耦只能在比例阻尼的假設(shè)下實(shí)現(xiàn)[8]。

        目前數(shù)值代數(shù)領(lǐng)域利用Lancaster解耦來研究二階系統(tǒng)的解耦問題,將三個矩陣同時對角化轉(zhuǎn)化為Lancaster擴(kuò)展系統(tǒng)塊陣同時對角化,已證明幾乎所有二階系統(tǒng)都可以通過這樣的方式實(shí)現(xiàn)解耦。文獻(xiàn)[8]詳細(xì)地研究和介紹了這一問題,并提出基于Sylvester方程的解耦變換求解方法。本文將利用基于Lancaster結(jié)構(gòu)的二階系統(tǒng)解耦方法建立新的力源識別模型,利用解耦變換通過理論推導(dǎo)建立系統(tǒng)受力相應(yīng)模型,假設(shè)力源在微小時間段內(nèi)成線性變化的前提下推導(dǎo)力源計算公式。

        2 基于系統(tǒng)解耦的力源識別模型建立

        2.1 基于Lancaster結(jié)構(gòu)的系統(tǒng)解耦

        顯然,當(dāng)M可逆時(1)式與(4)式對應(yīng)的特征值問題等價,

        左乘可逆矩陣Πl(fā)得

        若Πr、Πl(fā)可實(shí)現(xiàn)(4)式系數(shù)矩陣的塊陣同時對角化,便將(4)式轉(zhuǎn)化為

        其中:MD、CD和KD為對角矩陣。顯然,(5)式與(4)式對應(yīng)的特征值問題等價,若設(shè)代入

        (6)式即為(1)式所描述原始二階系統(tǒng)的非耦合形式。

        2.2 基于精細(xì)逐步積分的力源計算

        且MDr表示對角陣MD對角線上第r個元素,表示矩陣l11第r行第i列所對應(yīng)的元素)表示矩陣l21第r行第i列所對應(yīng)的元素。將(6)式降階為一階微分方程

        取等間距積分步長τ=tj+1-tj,利用tj時刻的值計算tj+1時刻,令s=tj+ξ,則上式可改寫為

        其中:[T (τ)]=eHτ為指數(shù)矩陣,利用精細(xì)積分方法計算[8-10]。

        假設(shè)在時間步長 (tj, tj+1)內(nèi)力r(t)呈線性變化,即

        其中:ξ∈ (tj, tj+1),{r1}和 {r2}在區(qū)間 (tj, tj+1)內(nèi)為時不變向量,且有

        則將(10)式代入(9)式,整理可得線性力源精細(xì)積分格式

        式中:Tij表示矩陣 [T (τ)]的第i行第j列,其余類推。

        2.3 基于線性力假設(shè)的力源識別公式

        假設(shè)在時間步長 (tj, tj+1)內(nèi),力源f(t)是線性變化的,即

        則有

        將上式一階導(dǎo)代入(7)式,再將(7)式代入(11)-(12)式,整理得

        將上式中第2個式子整理得出fi(tj+1)的表達(dá)式代入第1個式子,就可以得到fi( tj)的表達(dá)式,從而求解得力源的數(shù)值解。

        3 算 例

        例1針對文獻(xiàn)[6]中涉及的兩自由度無阻尼系統(tǒng)

        進(jìn)行模型仿真計算。首先利用文獻(xiàn)[8]中的方法構(gòu)造解耦后系統(tǒng)參數(shù)矩陣MD和KD,并利用基于Sylvester方程的解耦變換求解方法計算得解耦變換Πr、Πl(fā)。然后利用Πr、Πl(fā)獲得如(6)式所示的系統(tǒng)非耦合形式。最后在線性力源假設(shè)下,取時間步長△t=0.28,利用(13)式計算得力源估計值。表1給出了本文方法與文獻(xiàn)[6]中方法計算的力源估計值,由表中數(shù)據(jù)可知本文提出的基于系統(tǒng)解耦的力源識別方法精度較高,具體模型計算的絕對誤差如圖1所示。

        圖1 基于系統(tǒng)解耦力源識別絕對誤差Fig.1 The absolute error of load identification based on decouplingmethod

        表1 兩自由度系統(tǒng)力源識別計算結(jié)果Tab.1 Calculation result of dynam ic load identification for two DOF system

        例2本例將對一船舶振動系統(tǒng)進(jìn)行力源識別估算,且該力源為動態(tài)力,具體系統(tǒng)方程如下:

        很顯然,該系統(tǒng)為比例阻尼系統(tǒng),因此既可以利用模態(tài)分析法獲得系統(tǒng)非耦合形式,也可以利用本文方法進(jìn)行獲得系統(tǒng)解耦形式,下面將給出基于這兩種方法的力源識別模型計算情況。同樣選取時間步長τ=0.02,圖2(a)給出了本文方法計算的力源估計值和真實(shí)值曲線,可以看出估計值與真實(shí)值基本吻合,圖2(b)則給出了本文方法的誤差曲線;圖2(c)給出了基于模態(tài)分析法的力源估計值和真實(shí)值曲線,圖2(d)則給出了基于模態(tài)分析法的力源估計誤差曲線;從圖2(b)和圖2(d)的幅值比較來看,本文方法精度較高。

        圖2 (a) 基于系統(tǒng)解耦力源識別計算值與真值比較曲線Fig.2(a)The comparation between the calculated value and the real value of load identification based on decouplingmethod

        圖2 (b) 基于系統(tǒng)解耦力源識別絕對誤差曲線Fig.2(b)The absolute error curve of load identification based on decouplingmethod

        圖2 (c) 基于模態(tài)法力源識別計算值與真值比較曲線Fig.2(c)The comparation between the calculated value and the real value of load identification based onmodal identification

        圖2 (d) 基于模態(tài)法力源識別絕對誤差曲線Fig.2(d)The absolute error curve of load identification based onmodal identification

        例3本例將對文獻(xiàn)[8]中水池實(shí)驗獲得的船舶運(yùn)動系統(tǒng)數(shù)據(jù)進(jìn)行力源識別,所涉及船舶質(zhì)量為425 t,船長為60m,船舶浪向角為0°,其中共測得18 kns和27 kns航速下各20個遭遇頻率的艦船縱搖—升沉運(yùn)動數(shù)據(jù)。

        已知船舶運(yùn)動系統(tǒng)為非比例系統(tǒng),因此,無法利用模態(tài)分析法獲得系統(tǒng)的非耦合方程。利用本文提出的力源識別方法,選取迭代步長τ=0.02進(jìn)行解算,下面將給出兩組不同航速下不同遭遇頻率下的力源識別結(jié)果,圖3和圖4分別給出了計算值與真實(shí)值對比圖,以及相應(yīng)的絕對誤差曲線圖。實(shí)驗驗證了此方法對非比例阻尼系統(tǒng)同樣適用。

        (1)航速18 kns遭遇頻率Ome=0.889:

        (2)數(shù)據(jù)2:航速27 kns遭遇頻率Ome=0.945:

        圖3 (a)力源識別計算值與真值對比圖(Ome=0.889)Fig.3(a)The comparation between the calculated value and the real value of load identification(Ome=0.889)

        圖3 (b)力源識別絕對誤差曲線圖(Ome=0.889)Fig.3(b)The absolute error curve of load identification(Ome=0.889)

        圖4 (a)力源識別計算值與真值對比圖(Ome=0.945) Fig.4(a)The comparation between the calculated value and the real value of load identification(Ome=0.945)

        圖4 (b)力源識別絕對誤差曲線圖(Ome=0.945)Fig.4(b)The absolute error curve of load identification(Ome=0.945)

        4 結(jié) 論

        動態(tài)載荷的確定在結(jié)構(gòu)的分析與研究中意義重大,準(zhǔn)確地識別動態(tài)載荷是工程結(jié)構(gòu)可靠性與安全性的重要保證。然而在工程實(shí)行中,常常是很難對作用于結(jié)構(gòu)的外載荷作直接測量和計算。因此,對動態(tài)載荷的確定只能從可測的系統(tǒng)響應(yīng)信息來間接地確定。由于其廣泛的工程背景和明顯的應(yīng)用價值,使這一問題的理論和應(yīng)用成為研究的熱點(diǎn)。

        本文把現(xiàn)有一類載荷時域識別方法涉及到的解耦問題作為切入點(diǎn),將基于Lancaster結(jié)構(gòu)的二階系統(tǒng)解耦方法應(yīng)用于系統(tǒng)動力學(xué)方程解耦中,并利用解耦變換通過數(shù)學(xué)理論推導(dǎo)提出了基于系統(tǒng)解耦的載荷識別數(shù)學(xué)模型,最后采用精細(xì)時程積分法進(jìn)行模型計算。本文提出的模型一方面克服了模態(tài)分析法只適用于比例阻尼系統(tǒng)解耦的問題,另一方面從數(shù)值試驗來看,針對比例阻尼系統(tǒng)的力源識別精度要高于現(xiàn)有方法,而針對非比例阻尼系統(tǒng)該方法仍然有效。基于Lancaster結(jié)構(gòu)的二階系統(tǒng)解耦方法是數(shù)值代數(shù)領(lǐng)域的前沿研究之一,一直處于理論研究的階段,因此本文方法也是數(shù)學(xué)理論研究在工程應(yīng)用中的一個嘗試。

        參 考 文 獻(xiàn):

        [1]Desanghere G,Snoeys R.Indirect identification of excitation foreces bymodal coordinate transformation[C]//Proceedings of the 3rd IMAC.Florida,USA,1985:685-690.

        [2]Ory G,Glaser H,Hlzdeppe D.The reconstruction of force function based on aeroelasticity and structural dynamics[M]// Proceedings of 2nd International Symposium On Aeroelasticity and Structural Dynamics.Aschen,FRG,1985:164-168.

        [3]唐秀近.動態(tài)力識別的時域方法[J].大連工學(xué)院學(xué)報,1987,26(4):21-27. Tang Xiujin.Dynamic load identification by time domainmethod[J].Journal of DaLian Institute of Technology,1987,26 (4):21-27.

        [4]唐秀近.時域動態(tài)載荷識別的精度問題[J].大連理工大學(xué)學(xué)報,1990,30(1):31-37. Tang Xiujin.Precision problems of dynamic load identification in time domain[J].Journal of DaLian University of Technology,1990,30(1):31-37.

        [5]Chu Liangcheng,et al.Dynamic load identification in time domain[J].China Ocean Engineering,1990,5(3):279-286.

        [6]徐 倩,文祥榮,孫守光.結(jié)構(gòu)動態(tài)力源識別的精細(xì)逐步積分法[J].計算力學(xué)學(xué)報,2002,19(1):53-57. Xu Qian,Wen Xiangrong,Sun Shouguang.High precision direct integration scheme for structural dynamic load identification[J].Chinese Journal of ComputationalMechanics,2002,19(1):53-57.

        [7]張運(yùn)良,林 皋,王永學(xué),李志軍,李廣偉.一種改進(jìn)的動態(tài)力源時域識別方法[J].計算力學(xué)學(xué)報,2004,21(2):209-215. Zhang Yunliang,Lin Gao,Wang Yongxue,Li Zhijun,Li Guangwei.An improved method of dynamic load identification in time domain[J].Chinese Journal of ComputationalMechanics,2004,21(2):209-215.

        [8]王淑娟.基于Lancaster結(jié)構(gòu)的二階系統(tǒng)解耦算法研究及其應(yīng)用[D].哈爾濱:哈爾濱工程大學(xué),2009:1-145. Wang Shujuan.Quadratic system decoupling algorithm researches based on lancaster structure and its applications[D]. Harbin:Harbin Engineering University,2009:1-145.

        [9]王 靜,陳海波,王 靖.基于精細(xì)積分的沖擊力源時域識別方法研究[J].振動與沖擊,2013,32(20):81-85. Wang Jing,Chen Haibo,Wang Jing.Impulsive load identification in time domain based on precise time-integrationmethod [J].Journal of Vibration and Shock,2013,32(20):81-85.

        [10]鐘萬勰.結(jié)構(gòu)動力方程的精細(xì)時程積分法[J].大連理工大學(xué)報,1994,32(2):131-136. ZhongWanxie.On precise time-integrationmethod for structural dynamics[J].Journal of DaLian University of Technology,1994,32(2):131-136.

        [11]林家浩,沈為平,威廉斯FW.受演變隨機(jī)激勵結(jié)構(gòu)響應(yīng)的精細(xì)逐步積分法[J].大連理工大學(xué)學(xué)報,1995,35(5): 600-605. Lin Jiahao,Shen Weiping.High precision direct integration scheme for analysing structural nonstationary random responses[J].Journal of DaLian University of Technology,1995,35(5):600-605.

        [12]鐘萬勰,楊再石.連續(xù)時間LQ控制主要本征對的算法[J].應(yīng)用數(shù)學(xué)和力學(xué),1991,12(1):45-50. Zhong Wanxie,Yang Zaishi.On the computation of themain eigen-pairs of the continuous-time linear quadratic control problem[J].Applied Mathematics and Mechanics,1991,12(1):45-50.

        A load identification method based on decoup ling of Lancaster structure

        WANG Shu-juan,WANG Guo-qiao,SHEN Ji-hong,HANWei-min
        (School of Science,Harbin Engineering University,Harbin 150001,China)

        The way of the load identification in time domain which is based on themodal analysis that the system equation is transformed into a non-coupling form bymodalmatrix,is only applicable for the proportional damping or classical damping system.In this paper,based on structure-preserving decoupling method of second-order system,a new mathematicalmodel of load identification was established.The system-decouplingmodel was established based on the Lancaster structure and the decoupling transformation was obtained;then,themodel of uncoupled response system based on the decoupling transformation was established,and on the assumption of linear load identification in tiny-time,the formula of load identification was deduced;finally,themodelwas calculated based on high precision direct integration scheme, and the history of force identification was obtained by structural dynamic response.The experiments verify that the proposed method could improve the accuracy of recognition,and is effective for the non-proportional damping system.And it is the focus of the paper.

        Lancaster structure;decoupling;high precision direct integration scheme;linear assumption

        TB535

        :Adoi:10.3969/j.issn.1007-7294.2017.01.007

        2016-09-09

        國家自然科學(xué)基金資助項目(11226324);黑龍江省基金資助項目(A201407)

        王淑娟(1982-),女,碩士生導(dǎo)師,E-mail:wangshujuan@hrbeu.edu.cn;

        王國巧(1990-),女,碩士研究生,E-mail:Wgq19900424@163.com;

        沈繼紅(1966-),男,教授,博士生導(dǎo)師;

        韓偉民(1990-);男,碩士研究生。

        1007-7294(2017)01-0052-09

        猜你喜歡
        力源二階時域
        “童心向黨”征集作品展示
        未來教育家(2021年9期)2021-12-24 08:24:22
        Asymmetric coherent rainbows induced by liquid convection?
        一類二階迭代泛函微分方程的周期解
        一種光傳送網(wǎng)的建模及其價值評估
        軟件(2020年3期)2020-04-20 01:45:48
        一類二階中立隨機(jī)偏微分方程的吸引集和擬不變集
        基于時域信號的三電平逆變器復(fù)合故障診斷
        二階線性微分方程的解法
        一類二階中立隨機(jī)偏微分方程的吸引集和擬不變集
        包力源、鐘琦翔作品
        基于極大似然準(zhǔn)則與滾動時域估計的自適應(yīng)UKF算法
        午夜福利电影| 国产精品亚洲一区二区三区在线| av网站免费在线浏览| 欧美性白人极品1819hd| 内射合集对白在线| 亚洲欧美日韩一区二区三区在线| 91最新免费观看在线| 538在线视频| 精品熟女av中文字幕| 亚洲精品国产第一综合色吧| 久久蜜桃一区二区三区| 亚洲熟女乱一区二区三区| av人摸人人人澡人人超碰下载| 农村欧美丰满熟妇xxxx| 亚洲性无码av在线| 久久中文字幕久久久久91| 中文字幕人妻av四季| 午夜人妻久久久久久久久| 日韩人妻无码精品久久久不卡| 国产精品视频一区二区三区四| 日韩人妻无码精品系列专区无遮| 97人妻精品一区二区三区免费| 欧美午夜理伦三级在线观看| 国产特级毛片aaaaaa| 精品无码AV无码免费专区| 国产黄色污一区二区三区| 日韩人妻久久中文字幕| 无码人妻一区二区三区兔费| 人妻在线日韩免费视频| 亚洲中文字幕久爱亚洲伊人 | 青青草视频在线播放81| av手机在线观看不卡| 国产一区二区三区在线电影| 3344永久在线观看视频| 99久久久精品国产性黑人| 亚洲无毛成人在线视频| 乱人伦中文视频在线| 欧美粗大无套gay| 亚洲免费无毛av一区二区三区| 一级一片内射视频网址| 亚洲av无码无限在线观看|