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

基于近完全時滯補償的隱式實時混合試驗方法

2015-05-09 01:35:24劉進進許國山
振動工程學報 2015年3期
關鍵詞:結構方法

王 貞, 劉進進, 吳 斌, 許國山

(1.哈爾濱工業大學結構工程災變與控制教育部重點實驗室, 黑龍江 哈爾濱 150090; 2.哈爾濱工業大學土木工程學院, 黑龍江 哈爾濱 150090; 3.蘇州設計研究院股份有限公司, 江蘇 蘇州 215021)

基于近完全時滯補償的隱式實時混合試驗方法

王 貞1,2, 劉進進3, 吳 斌1,2, 許國山1,2

(1.哈爾濱工業大學結構工程災變與控制教育部重點實驗室, 黑龍江 哈爾濱 150090; 2.哈爾濱工業大學土木工程學院, 黑龍江 哈爾濱 150090; 3.蘇州設計研究院股份有限公司, 江蘇 蘇州 215021)

實時混合試驗是一種新型結構抗震混合試驗方法。隱式逐步積分算法雖然具有穩定性好的特點,但在實時混合試驗中實施困難。在分析兩種多自由度隱式實時混合試驗方法的基礎上,結合近完全時滯補償方法,提出了一種新型隱式實時混合試驗方法,并分析了該方法的性能。數值模擬表明,該方法具有較高的收斂速度和計算精度,能夠同時考慮時滯補償,能滿足自由度數目較多的多自由度結構實時混合試驗的要求。

實時混合試驗; 隱式逐步積分算法; 時滯補償; 迭代

引 言

傳統的抗震試驗方法有三類,即擬靜力、擬動力和振動臺。近些年,擬動力試驗方法得到了快速發展,出現了子結構化、網絡化和實時化的趨勢,分別產生了子結構擬動力試驗方法、網絡分布式擬動力試驗和實時混合試驗(有時也稱實時子結構試驗或者動力子結構試驗等[1-4])。實時混合試驗方法把原型結構分割為試驗子結構和數值子結構,分別通過真實試驗和實時計算進行模擬,二者之間的相互作用通過液壓伺服加載系統的實時加載保證。試驗子結構往往是具有很強非線性或者力學性能時間相關的構件或子結構,數值子結構是原型結構的其余部分。實時混合試驗由于實時執行加載和數值仿真,能反應試件的速度/加速度相關力學性能,彌補了傳統擬動力試驗方法的不足,近些年得到了快速發展。為了保證實時加載[5-6]和實時計算[7],尤其前者,使得實時混合試驗具有更大挑戰。

實時混合試驗需要在線實時求解結構的運動方程,對逐步積分算法提出了較高要求。常用的顯式算法,如中心差分法,由于其條件穩定性而可能在保證計算實時性時無法滿足對積分步長的限制條件。目前已經提出的基于隱式的實時混合試驗方法,主要有兩種,即吳斌等人提出的等效力控制方法[8-9]和Shing等人提出的實時迭代混合試驗方法[10]。此兩種方法取得了一定的成功,但也存在一定的局限性。本文在分析它們局限性的基礎上,結合作者提出的時滯補償方法[11-12],提出了一種新的隱式實時混合試驗方法,研究表明該方法對于多自由度結構實時混合試驗具有較明顯的優勢。

1 現有隱式實時混合試驗方法

實時混合試驗離散運動方程和基于常加速度假定的位移、速度表達式為:

MNai+1+CNvi+1+RN(di+1)+RE(ai+1,vi+1,di+1)=

Fi+1

(1)

(2)

(3)

式中 角標N表示數值子結構,角標E表示試驗子結構;M為結構質量矩陣;C為結構阻尼矩陣;R為子結構恢復力向量;F為外部荷載向量;d,v,a分別為結構位移、速度、加速度向量;Δt為逐步積分時間間隔。

對式(2)和式(3)進行整理,用di+1和上一步結構響應來表示vi+1和ai+1,并代入結構的運動方程(1),可得

RN(di+1)+KPDdi+1+RE(ai+1,vi+1,di+1)=FEQ,i+1

(4)

其中

(5)

(6)

式(4)是一個關于di+1的非線性方程,在采用牛頓迭代法求解時,有:

(7)

(8)

