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

橋梁混凝土收縮徐變效應(yīng)時變過程綜合仿真

2021-03-16 07:22:38王龍飛
公路交通科技 2021年2期
關(guān)鍵詞:有限元效應(yīng)混凝土

王龍飛

(1.甘肅省交通規(guī)劃勘察設(shè)計院股份有限公司,甘肅 蘭州 730030; 2.長安大學(xué) 公路學(xué)院,陜西 西安 710064)

0 引言

橋梁混凝土結(jié)構(gòu)往往分階段施工,為了縮短工期,混凝土加載齡期一般較早,所以結(jié)構(gòu)收縮徐變效應(yīng)較大。對于分階段施工的超靜定混凝土橋梁,收縮徐變會不斷影響結(jié)構(gòu)的內(nèi)力和變形,對于鋼混組合橋梁,還會導(dǎo)致截面應(yīng)力不斷重分布[1],所以規(guī)范要求必須考慮橋梁混凝結(jié)構(gòu)收縮徐變效應(yīng)[2]。由于徐變與結(jié)構(gòu)混凝土應(yīng)力有關(guān),所以收縮徐變引起的應(yīng)力不斷變化會導(dǎo)致混凝土部分應(yīng)力的加載齡期在不斷改變,因此,要精確計算收縮徐變效應(yīng)就應(yīng)考慮其相互作用的整個時間變化過程。

收縮徐變是混凝土橋梁研究和設(shè)計的重要內(nèi)容,所以,長期以來,許多學(xué)者對收縮徐變效應(yīng)進(jìn)行了多方面的研究。薛偉辰等[3]通過對6根預(yù)應(yīng)力混凝土模型梁進(jìn)行長期試驗,系統(tǒng)研究了結(jié)構(gòu)的徐變性能,并提出了徐變變形設(shè)計計算建議公式。Masovic等[4]設(shè)計了先簡支后連續(xù)的鋼筋混凝土試驗梁,并對其收縮、徐變和裂縫發(fā)展等特性進(jìn)行了4年以上的持續(xù)研究,獲得了結(jié)構(gòu)的長期性能和內(nèi)力簡化計算方法。劉沐宇等[5-6]研究了混凝土橋梁收縮徐變的時變性、不確定性以及有限元計算方法。項貽強(qiáng)等[7-9]對節(jié)段施工的連續(xù)剛構(gòu)橋混凝土收縮徐變的時變特性進(jìn)行了研究,并提出了相應(yīng)的模擬計算方法。隨著鋼混組合結(jié)構(gòu)的優(yōu)勢不斷為人們所認(rèn)識和接受,工程應(yīng)用在不斷增多,所以,對各種組合橋梁的受力性能和收縮徐變特性的研究也在快速發(fā)展。王春生等[10-12]對高強(qiáng)鋼組合梁、管翼緣組合梁、波形鋼腹板組合梁進(jìn)行了試驗研究,獲得了相應(yīng)的受力特性。樊健生等[13-14]對考慮收縮、徐變和裂縫影響的鋼混組合梁長期性能進(jìn)行了理論、試驗和計算方法研究,提出了基于這些因素的計算模型,并得到試驗驗證。Al-deen等[15-16]通過研究鋼混組合試驗梁的長期性能,提出了考慮收縮徐變、材料特性和剪力鍵的結(jié)構(gòu)計算方法、設(shè)計準(zhǔn)則和長期變形評價方法。劉亞茹等[17-18]研究了鋼混組合梁的長期變形、收縮徐變和剪力滑移等效應(yīng),獲得結(jié)構(gòu)相應(yīng)特點并提出了有限元模擬計算方法。總體來看,對橋梁混凝土結(jié)構(gòu)的收縮徐變效應(yīng)等長期性能的研究一直在廣泛開展,對其機(jī)理的認(rèn)識不斷加深,計算方法也在持續(xù)改進(jìn),獲得了許多重要成果,指導(dǎo)了大批混凝土結(jié)構(gòu)的設(shè)計和施工,保證了結(jié)構(gòu)的安全和質(zhì)量。但從結(jié)構(gòu)設(shè)計工作的角度來看,目前很多收縮徐變效應(yīng)分析方法比較復(fù)雜,不易為設(shè)計人員所掌握,常用的有限元分析軟件,要么不能直接計算收縮徐變效應(yīng),要么不能考慮時變過程,要么不便與規(guī)范計算方法相銜接,因此,普通設(shè)計人員難以對結(jié)構(gòu)進(jìn)行精確的收縮徐變效應(yīng)分析并嚴(yán)格按規(guī)范設(shè)計,影響了混凝土結(jié)構(gòu)在這方面分析的有效性,需要有更廣泛實用的計算方法。

