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

空間相關多點地震位移輸入峰值控制對結構非線性響應的影響

2022-01-12 13:58:50陳科旭俞瑞芳孫平寬
振動工程學報 2021年6期
關鍵詞:結構

陳科旭,俞瑞芳,孫平寬

(1.中國地震局地球物理研究所,北京100081;2.中國公路工程咨詢集團有限公司,北京100089)

引言

合理的地震動輸入是結構進行抗震分析的前提條件[1]。地震動受傳播路徑、距離及場地條件等因素的影響,在時間和空間上都具有復雜的變化,對于平面尺寸較小的建筑物,地震動在空間不同觀測點的差異變化影響通常可以忽略不計,在進行結構抗震分析時,各個支撐點處可以采用同樣的地震動輸入;但是對于大尺度空間結構,如長大橋梁、隧道、渡槽,以及管線、輸電塔等生命線工程,地震動的空間變化將對其產生重要影響[2],因此大尺度空間結構采用多點地震動輸入更加符合實際的輸入模式[3-4]。結構抗震設計理論從一致地震動輸入發展到多點地震動輸入是一個巨大的飛躍。已有研究表明,與多點激勵下的結構響應相比,一致激勵往往明顯高估或低估結構的某些響應[5-8],這些研究結果的差異除了和結構特性相關外,很大程度上在于對地震動輸入的空間變化或非平穩特性描述不同[8]。基于工程實踐的需要,能夠表征地震動空間變化的相干函數模型及空間相關多點地震動擬合方法得到了很大的發展[9-13]。目前用于工程實踐的相關多點地震動合成方法都是基于Hao等提出的三角級數法進行的[14],該方法假定地震動是一個平穩的隨機過程,用相干函數來描述地震動的空間變化。由于相干函數是基于平穩隨機過程假定給出的,其無法描述每個測點地震動強度和頻率隨時間的變化,因此工程中常采用具有統計參數的強度包線函數來反映地震動強度的非平穩特性。對頻率的非平穩特性,可采用分段合成再疊加的方法,或通過地震加速度相位差譜分布來近似描述[15],或直接采用時-頻包線函數來近似模擬地震動頻率非平穩特性[16]。

目前,在通用有限元分析軟件中,對大跨結構采用多點輸入模式進行地震響應分析時,一般采用位移輸入法和加速度輸入法,其中位移輸入法是直接在結構各基底輸入地震位移時程,加速度輸入法則是通過特殊的建模方式使結構基底的輸入等效為各點的加速度時程[17],其實質是在上部結構上施加等效的地震荷載。現有的空間相關多點地震動擬合方法中,一般都是以加速度反應譜和加速度峰值為控制目標擬合得到結構不同支撐點的加速度時程,此時相應于不同支撐點處位移時程的峰值則是不同的;再者,由于地震動擬合中初始加速度時程采用隨機相位,即使控制加速度時程的峰值相同,其相應的位移時程的峰值也有較大的離散性[18]。由于在場址地震動參數的估計中,目前的研究成果很難得到可應用于工程實踐的位移峰值的估計方法[19],因此采用加速度峰值作為擬合目標進行地震動模擬是工程實踐中必然的選擇。此時,面臨的問題是采用位移輸入法進行結構非線性響應分析時,是否要控制位移時程的峰值?如何控制位移峰值?不同的控制條件會對結構響應產生怎樣的影響?為了回答這些問題,本文首先基于不同的研究目標設計了地震動輸入方案,模擬得到了能夠表征地震動強度或頻率特性的空間相關多點地震動,并且控制不同支撐點有相同的加速度峰值或位移峰值;然后采用多點地震動加速度/位移輸入,對三跨連續梁橋進行了非線性時程響應分析,討論不同的地震動輸入特性和方案對結構響應的影響。本文的分析結果可為空間相關多點地震動人工擬合中關鍵控制參數的設置、地震動輸入方案的選擇提供理論依據。

