?

基于纖維模型的RC結構的鋼筋本構關系研究

2016-01-06 14:58張耀庭趙璧歸杜曉菊盧杰志??
湖南大學學報·自然科學版 2015年9期

張耀庭 趙璧歸 杜曉菊 盧杰志 ??

摘要:以基于纖維模型的OpenSees軟件為平臺,首先,選用常用的鋼筋本構(Steel01, Steel02)和混凝土本構(Mander, Hoshikuma),對PEER數據庫中的兩個鋼筋混凝土柱的擬靜力試驗進行了數值建模與對比分析;再采用OpenSees軟件中Reinforcing Steel鋼筋本構模型,對鋼筋單調拉壓試驗、循環加載試驗、疲勞試驗以及鋼筋混凝土柱擬靜力試驗進行了模擬與對比分析.結果表明:現有的混凝土本構模型已相對成熟,并能較好地滿足非線性數值模擬的精度要求,混凝土本構模型的改變對計算結果的影響不大,而鋼筋本構模型的改變則會對分析結果產生較大的影響;Reinforcing Steel模型有效地提高了非線性模擬分析的精度,但也存在著對鋼筋滯回曲線捏攏效應體現不足、分析結果偏于保守、計算過程較為復雜且較難收斂等缺點,尚需要對鋼筋本構模型進行進一步的改進與完善.

關鍵詞:非線性分析; 纖維模型; 本構模型; Reinforcing Steel模型; 擬靜力試驗

中圖分類號:TU375.3 文獻標識碼:A

在罕遇地震作用下,鋼筋混凝土結構往往會進入到彈塑性甚至是倒塌破壞階段.鋼筋混凝土結構抗震分析所采用的技術手段,主要包括結構試驗與數值分析[1].由于經濟及技術條件的限制,結構試驗往往只能用作結構抗震設計與分析中的輔助工具,因而,數值分析成為工程抗震中的關鍵技術措施.

鋼筋混凝土結構的非線性有限元分析已經發展到基于精細化纖維模型的分析階段,研究表明,影響纖維模型分析結果的準確性及可靠性的因素很多,材料的本構模型[2-6]即為其主要影響因素之一,這包括對現有材料本構模型的完善性與適用性的了解,以及建模時所選用本構模型的合理性.針對以上問題,本文以OpenSees軟件[7]為分析平臺,首先通過對比鋼筋混凝土結構非線性分析結果對鋼筋和混凝土材料的本構模型的敏感性差異,對各種常用材料本構關系的適用性進行研究.然后,通過對鋼筋力學性能試驗及鋼筋混凝土柱擬靜力試驗的模擬,探討OpenSees軟件中現有鋼筋本構模型的優缺點,并對其計算精度和適用性做出客觀評價,以期為同類工程問題提供參考.

1鋼筋、混凝土本構關系對結構彈塑性分析

的敏感性

現階段已有很多材料本構模型被提出并被導入OpenSees軟件中.其中,Steel01和Steel02鋼筋模型,以及Mander[2]和Hoshikuma[3]混凝土側限模型在結構構件響應模擬中最為常用.本章以OpenSees軟件為分析平臺,選用上述模型對PEER試驗數據庫中的2個鋼筋混凝土柱的擬靜力試驗[8-9]進行數值模擬.柱子的基本信息見表1,試驗加載制度為位移幅值遞增的循環加載.建模時,選用基于位移的非線性梁柱單元模型(即Nonlinear Beam Column單元模型)來模擬試驗柱,鋼筋的本構模型分別選用Steel01和Steel02模型,保護層混凝土采用Concrete01模型,約束混凝土的本構模型分別選用Mander和Hoshikuma模型,試驗結果與模擬結果的力位移滯回曲線如圖1和圖2所示.

由圖1及圖2可以看出:

1)模擬結果的力位移關系曲線,在形狀、力位移關系等方面均和試驗結果的滯回曲線相接近,反映出本文采用的纖維模型具有較好的適用性與合理性.

2)圖1(b)和(c)中的滯回曲線形狀十分接近,僅僅是柱頂最大水平承載力有微細差別,表明只改變約束混凝土的本構模型,不會對模擬結果產生很大的影響,這也說明現有的混凝土本構關系已發展得相對成熟,能滿足基于纖維模型的數值分析要求.

