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

隨機激勵下滯遲系統的穩態響應閉合解1)

2017-07-03 14:59:59陳林聰孫建橋
力學學報 2017年3期
關鍵詞:系統

劉 俊 陳林聰,2) 孫建橋

?(華僑大學土木工程學院,福建廈門361021)?(加州大學Merced分校工程學院,美國加利福尼亞州95343)

-動力學與控制

隨機激勵下滯遲系統的穩態響應閉合解1)

劉 俊?陳林聰?,2)孫建橋?

?(華僑大學土木工程學院,福建廈門361021)?(加州大學Merced分校工程學院,美國加利福尼亞州95343)

滯遲系統屬于一類典型的強非線性系統,滯遲力不僅取決于系統的瞬時變形,還與變形歷程有關.雖然滯遲系統的隨機振動問題已被廣泛研究,但至今尚未得到滯遲系統隨機響應概率密度函數的精確閉合解.本文運用迭代加權殘值法獲得了高斯白噪聲激勵下Bouc-Wen滯遲系統穩態響應概率密度函數的近似閉合解.首先,運用等效線性化法求出系統的穩態高斯概率密度函數;然后以此構造權函數,應用加權殘值法求得了系統指數多項式形式的非高斯概率密度函數;最后引入迭代的過程,逐步優化權函數,提高計算所得結果的精度.以隨機地震激勵下鋼纖維陶粒混凝土結構的穩態響應作為算例,其中Bouc-Wen模型的參數是基于擬靜力學試驗數據,并應用最小二乘法辨識獲得.與Monte Carlo模擬結果相比,等效線性化法得到的結果精度較差;由加權殘值法得到的結果能夠表現出非線性特征,但其精度依然無法令人滿意;采用迭代加權殘值法得到的近似閉合解與Monte Carlo模擬的結果吻合非常好;對于較強隨機激勵情形,采用漸進迭代加權殘值法具有較高的求解效率,所獲得的理論解析解具有較高的精度.結果表明,所獲得的近似閉合解不僅對于土木工程領域具有重要的實際應用價值,而且還可作為檢驗其他非線性系統隨機響應預測方法的精度的標準.

滯遲系統,穩態響應,閉合解,迭代加權殘值法,鋼纖維陶?;炷两Y構

引言

滯遲現象如彈塑性[13]、鐵電性[4]、形狀記憶合金材料[5]等常出現在科學研究和工程實際的不同領域,強烈震動載荷下的結構系統通常表現出滯遲現象[67].學術界提出了許多模型描述這些滯遲關系,如雙線性模型[89],Ramberg-Osgood模型[10],Bouc-Wen模型[67],Ozdemir模型[11],Masing模型[12],Duhem 模型[13],Preisach模型[14]等,其中Bouc-Wen模型是較為通用的一種.

滯遲力不僅與系統當前的狀態有關,而且還與系統過去的狀態有關.因此,滯遲動力學系統屬于一類典型的強非線性系統.近年來,滯遲系統的隨機振動問題已被廣泛研究,但至今尚未獲得精確的閉合解,因此,學術界提出了各種近似解法.如,Iwan[15]采用基于克雷洛夫-包哥留波夫假設的等效線性化法研究了雙線性滯遲系統的隨機響應.Roberts[16-17]分別采用標準隨機平均法及能量包線隨機平均法研究了雙線性滯遲系統的隨機響應.Bouc[7]釆用FPK方程法研究了滯遲系統的隨機響應.Wen[7,18]采用線性化方法研究了滯遲系統的隨機響應.Zhu等[19]采用能量包線隨機平均法研究分布彈塑性元件為滯遲恢復力模型的滯遲系統的隨機響應.Lin等[20]采用能量包線隨機平均法研究了Bouc-Wen滯遲系統的隨機響應.Ying等[21]采用能量包線隨機平均法研究了Duhem滯遲系統的隨機響應.對具有非局部記憶特性的Preisach滯遲系統,Mayergoyz和Korman[22]研究了隨機輸入下 Preisach系統的平均輸出,Ni等[23]基于方差和非局部記憶滯遲本構模型的切換概率分析,近似得到了Preisach滯遲系統的非線性隨機動力響應的二階統計量.Spanos等[24]和Wang等[25]分別采用標準隨機平均和能量包線隨機平均法研究了 Preisach滯遲系統的隨機響應.最近,Jin等[26]運用隨機平均法和伽遼金法獲得了隨機地震激勵下滯遲系統的近似瞬態響應.