1 大跨結構在多點地震加速度/位移輸入下響應的計算方法

大跨空間結構在進行地震響應分析時,如果假定結構為集中質量體系,在結構基底各支撐點處施加不同的地震動輸入,采用體系外的固定坐標系,用總位移表示各節點的運動,則結構的動力方程可以寫為[20]

式中M,C,K為質量、阻尼和剛度矩陣,下標ss,bb,sb分別表示上部結構、支撐點和二者耦合項;分別表示上部結構的加速度、速度和位移響應分別表示下部支撐點處的地震加速度、速度和位移時程,一般為已知量;Rb為下部支撐節點上的反力。

展開式(1)中的第一行可得

可以看出,如果已知地震位移激勵Ub及相應的速度U?b項,則求解方程(2)可得到上部結構的絕對位移響應Us。方程(2)是目前通用結構有限元分析軟件中位移輸入法(DM)[21]的基本方程,計算時首先釋放支承處加載方向自由度,然后直接輸入地震位移時程。相較于理論上的相對位移法,直接位移輸入法在計算時不需預先提取質量、剛度、阻尼矩陣,這將大大簡化前處理流程。

如果展開式(1)的第二行,則可得到

式中下部支撐節點上的反力Rb可以表示為

然而,基于方程(5)的LMM法計算得到的一般是結構響應的近似值,和位移輸入法得到的精確值相比,LMM法計算得到的結構響應的某些量會產生較大的誤差。如圖1所示為一個三跨橋梁的橋墩扭矩和滑動支座位移,采用LMM法計算,其結構產生的最大誤差達到了66.08%。現分析該誤差產生的原因,若計算中結構阻尼假定是Rayleigh阻尼,即

代入方程(5)可得

將式(6)和(8)代入式(5)中,則可以得到地震加速度U?g和支座處的實際響應加速度U?b之間更加精確的近似,文獻[23]討論了這種修正方法,并將這種修正了擬輸入加速度的方法稱為修正后大質量法(MLMM)。為了說明修正后大質量法的計算效果,將上述的3跨橋梁按照MLMM法計算得到的結果與位移輸入法的計算結果進行了對比,結果如圖1所示。可以看出,同樣的輸出量,采用MLMM法的誤差大約為1%-2%。因此,在本文的研究中,將采用MLMM法和DM法進行多點加速度/位移輸入下大跨結構的非線性響應分析。

圖1 不同輸入方法的結果對比Fig.1 Comparison of results of different input methods

2 空間相關多點地震動模擬及輸入方案

為了研究多點地震動輸入下,控制不同支撐點的加速度/位移峰值對結構非線性響應的影響,本節將設計大跨結構的地震動輸入方案,并模擬滿足預定控制目標的空間相關多點地震動。

2.1 結構模型

本文設計的橋梁模型如圖2所示,該模型為一座全長120 m的3跨連續梁混凝土直線橋(跨度為30 m+60 m+30 m)。橋的主梁為預應力箱型梁,橋墩為矩形獨柱式混凝土墩,墩高均為10 m,墩頂分別布置固定盆式支座、單/雙向活動盆式支座,主梁兩端為連接路堤的橋臺。主梁采用C50混凝土,橋墩采用C40混凝土。采用通用有限元軟件CSIBridge進行橋梁建模與分析,主梁、橋墩均采用框架單元,墩底固接;雙向活動支座和單向活動支座的活動方向采用橡膠隔震單元,固定支座采用線性彈簧單元,釋放其他方向自由度。在有限元模型中,主梁屬性設置為彈性;橋墩塑性鉸的恢復力模型選用Takeda模型;活動支座的恢復力模型選用雙線性模型,基本參數可由《公路橋梁抗震設計細則》[24]計算得到,如屈服強度為18.12 kN,屈服位移為0.005 m,屈服后剛度比設置為0。計算模型采用Rayleigh阻尼,阻尼比為5%,其中質量阻尼系數α=0.7778,剛度阻尼系數β=2.187×10-3;模型前三階自振周期分別為0.71,0.50和0.47 s。