(1)迭代收斂速度的波動導致作動器運動不平穩,容易對作動器造成損傷且得不到試件的真實力反應;

(2)不易在線獲得結構剛度,矩陣求逆耗時較大。

針對隱式實時混合試驗存在的問題,Shing等人[10]提出了以實時迭代、插值發送命令等為特征的方法,如圖1所示。主要內容包括:

(9)

(3)考慮到固定迭代次數及作動器時滯等因素帶來的誤差,在每積分步最后一個迭代步進行修正。取積分步末點處的位移di+1為

(10)

試驗子結構的恢復力修正為

(11)

Shing方法取得了一定的成功,但仍然存在如下問題:1)通過對迭代位移進行插值發送命令,雖然解決了作動器加載速度突變的問題,但改變了下一次迭代計算的起點,不僅使得迭代過程變得復雜,收斂速度也大幅下降,最后一迭代步可能仍然存在較大誤差。2)在積分步長內每個采樣時刻點迭代一次,意味著必須在一個控制采樣步長內完成迭代計算、位移命令的內插并發送給作動器,當結構自由度較多時,很可能無法滿足實時性的要求。3)即使迭代獲得精度很高的位移,位移命令進行加載時不可避免地涉及系統時滯的影響,引入新的誤差。

圖1 位移命令插值示意圖[10]Fig.1 Schematic for displacement command interpolation[10]

等效力控制方法[8-9]把式(4)等號左側看作是3個并聯元件,采用力反饋控制環取代了數值迭代來求解該方程,其原理如圖2所示。圖中CF稱為力-位移轉化矩陣,功能是把力信號轉化成位移信號,實現作動器的位移控制加載。采用該方法執行實時混合試驗的流程是,首先得到第i+1步的等效力,內插得到積分步內的等效力命令,并在對應的時刻發送給控制回路,作動器受相應的位移命令驅動;在積分步長結束時,把試件的反力反饋到式(4),再次求解位移,以降低控制誤差造成的影響。

可見,等效力控制方法創造性地提出采用控制的思路求解非線性方程,由于引入了控制的思想,使其具有鮮明的個性和突出的特點。分析表明,采用積分控制器并且階躍發送命令的等效力控制方法類似于傳統的迭代方法。合理設計的控制器可以較好地考慮等效力命令插值帶來的影響,而數值迭代則很難實現這一點。此外,由于采用了力反饋控制環,等效力控制方法可以同時考慮到系統作動性能和迭代過程,在一定程度上減小了系統時滯。不過,該方法存在以下缺點:1)內插發送等效力命令,雖然能保證作動器運動的連續性,但會導致系統的控制性能變差。2)雖然合理設計的等效力控制器可以有效減小系統時滯,但是有時仍需要其他時滯補償方法來補償系統時滯。3)積分步內固定計算次數,要求數值子結構靜反力和等效力偏差的求解等過程必須在一個控制采樣步長內完成,尤其是最后一個計算步內,還要完成位移的修正計算。當數值子結構自由度較多時,矩陣運算復雜耗時,很難滿足實時混合試驗實時性的要求。

2 基于近完全時滯補償方法的隱式實時混合試驗方法

縱觀前述兩種隱式實時混合試驗方法,它們的主要缺點是:1)計算耗時:每一計算步復雜耗時的矩陣計算必須在一個控制采樣步長內完成,限制了數值子結構的復雜程度;2)計算精度:每一采樣步的計算以當前所實現的位移為基礎,計算收斂速度有很大程度的下降,在固定次數計算后可能仍然存在較大誤差,該誤差在試驗過程中傳播、累積,對試驗結果可能產生較大影響;3)系統時滯:即使一個積分步長內的計算精度很高,將位移命令發送給作動器進行加載時,不可避免地涉及到時滯,因此會受到時滯補償效果的影響。

鑒于此,提出以近完全時滯補償方法為基礎的隱式實時混合試驗方法。近完全時滯補償方法[11,12]的主要特點是借助過預測結構位移響應,使得試驗子結構的位移超前于數值子結構,只需在已實現的實測位移數據中搜索與期望位移最接近的位移值,同一時刻測得的試件子結構反力即為試驗子結構與期望位移對應的反力。