然而,上述方法都存在一些不足,如隨機平均法僅限于能量耗散很小且弱激勵情形;等效線性化法能得到較準確的均方速度和均方位移,但是這種方法只局限于高斯統計的情形,在參激情況下常被認為是不充足或不合適的;非高斯閉合法在尾部通常會獲得負值概率,特別是參激情形.因此,還需進一步開展滯遲非線性系統的隨機響應預測研究.

Er[2729]提出了一種指數多項式閉合法求解穩態FPK方程,但由于采用了高斯概率密度函數構造權函數,因此指數多項式閉合解法的應用范圍受到了很大限制.Di Paola和Sof[30]改進了指數多項式閉合法,提出采用一種簡單有效的迭代過程以提高近似解的精確性,但該方法沒有從根本上解決指數閉合解的局限性問題.最近,文獻[31]提出了一種求解FPK方程穩態解的新方法——迭代加權殘值法,該方法的核心思想是逐步優化權函數,目前已被成功應用于求解多種復雜單自由度強非線性系統[3132].

本文將迭代加權殘值法進一步應用于Bouc-Wen滯遲系統,構造其穩態概率密度函數.首先,將穩態FPK方程的解設為指數多項式,然后借助加權殘值法確定假設解中的待定系數;引入迭代方法,逐步優化權函數,提高加權殘值法的精度,得到系統穩態響應的近似閉合解.作為應用,本文獲得了高斯白噪聲激勵下鋼纖維陶?;炷量蚣芙Y構的近似穩態響應概率密度函數閉合解.研究表明,所獲得的閉合解與Monte Carlo模擬結果吻合得較好.

1 迭代加權殘值法

迭代加權殘值法主要由兩個步驟組成.首先,應用加權殘值法獲得系統的近似穩態概率密度響應;然后,引入迭代過程,逐步優化權函數,提高加權殘值法的精度,最終得到具有較高精度的近似閉合穩態解.

1.1 加權殘值法

考慮外激高斯白噪聲激勵下的 Bouc-Wen系統,其運動方程如下

式中,X,Y與Z分別表示系統位移、速度與滯遲力;ξ為黏滯阻尼率;α∈(0,1)是屈服前后的剛度比;λ,β,γ與n為滯回環參數,β與γ控制滯回環的形狀,λ控制滯回環的幅值,n控制滯回環曲線的光滑性;W(t)是強度為2Dδ(t)的高斯白噪聲.

與系統(1)相應的穩態FPK方程為

式中,p=p(x,y,z).目前,式(2)尚未獲得精確解析解.

構造如下形式的近似解

其中,C0為歸一化常數,φ(x,y,z,cijk)為關于狀態變量的n階多項式,可表示為

其中,cijk為待求系數.ˉp(x,y,z)的存在條件為

其中,x=Rsinθcosφ,y=Rsinθsinφ,z=Rcosθ.將方程(3)代入式(2)中,得殘差

由于ˉp(x,y,z)只是p(x,y,z)的近似值.因此,殘差r(x,y,z,cijk)通常不為零.根據加權殘值法,引入一組權函數 Mijk(x,y,z),使其與殘差r(x,y,z,cijk)的乘積在域內的積分為零,即

權函數可取為[27-29]