本研究擬結(jié)合混凝土結(jié)構(gòu)收縮徐變的特點和有限元法計算原理,利用APDL語言改進(jìn)通用有限元軟件ANSYS中一些計算功能,形成可以同時考慮混凝土收縮和徐變在整個時變過程中相互作用效應(yīng)的仿真方法。由于ANSYS應(yīng)用非常廣泛,時變過程綜合仿真分析方法使用簡便,且可很好地與各種規(guī)范和試驗研究相銜接,具有較強(qiáng)的實用性,便于推廣應(yīng)用。

1 增量法有限元計算原理

徐變是混凝土在應(yīng)力作用下應(yīng)變不斷變化的現(xiàn)象,當(dāng)應(yīng)力不大(不超過混凝土強(qiáng)度50%)時,具有對應(yīng)力線性疊加的特點,根據(jù)試驗研究及工程經(jīng)驗,已提出了多種基于老化理論和先天理論的徐變模型。收縮是混凝土隨著時間的增長而出現(xiàn)不斷增大的壓縮應(yīng)變,與結(jié)構(gòu)應(yīng)力無關(guān),其變化規(guī)律可認(rèn)為與徐變相似。混凝土規(guī)范中一般會給出收縮徐變的經(jīng)驗?zāi)P汀?/p>

由于收縮徐變效應(yīng)隨時間不斷變化,對鋼混組合結(jié)構(gòu)和超靜定結(jié)構(gòu),會導(dǎo)致截面應(yīng)力重分布和變化的次內(nèi)力,非常復(fù)雜,其理論方程[19]如下:

(1)

式中,ε為含收縮徐變總應(yīng)變;σ為應(yīng)力;E為彈性模量;φ為徐變系數(shù);εs為收縮應(yīng)變;τ0為初始加載齡期;τ為加載齡期;t為計算齡期。

顯然,該方程并不能直接積分求解,引入老化系數(shù)(即為狄辛格法),再通過微分方法或代數(shù)方法可求解。但是實際工程結(jié)構(gòu)往往很復(fù)雜,有很多施工階段和體系轉(zhuǎn)換,采用狄辛格法計算非常困難,甚至不可能完成。隨著計算機(jī)和有限元技術(shù)的發(fā)展,利用有限元法對結(jié)構(gòu)進(jìn)行分析越來越普遍,有限元法不僅可較好地解決復(fù)雜計算問題,也容易掌握,實際應(yīng)用極為方便。參考相關(guān)的有限元法文獻(xiàn)[20],易得到適合收縮徐變效應(yīng)分析的增量法有限元計算原理。

假設(shè)結(jié)構(gòu)處于線彈性狀態(tài),僅考慮彈性內(nèi)力、徐變和收縮效應(yīng)。在某一時段tn內(nèi),單元應(yīng)變增量為:

Δε(tn)=Δεe(tn)+Δεc(tn)+Δεs(tn),

(2)

式中,Δεe(tn)為單元內(nèi)力變化引起的彈性應(yīng)變增量;Δεc(tn)為單元徐變應(yīng)變增量;Δεs(tn)為單元收縮應(yīng)變增量,則應(yīng)力增量為:

{Δσ(tn)}=De{Δεe(tn)},

(3)

式中,De為彈性矩陣,可得時段tn的切線矩陣:

(4)

于是單元外力、徐變和收縮效應(yīng)節(jié)點荷載增量分別為:

(5)

(6)

(7)

式中,B為應(yīng)變矩陣,單元切線剛度可表示為:

(8)

將各單元節(jié)點荷載增量集合起來,形成結(jié)構(gòu)總體平衡方程為:

K{Δδ(tn)}={Δp(tn)}L+

{Δp(tn)}C+{Δp(tn)}S,

(9)

式中,{Δδ(tn)}為節(jié)點位移增量。這樣,可求出節(jié)點位移增量,再計算應(yīng)變增量,再由式(2)、(3)得到應(yīng)力增量,然后各時段累加,即可得到所有時段的整體應(yīng)力:

(10)

從整個過程可以看出,有限元法思路簡潔明了,容易理解,一個時段接一個時段地計算,符合結(jié)構(gòu)整個受力時間演變過程,極便于對考慮施工階段、體系轉(zhuǎn)換和收縮徐變的結(jié)構(gòu)進(jìn)行時變過程仿真分析。

2 基于ANSYS的綜合仿真方法

ANSYS是國際著名的大型通用有限元軟件,獲得世界眾多專業(yè)技術(shù)協(xié)會認(rèn)可,廣泛應(yīng)用在各工程行業(yè)的科研、開發(fā)、設(shè)計、施工等領(lǐng)域,同時,該軟件提供了APDL,UPFs等多種二次開發(fā)工具,可依需要高效擴(kuò)充有限元計算功能。但是目前ANSYS有限元軟件并沒提供專門的混凝土收縮徐變效應(yīng)計算模塊,不能直接分析,但它提供了具有時變特性的金屬蠕變計算準(zhǔn)則和溫度效應(yīng)計算方法,通過二次開發(fā)后,可對混凝土收縮徐變效應(yīng)進(jìn)行分析。

2.1 徐變效應(yīng)仿真計算

可利用ANSYS中的金屬蠕變計算準(zhǔn)則來計算徐變效應(yīng),選取顯式分析中C6=0時的計算準(zhǔn)則:

(11)

由上文的增量法有限元計算原理可知,只要確定了tn時段徐變應(yīng)變增量就可以求出相應(yīng)效應(yīng)。假設(shè)該時段有效應(yīng)力不變,根據(jù)徐變系數(shù)定義,可得:

(12)

(13)

按照規(guī)范,通常混凝土初始加載齡期確定后,徐變系數(shù)對時間的變化模型是確定的,所以其任意時刻的變化率也是可知的,所以利用式(12)、式(13),有限元法就可以進(jìn)行徐變效應(yīng)計算。值得注意的是,徐變系數(shù)模型確定后,按初始加載齡期可準(zhǔn)確計算后續(xù)加載的徐變系數(shù)變化率,不受后期加載大小影響,只需確定后續(xù)加載相對初始加載時間間隔,因此,本方法對具有施工階段、體系轉(zhuǎn)換、組合梁和超靜定結(jié)構(gòu)等有復(fù)雜時變過程的結(jié)構(gòu)進(jìn)行徐變效應(yīng)分析,具有很大的優(yōu)勢。

2.2 收縮效應(yīng)仿真計算

收縮應(yīng)變與應(yīng)力無關(guān),隨著時間的增長而不斷增大,這顯然與結(jié)構(gòu)持續(xù)降溫相似,所以在ANSYS中可采用降溫計算的方法來進(jìn)行收縮效應(yīng)分析。

首先通過混凝土材料、結(jié)構(gòu)尺寸和使用環(huán)境確定混凝土結(jié)構(gòu)的收縮時變模型,利用tn時段的起、終點時間求出該時段的收縮應(yīng)變,計算出相應(yīng)的降溫值,再施加于混凝土結(jié)構(gòu)上,就可以分析收縮效應(yīng)。tn時段計算式如下:

(14)

(15)

2.3 收縮徐變效應(yīng)綜合仿真方法

根據(jù)以上的收縮、徐變效應(yīng)計算方法,可以看出,收縮徐變效應(yīng)均是每個時段的函數(shù),確定了時段,即可確定相應(yīng)的效應(yīng)。ANSYS具有荷載步計算功能,每一個荷載步相當(dāng)于一個時段,利用其強(qiáng)大的二次開發(fā)功能,將每一時段的收縮徐變作用輸入每一荷載步中,然后一步一步進(jìn)行計算,就可以得到任一時間的含收縮、徐變的綜合效應(yīng),綜合仿真方法整個計算過程與結(jié)構(gòu)實際受力時間演化過程極為相似。

圖1 時變過程綜合仿真流程Fig.1 Flowchart of comprehensive simulation during time-varying process

利用ANSYS中的APDL語言將整個計算過程編制成命令流,由計算機(jī)自動進(jìn)行,非常方便,具體流程如圖1所示,關(guān)鍵步驟說明如下:

(1)時段劃分。由于收縮徐變均先期影響較大,所以時間分段時,前期較小,后期可較大。對于有施工階段和體系轉(zhuǎn)換的結(jié)構(gòu),應(yīng)將其設(shè)置成瞬時時段,在ANSYS中建立數(shù)組輸入所有時間節(jié)點,便于計算時調(diào)取。

(2)收縮作用輸入。利用收縮應(yīng)變時變模型計算各時段的收縮應(yīng)變增量,轉(zhuǎn)化為各時段的降溫量,編制成相應(yīng)命令流。

(3)徐變作用輸入。按徐變系數(shù)時變模型計算各時段的徐變系數(shù)增量,形成各時段的C1值,編制成相應(yīng)命令流。

(4)對收縮徐變作用的輸入,應(yīng)盡量統(tǒng)一為時段的連續(xù)函數(shù),充分利用ANSYS 中的循環(huán)和判斷命令,編制成系統(tǒng)的命令流,讓程序自動計算。

3 與規(guī)范的銜接

綜合仿真方法可很方便地與各種設(shè)計規(guī)范相銜接,對結(jié)構(gòu)的收縮徐變效應(yīng)進(jìn)行模擬計算,下面以公路橋梁規(guī)范[21]為例進(jìn)行銜接應(yīng)用。

3.1 徐變效應(yīng)

文獻(xiàn)[21]提供的公式如下:

φ(t,t0)=φ0βc(t-t0),

(16)

φ0=φRHβ(fcm)β(t0),

(17)

(18)

式中,φ0為名義徐變系數(shù),結(jié)構(gòu)確定后,式(17)的計算參數(shù)均可以確定,是常數(shù),同時βH也為常數(shù);t0為加載時的混凝土齡期,t1=1 d。fcm為混凝土平均立方體抗壓強(qiáng)度。

顯然,規(guī)范公式是對時間的連續(xù)函數(shù),并具有連續(xù)導(dǎo)數(shù),也就是徐變系數(shù)變化率為:

(19)

將其代入式(12)即可直接計算徐變效應(yīng)。

3.2 收縮效應(yīng)

文獻(xiàn)[21]提供了如下公式:

εcs(t,ts)=εcs0βs(t-ts),

(20)

εcs0=εs(fcm)βRH,

(21)

(22)

式中,εcs0為名義收縮系數(shù),結(jié)構(gòu)確定后,式(21)的計算參數(shù)均可以確定,是常數(shù),同時βSH也為常數(shù);ts為收縮開始時的混凝土齡期,t1=1 d。

同樣,規(guī)范公式具有對時間的連續(xù)導(dǎo)數(shù),也就是收縮應(yīng)變率為:

(23)