在進行隱式實時混合試驗非線性方程的迭代求解時,若采用近完全時滯補償方法進行時滯補償,由于試驗子結構的位移響應超前于數值子結構,迭代計算得到位移后,無需發送給作動器進行加載,而是直接從當前完成的實測位移數據中搜索得到與迭代位移近似相等的試驗子結構位移值,同一時刻的試驗子結構實測反力即為與當前迭代位移相對應的反力值,將其直接帶入下一迭代步進行計算即可。圖3是基于近完全時滯補償方法[11-12]的隱式實時混合試驗方法原理示意圖。下面以試驗子結構為彈簧試件的單自由度結構體系為基礎,對照圖3詳細介紹基于近完全時滯補償方法的隱式實時混合試驗方法的流程(后文統一稱之為“新方法”):

圖2 基于等效力控制方法的實時混合試驗原理Fig.2 Principle of real-time hybrid testing based on equivalent force control method

圖3 基于近完全時滯補償方法的隱式實時混合試驗原理Fig.3 Principle of implicit real-time hybrid testing based on nearly-complete delay compensation method

1)位移過預測:ti時刻預測ti+1+τp時刻的位移,如采用顯式Newmark表達式,有dti+1+τp=di+vi(τp+Δt) +0.5ai(τp+Δt)2,τp=τa+Δτ,τa為預估的作動器最大時滯,Δτ為過預測時間跨度;

2)內插生成作動器位移命令:在dti+1+τp與ti時刻的作動器位移命令dti+τp之間,內插得到作動器位移命令dc(ti~ti+1),并發送給作動器。

3)采集數據:采集試驗子結構位移響應de(ti~ti+1)及試件反力RE(ti~ti+1)。

4)迭代求解非線性方程:

①每積分步第一次迭代以上一步的位移響應為初值,其余次迭代以上次迭代結果為初值;

④進行多次迭代,直到收斂;更新位移預測所需的位移di+1、速度vi+1、加速度ai+1;當試驗時間達到ti+1,開始下一個積分步的計算。

可見,基于近完全時滯補償方法的隱式實時混合試驗方法具有如下特點:1)采用簡化牛頓迭代,僅需試驗前計算逆矩陣;2)第k+1迭代步的計算完全基于第k迭代步的位移及其響應進行,實現了標準的簡化牛頓迭代,其迭代收斂速度要優于Shing方法和等效力控制方法,計算效率更高;3)積分步長內無固定迭代次數,位移迭代達到收斂精度即可,而收斂迭代次數預期少于前兩種方法;4)在積分步長內作動器的加載不依賴于當前迭代位移值,只要在一個積分步長內完成迭代收斂計算即可,對單步計算時間的限制比前兩種方法寬松,能適應自由度數目更多的結構;5)作動器的加載基于預測位移內插生成位移命令進行,不會出現加載速度失真現象,并能保證試件加載的實時性;6)在迭代求解非線性方程的同時已經考慮了作動器加載過程中的時滯,并進行了時滯補償;7)該方法的不足是,采取了位移過預測方法,預測的時間跨度更大,帶來了更大的預測誤差。該誤差可以通過預測精度更好的位移預測方法并通過回搜過程來減小。此外,還可通過其他方法先降低系統時滯,以減少預測跨度過大帶來的預測誤差。

3 數值模擬

3.1 單積分步計算精度

數值分析表明,對于不考慮時滯的小步長分析,三種方法都具有較好的精度。考慮到多自由度分析中存在高頻成分、實時計算要求采用較大步長,采取如下參數:采樣步長δt=0.004 s,積分步長Δt=0.04 s;采用積分等效力控制器,控制參數KI=400;結構參數為:KE=2×103N/m,MN=1 kg,KN=2×103N/m。三種方法單積分步位移計算結果如圖4所示。可見,在積分步長較大或結構頻率較高時,即使在最后一個計算步對積分點位移進行了修正,Shing方法和EFCM方法的收斂精度也不是很理想。

圖4 單積分步計算精度比較Fig.4 Computational accuracy comparisons in single integration step

時滯補償是實時混合試驗最關鍵的問題,鑒于時滯可能發生波動,采用常規補償方法很難實現實測的完全補償,接下來的分析中假定存在約1 ms的欠補償時滯。作動器模型取為一階傳遞函數TA(s)=700/(700+s),對于頻率為1 Hz的信號,該系統時滯約為1 ms。取積分步長Δt=0.01 s,不同方法的單積分步計算精度如圖5所示。可見,時滯得不到完全補償時,Shing方法和EFCM的單積分步計算精度低于新方法。