其中,pm(x,y,z)是由等效線性化法或隨機平均法等方法獲得的系統概率密度函數.值得注意的是,在數值計算時方程(7)中的積分域通常為有限積分域,但通常難以確定.為此,運用Monte Carlo模擬法粗略估計積分域Ωs.方程(7)可近似表示為

其中,0<i+j+k≤l,l為解的階數.將式(10)進行數值積分,即可得到一組關于cijk的二次非線性代數方程組.采用牛頓法求解該方程組,可獲得系統的近似穩態響應概率密度.

1.2 迭代過程

滯遲系統屬于強非線性系統,僅使用一次加權殘值法可能得不到具有足夠精度的解.因此,采用文獻[30-31]中的迭代技術來提高解的精度.令(k)為迭代k次后得到的近似穩態概率密度函數,取代式(9)中的pm,再應用加權殘值法求出下一個近似穩態概率密度函數ˉp(k+1).重復此步驟直到滿足以下收斂條件

其中,ε0為預設誤差,N1,N2與N3為常數,分別表示狀態空間離散的數目,Δx,Δy和Δz表示離散積分步長.

其中,pR(x,y,z)表示Monte Carlo模擬得到的穩態概率密度函數.

需要指出的是,當pm與真實解相差較大時,上述迭代過程可能不會收斂,特別是當系統為強非線性的情況.文獻[31]提出了一種漸近迭代的方案,即先應用加權殘值法獲得系統弱非線性情形時或弱阻尼弱激勵情形時的近似穩態概率密度函數,然后以此構造權函數,再在非線性參數空間漸近迭代.例如,針對滯遲系統(1),為了獲得激勵強度D=0.2情形時的近似穩態響應概率密度函數,首先采用迭代加權殘值法獲得弱隨機激勵D=0.1情形時的近似穩態概率函數pm1;然后將pm1代入式(9),以此構造新的權函數,再次利用迭代加權殘值法,求得隨機激勵D=0.15情形時的近似穩態概率函數 pm2;最后將pm2代入式(9),構造新的權函數,利用迭代加權殘值法求解D=0.2時的概率密度函數.這種漸近迭代的方法可有效避免選取的初始值 pm與真實解相差較大時導致的不收斂情況,提高了求解的效率.

2 算例

考慮如圖1所示的水平隨機地震作用下單自由度鋼纖維陶粒混凝土框架滯遲系統,其平衡條件為

其中,x(t)是集中質量的相對位移;xt(t)是集中質量的總位移;g(x,)是恢復力.

圖1 地震地面運動下單自由度鋼釬維陶?;炷量蚣芙Y構Fig.1 Steel fibe reinforced ceramsite concrete(SFRCC)single degree of freedom frame subjected to earthquake ground motion

由圖1中可得

式中,xg和分別表示地面的位移和加速度,將方程(15)代入方程(13)可得

其中,F(t)為圖2所示導致水平地面加速度的等效載荷.該框架結構被簡化為如圖2所示的等效系統.

圖2 等效系統Fig.2 Equivalent system

將方程(17)無量綱化,并把等效載荷理想化為獨立的高斯白噪聲.該系統的動力學方程化為

其中,x為無量綱化的位移,ξ為黏滯阻尼率,α是屈服前后的剛度比,k1為常數,W1(t)是強度為2D1的高斯白噪聲,z為恢復力中依賴于時間的滯遲力,由以下公式表示

本文基于鋼纖維陶?;炷林芩降椭芊磸图虞d試驗所得的數據,經最小二乘法辨識獲得滯遲系統的參數:n= 1,λ=1.1037,β=0.202和γ=0.3147.理論曲線和試驗曲線對比如圖3所示,二者吻合較好.系統的其他參數為:ξ=0.025,α=0.2,D=0.1,ε0=0.001.

圖3 鋼纖維陶?;炷林耆苄噪A段試驗曲線和理論曲線(實線表示理論結果,虛線表示試驗結果)Fig.3 Theoretical curve and practical curve of SFRCC column in fully plastic stage(the solid line denotes the theoretical curve and the dotted line denotes the practical curve)