3)圖1(c)(d)和圖2(b)(c)中的滯回曲線形狀差別較大,表明:僅改變鋼筋的本構模型時,力位移關系曲線將發生較大的變化.結合2)同時可以證明:大位移循環荷載作用下混凝土已發生脆性破壞,而鋼筋仍然可為結構構件耗散能量提供作用.因此,相對于混凝土本構模型,鋼筋混凝土結構非線性分析結果對鋼筋本構模型的變化更敏感.

4)圖1(b)(c)的曲線不夠平滑,不能充分體現出柱試驗滯回曲線的捏攏效應和大位移加載下的強度退化效應;圖2(b)不能體現出位移超過40 mm后柱的彈塑性響應;圖1(d)和圖2(c)中的曲線較為圓滑,有明顯的捏攏性,與試驗結果更加吻合.表明:Steel02鋼筋模型的適用性和模擬精度均優于Steel01模型.

5)圖1(b)(c)中的模擬曲線完全沒有顯示出大位移下柱頂水平承載力的退化效應,且柱頂水平承載力遠大于試驗值;圖1(d)和圖2(c)中,在位移較大時,柱頂水平承載力表現出了一定的退化效應,但不夠明顯,且柱頂水平承載力遠大于試驗值.這表明:采用Steel02鋼筋模型尚不能很好地模擬出試驗柱在較大位移幅值下的強度退化效應,其模擬精確度有待進一步提高.

綜上所述,在鋼筋混凝土結構基于纖維模型的彈塑性分析中,由于混凝土在大位移循環荷載作用下會較早發生脆性破壞,因而分析過程對混凝土模型的精度要求不高,現有的混凝土材料的本構關系已經能滿足結構彈塑性分析的計算要求;而鋼筋屈服后仍然可為結構構件抵抗外力作用,塑性分析過程也因此對鋼筋本構有更高的精度要求.上述模型中Steel01模型應力應變曲線呈折線形,未考慮鋼筋材料的包晶格效應;Steel02模型應力應變曲線雖呈光滑曲線形,適當考慮了包晶格效應,但參數定義過于簡略,未考慮鋼筋受力過程中橫截面積的變化、受壓屈曲效應以及疲勞損傷等因素.因此,2種鋼筋模型的模擬結果均未明顯體現出柱子的強度退化效應,需進一步完善與改進.2009年,針對Steel01和Steel02存在的缺點,文獻[10]提出了一個改進的鋼筋本構模型——Reinforcing Steel模型,并將其嵌入到OpenSees軟件中,由于該模型開發的時間不長,其適用性尚有待研究.對此,本文將以Reinforcing Steel模型為重點,探討其適用性.

2Reinforcing Steel模型簡介

鋼筋模型Reinforcing Steel[7]和Steel01,Steel02模型相比有如下特點:

1)參數定義復雜.Reinforcing Steel鋼筋模型除了定義鋼筋基本力學性能參數外,還需補充定義鋼筋長度與直徑比值Isr,屈服曲線形狀參數alpha,疲勞參數Cf和Cd等.只有在所有參數都得到合理賦值的情況下,才能準確模擬鋼筋的各項力學性能.

2)使用自然坐標描述鋼筋的應力應變關系.該模型考慮了鋼筋受力過程中橫截面積的變化,它可將用戶輸入的工程坐標下的應變εs,應力fs,以及彈性模量Es分別轉化成自然坐標下的值ε′s,fs′和E′s.坐標轉換式如式(1)~(3)所示.

εs′=ln1+εs, (1)

fs′=fs1+εs,(2)

Es′=Es+fsεs+1εs+12.(3)

坐標轉換前后的鋼筋單調受拉的應力應變曲線如圖3所示.

3)以受拉時的應力應變曲線為基礎,定義了一套復雜而完整的加、卸載流動法則用于模擬鋼筋的循環受力性能,并充分考慮了鋼筋的包晶格效應.這一套循環法則的原型是Chang和Mander[11]提出、由Kunnath[10]修正的鋼筋模型,具體的循環法則和對應的循環應力應變表達式,參見文獻[8-9].

4)能夠反映鋼筋的屈曲效應和疲勞效應.Reinforcing Steel鋼筋模型融合了Gomes和Appleton[12]提出的鋼筋屈曲模型、Dhakal和Maekawa[13]提出的鋼筋屈曲模型以及Coffin [14-15]和Manson [16]提出的鋼筋疲勞退化模型.