圖2 橋梁模型(單位:m)Fig.2 Bridge model(Unit:m)

2.2 地震動輸入方案及人工模擬

依據圖2所示的橋梁模型,進行0#橋臺(支點A)、1#橋墩(支點B)、2#橋墩(支點C)及3#橋臺(支點D)橫向的地震動擬合。為了綜合考慮空間相關多點地震動的強度(頻率)非平穩特性、不同峰值加速度/位移控制條件對橋梁結構非線性響應的影響,本文設計了兩類(六組)多點地震動輸入方案,如表1所示。

第一類地震動模擬時,引入具有統計參數的時-頻包線函數[16],同時考慮加速度時程強度和頻率非平穩特性,并分為三個方案,其中方案1-1控制4個支撐點的加速度峰值相同(相應的位移時程峰值隨機),方案1-2按照4個支撐點位移時程峰值的最大值控制4條位移時程的峰值,方案1-3按照4個支撐點位移時程峰值的最小值控制4條位移時程的峰值。

第二類地震動模擬時,采用強度包線函數[25],僅考慮地震動強度非平穩特性,同樣按照控制不同支撐點的加速度峰值相同或位移峰值相同設計了三個模擬方案,即方案2-1、方案2-2和方案2-3。

現按照表1所示的地震動模擬方案進行橋梁4個支撐點的地震動模擬。地震動模擬時,目標反應譜采用《建筑抗震設計規范》(GB 50011-2010)中的反應譜[26],其中地震動加速度峰值為0.102g,特征周期為0.4 s,曲線下降段的衰減指數為0.9,結構的阻尼比為0.05。地震動擬合時,取80個頻率控制點確定目標反應譜,且對目標譜的允許擬合誤差為5%。初始種子時程選取美國San Simeon地震P06臺站的東西向(EW)地震加速度記錄,并采用了基于San Simon地震臺陣記錄得到的相干模型[27]。

表1 多點地震動模擬方案Tab.1 Simulation schemes of multi-point ground motions

(1)第一類地震動模擬:表征地震動強度和頻率非平穩特性。

方案1-1:根據種子時程的實際場地條件,選取主頻率參數構建具有統計意義的時-頻包線函數[28],則對應于采樣頻率fk處的時頻聯合分布函數為

式中主頻率Fp(t)定義為一系列采樣時間點t1,t2,…,tn對應于時頻譜的最大幅值的頻率值;強度包線采用三段式包線函數E(t),如下式所示

按照種子時程的5%到75%的Arais強度定義三段式包線函數,即t1=5 s,t2=13 s,c=0.347。則基于主頻率構造的時-頻包線函數可由下式表示

按照文獻[28]建立的空間相關多點地震動擬合方法,擬合得到匹配反應譜、加速度峰值的4個支撐點地震加速度時程如圖3所示。4條加速度時程的峰值均為0.102g,對目標反應譜的擬合誤差均小于5%,且對地震動空間變化模擬效果較好。

圖3 空間相關多點加速度時程(方案1-1)Fig.3 Spatial correlation multi-point acceleration time history(Scheme 1-1)

方案1-2:對方案1-1擬合得到的4條加速度時程進行兩次積分得到相應的4條位移時程,其位移峰值分別為0.104,0.135,0.0854和0.107 m,以位移峰值的最大值(0.135 m)作為位移峰值目標對4條位移時程進行調整,最終得到4條峰值均為0.135 m的位移時程,如圖4所示。該組時程的其他控制條件,如反應譜、頻率和強度非平穩特性、空間相關性等,均與方案1-1相同。

圖4 空間相關多點位移時程(方案1-2)Fig.4 Spatial correlation multi-point displacement time history(Scheme 1-2)