最終可得如下形式的運動方程

本文首先運用等效線性化法獲得了系統的近似穩態響應概率密度函數,表達式如下

然后以式(21)構造權函數,運用加權殘值法獲得了新的近似穩態響應概率密度函數,即

根據誤差公式(11),式(21)與式(22)的誤差分別為0.163和0.249.顯然,式(22)的精度低于式(21).因此,為了進一步提高精度,現以解(22)構造權函數,經過1次迭代,獲得了誤差為0.035的解

圖4給出了關于p(x,z)的穩態邊緣概率密度函數.由圖4和圖5可知,等效線性化的結果與Monte Carlo模擬結果相差甚遠;加權殘值法的計算結果精度尚未令人滿意,但已經表現出系統的一些非線性特征;迭代加權殘值法得到結果與Monte Carlo模擬的結果吻合較好.

另外,以式(23)構造權函數,采用漸進迭代加權殘值法,經過3次迭代獲得了較強隨機激勵D=0.2情形時系統的近似穩態概率密度函數閉合解

圖4 D=0.1情形時關于p(x,z)的穩態邊緣概率密度函數Fig.4 The steady-state marginal probability density function p(x,z)in case of D=0.1

圖5 D=0.1情形時的穩態邊緣概率密度函數p1(x)和p2(z)(?,?表示Monte Carlo模擬數據)Fig.5 The steady-state marginal probability density function p1(x)and p2(z)in case of D=0.1(?,? represent the Monte Carlo simulation data)

式(24)的精度為0.040.圖6給出D=0.2情形時的穩態邊緣概率密度函數.圖6(a)與6(b)表示關于p(x,z)的穩態邊緣概率密度函數.圖6(c)與圖6(d)分別表示關于p1(x)及p2(z)的穩態邊緣概率密度及均方差.符號(?,?)表示樣本數為40000的Monte Carlo模擬結果.由圖6可知,理論解析解與模擬結果在較強隨機激勵情形時仍吻合得較好.

圖6 D=0.2情形時穩態邊緣概率密度函數(?,?表示Monte Carlo模擬結果)Fig.6 The steady-state marginal probability density function in case of D=0.2(?,? represent results from Monte Carlo simulation)

3 結論

本文首次獲得Bouc-Wen滯遲系統的穩態響應概率密度函數的近似閉合解.首先,利用等效線性化法獲得的高斯概率密度函數構造權函數,然后應用加權殘值法獲得系統指數多項式形式的非高斯概率密度函數,最后引入迭代的辦法提高了加權殘值法計算所得結果的精度.對較強隨機激勵情形,以弱隨機激勵情形時的穩態響應概率密度函數構造權函,再在隨機激勵的參數空間內漸近迭代.作為算例,本文研究了隨機激勵下鋼纖維陶?;炷两Y構的穩態響應,其中Bouc-Wen系統的參數是基于擬靜力學試驗由最小二乘法辨識獲得.對比Monte Carlo模擬結果,等效線性化的結果精度較差,加權殘值法的結果雖表現出了系統的非線性特征但精度尚未令人滿意,迭代加權殘值法得到結果與Monte Carlo模擬的結果吻合得較好.對于較強隨機激勵情形,基于漸近迭代加權殘值法具有較高的效率,且所得結果也具有較好的精度.本文所獲得的閉合解,不僅對于土木工程領域具有重要的實際應用價值,還可作為一個標準,用來檢驗其他非線性系統隨機響應預測方法的精度.

1萬征,姚仰平,孟達.復雜加載下混凝土的彈塑性本構模型.力學學報,2016,48(5):1159-1171(Wan Zheng,Yao Yangping,MengDa.An elastoplastic constitutive model of concrete of complicated load.Chinese Journal of Theoretical and Applied Mechanics,2016,48(5):1159-1171(in Chinese))