將其代入式(15),即可求出各時段的降溫值,從而計算結(jié)構(gòu)收縮效應(yīng)。混凝土收縮徐變效應(yīng)綜合計算時,將兩種作用同時輸入。當(dāng)進(jìn)行施工階段、體系轉(zhuǎn)換的結(jié)構(gòu)計算時,應(yīng)注意各階段徐變、收縮齡期時間的協(xié)調(diào)。

4 算例驗證

4.1 規(guī)范方法算例

為了驗證綜合仿真方法的正確性,擬設(shè)計算例。由于要考慮整個時變過程,也為了便于理論計算,算例設(shè)計成方形受壓柱,截面邊長2 m,柱高10 m,采用C50混凝土,底部固結(jié)。分三階段向柱端施加均布荷載,每階段均增加5 MPa,不計方柱重力荷載。初期加載齡期14 d,二期為44 d,三期為134 d,整個計算時間長度360 d。采用3個工況計算收縮徐變效應(yīng),分別為階段荷載+收縮作用、階段荷載+徐變作用、階段荷載+徐變作用+收縮作用。

用理論方法和綜合仿真方法計算各工況的收縮、徐變和荷載效應(yīng)在計算時間長度的變化過程,并進(jìn)行比較。理論和綜合仿真計算均采用橋梁規(guī)范[21]提供的收縮徐變模型。環(huán)境年平均相對濕度取55%,收縮開始齡期3 d,徐變開始加載齡期為14 d。理論計算按天進(jìn)行,分別形成3個工況整個時間過程的收縮徐變效應(yīng)曲線。模擬計算利用ANSYS有限元軟件進(jìn)行,采用SOLID65單元模擬方柱實體(見圖2),施加階段面荷載,整個計算時間分為71個荷載步,按上文方法計算各參數(shù),利用APDL語言編制整個時變過程計算命令流。

圖2 混凝土柱ANSYS有限元模型Fig.2 Finite element model of concrete column in ANSYS

經(jīng)計算,在ANSYS中可得到柱頂位移時間變化曲線,其中工況3的位移時變曲線見圖3,該圖為軟件計算結(jié)果截圖,橫坐標(biāo)為持荷時間(d),縱坐標(biāo)為頂點豎向位移(cm),可以看出模擬計算結(jié)果具有明顯的階段加載特點,同時整個變化過程中,初期收縮徐變效應(yīng)變化較快,后期變化較慢,與混凝土收縮徐變效應(yīng)特性相符。

圖3 工況3頂點位移時變曲線Fig.3 Curve of displacement vs. time in case 3

圖4 各工況頂點位移時變過程對比(模擬/理論)Fig.4 Comparison of displacement of top point during time-varying process in different cases (simulating/theory)

圖4為3個工況理論計算與模擬計算的柱頂位移時間變化曲線,可以看出,3個工況階段加載效應(yīng)和收縮徐變效應(yīng)均很明顯,模擬計算與理論計算偏差在整個時變過程中前期略大于后期。由表1可以看出,前期最大偏差不超過3%,后期不大于0.6%,整個時變過程二者吻合很好。

從算例的各工況模擬計算時變過程結(jié)果對比可以看出,綜合仿真方法計算過程及精度均較好,可以計算具有復(fù)雜施工階段的收縮徐變時變過程綜合效應(yīng),容易銜接規(guī)范,使用簡便,很適合于設(shè)計分析。

表1 各工況及混凝土齡期頂點位移對比(模擬/理論)Tab.1 Comparisons of displacement of top point in different cases and ages simulation/theory (simulating/theoretical)

4.2 試驗研究算例

文獻(xiàn)[22]中對4片無黏結(jié)預(yù)應(yīng)力縮尺模型梁進(jìn)行了為期1年的徐變試驗,其中1#梁為普通混凝土,其余均為不同參量的粉煤灰高性能混凝土。在此取1#梁與粉煤灰參量40%的4#梁進(jìn)行模擬計算驗證。