方案1-3:以位移峰值的最小值(0.0854 m)作為位移峰值目標對4條位移時程進行調整,得到4條位移峰值均為0.0854 m的位移時程,如圖5所示。同樣,該組時程的其他控制條件與方案1-1相同。

圖5 空間相關多點位移時程(方案1-3)Fig.5 Spatial coherention multi-point displacement time history(Scheme 1-3)

(2)第二類地震動模擬:僅考慮了地震動強度非平穩特性。

方案2-1:采用三段式強度包絡函數E(t)在時域內進行初始地震動調整,強度包線參數同方案1-1。按照預定的目標反應譜、加速度峰值、空間相關參數等,擬合得到4個支點的地震加速度時程,如圖6所示。模擬得到的4條加速度時程的峰值均為0.102g,對目標反應譜的擬合誤差均小于5%,且對地震動空間變化模擬效果較好。

圖6 空間相關多點加速度時程(方案2-1)Fig.6 Spatial coherention multi-point acceleration time history(Scheme 2-1)

方案2-2:對方案2-1擬合得到的4條加速度時程進行兩次積分得到相應的4條位移時程,其位移峰值分別為0.0732,0.106,0.0776和0.0955 m,以位移峰值的最大值(0.106 m)作為位移峰值目標對4條位移時程進行調整,最終得到4條位移峰值均為0.1060 m的位移時程,如圖7所示。該組時程的其他控制條件,如反應譜、強度非平穩特性、空間相關性等,均與方案2-1相同。

圖7 空間相關多點位移時程(方案2-2)Fig.7 Spatial coherention multi-point displacement time history(Scheme 2-2)

方案2-3:以位移峰值的最小值(0.0732 m)作為位移峰值目標,對4條位移時程進行調整,最終得到4條位移峰值均為0.0732 m的位移時程,如圖8所示。該組時程的其他控制條件(反應譜、強度非平穩特性、空間相關性)與方案2-1相同。

圖8 空間相關多點位移時程(方案2-3)Fig.8 Spatial coherention multi-point displacement time history(Scheme 2-3)

3 關鍵擬合參數對結構非線性響應的影響

將擬合得到的6組空間相關多點地震動,采用多點輸入的方式對圖2所示的結構進行非線性反應分析,可得到橋墩和主梁的內力、位移等響應參數。表2給出了6組輸入方案計算得到的橋墩的剪力和扭矩;表3為3跨主梁的跨中彎矩;表4為橋墩和主梁的橫向位移。本節將分析采用不同輸入方案時結構響應的變化,來說明地震動特性及不同加速度/位移峰值控制條件對結構響應的影響。

表3 不同輸入方案主梁跨中彎矩對比/(kN·m)Tab.3 Comparison of mid-span bending moment of main girders with different input schemes/(kN·m)

表4 不同輸入方案橋梁橫向位移/mTab.4 Lateral displacements of bridge with different input schemes/m

橋墩屈服時對應的等效屈服剪力為1560 kN,超過這個值橋墩將進入非線性,表2中的數據表明,2#橋墩在方案1-2和方案2-2中均已率先屈服,進入了非線性,恢復力曲線如圖9所示。表明控制位移峰值的大小將會影響橋梁結構進入非線性程度的強弱,當按照最大位移峰值進行控制時,橋墩率先進入了非線性,滯回環的耗能面積開始增加,相應的橋墩內力及位移響應量增大。圖10給出了3#橋臺的支座在第一類輸入方案下的恢復力曲線,可以看出在3種輸入方案下支座均已進入非線性,并且支座在方案1-2中的滯回耗能面積和位移峰值要明顯大于方案1-1和方案1-3的結果;圖11進一步給出了兩類不同輸入方案下支座滯回耗能面積的對比,可以看出,當按照最大位移峰值進行控制時(方案1-2和方案2-2),支座的滯回耗能面積更大,因而結構進入的非線性程度更深,相應的橋梁橫向位移也將增大。