2李龍彪.纖維增強陶瓷基復合材料疲勞遲滯回線模型研究.力學學報,2014,46(5):710-729(Li Longbiao.Investigation on fatigue hysteresis loops models of fibre-reinforce ceramic-matrix composites.Chinese Journal of Theoretical and Applied Mechanics,2014,46(5):710-729(in Chinese))

3郭洪寶,賈普榮,王波等.基于遲滯行為的2D-SiC/SiC復合材料組份力學性能分析.力學學報,2015,47(2):260-269(Guo Hongbao,Jia Purong,Wang Bo et al.Study on constituent properties of a 2D-SiC/SiC composite by hysteresis measurements.Chinese Journal of Theoretical and Applied Mechanics,2015,47(2):260-269(in Chinese))

4萬強,陳常青,沈亞鵬.鐵電陶瓷PZT53復雜力電耦合行為的實驗研究.力學學報,2005,37(4):413-420(Wan Qiang,Chen Changqing,Shen Yapeng.An experimental investigation into the complex electromechanical behavior of PZT53.Chinese Journal of Theoretical and Applied Mechanics,2005,37(4):413-420(in Chinese))

5 Mayergoyz I.Mathematical models of hysteresis.IEEE Transactions on Magnetic,1986,22(5):603-608

6 Wen YK.Methods of random vibrationfor inelastic structures.Applied Mechanics Review,1989,42(2):39-52

7 Bouc R.Forced vibration of mechanical systems with hysteresis//Proceedings of the Fourth Conference on Non-Linear Oscillation,Prague,Czechoslovakia,1967

8 Caughey T.Random excitation of a system with bilinear hysteresis.ASME Journal of Applied Mechanics,1960,27(4):649-652

9 Suzuki Y,Minai R.Application of stochastic di ff erential equations to seismic reliability analysis of hysteretic structures.Probabilistic Engineering Mechanics,1988,3(1):43-52

10 Jennings PC.Periodic response of a general yielding structure.Journal of Engineering Mechanics-ASCE Division,1964,90(2):131-166

11 Bhatti M,Pister K.A dual criteria approach for optimal-design for earthquake-resistant structural systems.Earthquake Engineering Structure Design,1981,9(6):557-572

12 Beck JL,Jayakumar P.Class of Masing models for plastic hysteresis in structures//Proceedings of 14th ASCE Structures Congress,1996:1083-1090

13 Visintin A.Di ff erential Models of HysteresisⅢ.Springer Science&Business Media,2013

14 Mayergoyz I.Mathematical models of hysteresis.IEEE Transactions on Magnetic,1986,22(5):603-608

15 Iwan WD.A distributed-element model for hysteresis and its steadystate dynamic response.ASME Journal of Applied Mechanics,1966,33:893

16 Roberts JB.The response of an oscillator with bilinear hysteresis to stationary random excitation.ASME Journal of Applied Mechanics,1978,45(4):923-928

17 Roberts JB.Application of averaging methods to randomly excited hysteretic systems//Nonlinear Stochastic Dynamic Engineering Systems,Berlin,Springer,1988

18 Wen YK.Equivalent linearization for hysteretic systems under random excitation.ASME Journal of Applied Mechanics,1980,47(1):150-154

19 Zhu WQ,Lei Y.Stochastic averaging of energy envelope of bilinear hysteretic systems//Nonlinear Stochastic Dynamic Engineering Systems,Berlin:Springer,1988

20 LinYK,Cai GQ.Probabilistic Structural Dynamics:AdvancedTheory and Applications.New York:McGraw-Hill,1995

21 ZG Ying,WQ Zhu,Ni YQ,et al.Stochastic averaging of Duhem hysteretic systems.Journal of Sound and Vibration,2002,254(1):91-104

22 Mayergoyz ID,Korman CE.The Preisach model with stochastic input as a model for aftere ff ect.Journal of Applied Physics,1994,75(10):5478-5480