模型梁依據(jù)鐵路標(biāo)準(zhǔn)圖中的32 m跨簡支梁制作,縮尺比例為1∶5,試驗梁長6.52 m,模擬計算梁取6.6 m,無黏結(jié)預(yù)應(yīng)力筋根據(jù)標(biāo)準(zhǔn)圖布置采用3根φ15.2鋼絞線,半跨模型梁立面如圖5所示,梁端及跨中截面如圖6所示。

圖5 模型梁立面示意圖(單位:mm)Fig.5 Schematic diagram of elevation of model beams (unit: mm)

圖6 模型梁截面示意圖(單位:mm)Fig.6 Schematic diagram of cross-section of model beams (unit: mm)

1#梁加載齡期為36 d,此時彈性模量為48.3 GPa;4#梁加載齡期為33 d,此時彈性模量為51.9 GPa。通過百分表測量的試驗梁跨中撓度時間變化,計算徐變系數(shù)實測值見圖7。文獻(xiàn)[22]中給出的徐變系數(shù)擬合公式如下:

(24)

式中,φ(∞,t0)為徐變系數(shù)終級值,1#梁為2.24,4#梁為1.76。兩片模型梁的徐變系數(shù)擬合值時程曲線見圖7,前期擬合值與實測值相差較大,后期較為接近。

圖7 模型梁徐變系數(shù)時程圖(擬合/實測)Fig.7 Curves of creep coefficient of model beams vs. time (fitting/test)

由于文獻(xiàn)[22]根據(jù)試驗提出了1#梁和4#梁的徐變系數(shù)擬合公式,該公式為時間的連續(xù)函數(shù),且具有連續(xù)導(dǎo)數(shù),導(dǎo)數(shù)見式(25),所以適合采用綜合仿真方法對結(jié)構(gòu)在整個時變過程中的徐變效應(yīng)進(jìn)行分析。

(25)

計算時,混凝土均采用SOLID65單元模擬,材料特性采用試驗值,由于試驗徐變系數(shù)已考慮鋼筋影響,所以模擬不計鋼筋,無黏結(jié)預(yù)應(yīng)力筋采用LINK10單元模擬,通過耦合設(shè)置來釋放預(yù)應(yīng)力鋼絞線切向約束,利用初應(yīng)力施加預(yù)應(yīng)力,計算模型見圖8。自初始加載開始,試驗持荷時間為365 d,分為64個荷載步,通過綜合仿真方法編制整個時變過程計算命令流。

圖8 模型梁ANSYS有限元模型Fig.8 Finite element model of model beam in ANSYS

利用文獻(xiàn)[22]提出的擬合徐變系數(shù),在無黏結(jié)預(yù)應(yīng)力作用下,利用本研究有限元綜合仿真方法計算的1#梁和4#梁跨中上拱度在1年內(nèi)的時變過程曲線如圖9所示,可以看出變化趨勢與徐變系數(shù)時間變化相似,明確反映了徐變效應(yīng)的規(guī)律。與兩片試驗?zāi)P土嚎缰猩瞎岸鹊膶崪y值進(jìn)行比較,可以看出,在初期模擬值與實測值相差較大,終期二者吻合較好,整個變化趨勢相似。因為擬合徐變系數(shù)初期擬合值與實測值相差較大,所以前期效應(yīng)也相差較大,如果根據(jù)實測值變化情況采用分段函數(shù)來擬合,徐變效應(yīng)在整個時變過程的吻合情況將會更好。

圖9 模型梁跨中上拱度時程圖(模擬/實測)Fig.9 Curves of upward camber in middle span of model beams vs. time (simulating/test)

從以上兩個算例可以看出,只要規(guī)范和試驗研究所提出的結(jié)構(gòu)混凝土收縮徐變模型具有對時間的函數(shù)關(guān)系,就可以利用綜合仿真方法計算結(jié)構(gòu)在整個時變過程的收縮徐變效應(yīng),因此,綜合仿真方法與規(guī)范和試驗研究均銜接方便,可廣泛地應(yīng)用在結(jié)構(gòu)試驗、設(shè)計和監(jiān)測等領(lǐng)域,具有良好的實用性。