圖9 2#橋墩的恢復力曲線Fig.9 Restoring force curves of 2# pier

圖10 不同方案下3#橋臺支座恢復力曲線對比Fig.10 Comparison of restoring force curves of bearings in 3# abutment under different schemes

圖11 3#橋臺支座的非線性耗能對比Fig.11 Comparison of nonlinear energy dissipation of bearings in 3# abutment

表2 不同輸入方案橋墩的剪力和扭矩對比Tab.2 Comparison of shear force and torque of bridge piers with different input schemes

因此以下將分別對橋梁的內力和位移響應情況進行討論,分析空間相關多點地震動的位移峰值控制對結構非線性響應的影響。

3.1 控制不同的位移峰值對結構響應的影響

3.1.1 對內力的影響

基于表2和表3中的結果,首先計算了結構在方案1-2和方案1-3輸入下的內力響應相對于方案1-1輸入下結構響應的相對誤差,其結果如圖12(a)所示。分析誤差變化可以看出,當按照最大的位移峰值進行位移時程的控制時(方案1-2),橋墩的剪力和扭矩、主梁跨中彎矩大都會得到相對較大的結果,相應于各參數的最大增幅分別約為25%,21%和50%;若按照最小的位移峰值進行位移時程的控制時(方案1-3),無論是橋墩的剪力和扭矩,還是主梁的跨中彎矩,與方案1-1的計算結果相比,減小的幅度較大,相對于各參數最大降低幅度分別約為37%,26%和28%。在方案2-2和方案2-3輸入下的結構響應相對于方案2-1輸入下計算結果的相對誤差如圖12(b)所示。同樣地,當按照最大位移進行控制時(方案2-2),也得到相對較大的結構響應,各參數的最大增幅分別為17%,25%和38%;若按照最小的峰值位移進行控制,計算得到的橋墩剪力和扭矩、跨中主梁彎矩的最大減小幅度分別為30%,18%和20%。

圖12 不同地震動輸入方案橋梁內力變化Fig.12 Internal force changes of bridge in different ground motion input schemes

由以上分析可以看出,按照位移輸入法進行結構分析時,控制位移峰值的大小對結構內力響應有較大的影響,但是這種影響似乎與位移時程的形狀無明顯的相關性,如圖13所示的主梁跨中彎矩的變化來看,雖然方案1-2和方案2-2的位移時程來自不同的地震動特性分組,且形狀不同,但由于位移峰值(PGD,Peak Ground Motion Displacement)較大,都得到了相對較大的計算結果,而且大部分的內力結果都是隨著峰值位移的減小逐漸減小。分析表2中的結構響應,這種變化也同樣出現在橋墩的扭矩和剪力的結果中,進一步說明位移輸入法中各支撐點處地震動的位移峰值與結構的內力響應之間呈現正相關。由于在以加速度反應譜和加速度峰值為目標的地震動擬合中,不同支撐點處位移的峰值離散性比較大,所以為了保證結構的安全,地震動擬合中應控制位移峰值的最小取值。

圖13 地震位移輸入方案下主梁跨中彎矩情況Fig.13 The mid-span bending moment of the main girders under the seismic displacement input schemes

3.1.2 對橫向位移的影響

分析表4的數據,發現采用加速度輸入(方案1-1)和位移輸入(方案1-2)時,計算得到的2#墩頂的位移響應出現了不同的變化規律,即當采用加速度輸入時,其位移響應比1#墩頂位移減少約35%,但當采用位移輸入時,兩個墩頂得到位移響應相當。產生這種情況的原因可能有兩個,其一是加速度輸入時對應于2#橋墩的位移峰值是4個支點中的最小值,故方案1-2中對其調整的幅度最大;其二是2#墩頂設置了縱向滑動支座,相對于1#墩頂的固定支座,滑動支座可能對位移輸入的峰值變化更加敏感。第二類方案中的結果也與之類似。