由前述分析可知,數值子結構自由度較多或結構頻率較高時,Shing方法和EFCM的計算精度會有所下降;實際試驗的欠補償時滯也會使得兩種方法計算精度下降。即使進行了位移修正,計算精度也不如所提出的基于簡化牛頓迭代的新方法。

圖5 時滯欠補償時單積分步位移計算精度比較Fig.5 Computational accuracy comparisons in single integration step with under-compensated delay

3.2 計算耗時

本小節比較應用于多自由度結構實時混合試驗時Shing方法、EFCM和新方法的計算耗時情況。取如圖6所示的5層剪切型結構為分析對象,層間剛度為ki=40×103N/m,集中質量為mi=100 kg(i=1~5),結構自振頻率為0.91,2.64,4.17,5.36,6.11 Hz,結構無阻尼。選底層柱為試驗子結構,剛度取kE=k1/2。采用理想作動器模型進行數值模擬,積分步長取Δt=0.01 s,每個積分步計算10次。新方法的計算精度取10-7m。計算平臺參數如表1所示。

圖6 5自由度結構計算示意圖Fig.6 Computational schematic of 5DOF structure

單積分步內的各迭代步耗時如圖7所示。可見,3種方法單迭代步的耗時并無太大差別,均約為1×10-5s。3種方法的第一迭代步耗時較其余步更大。EFCM在最后一個迭代步要完成位移命令計算和位移修正計算,計算耗時增加。新方法由于收斂速度較快,只進行3次迭代計算即達到精度要求,完成了該積分步內方程求解。在實時混合試驗中,Shing方法和EFCM方法要求采樣步長大于單步迭代最大耗時,新方法的迭代次數遠小于它們,因此新方法大大提高了計算效率,更適用于多自由度結構實時混合試驗。進一步的時程耗時分析見文[13]。需要說明的是,計算耗時統計是在性能較高的計算機上完成,各計算步耗時均較小。實際試驗中,控制板不僅要完成數值子結構的計算、時滯補償計算,還要完成液壓伺服作動器的控制計算,多任務的同步實時處理,使得所能完成的實時混合試驗的自由度數目大大受限。

圖7 單積分步計算耗時Fig.7 Time consumed for computation in single integration step

表1 計算平臺參數

3.3 時程分析

為了進一步檢驗3種方法的性能,分別完成了無阻尼自由振動、無阻尼地震激勵、非線性自由振動、非線性地震激勵等多種工況的數值模型。限于篇幅,此處僅介紹最后一種工況的數值模擬。

分析模型為圖6所示的5層剪切型結構,假設試驗子結構具有雙線性恢復力關系,起始剛度為kE= 20×103N/m,屈服位移為0.008 m,屈服后剛度為起始剛度的0.115倍。首層自由度初始位移為0,初始速度為0;其他自由度的初始狀態為0。外荷載激勵取El Centro(NS,1940)地震記錄,峰值加速度調為0.2g。

考慮未能完全補償時滯的影響,即假設欠補償時滯約為1 ms,取作動器模型為TA(s)=700/(700+s)。等效力控制器參數不變。首層自由度位移時程及恢復力-位移關系見圖8和9所示。可見新方法的計算精度要好于其他兩種方法。

圖8 位移時程Fig.8 Displacement time histories

圖9 試件恢復力-位移曲線Fig.9 Restoring force-displacement relation of specimen

4 結 論

分析了Shing方法和等效力控制方法(EFCM)在實時混合試驗中可能存在的困難,提出了基于近完全時滯補償方法的隱式實時混合試驗方法,詳細闡述了其原理和實現流程,并從理論上分析了其特點。最后,通過數值模擬比較了三種方法在收斂速度和計算精度等方面的性能,驗證了理論分析結果。主要結論如下:

1)Shing方法和EFCM應用于多自由度結構實時混合試驗時,可能會由于自由度數目的增多,使得運算耗時增加,難以滿足實時性的要求;

2)當積分步長增大或試驗子結構的初始剛度較大,以及試驗中可能存在時滯欠補償情況時,Shing方法和EFCM的計算精度會有所下降;