有關Reinforcing Steel鋼筋模型的使用方法和各個模塊參數意義,可參考文獻[7-16].Reinforcing Steel模型能考慮鋼筋的屈曲和疲勞效應,在理論上比Steel01和Steel02模型更加完善.

3鋼筋力學性能及柱擬靜力試驗的模擬分析

3.1鋼筋單調拉壓響應的模擬

本文選用Dodd(1992)[17]完成的鋼筋單調受壓試驗進行模擬與分析.鋼筋基本參數:fy =295.4 MPa,fu=479.5 MPa,Es=210 000 MPa,Esh=3 367.6 MPa,esh=0.016,eult=0.207.由于該試件長細比很小,試驗過程中沒有發生屈曲,故模擬過程中

忽略受壓屈曲效應的影響.本文采用Reinforcing Steel模型同時進行了鋼筋的單調拉壓模擬.試驗結果的應力應變曲線、模擬結果的應力應變曲線以及轉換坐標后的曲線(拉為正,壓為負)如圖4所示.

由圖4可以發現:

1)單調受壓的模擬結果曲線和試驗結果曲線在形狀、屈服強度、屈服平穩段長度、應力應變關系等方面吻合很好.

2)圖4(c)的單調拉壓模擬曲線中,鋼筋進入強化階段后,在應變絕對值相等時,拉應力比壓應力絕對值小.由于試驗過程中忽略了鋼筋橫截面積和原始長度變化的影響,因此,這符合鋼筋受壓試驗對壓應力的測量值會偏大而受拉試驗對拉應力的測量值會偏小的事實.表明:Reinforcing Steel模型考慮了鋼筋受力過程中橫截面積的變化.

3)圖4(c)中對模擬曲線進行了工程坐標至自然坐標的轉換,轉換結果是考慮了鋼筋橫截面積變化的應力應變曲線,該曲線關于原點對稱,符合鋼筋實際受拉與受壓特性之間的關系.

因此,模擬結果表明:Reinforcing Steel模型對于鋼筋單調受力性能的模擬能很好地吻合試驗結果;采用自然坐標系下的應力應變關系能夠有效提升模型精度.

3.2鋼筋拉壓循環響應的模擬

本文采用Reinforcing Steel模型對鄭江峰[18]完成的鋼筋拉壓循環加載試驗進行模擬與分析.2根鋼筋試件的編號分別為LY1和LY2,鋼筋型號均為HRB335.

基本參數:fy =344.7 MPa,fu=492.8 MPa,Es= 191 922 MPa,Esh=3 838.4 MPa,esh=0.028,eult=0.18.

屈曲參數:選用Dhakal和Maekawa屈曲模型.對于LY1,取Isr=8,alpha=1.0;對于LY2,取Isr=10,alpha=1.0.

不考慮疲勞效應的影響.位移幅值分別取1 mm, 2 mm, 3 mm, 4 mm(LY2僅取前3個幅值).模擬結果的應力應變曲線和試驗曲線分別如圖5和圖6所示.

對比圖5及圖6的2組曲線可以發現:

1)當應變幅值較小時,模擬曲線較好地反映了試驗鋼筋的受壓屈曲性能,模擬曲線在形狀、應力應變關系、屈曲發生位置、屈曲后曲線走向上都和試驗曲線具有較好的一致性;當應變幅值較大時,模擬曲線的捏攏效果差,受壓后的卸載剛度與試驗相差較大.這表明:采用Reinforcing Steel模型模擬大變形加載下鋼筋的響應時,還需要進一步提高其精度.

2)試件模擬破壞的發生時間比試驗結果提前.對于試件LY1,模擬結果顯示在應變達到0.04后,試件的剛度和強度明顯減小,預示試件即將破壞,但試驗結果顯示在加載應變幅值達到0.06后試件才逐漸進入破壞階段;試件LY2進入破壞階段時的應變模擬值也比試驗值小.這表明:Reinforcing Steel鋼筋模型對鋼筋承受循環荷載時體現出的延性性能的估計偏于保守,對應力應變曲線捏攏性體現不足,其本構模型曲線的凹凸性和實際曲線有所相悖,尚需進一步修正.