23 NiYQ,YingZG,KoJM.RandomresponseanalysisofPreisachhysteretic systems with symmetric weight distribution.ASME Journal of Applied Mechanics,2002,69(2):171-178

24 Spanos PD,Cacciola P,Muscolino G.Stochastic averaging of Preisach hysteretic systems.ASCE Journal of Engineering Mechanics,2004,130(11):1257-1267

25 Wang Y,Ying ZQ,Zhu WQ.Stochastic averaging of energy envelope of Preisach hysteretic systems.Journal of Sound and Vibration,2009,321(3-5):976-993

26 Jin XL,Wang Y,Huang ZL.Analysis and control for transient responses of seismic-excited hysteretic structures.Soil Dynamics and Earthquake Engineering,2015,73(6):58-65

27 Er GK.Exponential closure method for some randomly excited non-linear systems.International Journal of Nonlinear Mechanics,2000,35(1):69-78

28 Er GK.A consistent method for the solution to reduced FPK equation in statistical mechanics.Physica A,1999,262(1):118-128

29 Er GK.An improved closure method for analysis of nonlinear stochastic systems.Nonlinear Dynamics,1998,17(3):285-297

30 Di Paola M,Sof A.Approximate solution of the Fokker-Planck-Kolmogorov equation.Probabilistic Engineering Mechanics,2002,17(4):369-384

31 Chen LC,Sun JQ.The closed-form solution of the reduced Fokker-Planck-Kolmogorov equation for nonlinear systems.Communications in Nonlinear Science&Numerical Simulation,2016,41(12):1-10

32 Chen LC,Sun JQ,The closed-form steady-state probability density function of van der pol oscillator under random excitations.Journal of Applied Nonlinear Dynamics,2016,5(4):495-502

THE CLOSED-FORM SOLUTION OF STEADY STATE RESPONSE OF HYSTERETIC SYSTEM UNDER STOCHASTIC EXCITATION1)

Liu Jun?Chen Lincong?,2)Sun Jian-Qiao??(College of civil Engineering,Huaqiao University,Xiamen 361021,Fujian,China)?(School of Engineering University of California Merced,CA 95343,USA)

The hysteretic system is one of the typical strongly nonlinear systems.Hysteretic force depends not only on the instantaneous deformation but also on the past history of deformation.In the last few decades,random vibration of hysteretic system has been studied extensively,but no closed-form solution of random response of hysteretic systems is available so far.In this paper,the newly developed nonlinear random vibration scheme called iterative method of weighted residuals is explored to obtain the closed-form solution of steady-state probability density function(PDF)of the Bouc-Wen hysteretic system under Gaussian white noise excitation.First,a Gaussian PDF is obtained with equivalent linearization technique,which is used as a weighting function.Then,the method of weighted residuals is utilized to determine the non-Gaussian PDF of exponential polynomial type.Finally,an iterative procedure is introduced to improvethe accuracy of the solutions obtained from the method of weighted residuals.As an illustrative example,the steadystate stochastic response of the steel fibe reinforced ceramsite concrete column under random excitation is studied,in which the hysteretic parameters associated with Bouc-Wen hysteretic model is identifie from the pseudo-static test by using the method of least square.Compared to the Monte Carlo results,the accuracy of results obtained from equivalent linearization method is poor.The results obtained from weight residue method can show the nonlinearity of Bouc-Wen systems,but its accuracy is still unsatisfactory.The iterative method of weight residuals can lead to results with higher accuracy.In the case of stronger random excitation,the progressive iterative method of weighted residuals has high efficiency.The obtained solutions agree well with the Monte Carlo simulation data.The proposed closed-form solution of PDF of Bouc-Wen hysteretic system not only is significan to the civil engineering,but also can be a benchmark to examine the accuracy of solutions obtained by other methods.

hysteretic system,steady-state response,closed-form solution,iterative method of weight residuals,steel fibe reinforced ceramsite concrete structure