圖14給出了采用不同地震動輸入方案時,橋梁關鍵位置橫向位移相對于方案1-1和方案2-1計算結果的相對誤差變化。從圖14(a)可以看出,當控制位移輸入峰值為最大值時(方案1-2),與方案1-1的結果相比,2#墩頂的位移增幅達到56%,主梁右端的橫向位移增幅也達到47%;此外,如果按照較小的峰值位移進行4個支點的位移控制,一般會得到相對較小的位移結果,減小幅值達66%。圖14(b)所示的第二類地震動輸入方案也計算得到了類似的結果。綜合分析兩類(6組)地震動輸入下結構的位移響應可以看出,控制不同的位移峰值條件得到的位移響應的變化和結構內力響應變化基本一致,這進一步說明地震位移輸入峰值的控制對結構的位移響應有較大的影響,當采用位移輸入法進行計算時,應該合理控制輸入時程的位移峰值。

圖14 不同地震動輸入方案橋梁橫向位移變化Fig.14 Lateral displacement changes of bridge in different ground motion input schemes

3.2 地震動頻率非平穩特性對結構響應的影響

本文在4個支撐點的相關加速度時程模擬中,采用時-頻包線函數模擬地震動頻率和強度非平穩特性(方案1-1),采用強度包線函數模擬地震動強度非平穩特性(方案2-1)。雖然方案1-2、方案1-3、方案2-2和方案2-3是加速度時程的位移時程,但擬合方案中按不同的位移峰值進行了調整,因此本節僅比較控制了相同加速度峰值的方案1-1和方案2-1的計算結果,來說明頻率非平穩特性對結構響應的影響。

3.2.1 對內力的影響

為了直觀地分析地震動頻率非平穩特性對結構內力響應的影響,圖15給出了采用方案1-1和方案1-2地震動輸入時計算得到的橋墩剪力、扭矩及主梁的跨中彎矩。可以看出,若考慮輸入地震動的頻率非平穩特性,1#橋墩的剪力和扭矩、2#橋墩的扭矩都得到相對較大的值,與僅考慮地震動強度非平穩特性(方案2-1)的計算結果相比最大的增幅達到了63%。比較3跨主梁的跨中彎矩可以發現,雖然方案1-1計算得到的第1和第3跨的跨中彎矩與方案2-1的結果相比略有增加,但增加的幅度不大。由以上的分析可以看出,地震動輸入特性對主梁的跨中彎矩、橋墩剪力的影響不明顯,但對橋墩扭矩的影響比較顯著。總體來看,若地震動擬合中不考慮頻率的非平穩特性,則存在低估結構內力響應的風險。

圖15 方案1-1和方案2-1計算得到的橋梁內力Fig.15 Internal forces of bridge in schemes 1-1 and 2-1

3.2.2 對橫向位移的影響

圖16直觀地給出了兩種輸入方案下,計算得到的橋梁4個關鍵位置的橫向位移。相較于內力的變化,地震動特性變化對位移響應的影響比較明顯,采用考慮頻率非平穩特性的加速度輸入都會得到相對較大的位移響應,與方案2-1輸入時計算得到的位移響應相比,最大的增幅達到了50%。這些結果進一步說明,在空間相關多點地震動擬合中,應該合理地考慮地震動的頻率非平穩變化特性。

圖16 方案1-1和方案2-1計算得到的橋梁橫向位移Fig.16 Lateral displacements of bridge in schemes 1-1 and 2-1

4 結論

本文通過引入不同的包線函數和控制不同的加速度/位移峰值目標,合成了能夠反映空間相關地震動非平穩特性和峰值控制的兩類(6組)多點地震動,對大跨橋梁進行非線性時程分析,得出以下主要結論:

(1)控制地震位移時程的峰值對結構內力和位移響應有較大的影響,輸入位移時程峰值越大,結構進入非線性的程度越深,得到的結構響應越大,且與位移時程的形狀無明顯相關性;

(2)以加速度反應譜和加速度峰值為目標的空間相關地震動擬合中,各支撐點處位移時程的峰值離散性比較大,為了保證結構的安全,空間相關多點地震動擬合中應控制位移峰值的最小取值;

(3)地震動的頻率非平穩特性對橋梁的位移響應、橋墩扭矩影響明顯,對主梁內力、橋墩剪力影響較小;但如果僅考慮地震動強度非平穩特性,則存在低估結構響應的風險。

因此,在實際工程應用中計算大跨橋梁結構的非線性時程響應時,不僅需要合理地描述地震動的空間變化特性,還需要考慮地震動頻率的非平穩特性,同時控制位移輸入的最小峰值,從而合理估計結構的地震響應,進而保證結構抗震設計的安全性。

猜你喜歡
結構
DNA結構的發現
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
循環結構謹防“死循環”
論《日出》的結構
縱向結構
縱向結構
我國社會結構的重建
人間(2015年21期)2015-03-11 15:23:21
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 国产精品永久在线| 美女被操黄色视频网站| 亚洲视屏在线观看| 成人永久免费A∨一级在线播放| 国产三级精品三级在线观看| 午夜福利无码一区二区| 亚洲精品无码人妻无码| 亚洲日本www| 最新国产午夜精品视频成人| 国产精品亚洲一区二区在线观看| 97亚洲色综久久精品| 有专无码视频| 潮喷在线无码白浆| 欧美成人一区午夜福利在线| 呦女亚洲一区精品| 午夜日本永久乱码免费播放片| 亚洲最大在线观看| 99re这里只有国产中文精品国产精品 | 久青草免费在线视频| 99久久精彩视频| 国产国模一区二区三区四区| 色偷偷男人的天堂亚洲av| 亚洲第一精品福利| 一级毛片在线播放免费观看| 2021无码专区人妻系列日韩| 国产 在线视频无码| 99热这里都是国产精品| 伊人久久综在合线亚洲2019| 日韩国产一区二区三区无码| 美女一级毛片无遮挡内谢| 国产成人福利在线视老湿机| 国产欧美日韩另类精彩视频| 亚洲日本中文字幕乱码中文| 四虎成人精品| 亚洲综合亚洲国产尤物| 亚洲乱亚洲乱妇24p| 综合五月天网| 欧美亚洲一区二区三区导航| 永久成人无码激情视频免费| 91色爱欧美精品www| 欧美在线精品一区二区三区| 国产幂在线无码精品| 免费国产无遮挡又黄又爽| 在线观看精品国产入口| 国模视频一区二区| 免费在线成人网| 国产精品综合久久久| 欧美国产综合色视频| 国产精品无码制服丝袜| 九色综合视频网| 五月婷婷综合网| 精品国产中文一级毛片在线看 | 99无码中文字幕视频| 亚洲三级影院| 欧美成人怡春院在线激情| 亚洲成人网在线观看| 亚洲美女一区| 免费jizz在线播放| 国产一级毛片yw| 国产极品粉嫩小泬免费看| 免费高清a毛片| 免费欧美一级| 午夜人性色福利无码视频在线观看 | 亚洲综合久久成人AV| 青青草原偷拍视频| 57pao国产成视频免费播放| 中文字幕1区2区| 极品国产在线| 亚洲中文字幕久久精品无码一区| 97se亚洲综合| 亚洲第一香蕉视频| 再看日本中文字幕在线观看| 她的性爱视频| www.日韩三级| 日本黄网在线观看| 精品撒尿视频一区二区三区| 国产精品hd在线播放| 国产精品黑色丝袜的老师| 午夜久久影院| 国产成人一区免费观看| 欧美一级高清片久久99| 亚洲人成影院在线观看|