5 結(jié)論

(1)在混凝土結(jié)構(gòu)的整個時變過程中,綜合仿真方法計算的收縮、徐變效應(yīng)與理論計算結(jié)果和試驗徐變效應(yīng)吻合較好,且結(jié)果精確度較高,說明方法完全可行。

(2)綜合仿真方法將復(fù)雜的收縮徐變和施工階段過程分成一系列的時段,作用過程轉(zhuǎn)化為時間函數(shù),適合有限元法中的荷載步分析,并容易編制成命令流,讓計算機(jī)自動計算,使用極簡便。

(3)只要收縮徐變模型是時間的函數(shù)或分段函數(shù),均可采用綜合仿真方法,尤其適合收縮徐變模型對時間有連續(xù)導(dǎo)數(shù)的情況,容易與規(guī)范及試驗研究相銜接,所以方法具有較廣泛的實用性。

猜你喜歡
有限元效應(yīng)混凝土
混凝土試驗之家
關(guān)于不同聚合物對混凝土修復(fù)的研究
鈾對大型溞的急性毒性效應(yīng)
懶馬效應(yīng)
混凝土預(yù)制塊模板在堆石混凝土壩中的應(yīng)用
混凝土,了不起
應(yīng)變效應(yīng)及其應(yīng)用
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 亚洲国产综合精品中文第一| 91丝袜乱伦| 五月婷婷亚洲综合| 91美女视频在线| 国产午夜福利亚洲第一| 日韩经典精品无码一区二区| 成人字幕网视频在线观看| 国产一级毛片高清完整视频版| 一本大道香蕉久中文在线播放| 免费福利视频网站| 人妻丰满熟妇啪啪| 久久久久国产一级毛片高清板| 国产黄在线观看| 国产在线第二页| 国产精品视频导航| 日韩东京热无码人妻| 在线观看av永久| 欧美色99| 亚洲伊人久久精品影院| 欧美日韩久久综合| 欧美19综合中文字幕| 午夜福利免费视频| 亚洲AⅤ永久无码精品毛片| 99精品国产高清一区二区| 亚洲swag精品自拍一区| 波多野结衣一二三| 一级成人欧美一区在线观看 | 亚洲一区黄色| 在线观看无码av免费不卡网站| 中日无码在线观看| 国产黄网站在线观看| 久久人人爽人人爽人人片aV东京热 | 国产91精品调教在线播放| 一区二区三区高清视频国产女人| 国产成+人+综合+亚洲欧美| 日韩人妻精品一区| 一级一毛片a级毛片| 亚洲男人的天堂视频| 激情综合网激情综合| 欧美成人综合在线| 自拍偷拍欧美| 无码在线激情片| 国产不卡网| av在线5g无码天天| 国产拍在线| 国产国产人免费视频成18| 亚洲69视频| 波多野结衣一区二区三区四区视频 | 亚洲男人在线天堂| 男女男免费视频网站国产| 精品久久久久无码| 日韩欧美中文| 色国产视频| 国产制服丝袜无码视频| 久久婷婷六月| 久久国产乱子伦视频无卡顿| 日韩毛片在线视频| 午夜日本永久乱码免费播放片| 亚洲成在线观看 | 国产午夜看片| 熟女日韩精品2区| 一级黄色片网| 亚洲美女操| 国产福利免费视频| 久久免费看片| 亚洲成人精品在线| 91无码网站| 亚洲av日韩综合一区尤物| 一本大道无码日韩精品影视| 天堂网亚洲系列亚洲系列| 日韩高清中文字幕| 国产理论最新国产精品视频| 萌白酱国产一区二区| 91久久偷偷做嫩草影院免费看 | 亚洲精品在线观看91| 国产在线观看成人91| 日韩黄色大片免费看| 国产在线高清一级毛片| 999精品在线视频| 亚洲成人福利网站| 日韩在线中文| 综合五月天网|