O324

:A

10.6052/0459-1879-17-003

2017–01–11 收稿,2017–04–11 錄用,2017–04–11 網絡版發表.

1)國家自然科學基金(11172197,11332008,11572215,11672111,51608211)、福建省自然科學基金(2013J05080)、福建省高校杰出青年科研人才培育計劃和華僑大學優秀青年科技創新人才(ZQN-YX307)資助項目.

2)陳林聰,副教授,主要研究方向:非線性隨機振動與控制.E-mail:lincongchen@hqu.edu.cn

劉俊,陳林聰,孫建橋.隨機激勵下滯遲系統的穩態響應閉合解.力學學報,2017,49(3):685-692

Liu Jun,Chen Lincong,Sun Jian-Qiao.The closed-form solution of steady state response of hysteretic system under stochastic excita-tion.Chinese Journal of Theoretical and Applied Mechanics,2017,49(3):685-692

猜你喜歡
系統
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
基于PowerPC+FPGA顯示系統
基于UG的發射箱自動化虛擬裝配系統開發
半沸制皂系統(下)
FAO系統特有功能分析及互聯互通探討
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
一德系統 德行天下
PLC在多段調速系統中的應用
主站蜘蛛池模板: 国产精品视频第一专区| 亚洲最大综合网| 激情无码字幕综合| av大片在线无码免费| 制服丝袜 91视频| 日韩av无码精品专区| 免费一级毛片完整版在线看| 国内精品小视频在线| 亚洲欧美在线综合一区二区三区| 免费在线看黄网址| 性视频久久| jizz国产在线| 91精品小视频| 91偷拍一区| 在线无码九区| 2021国产精品自产拍在线| 婷婷中文在线| 国产午夜无码专区喷水| 天堂网亚洲综合在线| 亚洲中文字幕手机在线第一页| 中文字幕第1页在线播| 黄色免费在线网址| 欧美精品一区在线看| 特级毛片8级毛片免费观看| 国产精品无码AV中文| 成人福利在线看| 欧美日韩国产精品综合| 五月天香蕉视频国产亚| 国产精选小视频在线观看| 日韩黄色大片免费看| 久久黄色影院| 欧美69视频在线| 欧美成人a∨视频免费观看 | 成人午夜福利视频| 国产h视频免费观看| 99热亚洲精品6码| 国产精品亚洲αv天堂无码| 亚洲av无码牛牛影视在线二区| 亚洲欧美一区二区三区麻豆| 久久国产V一级毛多内射| 亚洲视频无码| 91丨九色丨首页在线播放| 色婷婷久久| 免费无码AV片在线观看国产| 国产区在线看| 国产精品欧美日本韩免费一区二区三区不卡| 国产福利拍拍拍| 日韩高清无码免费| 国产大片黄在线观看| 日韩欧美国产另类| 亚洲一区二区三区麻豆| 久久精品人人做人人爽97| 欧美精品在线看| 精品视频福利| 在线va视频| 波多野结衣无码中文字幕在线观看一区二区 | 在线色国产| 2020久久国产综合精品swag| 亚洲欧美综合在线观看| 久久狠狠色噜噜狠狠狠狠97视色 | 国产精品伦视频观看免费| 日韩一区精品视频一区二区| jijzzizz老师出水喷水喷出| 亚洲第一成年网| 99精品视频九九精品| 91免费国产在线观看尤物| 宅男噜噜噜66国产在线观看| 在线精品视频成人网| 亚洲婷婷丁香| 亚洲成人在线免费| 欧美黄色网站在线看| 国产经典在线观看一区| 四虎影视库国产精品一区| 成人午夜天| 视频国产精品丝袜第一页| h视频在线观看网站| 亚洲热线99精品视频| 午夜三级在线| AV在线天堂进入| 欧美一级在线看| 久久综合婷婷| 日本国产一区在线观看|