3)基于近完全時滯補償方法的隱式實時混合試驗方法,具有較高的計算精度和收斂速度,能夠同時考慮時滯補償,能滿足較多自由度數目的多自由度結構實時混合試驗的要求。

[1] Nakashima M, Kato H, Takaoka E. Development of real-time pseudo-dynamic testing[J]. Earthquake Engineering & Structure Dynamics, 1992,21(1):79—92.

[2] 吳斌,王倩穎.實時子結構實驗的研究進展[J].實驗力學,2007,22(5):1—9.

Wu Bin, Wang Qianying. Development of real-time substructure testing[J]. Journal of Experimental Mechanics, 2007,22(5):1—9.

[3] 遲福東,王進廷,金峰,等.土-結構動力相互作用的實時耦聯動力試驗的時滯穩定性[J].工程力學,2012,29(8):1—7.

Chi Fudong, Wang Jinting, Jin Feng, et al. Delay-dependent stability of real-time dynamic hybrid testing for soil-structure interaction analysis[J]. Engineering Mechanics, 2012,29(8):1—7.

[4] 蔡新江,田石柱.振動臺試驗方法的研究進展[J].結構工程師,2011,27(增刊):42—46.

Cai Xinjiang, Tian Shizhu. Research advances of shaking tabel testing method[J]. Structural Engineers, 2011,27(Sup):42—46.

[5] Horiuchi T, Inoue M, Konno T, et al. Namita. real-time hybrid experimental system with actuator delay compensation and its application to a piping system with energy absorber[J]. Earthquake Engineering & Structural Dynamics, 1999,28(10):1 121—1 141.

[6] 王貞,吳斌.基于最小二乘法的時滯實時在線估計方法[J].振動工程學報,2009,22(6):625—631.

Wang Zhen, Wu Bin. A real-time approach to delay estimation based on the Least-Square algorithm[J]. Journal of Vibration Engineering, 2009,22(6):625—631.

[7] 賈傳果,李英民,劉立平,等.Rosenbrock耦合積分方法及其收斂性分析[J].重慶大學學報,2013,36(2):127—133.

Jia Chuanguo, Li Yingmin, Liu Liping, et al. Rosenbrock-based coupled time integration method and its convergence analysis[J]. Journal of Chongqiong University, 2013,36(2):127—133.

[8] Wu B, Wang Q, Shing P B, et al. Equivalent force control method for generalized real-time substructure testing with implicit integration[J]. Earthquake Engineering & Structural Dynamics, 2007,36:1 127—1 149.

[9] 許國山,吳斌.采用等效力控制方法的非線性結構實時子結構試驗[J].振動工程學報,2010,23(3):237—242.

Xu Guoshan, Wu Bin. Real-time substructure testing of nonlinear structures with equivalent force control method[J]. Journal of Vibration Engineering, 2010,23(3):237—242.

[10]Jung R, Shing P, Stauffer E, et al. Performance of a real-time pseudo-dynamic test system considering nonlinear structural response[J]. Earthquake Engineering & Structural Dynamics, 2007,36(12):1 785—1 809.

[11]Bin Wu, Zhen Wang, Oreste S. Bursi. Actuator Dynamics Compensation Based on Upper Bound Delay for Real-time Hybrid Simulation[J]. Earthquake Engineering & Structural Dynamics, 2013,42(12):1 749—1 765.

[12]王貞,劉進進,吳斌.實時混合試驗近完全時滯補償方法的參數確定[J]. 工程力學,2014,31(10):158—166.

Wang Zhen, Liu Jinjin, Wu Bin. Parameter determination of Nearly-complete compensation scheme for time delay in real-time hybrid simulation[J]. Engineering Mechanics, 2014,31(10):158—166.

[13]劉進進.近完全時滯補償方法及其在實時混合試驗中的應用[D].哈爾濱:哈爾濱工業大學, 2013.

Liu Jinjin. Nearly complete compensation scheme for time delay and its application to real-time hybrid testing[D]. Harbin: Harbin Institute of Technology, 2013.

Implicit method for real-time hybrid testing based on nearly-complete delay compensation scheme

WANGZhen1,2,LIUJin-jin3,WUBin1,2,XUGuo-shan1,2

(1.Key Lab of Structures Dynamic Behavior & Control (Harbin Institute of Technology), Ministry of Education,Harbin 150090, China; 2.School of Civil Engineering, Harbin Institute of Technology, Harbin 150090, China;3.Suzhou Institute of Architectural Design Co., LTD., Suzhou 215021, China)

