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

基于自適應模型的在役斜拉橋時變力學性能追蹤

2023-03-20 02:21:12孫華懷陳惟珍
振動與沖擊 2023年5期
關鍵詞:有限元模型

孫華懷, 徐 俊, 陳惟珍

(1.揚州大學 建筑科學與工程學院,江蘇 揚州 225127; 2.同濟大學 橋梁工程系,上海 200092)

由于結構的復雜性,斜拉橋數值模型往往與實際結構存在一定的差異。在整個使用壽命期內,混凝土收縮、徐變和環境溫度效應又始終影響著斜拉橋的力學性能[1]。故在役斜拉橋力學性能的數值評估結果往往與周期性現場檢測不符。因此,如何根據周期性現場檢測準確地追蹤在役斜拉橋時變力學性能極具挑戰。

以往文獻表明,數值模型與實際結構通常存在3個方面的差異:模型結構誤差,模型階次誤差和模型參數誤差[2]。模型參數誤差來自不精確的物理參數,也是斜拉橋數值模型與實際結構差異的主要來源[3]。為了消除模型參數誤差,有限元模型修正技術已成為一種應用廣泛且有效的橋梁結構分析工具[4]。有限元模型修正技術可分為直接法和迭代法兩類。直接法根據測量數據重構結構整體質量和剛度矩陣。迭代法對模型中的物理參數進行修正。目前,使用最廣的迭代法是基于靈敏度分析的有限元模型修正[5]。它通常被描述為一個優化問題,目標是將結構實測特性和模型計算結果間的差異最小化[6-8]。然而,上述研究通過多個誤差函數的加權求和建立一個目標函數,并作為單目標優化問題。在此過程中,假設的權重因子代表不同目標之間的妥協,不能考慮個體目標間的沖突關系[9]。考慮個體目標間的沖突,多目標優化法隨之被采用對橋梁有限元模型進行修正[10-11]。但他們的研究采用傳統的基于梯度的優化算法求解派生的多目標優化問題。這些算法通常獲得多目標優化問題的局部最優解,而不是全局最優解[12]。近年來,人工智能技術并被證明是求解優化問題的替代方法。因此,橋梁結構有限元模型修正問題的求解開始采用多目標優化進化算法,包括遺傳算法[13-14],神經網絡算法[15],粒子群優化算法[16]等。然而,上述研究僅考慮了斜拉橋模型的參數誤差,卻忽略服役階段環境溫度和混凝土時變效應。因此,斜拉橋力學性能的評估結果無法與周期性現場檢測始終保持一致,也無法解釋實測中斜拉橋力學性能的時變現象。

為了分析混凝土結構時變效應,研究人員提出了幾種計算方法,如按齡期調整的等效彈性模量法、徐變率法和逐步數值法[17]。然而,僅少量的研究涉及考慮混凝土時變效應的斜拉橋結構分析。起初,研究重點是分析混凝土收縮、徐變對斜拉橋靜力性能的時變效應,如應力、撓度和索力[18-20]。研究結果表明,在斜拉橋靜力性能的評估中,需要考慮混凝土的時變效應。但他們的研究主要集中于斜拉橋靜力性能。故Si等[21]考慮混凝土老化、徐變和收縮的影響研究了斜拉橋的動力性能。由于混凝土時變效應,斜拉橋自振頻率隨時間逐漸增加。混凝土時變效應在斜拉橋結構動力性能中同樣起著重要作用。因此,研究者提出在斜拉橋設計階段考慮混凝土時變效應。Lozano-Galant等[22]研究了混凝土徐變和收縮對斜拉橋服役階段應力的影響,并提出了在斜拉橋設計階段將混凝土時變效應計入目標服役階段的方法。此外,為了達到理想的最終狀態,Martins等[23]提出了考慮混凝土時變效應的斜拉橋索力的優化算法。然而,他們的研究都忽略了數值模型固有的參數誤差。同時,鮮有研究考慮混凝土時變效應對在役斜拉橋時變力學性能進行數值追蹤。