3.3鋼筋屈曲響應模擬

本文選取Bae[19]完成的鋼筋單調壓縮試驗,來校驗Reinforcing Steel鋼筋模型屈曲響應的模擬精度.鋼筋試件的基本參數:fy = 437.0 MPa,fu = 728.0 MPa,Es = 198 600 MPa,Esh = 4 000 MPa,esh = 0.009 2, eult = 0.147.屈曲參數:選用Dhakal和Maekawa屈曲模型,形狀參數alpha取1.0,長細比Isr分別為4, 6, 8, 12.

加載制度為單調壓縮,最大壓應變為0.15.模擬結果曲線和試驗結果曲線如圖7所示.

由圖7的計算結果可以看出:

1)模擬曲線在形狀上和試驗曲線基本一致,可大致分為彈性階段、屈服平穩段、強化上升段、屈服下降段和受壓破壞平穩段等部分,這與鋼筋的實際受壓性能相吻合.表明:Reinforcing Steel鋼筋模型可以較為精確地體現出鋼筋的受壓屈曲性能.

2)模擬曲線在峰值應力、屈曲發生位置、破壞發生位置、屈曲后強度下降速度等方面與試驗曲線有所不同,特別是對于長細比較?。↙/D<6)的試件,模擬曲線的峰值應力偏大,屈曲發生時的應變偏小,屈曲后強度下降快,而試驗曲線屈曲效應尚不明顯;對于長細比較大(L/D>8)的試件,模擬曲線由幾段曲線組成,有明顯的分界點和破壞發生點,破壞發生后曲線的斜率立即減小至零,而試驗曲線則比較平滑,斜率逐漸減小,沒有明顯的破壞發生點.這是由于Reinforcing Steel鋼筋模型的屈曲模塊把屈曲點和破壞點明確定義為長細比L/D的函數,同時把退化曲線簡化為一段斜直線的緣故.實際的鋼筋屈曲點發生位置和屈曲后應力應變曲線形狀受制于多種因素,具有較強的隨機性.

3)Reinforcing Steel鋼筋模型對長細比較大(L/D>8)的試件的屈曲響應模擬精度比長細比較小的試件的模擬精度更高;模擬長細比較小的鋼筋的響應時,忽略屈曲效應的模擬效果可能更好.

3.4鋼筋疲勞破壞效應的模擬

本文選取Brown 和 Kunnath[20]完成的鋼筋的等應變幅值循環拉壓加載試驗進行疲勞破壞模擬.基本參數:fy = 538.0 MPa,fu = 687.0 MPa,Es = 200 000 MPa, Esh = 3 000 MPa,esh = 0.012,eult = 0.116.疲勞參數:根據文獻[20]試驗數據處理結果,alpha=0.44,Cf=0.16,Cd=0.35.若考慮屈曲效應,則取長細比Isr=6.0,加載應變幅值為0.025.考慮屈曲、不考慮屈曲的模擬結果以及試驗結果如圖8所示.

由圖8可以看出:

1)模擬結果曲線在形狀、應力應變關系、最大拉壓應力、強度退化速率等方面和試驗結果曲線吻合度很高,這說明Reinforcing Steel模型模擬鋼筋疲勞損失效應的精度比較高;對于鋼筋破壞所需半周期數,試驗結果是28個半周期,模擬結果是21個半周期,比試驗結果小,表明疲勞周期參數Cf取值偏小,Reinforcing Steel鋼筋模型的精準度還需進一步提高.

2)在圖8(b)(c)中,考慮屈曲后的模擬曲線與試驗曲線更吻合,每個循環周期的最大壓應力值和試驗值的誤差更??;未考慮屈曲的模擬曲線的最大壓應力值偏大,但對于鋼筋破壞所需半周期數,其結果和考慮屈曲的模擬估計值相等.表明:考慮鋼筋屈曲效應會對Reinforcing Steel鋼筋模型的疲勞應力應變曲線模擬有利,但對其疲勞破壞周期數的評估結果沒有影響.即:Reinforcing Steel鋼筋模型的疲勞模塊和屈曲模塊并不相互耦聯.

3.5鋼筋混凝土柱擬靜力試驗的模擬與鋼筋模型

的比較分析