Real-time Hybrid Testing (RHT) is a novel hybrid testing method for evaluating structural seismic performance. Even though implicit time-stepping approaches exhibit better stability, their applications to RHT are limited due to their difficulties. This paper firstly analyses limitations of two existing implicit RHT schemes for Multi-Degree of Freedom (MDoF) structures, and then proposes a new implicit RHT based on the nearly-complete delay compensation developed by the authors formerly. Theoretical performance analysis and numerical simulations are performed following that. Numerical simulations show that the proposed scheme is endowed with better convergence speed and accuracy and incorporates delay compensation, which indicate that it is more suitable to RHT for MDoF systems.

real-time hybrid testing; implicit time-stepping method; delay compensation; iteration

2013-12-01;

2014-07-14

國家自然科學基金國際合作項目(51161120360)和重點項目(91315301-09);中國博士后科學基金面上資助項目(2013M531046);中央高校基本科研業務費專項資金資助項目(HIT.NSRIF.2014102);哈爾濱工業大學海外博士科研啟動項目

TU317; P315.9

A

1004-4523(2015)03-0418-07

10.16385/j.cnki.issn.1004-4523.2015.03.011

王貞(1983—), 男, 博士,講師。電話: 15004548396; E-mail: zhenwang@hit.edu.cn

猜你喜歡
結構方法
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
學習方法
論《日出》的結構
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 亚洲国产精品人久久电影| 国产成人高清在线精品| 国产成人8x视频一区二区| 精品久久香蕉国产线看观看gif| 亚洲福利一区二区三区| 国产精品成人免费综合| 成人国产精品一级毛片天堂| 国产aⅴ无码专区亚洲av综合网| 亚洲欧美极品| 亚洲av日韩综合一区尤物| 国产真实乱子伦精品视手机观看| 尤物午夜福利视频| 国产成人亚洲无码淙合青草| 尤物在线观看乱码| 亚洲一区波多野结衣二区三区| 久久毛片网| 国产第一页免费浮力影院| 久久黄色视频影| 亚洲色无码专线精品观看| 免费av一区二区三区在线| 国产无码精品在线播放| 91欧洲国产日韩在线人成| 亚洲精品第一页不卡| 国产精品视频第一专区| 精品久久久久无码| 国产亚洲视频中文字幕视频 | 无码有码中文字幕| 国产黄在线免费观看| 国产成人乱无码视频| 亚洲国产系列| 国产精品网曝门免费视频| 欧美成人一区午夜福利在线| 国产成人高清亚洲一区久久| 日韩毛片免费| 亚洲综合第一区| 五月婷婷精品| 在线播放91| 天天躁夜夜躁狠狠躁躁88| 亚洲中文制服丝袜欧美精品| 性色生活片在线观看| 久青草免费在线视频| 丰满人妻久久中文字幕| а∨天堂一区中文字幕| 亚洲全网成人资源在线观看| 亚洲综合色区在线播放2019| 国产欧美一区二区三区视频在线观看| 色噜噜狠狠色综合网图区| 在线人成精品免费视频| 欧美一区二区人人喊爽| 国产欧美视频综合二区| 国产经典免费播放视频| 欧美特黄一级大黄录像| 亚洲天堂日本| 日本午夜视频在线观看| 本亚洲精品网站| 色婷婷天天综合在线| 欧洲日本亚洲中文字幕| 99久久精品无码专区免费| 精品超清无码视频在线观看| 五月六月伊人狠狠丁香网| 色九九视频| 成人亚洲天堂| 日本三区视频| 亚洲综合中文字幕国产精品欧美| 色婷婷在线播放| 久久精品娱乐亚洲领先| 日韩小视频在线观看| 99热这里只有免费国产精品 | 免费A级毛片无码无遮挡| 尤物精品国产福利网站| 免费在线播放毛片| 中日韩欧亚无码视频| 亚洲日本韩在线观看| 亚洲 欧美 偷自乱 图片| 凹凸国产分类在线观看| 亚洲精品中文字幕无乱码| 四虎影视国产精品| 亚洲中文字幕日产无码2021| 91视频99| 亚洲欧洲日产国产无码AV| 99久久国产自偷自偷免费一区| 日本成人一区|