本研究分析服役階段斜拉橋自適應模型的建立流程。根據斜拉橋設計資料,建立斜拉橋成橋時初始有限元模型。采用多目標優化法修正斜拉橋成橋時有限元模型。基于修正后的成橋有限元模型,考慮混凝土收縮、徐變及環境溫度效應,建立斜拉橋服役階段自適應模型,數值追蹤在役斜拉橋長期荷載作用下的時變力學性能。同時,設計一座在役斜拉橋的周期性現場檢測方案,對服役前2年該橋實際力學性能進行定期現場檢測。

1 在役斜拉橋自適應模型

斜拉橋服役階段,影響斜拉橋力學性能的因素較多,且各影響因素的作用時間域存在差別。混凝土收縮、徐變主要作用于構件開始受力后的三年內,拉索的退化主要發生在斜拉橋服役接近二十年的數年內,而環境溫度效應則始終存在。為了有效分析斜拉橋時變力學性能,準確把握斜拉橋在服役特定時刻的結構狀態,可建立斜拉橋自適應模型。斜拉橋自適應模型是考慮時變效應的有限元模型,通過已記錄的結構歷史行為、結構狀態及環境信息,沿時間正向,從成橋狀態開始,一直追蹤到最近的橋梁檢測時間節點處。首先,根據斜拉橋的檔案信息資料,建立成橋時有限元模型,作為斜拉橋真實結構的初始模型。基于斜拉橋成橋時現場檢測對斜拉橋初始有限元模型進行修正。基于修正的成橋模型,考慮服役階段橋梁結構信息得到斜拉橋自適應模型,實現對各服役時間節點結構的真實模擬,如圖1所示。

圖1 斜拉橋自適應模型的建立Fig.1 Establishing adaptive models of cable-stayed bridges

2 在役斜拉橋時變力學性能數值追蹤

2.1 基于靈敏度分析的斜拉橋有限元模型修正

基于靈敏度分析的有限元模型修正能夠為斜拉橋建立準確的數值模型。為了保證數值結果與現場檢測一致,通常建立多個目標函數,使得斜拉橋有限元模型修正成為一個多目標優化問題。因此,選取修正參數、建立目標函數和采用穩健的優化算法是三個關鍵步驟。

(1)

以斜拉橋成橋時全橋索力數值結果與現場實測值的差異為基礎,構造標準化目標函數gi(X)為

(2)

minimizeG(X)={g1(X),g2(X)}T

s.t.Ω={X|Xl≤X≤Xu}

(3)

式中,Xl、Xu和Ω分別是修正參數向量的下限、上限和可行域。

遺傳算法是一種高效的隨機搜索方法,并用于求解有限元模型修正問題的全局最優解。在求解時,遺傳算法從多個可行解Xi(i=1 toNs,i∈N+)出發。這多個可行解的集合稱為一個種群,記為pop(k),k表示當前迭代步。種群中每個個體代表著多目標優化問題的一個可行解Xi。整個算法循環運行選擇、交叉和變異三個遺傳算子,直至適應度函數滿足收斂條件,如圖2所示。在求解式3的多目標優化問題時,遺傳算法種群規模Ns設置為50。將交叉比Pc設為1.0,為下一種群提供足夠的變異性。求解得到由m個非劣解組成的Pareto前沿Ω*={X1,X2,…,Xm}。從多目標優化問題的Pareto前沿Ω*中,根據距原點的歐氏距離為修正參數向量X選擇一個最佳妥協解Xb。采用修正后的模型對斜拉橋成橋時力學性能進行數值評估。

圖2 遺傳算法求解基于靈敏度分析的有限元模型修正Fig.2 Genetic algorithm for sensitivity-based finite element model updating

2.2 斜拉橋時變力學性能有限元逐步分析

(4)

表A.1 結構構件徐變和收縮模型中各參數細節Tab.A.1 Details of parameters in the creep and shrinkage model of structural members

(5)

(6)

(7)

在數值分析時,斜拉橋的混凝土構件一般采用空間梁單元模擬。在第i個時間間隔Δti內,混凝土單元ele(ele=1 toNe, ele∈N+)的應變增量列陣Δεele(i)由四個部分組成

(8)

(9)

(10)

3 海河大橋結構及周期性現場檢測