選取陸新征和文獻[8]所完成的鋼筋混凝土柱的擬靜力試驗進行再次模擬驗證.陸新征試驗柱的基本信息見表2,文獻[8]試驗柱信息仍如表1所示.對2個試驗柱進行建模時,保護層混凝土均采用Concrete01模型,核心混凝土采用Mander模型;為便于對比,2個柱鋼筋模型均分別選用Steel02模型和Reinforcing Steel模型.相應的材料參數取值見表3,兩根柱試驗模擬結果分別見圖9和圖10.

由圖9和圖10均可以看出:

1)在位移較小時,Reinforcing Steel模型比Steel02模型能更精確地體現循環加載條件下鋼筋混凝土柱的力位移響應,特別是對于強度退化速度和水平峰值承載力,采用Reinforcing Steel鋼筋模型的模擬結果更加接近試驗結果.表明:Reinforcing Steel鋼筋模型的模擬精確度比Steel02模型更高.

2)采用Reinforcing Steel鋼筋模型時,適當提高疲勞參數Cf和Cd的取值,可以提高較大變形下鋼筋混凝土柱響應的模擬精度,OpenSees軟件對疲勞參數的建議值并不一定準確.對國產鋼筋來講,軟件建議值尚顯保守.

3)Reinforcing Steel鋼筋模型依賴疲勞參數取值,難以實現對柱臨近破壞時的響應模擬,其計算收斂性能比Steel02模型稍差.

4結論與展望

本文以OpenSees軟件為平臺,對鋼筋力學性能及鋼筋混凝土柱擬靜力試驗進行了數值模擬與分析,重點討論了鋼筋混凝土結構的非線性分析結果對現有鋼筋和混凝土本構模型的敏感性,及以Reinforcing Steel為主的鋼筋本構模型的優缺點與適用性,得到以下結論:

1)在鋼筋混凝土結構的非線性分析中,與已有的鋼筋本構模型相比,現有的混凝土本構模型已經相對成熟,混凝土本構模型的改變對模擬結果影響不大,而選用不同的鋼筋本構模型,對分析結果則會產生較大影響,提高鋼筋混凝土結構彈塑性分析結果的精度,應該從改進鋼筋本構模型的精度入手.

2)OpenSees中的鋼筋本構Steel01和Steel02模型,其參數定義較為簡潔,用于模擬鋼筋響應時計算量較小.這2種模型沒有考慮屈曲、疲勞、強度退化等效應的影響,模擬精度較低.特別是Steel01模型,遠遠不能達到現階段鋼筋混凝土結構非線性分析的精度要求.

3)Reinforcing Steel鋼筋模型較為準確地考慮了鋼筋受壓屈曲效應和承受循環荷載下的疲勞效應、強度退化效應,它能夠較為準確地模擬鋼筋單調加載、循環加載和疲勞加載試驗結果;在模擬循環荷載作用下的結構或構件的響應時,其精度比Steel01和Steel02模型更高.

4)Reinforcing Steel鋼筋模型對鋼筋性能、鋼筋混凝土結構或構件的彈塑性響應的模擬偏于保守,Reinforcing Steel鋼筋模型還存在著一些值得改進的地方:鋼筋循環加載響應的模擬曲線的捏攏性不充分;疲勞參數(Cf,Cd和alpha)取值還需進一步確定化;模擬鋼筋屈曲效應時,對鋼筋屈曲發生點應變、破壞發生點應變的估計不夠準確,長細比較小的鋼筋受壓屈曲的極限應力偏于保守,受壓曲線也不夠平滑;對變形較大、瀕于破壞的結構或構件響應的模擬精確度還需進一步提高;計算量較大,且計算過程較難收斂.

5)本文驗證材料本構模型時,大多選擇了國外的試驗,這些既有的鋼筋本構關系對我國鋼筋及鋼筋混凝土結構構件的適用性有待進一步研究.另外,鋼筋與混凝土之間的黏結滑移也會影響到結構構件響應,而Reinforcing Steel模型還未能準確考慮滑移效應.如何配合各種鋼筋模型準確模擬鋼筋混凝土黏結滑移效應,以期進一步提高模擬精度,還需更深入細致的研究工作.

參考文獻

[1]YEONGAE H. Framework for damagebased probabilistic seismic performance evaluation of reinforced concrete frames [D]. San Francisco: School of Civil and Environmental Engineering, University of California, 2009.