海河大橋(新橋)是一座鋼-混凝土混合梁獨塔斜拉橋。該橋全長490 m,主跨為310 m,邊跨為180 m,全橋跨徑組合為(310+2×50+2×40) m,如圖3所示。主跨梁體采用鋼箱梁與預應力混凝土箱梁混合結構,邊跨采用預應力混凝土箱梁結構。混凝土箱梁為預應力結構,縱向預應力鋼束采用7Φs15.20和9Φs15.20高強度低松弛預應力鋼絞線。全橋上、下游共74根斜拉索,斜拉索采用直徑Φ7半平行鋼絲索。橋塔為鋼筋混凝土雙柱變截面結構。

圖3 海河大橋(新橋)結構及周期性現場檢測方案(m)Fig.3 ConFiguration and periodic on-site inspections of Haihe Bridge (the new one) (m)

對海河大橋(新橋)主梁撓度和全橋索力進行定期檢測。根據拉索縱向位置,主跨側拉索編號為C1~C18,橋塔處0號拉索編號為C0,邊跨側拉索對應編號為C1′~C18′(圖3)。橋面高程測點布置在橋面兩端如圖3所示,主橋上、下游各35個高程測點。高程測點編號從主跨至邊跨為“M索號-側”,上游側記為“1”,下游側記為“2”。斜拉索的索力采用無線索力測試分析系統(TST5927)進行檢測。TST5927系統包括數據采集模塊和數據分析模塊(圖3)。模塊間的信號傳遞采用無線傳輸技術(WiFi)。TST5927系統的工作頻率為0~50 Hz。數據采集模塊配有16位A/D轉換器,最大采樣頻率為200 Hz。橋面高程采用全站儀和精密電子水準儀進行測量。竣工后,對服役2年內該橋恒載作用下橋面高程、索力進行定期現場檢測,具體檢測時間及環境條件如表1所示。

表1 周期性現場檢測具體時刻和環境條件細節Tab.1 Details of detection time and environmental conditions of periodic field measurements

4 海河大橋數值模擬和有限元模型修正

4.1 海河大橋(新橋)數值模擬