[2]MANDER J B, PRIESTLEY M J N, PARK R. Theoretical stress strain model for confined concrete[J]. Journal of Structural Engineering, 1988, 114(8): 1804-1826.

[3]HOSHIKUMA J, KAWASHIMA K, NAGAYA K, et al. Stressstrain model for confined reinforced concrete in bridge piers[J]. Journal of Structural Engineering, 1997, 123(5): 624-633.

[4]MENEGOTTO M, PINTO P E. Method of analysis for cyclically loaded RC plane frames including changes in geometry and nonelastic behaviour of elements under combined normal force and bending[C]//Proceedings of Symposium on the Resistance and Ultimate Deformability of Structures Acted on by Well Defined Repeated Loads. Zurich, Switzerland: International Association for Bridge and Structural Engineering, 1973: 15-22.

[5]張耀庭,盧怡思,杜曉菊,等. 柱端彎矩增大系數對PC框架結構抗震性能影響的研究[J]. 湖南大學學報:自然科學版,2014 ,41(11):37-47.

ZHANG Yaoting, LU Yisi, DU Xiaoju, et al. Study on the effects of moment magnifying coefficients at column ends on seismic capacity of prestressed concrete frame[J]. Journal of Hunan University: Natural Sciences, 2014, 41(11):37-47.(In Chinese)

[6]任亮,方志,王誠. 基于截面纖維模型的RPC箱型橋墩抗震性能分析[J]. 湖南大學學報:自然科學版,2013,40(8):16-21.

REN Liang, FANG Zhi, WANG Cheng. Seismic behavior analysis of PRC box piers based on the fiber element model[J]. Journal of Hunan University: Natural Sciences, 2013,40(8):16-21.(In Chinese)

[7]MAZZONI S, MCKENNA F, SCOTT M H, et al. OpenSees command language manual[M]. Berkeley: PEER, University of California, 2007.

[8]MO Y L, WANG S J. Seismic behavior of RC column with various tie configurations[J]. Journal of Structural Engineering, 2000,126 (10): 1122-1130.

[9]OHNO T, NISHIOKA T. An experimental study on energy absorption capacity of columns in reinforced concrete structures[J]. Structural Engineering and Earthquake Engineering, 1984, 1(2): 137-147.

[10]KUNNATH S K, YEONGAE H, MOHLE J F. Nonlinear uniaxial material model for reinforcing steel bars[J]. Journal of Structural Engineering, 2009, 135(4): 335-343.

[11]CHANG G, MANDER J. Seismic energy based fatigue damage analysis of bridge columns: part I-evaluation of seismic capacity[R]. Buffalo:NCEER,1994.

[12]GOMES A, APPLETON J. Nonlinear cyclic stressstrain relationship of reinforcing bars including buckling[J]. Engineering Structures, 1997,19(10): 822-826.

[13]DHAKAL R, MAEKAWA K. Modeling for post yield buckling of reinforcement[J]. Journal of Structural Engineering, 2002,128(9): 1139-1147.

[14]COFFIN L F. A study of the effect of cyclic thermal stresses on a ductile metal[J]. American Society of Mechanical Engineers, 1954,76(1): 931-950.

[15]COFFIN L F. A note on low cycle fatigue laws[J]. J Mater, 1971, 6: 388-402.

[16]MANSON S S. Behavior of materials under conditions of thermal stress[R]. Ann Arbor, Michigan: Engineering Research Institute,University of Michigan, 1953.

[17]DODD L. The dynamic behavior of reinforcedconcrete bridge piers subjected to New Zealand seismicity[D]. Christchurch, New Zealand: University of Canterbury Christchurch, 1992.

[18]鄭江峰. 考慮屈曲的鋼筋滯回性能試驗研究[D].重慶:重慶大學土木工程學院, 2012.

ZHENG Jiangfeng. Experimental studies on cyclic behavior of reinforcing bars including buckling[D]. Chongqing: College of Civil Engineering, Chongqing University, 2012.(In Chinese)

[19]BAE S, MIESES A, BAYRAK O. Inelastic buckling of reinforcing bars[J]. Journal of Structural Engineering, 2005,131(2): 314-321.

[20]BROWN J, KUNNATH S K. Lowcycle fatigue failure of reinforcing steel bars[J]. ACI Materials Journal, 2004, 101(6): 457-466.

91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合