海河大橋(新橋)有限元模型包括橋塔、橋墩、主梁、斜拉索、預應力鋼束、橫隔板和支座等構件。橋塔和橋墩采用三維梁單元(BEAM188)模擬。混凝土主梁和鋼主梁也采用三維梁單元(BEAM188)模擬。混凝土主梁BEAM188單元縱向尺寸與混凝土主梁節段長度一致(圖4)。鋼主梁BEAM188單元縱向尺寸為鋼主梁節段長度的1(5(圖4)。采用三維梁單元(BEAM4)模擬混凝土箱梁中預應力鋼束。預應力效應通過對鋼束單元施加初應變實現。預應力鋼束單元通過剛度較大的三維梁單元(BEAM4)與主梁BEAM188單元的節點相連。混凝土箱梁和鋼箱梁的橫隔板模擬為138個集中質量單元(MASS21)。斜拉索采用三維桿單元(LINK10)模擬。拉索初始張拉力通過給拉索LINK單元施加初應變實現。拉索垂度采用等效彈性模量予以考慮。采用剛度較大的三維梁單元(BEAM4)連接拉索單元與主梁BEAM188單元。橋墩/主梁處的每個支座分別采用在X、Y和Z方向上三個獨立的一維彈簧單元(COMBIN14)模擬。將全橋二期恒載模擬為均布線荷載,施加在主梁的梁單元上。環境溫度變化模擬為溫度荷載,并施加于全橋相應構件。海河大橋(新橋)有限元模型如圖4所示。

圖4 海河大橋(新橋)有限元模型細節Fig.4 Details of finite element model of Haihe Bridge (the new one)

4.2 海河大橋(新橋)有限元模型修正

由于鋼箱梁內形狀復雜的加勁肋和輔助設施,很難使模擬的鋼箱梁截面(圖4)與設計截面完全一致。鋼箱梁節段的截面剛度和自重較初始設計值略有下降。因此,將鋼箱梁材料的彈性模量和密度(即x1和x2)作為修正參數。考慮構造鋼筋的影響,模型中混凝土橋塔(即x3和x4)、橋墩(即x5和x6)和主梁(即x7和x8)材料的彈性模量和密度也選為修正參數。混凝土箱梁和鋼箱梁的橫隔板的計算質量通常會與實際情況存在差異。因此,混凝土箱梁(即x9~x11)和鋼箱梁(即x12~x14)的橫隔板質量也需要修正。橋面附屬設施將對二期恒載產生一定的影響,故二期恒載(x15)也作為修正參數。已采用等效彈性模量考慮拉索垂度效應并計入護套重量,斜拉索的彈性模量和材料密度不作為修正參數。成橋時拉索的索力檢測存在誤差,斜拉索張拉力(即x16~x52)作為修正參數。主梁/橋墩的連接彈簧剛度的計算值往往支座的實際剛度不一致。模型中主梁/橋墩連接彈簧剛度(即x53~x64)也作為修正參數。

根據式(1)對預選的修正參數進行索力靈敏度分析。結果表明,修正參數x2(鋼箱梁材料密度)具有最大的無量綱的索力靈敏度,數值為0.41。這是因為鋼箱梁密度的變化會引起主跨各拉索索力的變化。修正參數x16~x52(拉索C18~C18′初始張拉力)的無量綱索力靈敏度差異很小,在0.04~0.07。全橋索力對主梁/橋墩連接彈簧剛度(即x53~x64)不敏感,無量綱索力靈敏度數值接近為0。為保證模型修正的效率,去除修正參數x5和x6(輔助墩材料特性)、x10(混凝土箱梁橫隔板)、x14(鋼箱梁橫隔板)、x53~x64(主梁/橋墩連接彈簧剛度)。最后,得到一個48維的待修正參數向量X。

以各修正參數的設計值作為有限元模型修正過程中的初始值(表2)。基于成橋時現場實測索力按式(2)構造目標函數。根據斜拉橋施工中的不確定性和工程結構數值分析經驗,每個修正參數設有±10%的變化范圍(表2),以保證模型修正過程中各修正參數的物理意義。采用遺傳算法求解,733代后得到含有22個非劣解的Pareto前沿Ω*={X1,X2,…,X22}。由距原點的歐氏距離,選擇Pareto前沿Ω*中非劣解X13的作為修正參數的最佳妥協解Xb。模型修正前,目標函數g1(X)和g2(X)初始值分別為6.87%和6.50%。修正后,目標函數g1(X)和g2(X)計算值分別下降為1.06%和1.23%,且非劣解X13將導致下游側拉索索力與實測值的相對差異略大一些。對于非劣解X13,部分參數的具體修正細節如表2。模型修正后,各構件材料特性和二期恒載的修正值與初始值的相對變化均在±10%范圍內。各構件材料特性和二期恒載的修正值均大于初始值。這是由于初始有限元模型中忽略形狀復雜的肋、混凝土構件的構造鋼筋和一些輔助設施,各構件材料特性和二期恒載的修正值補償了截面剛度和質量的損失。模型修正后,全橋拉索張拉力修正值與初始值的相對變化均在±10%范圍內(見附表A.2)。拉索C3張拉力的修正值擁有最大的相對變化,為10%。拉索C0張拉力的修正值擁有最小的相對變化(-4.7%)。同時,從附表A.2中可以發現,全橋大部分拉索張拉力的修正值大于初始值。

表2 修正參數初始值和修正值細節Tab.2 Initial and updated details of updating parameters

5 海河大橋(新橋)時變力學性能

基于斜拉橋自適應模型,數值評估海河大橋(新橋)服役1年和2年后時變力學性能。

5.1 服役1年后力學性能

以服役1年后現場實測索力為基準,圖5(a)和圖5(b)分別分析了上、下游索力初始有限元計算結果和自適應模型計算結果與實測值的相對差異。全橋索力初始有限元計算值普遍小于服役1年后索力現場實測值。上游側索力初始有限元計算結果與實測索力的相對差異在-15.31%~-4.16%。下游側索力初始有限元計算結果與實測索力的相對差異在-15.05%~-2.00%。基于自適應模型,上游側索力計算值與實測索力的相對差異縮小在-6.31%~2.48%之內如圖5(a)所示。下游側索力有限元計算值與實測索力的相對差異縮小在-6.31%~3.01%之內如圖5(b)所示。服役1年后,全橋拉索中上游側拉索C0的索力計算結果與實測索力相對差異最大,具體數值為-10.0%。

(a) 上游側

(b) 下游側圖5 服役1年后實測與計算索力相對差異Fig.5 Relative difference of measured and calculated cable forces after 1 year in service

圖6描述了服役1年后海河大橋(新橋)主梁撓度初始有限元模型結果、自適應模型數值結果及現場檢測結果。以橋面設計高程為參考基準值,計算了成橋狀態和服役1年后主梁撓度的實測結果。數值結果及現場檢測均表明,服役1年后,該橋主跨主梁撓度較大,而邊跨主梁由于輔助墩的作用撓度較小。根據周期性現場檢測,成橋至服役1年間,主跨主梁高程測點M13處撓度向下增加了0.063 m,橋塔處主梁高程測點M0處撓度向下增加了0.015 m。服役1年后,主跨主梁實測撓度最大值為-0.202 m(高程測點M13處)。主梁最大撓度初始有限元模型計算結果為0.098 m,位于M14測點。基于自適應模型,主跨主梁最大撓度數值計算值為-0.225 m,縱向位置位于橋面高程測點M12處。橋塔處主梁(高程測點M0處)撓度實測值和數值結果基本一致,分別為-0.022 m和-0.026 m。全橋主梁撓度實測值和自適應模型數值結果最大差異為-0.025 m。

圖6 服役1年后主梁撓度檢測結果與數值計算結果Fig.6 Measured and calculated girder deflections after 1 year in service

5.2 服役2年后力學性能

以服役2年后現場實測索力為基準,圖7(a)和圖7(b)分別分析了上、下游索力初始有限元計算結果和自適應模型計算結果與實測值的相對差異。全橋索力初始有限元計算值普遍小于服役2年后索力現場實測值。上游側索力初始有限元計算結果與實測索力的相對差異在-16.49%~-1.43%。下游側索力初始有限元計算結果與實測索力的相對差異在-17.36%~-1.12%。基于自適應模型,上游側索力計算值與實測索力的相對差異縮小在-9.02%~3.58%之內如圖7(a)所示。拉索C0索力計算結果與實測索力相對差異最大。下游側索力計算值與實測索力的相對差異縮小在-8.28%~3.94%之內如圖7(b)所示。服役2年后,全橋拉索中上游側拉索C0索力計算結果與實測索力相對差異最大(-9.02%)。

(a) 上游側

(b) 下游側圖7 服役2年后實測索力與計算索力相對差異Fig.7 Relative differences of measured and calculated cable forces after 2 years in service

圖8描述了服役2年后主梁撓度數值結果及現場檢測結果。以橋面設計高程為參考基準值,計算服役2年后主梁撓度的實測結果。同樣,服役2年后,該橋主跨主梁撓度較大,而邊跨主梁撓度較小。根據周期性現場檢測,成橋至服役2年間,主跨主梁高程測點M13處撓度向下增加了0.088 m,橋塔處主梁高程測點M0處撓度向下增加了0.021 m。服役2年后,主跨主梁實測撓度最大值為-0.230 m,縱向位置位于橋面高程測點M11處。主梁最大撓度初始有限元模型計算結果為0.098 m,位于M14測點。基于自適應模型,主跨主梁最大撓度數值計算值為-0.240 m,縱向位置位于橋面高程測點M12處。橋塔處主梁(高程測點M0處)撓度實測值和數值結果分別為-0.028 m和-0.026 m。全橋主梁撓度實測值和數值結果最大差異為-0.013 m(高程測點M12處)。

圖8 服役2年后主梁撓度檢測結果與數值計算結果Fig.8 Measured and calculated girder deflections after 2 years in service

服役1年和2年后,全橋索力數值計算結果與實測索力相對差異最大的均是上游側拉索C0。這可能是由于海河大橋(新橋)是漂浮體系,全橋拉索中橋塔處拉索C0的初始張拉力最大(見附表A.2)。服役1~2年內,全橋主梁撓度自適應模型計算結果與實測值間的最大差異從0.025 m降為0.013 m。這可能是由周期性現場檢測的隨機誤差和橋梁服役性能退化所共同導致。

表A.2 拉索張拉力初始值和修正值細節Tab.A.2 Details of initial and updated cable tension forces

6 結 論

基于服役階段斜拉橋自適應模型,數值評估了海河大橋(新橋)服役1年和2年后的時變力學性能,并與周期性現場檢測進行了對比,得到以下幾點結論:

(1) 服役1年和2年后,海河大橋(新橋)全橋索力初始有限元結果與實測值的相對差異分別達-15.31%和-17.36%,實測主梁最大撓度分別向下增加了0.063 m和0.088 m。初始有限元模型評估在役斜拉橋力學性能會存在超出10%的相對誤差,且無法追蹤斜拉橋時變力學性能。

(2) 基于自適應模型,海河大橋(新橋)全橋索力數值結果與周期性現場檢測值的相對差異始終保持在±10%之內。服役1年和2年后,主梁撓度數值結果和實測值最大差異分別為-0.025 m和-0.013 m。服役階段斜拉橋自適應模型能夠有效追蹤在役斜拉橋時變力學性能,且數值結果與實測值的相對差異保持在10%以內。

(3) 自適應模型的數值結果表明,服役階段海河大橋(新橋)主跨主梁撓度逐漸增大,邊跨主梁撓度變化較小。服役1年和2年后,主跨主梁最大撓度數值結果分別為-0.225 m和-0.240 m,橋塔處主梁撓度數值結果均為-0.026 m。

本文主要研究在役斜拉橋時變力學性能數值追蹤方法這一復雜的實際工程問題。海河大橋周期性現場檢測數據用于有限元模型修正及自適應模型驗證。實際上,現場檢測存在傳感器測試噪聲和檢測誤差,且與環境條件(風、溫度、交通荷載)和技術人員水平相關。由于目前索力實測數量有限,本研究未對索力實測數據的精度進行分析。后期將持續收集索力長期檢測數據,研究基于索力檢測不確定性的斜拉橋有限元模型修正方法。

附錄A

猜你喜歡
有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 国模视频一区二区| 亚洲中文字幕av无码区| 青青草原国产av福利网站| 国产一级二级在线观看| 国产69精品久久久久孕妇大杂乱| 激情乱人伦| 日韩第九页| 国产一区二区三区免费| jizz国产在线| 国产交换配偶在线视频| 国产白浆一区二区三区视频在线| 综合社区亚洲熟妇p| 亚洲午夜福利在线| 91啦中文字幕| 18禁不卡免费网站| 99热这里只有成人精品国产| 国产精品免费露脸视频| 色亚洲成人| 99在线视频精品| 亚洲一级毛片免费看| 欧美国产精品不卡在线观看| 国产精品永久在线| 玖玖精品在线| 大香网伊人久久综合网2020| 国产精品999在线| 亚洲av无码成人专区| 男女男免费视频网站国产| 影音先锋丝袜制服| 色综合五月| 国产一级α片| 国产国产人成免费视频77777| 欧美高清国产| 欧美国产日韩在线播放| 久久网欧美| 欧美在线视频不卡第一页| 老色鬼欧美精品| 不卡网亚洲无码| 国产精品无码翘臀在线看纯欲| 99精品福利视频| 毛片免费试看| 国产美女在线免费观看| 免费毛片a| 最新国产你懂的在线网址| 欧美精品另类| 国产高清又黄又嫩的免费视频网站| 国产一区二区人大臿蕉香蕉| 伊人精品视频免费在线| 国产一区二区影院| 99在线视频网站| 日韩久草视频| 少妇精品网站| 欧美成人区| 精品视频一区二区观看| 中文字幕欧美日韩高清| 久久精品这里只有精99品| 真实国产乱子伦高清| aaa国产一级毛片| 91无码人妻精品一区| 亚洲日本中文字幕乱码中文| 中文字幕在线看视频一区二区三区| 亚欧成人无码AV在线播放| 91蝌蚪视频在线观看| 国产精品网址在线观看你懂的| 久久午夜夜伦鲁鲁片不卡| 日本一本在线视频| 制服丝袜 91视频| 九九香蕉视频| 99视频在线免费| 无码免费试看| 日本久久免费| 亚洲第一色视频| 尤物国产在线| 久久不卡国产精品无码| 日日拍夜夜操| 欧美综合激情| www成人国产在线观看网站| 中文字幕丝袜一区二区| 天天躁夜夜躁狠狠躁躁88| 看国产毛片| 欧美一区二区啪啪| 91成人在线免费视频| 在线日本国产成人免费的|