• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于統(tǒng)一相場理論的早齡期混凝土化-熱-力多場耦合裂縫模擬與抗裂性能預(yù)測1)

    2021-05-31 01:34:36吳建營陳萬昕黃羽立
    力學(xué)學(xué)報(bào) 2021年5期
    關(guān)鍵詞:相場齡期水化

    吳建營 陳萬昕 黃羽立

    ?(華南理工大學(xué)亞熱帶建筑科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,廣州 510641)

    ?(華南理工大學(xué)土木工程系,廣州 510641)

    ??(清華大學(xué)土木工程系,北京 100084)

    引言

    混凝土是全世界范圍內(nèi)應(yīng)用最為廣泛的土木工程材料,在建筑結(jié)構(gòu)、橋梁隧道、能源設(shè)施(水電大壩、核電安全殼、液化天然氣或石油儲罐、風(fēng)電基礎(chǔ)等)等土木工程和基礎(chǔ)設(shè)施建設(shè)中發(fā)揮著重要作用.然而,一旦混凝土結(jié)構(gòu)在建造期出現(xiàn)裂縫,在服役期載荷和環(huán)境作用下這些早齡期裂縫將逐漸擴(kuò)展和演化,對結(jié)構(gòu)全生命周期的完整性、耐久性和安全性造成嚴(yán)重影響.因此,準(zhǔn)確預(yù)測早齡期混凝土的抗裂性能并由此控制混凝土結(jié)構(gòu)建造質(zhì)量,是預(yù)防結(jié)構(gòu)先天病患、保證結(jié)構(gòu)優(yōu)生健康的重要手段.

    混凝土澆筑養(yǎng)護(hù)過程中,水泥水化過程會釋放熱量使得混凝土升溫,熱傳導(dǎo)和熱對流等機(jī)制將導(dǎo)致混凝土溫度達(dá)到峰值后逐漸降低.這一過程中,與溫度有關(guān)的熱變形以及與水化度有關(guān)的自收縮變形會引發(fā)較大的拉應(yīng)力,而此時(shí)混凝土力學(xué)性能如強(qiáng)度等仍然處于較低水平,導(dǎo)致早齡期混凝土開裂并逐漸擴(kuò)展和演化[1].上述過程涉及化學(xué)、熱學(xué)、力學(xué)等復(fù)雜的多場耦合效應(yīng),已有理論模型和數(shù)值方法尚不能定量預(yù)測早齡期混凝土的抗裂性能并準(zhǔn)確描述早齡期裂縫的演化過程,這也成為混凝土結(jié)構(gòu)全壽命性能調(diào)控研究的主要難點(diǎn)問題之一.

    早齡期混凝土裂縫理論模型和數(shù)值模擬方法的相關(guān)研究可以追溯至20 世紀(jì)90 年代.根據(jù)所采用的混凝土裂縫模型,大致可以分為以下4 類:

    (1)經(jīng)驗(yàn)或塑性力學(xué)模型.早期研究大多采用塑性模型描述開裂后的早齡期混凝土性能,如文獻(xiàn)[2-5]等.

    (2) 彌散裂縫模型.這方面的代表性研究工作包括De Borst 等[6-8]提出的早齡期混凝土彌散裂縫模型.

    (3)損傷力學(xué)模型.隨著混凝土損傷力學(xué)的興起,若干學(xué)者[9-14]將部分經(jīng)典的損傷模型推廣至早齡期混凝土.

    (4)微細(xì)觀力學(xué)模型.為了考慮混凝土細(xì)觀結(jié)構(gòu)特別是粗骨料對早齡期混凝土抗裂性能的影響,部分學(xué)者發(fā)展了相應(yīng)的微細(xì)觀力學(xué)模型[15-17].

    上述工作為研究早齡期混凝土力學(xué)性能奠定了良好的基礎(chǔ).然而,這些研究大多采用較為簡單的裂縫模型,也未能考慮裂縫演化與水化反應(yīng)、熱量傳輸?shù)冗^程之間的相互耦合,自然也很難描述早齡期混凝土裂縫的演化過程并量化開裂后材料的力學(xué)性能.

    近年來,法國學(xué)者Nguyen 等[18-20]采用脆性斷裂相場模型描述開裂混凝土的力學(xué)行為,并考慮裂縫與水化反應(yīng)、熱量傳輸之間的耦合效應(yīng),能夠較好地描述早齡期的裂縫演化過程.脆性斷裂相場模型根植于Francfort 和Marigo[21]提出的線彈性斷裂能量變分原理,解決了Griffith 能量方法需要預(yù)設(shè)裂縫路徑的難題;由于將尖銳裂縫正則化為具有一定尺度的裂縫帶,并引入在[0,1]之間連續(xù)分布的裂縫相場及其梯度描述裂縫狀態(tài),脆性斷裂相場模型非常便于通過有限元、無網(wǎng)格、物質(zhì)點(diǎn)等多種數(shù)值方法加以實(shí)現(xiàn).可以嚴(yán)格證明[22]:當(dāng)裂縫尺度趨近于零時(shí),脆性斷裂相場模型Γ-收斂于Griffith 線彈性斷裂力學(xué).

    然而,脆性斷裂相場模型僅適用于脆性材料,難以合理反映具有明顯軟化段的混凝土損傷破壞行為.更糟糕的是,模型的定量分析結(jié)果存在嚴(yán)重的裂縫尺度敏感性問題:裂縫尺度越小,結(jié)構(gòu)起裂載荷和峰值載荷越大.因此,若將裂縫尺度視為數(shù)值參數(shù),則此類模型與Griffith 線彈性斷裂力學(xué)類似,同樣無法描述完好或僅有弱奇異性固體的裂縫起裂.迫不得已,研究人員只好將裂縫尺度視為材料屬性[23-25],這一處理方式仍然無法完全解決裂縫起裂問題[26-27].

    近年來,筆者將斷裂力學(xué)和經(jīng)典損傷力學(xué)有機(jī)結(jié)合,建立了同時(shí)適用于脆性斷裂和準(zhǔn)脆性破壞的統(tǒng)一相場理論[28].從該理論出發(fā),不僅經(jīng)典的脆性斷裂相場模型[23,29]可以作為特例給出,還在國際上首次提出了一類相場內(nèi)聚裂縫模型PF-CZM[30-32].理論分析證明:當(dāng)裂縫尺度趨近于零時(shí),PF-CZM 收斂為一類混合型破壞的內(nèi)聚裂縫模型[33-34],能夠準(zhǔn)確預(yù)測線性軟化以及指數(shù)軟化、雙曲軟化、混凝土科內(nèi)列森軟化[35]等非線性軟化曲線.該模型同時(shí)考慮了基于強(qiáng)度的裂縫起裂準(zhǔn)則、基于能量的裂縫擴(kuò)展準(zhǔn)則以及基于變分原理的裂縫擴(kuò)展方向判據(jù),僅需少量標(biāo)準(zhǔn)材料參數(shù)即可定量預(yù)測工程結(jié)構(gòu)的損傷破壞行為;特別是,當(dāng)裂縫尺度小于某一上限時(shí),其取值對脆性斷裂和準(zhǔn)脆性破壞的分析結(jié)果幾乎不產(chǎn)生影響.歸功于上述優(yōu)點(diǎn),統(tǒng)一相場理論和相場內(nèi)聚裂縫模型一經(jīng)提出,迅速得到了國內(nèi)外學(xué)者的廣泛認(rèn)可和跟蹤研究[36-41].

    基于上述考慮,在統(tǒng)一相場理論框架內(nèi),本工作擬進(jìn)一步考慮裂縫相場與水化反應(yīng)、熱量傳輸?shù)冗^程之間的雙向耦合影響,建立裂縫相場演化特征與混凝土水化度和溫度之間的定量聯(lián)系,提出混凝土化-熱-力多場耦合相場內(nèi)聚裂縫模型,發(fā)展相應(yīng)的多場有限元數(shù)值實(shí)現(xiàn)算法,并應(yīng)用于早齡期混凝土的裂縫模擬和抗裂性能預(yù)測.

    1 化?熱?力多場耦合相場內(nèi)聚裂縫模型

    圖1 固體內(nèi)嵌裂縫/界面及其幾何正則化[34]Fig.1 A cracking solid and its geometric regularization[34]

    相場模型中,尖銳裂縫/界面S 被彌散為如圖1所示的裂縫帶B ?Ω,其尺度記為b>0,并引入相場變量d(x) :B →[0,1]描述裂縫狀態(tài);裂縫帶B 的外邊界記為?B,其外法向矢量表示為nB.需要指出的是,在固體損傷破壞過程中,裂縫帶B 并非預(yù)先設(shè)定或保持固定,而是遵循自身特定本構(gòu)行為發(fā)生起裂、擴(kuò)展和演化等過程.類似于位移場,裂縫相場也可以施加某種強(qiáng)迫邊界條件:例如,對于彈性區(qū)域有d(x)=0;對于初始預(yù)設(shè)裂縫有d(x)=1 等.

    1.1 混凝土開裂過程

    首先介紹描述混凝土開裂行為的相場內(nèi)聚裂縫模型,包括裂縫相場-位移場耦合方程、應(yīng)力-應(yīng)變關(guān)系以及本構(gòu)特征函數(shù)等.需要指出的是,早齡期混凝土的變形和力學(xué)性能與齡期或水化程度、溫度等因素有關(guān),裂縫建模時(shí)必須合理地加以考慮.

    1.1.1 裂縫相場-位移場耦合方程

    統(tǒng)一相場理論[28,32]中,準(zhǔn)靜態(tài)條件下開裂固體的行為由如下耦合方程描述

    其中,b?為體積均布力,裂縫相場通量q和裂縫相場源Q分別表示為

    式中,應(yīng)力σ 和裂縫有效驅(qū)動(dòng)力H 由下節(jié)的本構(gòu)關(guān)系給出;尺度b>0 為表征裂縫帶寬度的數(shù)值參數(shù),其值越小,裂縫帶寬越小,越接近尖銳裂縫;裂縫幾何函數(shù)α(d) ∈[0,1]和能量退化函數(shù)ω(d) ∈[0,1]分別為裂縫相場d的遞增和遞減函數(shù),并滿足條件[28,43]

    這里,引入了歸一化參數(shù)cα:=,以保證固體完全破壞時(shí)單位裂縫面積的耗能為材料斷裂能屬性.

    需要指出,上述相場模型隸屬于斷裂力學(xué)或梯度損傷力學(xué)領(lǐng)域的Bourdin 等[29]變分相場斷裂模型,與物理學(xué)領(lǐng)域常用的Ginzburg-Landau 相變方程[44]、Cahn-Hillard[45]或Allan-Cahn[46]相分離、晶體生長或多相流方程有明顯區(qū)別,雖然后者也曾經(jīng)被美國學(xué)者Karma 等[47]用于固體裂縫擴(kuò)展分析.

    變分相場斷裂模型的理論基礎(chǔ)來源于線彈性斷裂力學(xué)中的Griffith 能量原理[48]以及后續(xù)Francfort–Marigo 提出的變分原理[21],并借鑒了圖像分割理論中Mumford-Shah 泛函[49]的Ambrosio-Tortorelli 橢圓正則化[50].因此,變分相場斷裂模型通常滿足Γ-收斂性質(zhì),即:當(dāng)正則化尺度參數(shù)b趨近于零時(shí),其近似解在能量泛函意義上收斂于原尖銳裂縫問題的精確解,具體可參見Braide[51]和Bourdin 等[22]中的證明和論述.

    除應(yīng)用范圍不同外,以上變分相場斷裂模型與Ginzburg-Landau 方程主要有以下區(qū)別.

    (1)Ginzburg-Landau 方程中的裂縫幾何函數(shù)通常采用雙勢阱函數(shù)α(d)=d2(1-d)2,該函數(shù)滿足性質(zhì)α(0)=α(1)=0,即裂縫相場(或序參數(shù))處于中間過渡狀態(tài)d∈(0,1)時(shí)的裂縫耗散比純相d=0 或d=1時(shí)更大,這一特點(diǎn)與材料相變、多相流或結(jié)構(gòu)拓?fù)鋬?yōu)化等問題較為契合.然而,雖然雙勢阱函數(shù)也滿足Γ-收斂特性,但該函數(shù)在[0,1]之間不具有單調(diào)遞增性,也無法區(qū)分d=0 和d=1 這兩種截然不同的裂縫狀態(tài).因此,若將其作為裂縫幾何函數(shù),將會導(dǎo)致裂縫面積和能量耗散隨裂縫擴(kuò)展呈現(xiàn)先增大后減小的趨勢,這與固體斷裂時(shí)Griffith 能量原理(即應(yīng)變能持續(xù)減小、裂縫表面能單調(diào)增大、最終總能量達(dá)到最小值)相矛盾.在斷裂力學(xué)領(lǐng)域被廣泛認(rèn)可的變分相場斷裂模型中,其裂縫幾何函數(shù)必須為裂縫相場(或序參數(shù))d的單調(diào)遞增函數(shù),且滿足條件(2a),這與裂縫擴(kuò)展的不可逆特征一致,其具體形式將在后文給出.

    (2)Ginzburg-Landau 方程通常會引入與裂縫相場(或序參數(shù))一階時(shí)間導(dǎo)數(shù)有關(guān)的松弛項(xiàng)或黏性項(xiàng),當(dāng)松弛模量趨近于零(或黏性系數(shù)趨近無窮大) 時(shí),方程退化為準(zhǔn)靜態(tài)式(1b)的形式.因此,此類模型先天具有一定的動(dòng)力阻尼效應(yīng).類似地,對于動(dòng)態(tài)裂縫擴(kuò)展問題,變分相場斷裂模型中的控制方程(1b)1(下標(biāo)1 指式(1b)中的第1 式,下同),其等號右邊也可以引入與裂縫相場一階時(shí)間導(dǎo)數(shù)(裂縫擴(kuò)展速率)有關(guān)的粘性阻尼項(xiàng)或(和) 與裂縫相場二階時(shí)間導(dǎo)數(shù)(裂縫擴(kuò)展加速度)有關(guān)的慣性項(xiàng),以描述裂縫動(dòng)態(tài)擴(kuò)展時(shí)的應(yīng)變率效應(yīng),具體參見作者的相關(guān)工作[52-53].這里僅考慮早齡期混凝土裂縫的準(zhǔn)靜態(tài)擴(kuò)展過程,因此并未引入阻尼項(xiàng)和慣性項(xiàng).

    需要指出的是,這里所采用的相場內(nèi)聚裂縫模型基于作者2017 年提出的統(tǒng)一相場理論[28],雖也屬于變分相場斷裂范疇,但與脆性斷裂變分相場模型有本質(zhì)不同:歸功于后文給出的本構(gòu)特征函數(shù),相場內(nèi)聚裂縫模型Γ-收斂于Barenblatt 內(nèi)聚裂縫模型[54],而非傳統(tǒng)的Griffith 線彈性斷裂理論[48].因此,相場內(nèi)聚裂縫模型不僅同時(shí)適用于脆性斷裂,也適用于粘聚斷裂,且模型預(yù)測給出的整體響應(yīng)(包括裂縫路徑和破壞模式等定性結(jié)果以及峰值載荷和軟化段等定量結(jié)果)與裂縫正則化尺度無關(guān).

    1.1.2 應(yīng)力-應(yīng)變本構(gòu)關(guān)系

    為了考慮水化過程對混凝土變形的影響,通??梢詫⒖倯?yīng)變張量? 分解為力學(xué)應(yīng)變?m、熱膨脹應(yīng)變?t和自收縮應(yīng)變?a三部分之和的形式,即

    為簡單起見,熱膨脹應(yīng)變?t和自收縮應(yīng)變?a均假定為二階各向同性張量,分別表示為

    式中,1為二階單位張量;ηt為熱膨脹系數(shù),表征與當(dāng)前溫度θ 和初始溫度θ0之差有關(guān)的熱膨脹應(yīng)變;MacAuley 括號定義為〈x〉:=max(x,0);當(dāng)水化度χ 超過閾值χ0(本文統(tǒng)一取為χ0=0.05) 后,混凝土出現(xiàn)自收縮變形,其大小與系數(shù)ηa有關(guān);ξ(χ)為截?cái)嗑€性函數(shù),表示為[55]

    其中,混凝土最終水化度χ∞與水灰比w/c之間服從上述指數(shù)關(guān)系[56].

    早齡期混凝土裂縫通常處于受拉張開狀態(tài),因此可采用如下各向同性損傷本構(gòu)關(guān)系

    式中G0為彈性剛度張量,E0為楊氏模量;為有效應(yīng)力張量的最大主應(yīng)力,ft為材料破壞強(qiáng)度.由于齡期的影響,這里楊氏模量和材料破壞強(qiáng)度均為水化度的遞增函數(shù),詳見1.2.2 節(jié).

    1.1.3 本構(gòu)特征函數(shù)

    為了反映混凝土開裂引起的應(yīng)變軟化行為,相場內(nèi)聚裂縫模型采用如下裂縫幾何函數(shù)和能量退化函數(shù)[28,30-31]

    對于線性軟化和混凝土常用的科內(nèi)列森軟化曲線[35],模型給出的軟化曲線如圖2 所示;其他如指數(shù)軟化、雙曲線軟化等也可以類似給出,見文獻(xiàn)[57].

    圖2 相場內(nèi)聚裂縫模型PF-CZM 給出的軟化曲線[28,32]Fig.2 Softening curves given by the PF-CZM[28,32]

    前期研究表明[58]:裂縫尺度b的取值越小,裂縫面積和相應(yīng)能量耗散的計(jì)算精度就越高;極限情況下,裂縫尺度b→0 時(shí),裂縫帶寬度也趨近于零,模型退化為一類混合型的內(nèi)聚裂縫模型[34].另一方面,模型給出的應(yīng)力-名義位移跳躍曲線σ(w)與裂縫尺度b無關(guān),而僅取決于彈性模量E0、破壞強(qiáng)度ft、斷裂能Gf以及軟化段參數(shù)(即初始斜率k0、裂縫極限位移wc)等材料屬性.因此,對于裂縫幾何函數(shù)和能量退化函數(shù)(8),統(tǒng)一相場理論退化為一類相場內(nèi)聚裂縫模型.于是,相場內(nèi)聚裂縫模型不僅適用于各類準(zhǔn)脆性破壞,還可通過線性軟化描述脆性斷裂[31],實(shí)現(xiàn)了固體材料結(jié)構(gòu)兩種典型破壞特征的統(tǒng)一反映;詳見文獻(xiàn)[28,31-33]中的討論.

    1.2 混凝土水化過程

    早齡期混凝土養(yǎng)護(hù)過程中,水化反應(yīng)會釋放熱量而導(dǎo)致材料升溫膨脹;隨后,傳導(dǎo)和對流等熱傳輸過程導(dǎo)致溫度達(dá)到最高值后逐漸降低并引發(fā)收縮變形.

    一方面,裂縫會對混凝土水化反應(yīng)和熱傳輸過程產(chǎn)生影響;另一方面,混凝土力學(xué)性能也會隨著齡期或水化度而發(fā)生變化,共同構(gòu)成了一類典型的化-熱-力多場耦合問題,如圖3 所示.

    圖3 早齡期混凝土開裂過程中的化-熱-力多場耦合Fig.3 Multi-physical coupling in the cracking process of ear-age concrete

    1.2.1 水化度場-溫度場耦合控制方程

    早齡期混凝土水化過程中,化-熱耦合問題由如下傅里葉定律和阿倫尼烏斯定律描述

    式中,熱流(通量)h和水化速率˙χ 分別表示為

    其中,ρc為單位體積熱容;Θ 為水化反應(yīng)所釋放的熱量,p∞為水化反應(yīng)的潛熱;熱傳導(dǎo)系數(shù)假設(shè)為二階各向同性張量且與水化度χ 無關(guān),即κ0=κ01,其中κ0為初始熱傳導(dǎo)系數(shù);邊界?Ωh上給定的熱流記為h?,邊界?Ωc上與空氣溫度θa有關(guān)的對流/輻射系數(shù)記為g?,剩余邊界?Ωθ:=?ΩΩhΩc的溫度為給定值θ?;活化能Ea表征熱生成速率,R=8.314 J/(K ·mol)為理想氣體常數(shù);對于混凝土,水化反應(yīng)演化函數(shù)φ(χ)一般可由如下多項(xiàng)式函數(shù)擬合給出[13]

    系數(shù)bi(i=0,1,2,···,6) 根據(jù)等溫量熱法試驗(yàn)數(shù)據(jù)擬合確定,如圖4 所示.

    圖4 典型混凝土水化反應(yīng)演化函數(shù)[13]Fig.4 Evolution function of the hydration reaction[13]

    由式(11b)可見,早齡期裂縫演化會影響混凝土水化反應(yīng)和熱量傳輸過程.

    1.2.2 混凝土齡期效應(yīng)

    隨著齡期增加,混凝土表現(xiàn)出明顯的齡期效應(yīng),即彈性模量E0、破壞強(qiáng)度ft和斷裂能Gf均逐漸增大,而泊松比ν0基本保持不變.

    為簡化起見,這里假設(shè)早齡期混凝土開始形成彈性模量、破壞強(qiáng)度和斷裂能的水化度閾值與自收縮應(yīng)變的水化度閾值χ0取為相同值,即考慮如下混凝土齡期效應(yīng)[4,59-60]

    式中,函數(shù)ξ(χ)由式(5)給出;E∞,Gf∞和ft∞分別為水化度達(dá)到飽和χ∞時(shí)的彈性模量、斷裂能和抗拉強(qiáng)度.

    上述混凝土力學(xué)性能的齡期效應(yīng)也可以采用不同的水化度閾值或其他函數(shù)形式,如文獻(xiàn)[19].此時(shí),混凝土的特征長度lch將與水化度χ 有關(guān);相應(yīng)地,相場內(nèi)聚裂縫模型中,參數(shù)a1也將隨之發(fā)生變化.另一方面,參數(shù)a2和a3僅與軟化曲線類型有關(guān),不受材料力學(xué)參數(shù)變化的影響.

    2 數(shù)值實(shí)現(xiàn)和求解算法

    為便于有限元數(shù)值實(shí)現(xiàn),通常將化-熱-力多場耦合控制方程(1)和式(11)表示為如下弱形式

    式中,δθ,δχ,δu,δd分別為任意容許的節(jié)點(diǎn)溫度、水化度、位移和裂縫相場.與溫度和位移條件類似,裂縫相場d(x) 也可以施加諸如d(x∈ΩB)=0 或d(x∈S0)=1 等強(qiáng)制邊界條件,這里S0為預(yù)制裂縫.

    上述弱形式的化-熱-力多場耦合問題可以采用多場有限元方法進(jìn)行空間離散,并選用合理的數(shù)值算法進(jìn)行求解.

    2.1 多場有限元數(shù)值實(shí)現(xiàn)

    在進(jìn)行有限元空間離散時(shí),可以將整個(gè)計(jì)算區(qū)域內(nèi)所有單元節(jié)點(diǎn)均賦予溫度、水化度、位移和裂縫相場自由度.為提高計(jì)算效率,也可以(不一定必要)將整個(gè)計(jì)算域Ω 分為兩部分[58],即:可能出現(xiàn)裂縫的相場子區(qū)域B 和剩余部分ΩB.對于前者,單元節(jié)點(diǎn)同時(shí)具有溫度、水化度、位移和裂縫相場等自由度,而后者的單元節(jié)點(diǎn)僅考慮溫度、水化度和位移等自由度.已有研究表明[28,30-31,57,61],對于三角形常應(yīng)變單元,為保證計(jì)算精度,相場子區(qū)域B 內(nèi)的單元大小一般取為he≤b/5;四節(jié)點(diǎn)雙線性單元,其大小可適當(dāng)放寬至he≤b/3.

    有限元空間離散后,溫度、水化度、裂縫相場和位移等場變量及其梯度分別表示為(數(shù)值實(shí)現(xiàn)時(shí),張量采用Voigt 標(biāo)記)

    方程組(16) 通常采用增量迭代法進(jìn)行求解,即將時(shí)間步劃分為N個(gè)子步.對于某典型子步[Tn,Tn+1],其時(shí)間增量為ΔT:=Tn+1-Tn且Tn時(shí)刻所有的狀態(tài)變量已知,在此基礎(chǔ)上采用迭代方法求解Tn+1時(shí)刻的狀態(tài)變量.

    2.2 控制方程的數(shù)值求解

    考慮早齡期混凝土開裂過程的時(shí)間尺度較水化作用小得多,上述多場耦合問題可以視為化-熱耦合以及裂縫相場-位移耦合兩個(gè)串行過程,分別進(jìn)行求解后判斷是否收斂;如不收斂,則重新進(jìn)行迭代.上述數(shù)值求解過程如流程圖5 所示.

    圖5 化-熱-力多場耦合相場內(nèi)聚裂縫模型數(shù)值求解流程圖Fig.5 Flow chart in solving chemo-thermo-mechanical coupled problems

    2.2.1 水化過程

    對于水化過程,采用后退歐拉法進(jìn)行時(shí)間離散

    下標(biāo)n標(biāo)識上一增量步結(jié)束后的狀態(tài)變量.

    水化過程涉及的裂縫相場均采用上一增量步結(jié)束后的dn.由此,化-熱耦合方程(16)1和(16)2可采用牛頓迭代方法進(jìn)行求解,其線性化形式表示為

    式中,子剛度矩陣分別表示為

    上述剛度矩陣不對稱,這是耦合問題的特點(diǎn)之一.

    2.2.2 開裂過程

    對于早齡期混凝土開裂和演化過程,有效應(yīng)力ˉσ需考慮熱膨脹應(yīng)變?t和自收縮應(yīng)變?a,這些狀態(tài)變量可通過水化過程求解得到的溫度和水化度計(jì)算給出.相應(yīng)的,裂縫相場-位移場耦合控制方程的求解與純力學(xué)問題基本相同.

    在求解耦合方程(16)3和(16)4時(shí),整體牛頓迭代或子問題交錯(cuò)迭代等算法存在收斂性差、計(jì)算速度慢等問題.近年來,筆者[62-63]首次將整體BFGS 擬牛頓迭代算法[64]應(yīng)用于這一耦合方程的求解,在保證算法收斂性的同時(shí),收斂速度和計(jì)算效率大幅度提升.這里對該算法做簡要介紹.

    在整體BFGS 擬牛頓算法中,對于非線性方程組f(z)=0,其狀態(tài)變量z按如下格式進(jìn)行更新

    整體BFGS 擬牛頓算法可以根據(jù)需要設(shè)置一定迭代次數(shù)(5~10 次)間隔后才對剛度矩陣?K進(jìn)行更新計(jì)算.因此,雖然所需迭代次數(shù)可能比牛頓迭代算法更多(后者不一定能收斂),但BFGS 算法計(jì)算過程較為簡單且無需每次迭代都更新剛度矩陣,計(jì)算效率通常更高.較之子問題交錯(cuò)迭代算法,整體BFGS擬牛頓算法的求解效率提高5~8 倍[62-63,66-67].

    3 應(yīng)用算例

    針對上述化-熱-力多場耦合相場內(nèi)聚裂縫模型,本節(jié)給出其在早齡期混凝土裂縫演化模擬和抗裂性能預(yù)測的驗(yàn)證算例.對于混凝土,采用科內(nèi)列森軟化曲線.這里,二維問題均假定為平面應(yīng)力狀態(tài),網(wǎng)格劃分采用四節(jié)點(diǎn)雙線性插值單元.

    由于水化反應(yīng)過程的控制方程(11)涉及時(shí)間導(dǎo)數(shù),在求解過程中按照式(18) 進(jìn)行時(shí)間離散并需要給定節(jié)點(diǎn)溫度和節(jié)點(diǎn)水化度的初始值.算例1 中,節(jié)點(diǎn)溫度的初始值取為室溫20°C,算例2 中,節(jié)點(diǎn)溫度的初始值取為室溫21°C;算例1 和算例2 中,節(jié)點(diǎn)水化度的初始值均取為0.

    3.1 帶方形孔早齡期混凝土圓板收縮開裂[18]

    首先考慮如圖6 所示的早齡期混凝土圓形板,其外直徑為0.5 m,方形孔的尺寸為0.4 m.材料參數(shù)由文獻(xiàn)[13,18]給出(詳見表1);所采用混凝土水化反應(yīng)演化函數(shù)如圖4 所示.由于模型具有對稱性,建模時(shí)僅考慮右上角1/4 模型,并分別在模型左側(cè)和底部施加X向和Y向位移約束.

    圖6 帶方形孔圓板收縮開裂:試件幾何尺寸、邊界和環(huán)境條件[18]Fig.6 A circular plate with a square hole:Specimen geometry,boundary and environment conditions[18]

    表1 帶方形孔早齡期混凝土圓形板模型參數(shù)取值Table 1 Model parameters for the cracking analysis of an early-age circular concrete?plate with a square hole

    首先考慮裂縫尺度對數(shù)值模擬結(jié)果的影響.模擬中,采用b=5 mm 和b=10 mm 兩種裂縫尺度,相應(yīng)的單元大小取為he=b/5.模擬給出典型時(shí)刻的裂縫相場云圖如圖7 所示.可以看出,裂縫尺度b的取值大小僅改變損傷帶的寬度,而對裂縫路徑和破壞模式幾乎沒有任何影響.

    圖7 帶方形孔圓板收縮開裂:不同時(shí)刻的裂縫相場云圖(上排:b=10 mm;下排: b=5 mm)Fig.7 A circular plate with a square hole:Predicted damage profiles at various instants(top: b=10 mm;bottom: b=5 mm)

    圖8 給出了兩種尺度參數(shù)情況下圓形板內(nèi)、外邊緣和中間位置處混凝土的溫度、水化度、彈性模量、破壞強(qiáng)度以及斷裂能隨養(yǎng)護(hù)齡期的演化曲線.同樣,相場尺度參數(shù)b不影響上述定量計(jì)算結(jié)果.

    圖8 帶方形孔圓板收縮開裂:不同位置處溫度、水化度和典型力學(xué)參數(shù)隨時(shí)間的演化曲線Fig.8 A circular plate with a square hole:Predicted evolution curves of the temperatures,hydration degrees and mechanical properties

    對于b=5 mm 和he=1 mm 情況,圖9 中給出了不同時(shí)刻的裂縫相場、溫度場和水化度場云圖.可以看出,早齡期混凝土的裂縫演化過程主要分為以下3個(gè)階段:

    圖9 帶方形孔圓板收縮開裂:不同時(shí)刻的裂縫相場、溫度場和水化度場云圖Fig.9 A circular plate with a square hole:Predicted evolution curves of the temperatures and hydration degrees at various instants

    (1)時(shí)間段[0 h,16 h]:水化加速和裂縫起裂階段.該階段混凝土水化度快速增大,水化反應(yīng)釋放的熱量超過周邊環(huán)境對流和熱傳導(dǎo)引起的熱量損失,導(dǎo)致混凝土溫度升高,方孔邊緣處混凝土的最高溫度可達(dá)35.3°C.熱膨脹變形引起方形孔頂點(diǎn)處應(yīng)力集中,裂縫起裂首先出現(xiàn)此處.

    (2)時(shí)間段[22 h,80 h]:水化減速和裂縫緩慢發(fā)展階段.該階段水化速率降低,釋放的熱量小于環(huán)境對流和熱傳導(dǎo)引起的熱量損失,導(dǎo)致混凝土溫度降低,材料逐漸從熱膨脹變形轉(zhuǎn)變?yōu)槭湛s變形.由于水化作用,混凝土獲得了70%的破壞強(qiáng)度,足以抵抗熱膨脹和自收縮引起的裂縫擴(kuò)展,該階段裂縫發(fā)展緩慢.

    (3) 時(shí)間段[80 h,+∞]:水化減速和裂縫迅速擴(kuò)展并貫穿階段.該階段,自收縮應(yīng)變進(jìn)一步增大,方形孔頂點(diǎn)位置處的拉應(yīng)力逐漸增加并超過破壞強(qiáng)度,致使裂縫迅速擴(kuò)展并最終貫穿整個(gè)試件,應(yīng)力得到釋放.最后,混凝土獲得了80%左右的破壞強(qiáng)度,其他部位不再有裂縫進(jìn)一步發(fā)展.

    可以看出,在不同階段,混凝土熱膨脹變形和自收縮變形兩種機(jī)制相互競爭,并對早齡期混凝土的應(yīng)力發(fā)展和抗裂性能產(chǎn)生顯著影響.

    與純力學(xué)載荷作用下的分析結(jié)果類似,化-熱-力多場耦合相場內(nèi)聚裂縫模型的數(shù)值模擬結(jié)果不受裂縫尺度和網(wǎng)格大小影響,能夠合理描述早齡期混凝土水化反應(yīng)和熱量傳輸耦合作用下復(fù)雜的開裂過程,特別是熱膨脹變形和自收縮變形這兩類相互競爭的物理機(jī)制得到了準(zhǔn)確反映.

    3.2 受約束圓環(huán)試驗(yàn)收縮開裂試驗(yàn)

    考慮如圖10 所示的受約束圓環(huán)收縮開裂試驗(yàn)[68].

    圖10 約束混凝土圓環(huán)試驗(yàn):試件幾何尺寸和環(huán)境條件[68]Fig.10 Restrained concrete ring test:Specimen geometry and environment conditions[68]

    混凝土圓環(huán)試件的厚度為75 mm,內(nèi)直徑為150 mm,外直徑為300 mm.試驗(yàn)中考慮了3.1 mm,9.5 mm 和19.0 mm 3 種厚度的約束鋼環(huán)(其中,3.1 mm 厚度圓環(huán)試件由于剛度較小、約束效應(yīng)較低,混凝土試件并未開裂,故這里僅考慮后兩種厚度情況),并通過貼在鋼環(huán)中部位置的4 個(gè)應(yīng)變片監(jiān)測鋼環(huán)應(yīng)變.所采用的材料(水灰比w/c=0.3 的砂漿)和模型參數(shù)見表2;混凝土水化反應(yīng)演化函數(shù)如圖11 所示.

    表2 受約束圓環(huán)試驗(yàn)收縮開裂試驗(yàn)參數(shù)取值Table 2 Model parameters for the cracking analysis of a constrained concrete ring

    圖11 約束混凝土圓環(huán)試驗(yàn):混凝土水化反應(yīng)演化函數(shù)Fig.11 Restrained concrete ring test:Evolution function of the hydration reaction

    類似地,為考察裂縫尺度和網(wǎng)格大小的影響,數(shù)值模擬采用兩種裂縫尺度b=2.5 mm 和b=5.0 mm,相應(yīng)的單元網(wǎng)格大小取為he=b/5.

    圖12 和圖13 分別給出了9.5 mm 和19.0 mm 兩種厚度試件不同時(shí)刻的裂縫相場云圖.可以看出,在開始階段,由于試件的環(huán)向應(yīng)力均勻,微裂縫沿內(nèi)邊緣均勻分布;隨著齡期增長,收縮應(yīng)變增大,由于隨機(jī)性的影響(這里體現(xiàn)為網(wǎng)格劃分引起的數(shù)值誤差),部分微裂縫較其他裂縫擴(kuò)展更為明顯,但混凝土強(qiáng)度也逐漸提高,導(dǎo)致該階段總體上裂縫擴(kuò)展緩慢;隨著齡期進(jìn)一步增長,收縮應(yīng)變顯著增大,導(dǎo)致環(huán)向某一位置出現(xiàn)控制性宏觀裂縫并迅速擴(kuò)展至圓環(huán)試件外邊緣.

    圖12 約束混凝土圓環(huán)試驗(yàn):約束混凝土圓環(huán)試驗(yàn):兩種相場尺度參數(shù)條件下不同時(shí)刻的裂縫相場云圖模擬結(jié)果(圓環(huán)厚度19 mm)Fig.12 Restrained concrete ring test:Predicted contours of the crack phase-field at various time instants for different length scale parameters(thickness of the steel ring:19 mm)

    圖13 約束混凝土圓環(huán)試驗(yàn):兩種相場尺度參數(shù)條件下不同時(shí)刻的裂縫相場云圖模擬結(jié)果(圓環(huán)厚度為9.5 mm)Fig.13 Restrained concrete ring test:Predicted contours of the crack phase-field at various time instants for different length scale parameters(thickness of the steel ring:9.5 mm)

    數(shù)值模擬給出的9.5 mm 和19.0 mm 兩種厚度圓鋼環(huán)應(yīng)變演化曲線與試驗(yàn)結(jié)果的對比如圖14 所示.可以看出,考慮化-熱-力多場耦合效應(yīng)的相場內(nèi)聚裂縫模型能夠反映兩種不同厚度鋼環(huán)試件的應(yīng)變演化規(guī)律,特別是能夠準(zhǔn)確地預(yù)測混凝土開裂后、環(huán)向拉應(yīng)力得到釋放導(dǎo)致鋼環(huán)所受壓應(yīng)力迅速降低、壓應(yīng)變出現(xiàn)陡降等典型特征.

    圖14 約束混凝土圓環(huán)試驗(yàn):圓環(huán)應(yīng)變演化曲線Fig.14 Restrained concrete ring test:predicted evolution curves of the compressive strains of the steel ring at various time instants

    同樣,與純力學(xué)載荷作用下的結(jié)果類似,無論是定性的裂縫演化過程、亦或是定量的圓環(huán)壓應(yīng)變演化曲線,裂縫尺度和網(wǎng)格大小對預(yù)測結(jié)果幾乎不產(chǎn)生影響,這一結(jié)果對模型應(yīng)用于實(shí)際混凝土結(jié)構(gòu)早齡期裂縫演化過程的準(zhǔn)確模擬、抗裂性能的定量預(yù)測至關(guān)重要.

    4 結(jié)論與展望

    針對早齡期混凝土的抗裂性能預(yù)測和裂縫擴(kuò)展全過程模擬,本工作在統(tǒng)一相場理論的理論框架內(nèi),采用相場內(nèi)聚裂縫模型描述裂縫相場-位移場耦合效應(yīng),進(jìn)一步考慮了溫度引起的混凝土熱膨脹變形、水化反應(yīng)導(dǎo)致的混凝土收縮和隨齡期變化的力學(xué)性能(彈性模量、抗拉強(qiáng)度和斷裂能等),以及裂縫對水化過程和傳熱過程的影響,建立了化-熱-力多場耦合相場內(nèi)聚裂縫模型.

    與純力學(xué)載荷作用下的相場內(nèi)聚裂縫模型一樣,該模型同時(shí)融入了基于強(qiáng)度的裂縫起裂準(zhǔn)則、基于能量的裂縫擴(kuò)展準(zhǔn)則以及基于變分原理的裂縫擴(kuò)展方向判據(jù),因此能夠描述早齡期混凝土水化作用和熱量傳輸引起的復(fù)雜裂縫起裂和演化過程.不同的是,刻畫裂縫起裂的混凝土抗拉強(qiáng)度并非固定值,而是隨齡期增長而逐漸增大并趨于最終強(qiáng)度;類似地,隨齡期變化的混凝土彈性模量和斷裂能,也影響著早齡期混凝土的力學(xué)性能和裂縫演化行為.因此,不同階段的水化作用和熱量傳輸導(dǎo)致混凝土變形在熱膨脹和自收縮之間相互競爭,并對混凝土的早齡期抗裂性能產(chǎn)生顯著影響.

    上述化-熱-力多場耦合相場內(nèi)聚裂縫模型可以采用有限元等數(shù)值方法加以實(shí)現(xiàn),相應(yīng)單元節(jié)點(diǎn)的自由度除位移和裂縫相場外,還包括溫度和水化度.考慮到水化過程和傳熱過程的時(shí)間尺度遠(yuǎn)比開裂過程大,將化-熱-力多場耦合作用作用下早齡期混凝土的裂縫演化過程視為化學(xué)-熱學(xué)耦合以及位移-裂縫相場耦合兩個(gè)串行過程并分別進(jìn)行求解.前者采用標(biāo)準(zhǔn)的牛頓迭代方法進(jìn)行求解,而后者則利用筆者近期發(fā)展的整體BFGS 擬牛頓算法,實(shí)現(xiàn)收斂性和計(jì)算效率的雙提升.

    數(shù)值算例表明,所提出的化-熱-力多場耦合相場內(nèi)聚裂縫模型不僅能夠描述早齡期混凝土水化反應(yīng)和熱量傳輸條件下復(fù)雜的裂縫演化過程,特別是熱膨脹變形引起的裂縫起裂以及自收縮變形導(dǎo)致的裂縫擴(kuò)展這兩類相互競爭的裂縫演化機(jī)制,而且數(shù)值模擬結(jié)果不受裂縫尺度和網(wǎng)格大小的影響,實(shí)現(xiàn)了早齡期混凝土抗裂性能的準(zhǔn)確預(yù)測.

    需要指出的是,上述工作忽略了徐變變形(包括基本徐變和瞬態(tài)徐變)[20,69]以及水分揮發(fā)引起的干縮裂縫[1,70]對混凝土早齡期力學(xué)性能的影響,這將在后續(xù)研究中進(jìn)一步加以考慮.

    猜你喜歡
    相場齡期水化
    基于聲發(fā)射技術(shù)的早齡期混凝土斷裂性能
    塑鋼纖維混凝土早齡期力學(xué)性能研究
    基于子單元光滑有限元的混凝土相場損傷模型研究
    橋梁錨碇基礎(chǔ)水化熱仿真模擬
    低溫早強(qiáng)低水化放熱水泥漿體系開發(fā)
    鑄件凝固微觀組織仿真程序開發(fā)
    基于相場理論的瀝青自愈合微觀進(jìn)程與機(jī)理研究進(jìn)展
    石油瀝青(2018年1期)2018-04-12 07:31:51
    基于COMSOL的相場模擬研究
    科技視界(2017年8期)2017-07-31 10:31:17
    水溶液中Al3+第三水化層靜態(tài)結(jié)構(gòu)特征和27Al-/17O-/1H-NMR特性的密度泛函理論研究
    早齡期混凝土對大跨連續(xù)剛構(gòu)橋撓度的影響
    四川水泥(2016年8期)2016-07-18 12:06:31
    中文亚洲av片在线观看爽| 成人午夜高清在线视频| 亚洲欧洲国产日韩| 色噜噜av男人的天堂激情| av国产免费在线观看| 欧美激情久久久久久爽电影| 亚洲欧美精品综合久久99| av视频在线观看入口| 久久精品久久久久久噜噜老黄 | 身体一侧抽搐| 97超碰精品成人国产| 亚洲国产色片| 青春草亚洲视频在线观看| 午夜福利高清视频| 免费一级毛片在线播放高清视频| 免费黄色在线免费观看| or卡值多少钱| 哪个播放器可以免费观看大片| 三级国产精品欧美在线观看| 免费一级毛片在线播放高清视频| 亚洲av免费高清在线观看| 99热网站在线观看| 国产精品99久久久久久久久| 亚洲五月天丁香| 搡老妇女老女人老熟妇| 久久国产乱子免费精品| 欧美日韩一区二区视频在线观看视频在线 | 日韩av在线免费看完整版不卡| 国产在线男女| 亚洲成人久久爱视频| 免费观看人在逋| av播播在线观看一区| 天美传媒精品一区二区| 久久久亚洲精品成人影院| 精品国产一区二区三区久久久樱花 | 欧美97在线视频| 精品无人区乱码1区二区| 青青草视频在线视频观看| 男人舔女人下体高潮全视频| 国产精品乱码一区二三区的特点| 尤物成人国产欧美一区二区三区| 亚洲四区av| 青春草亚洲视频在线观看| 日韩av在线大香蕉| 国产淫片久久久久久久久| 久久精品久久久久久久性| 免费看av在线观看网站| 干丝袜人妻中文字幕| 国产精品美女特级片免费视频播放器| 日韩高清综合在线| 亚洲人与动物交配视频| 国产精品一二三区在线看| 亚洲精品一区蜜桃| 蜜桃亚洲精品一区二区三区| 69人妻影院| 国产精品精品国产色婷婷| 亚洲18禁久久av| 亚洲av电影在线观看一区二区三区 | 亚洲成人精品中文字幕电影| 亚洲精品乱码久久久久久按摩| 国产私拍福利视频在线观看| 亚洲精品日韩av片在线观看| 亚洲激情五月婷婷啪啪| ponron亚洲| 成人亚洲欧美一区二区av| av在线观看视频网站免费| 免费不卡的大黄色大毛片视频在线观看 | 亚洲av成人精品一区久久| 国产亚洲精品久久久com| 欧美3d第一页| 成人一区二区视频在线观看| 欧美日韩综合久久久久久| av播播在线观看一区| 国产成人精品婷婷| 91精品伊人久久大香线蕉| 又爽又黄无遮挡网站| 欧美高清成人免费视频www| 日本黄大片高清| www.色视频.com| 一个人观看的视频www高清免费观看| 99在线视频只有这里精品首页| 国产久久久一区二区三区| 色5月婷婷丁香| 亚洲电影在线观看av| 日韩中字成人| 有码 亚洲区| 性色avwww在线观看| 色视频www国产| 日韩视频在线欧美| 全区人妻精品视频| 亚洲一级一片aⅴ在线观看| 免费大片18禁| 熟女人妻精品中文字幕| 久久久a久久爽久久v久久| 天天一区二区日本电影三级| 日韩一本色道免费dvd| 精品午夜福利在线看| 综合色丁香网| 国产淫片久久久久久久久| 亚洲在线自拍视频| 两性午夜刺激爽爽歪歪视频在线观看| 免费无遮挡裸体视频| 五月伊人婷婷丁香| 偷拍熟女少妇极品色| 中文乱码字字幕精品一区二区三区 | 我要搜黄色片| 水蜜桃什么品种好| 真实男女啪啪啪动态图| 日韩欧美在线乱码| 久久久午夜欧美精品| 搡女人真爽免费视频火全软件| 日本av手机在线免费观看| 精品人妻熟女av久视频| 男女那种视频在线观看| 国产精品不卡视频一区二区| 国产极品精品免费视频能看的| 日韩精品有码人妻一区| 有码 亚洲区| 一区二区三区乱码不卡18| 亚洲久久久久久中文字幕| 国产午夜精品论理片| 国产精品一区二区性色av| 草草在线视频免费看| 天天躁日日操中文字幕| 蜜桃亚洲精品一区二区三区| 一个人看视频在线观看www免费| 国产三级中文精品| 久久婷婷人人爽人人干人人爱| 国产精品日韩av在线免费观看| 亚洲国产色片| 亚洲久久久久久中文字幕| 欧美日韩国产亚洲二区| 午夜日本视频在线| 国产精品一区二区三区四区久久| 日本一二三区视频观看| 免费观看精品视频网站| 国产成人精品久久久久久| www.色视频.com| 99久久精品热视频| 天天躁日日操中文字幕| 高清在线视频一区二区三区 | 欧美日韩国产亚洲二区| 两个人的视频大全免费| 亚洲精品色激情综合| 97热精品久久久久久| 欧美高清性xxxxhd video| 大香蕉久久网| 亚洲av男天堂| 欧美性猛交╳xxx乱大交人| 国产人妻一区二区三区在| 在线免费观看不下载黄p国产| 能在线免费观看的黄片| 成年免费大片在线观看| 国产欧美日韩精品一区二区| 国产精品一区二区性色av| 亚洲av男天堂| 亚洲综合精品二区| 欧美潮喷喷水| 亚洲精品亚洲一区二区| 偷拍熟女少妇极品色| 特大巨黑吊av在线直播| 真实男女啪啪啪动态图| 免费观看在线日韩| 97人妻精品一区二区三区麻豆| 国产黄色小视频在线观看| 亚洲欧美精品综合久久99| 午夜免费男女啪啪视频观看| av在线观看视频网站免费| 两个人视频免费观看高清| 18禁动态无遮挡网站| 国产又色又爽无遮挡免| 大香蕉久久网| 亚洲精品aⅴ在线观看| 国产 一区精品| 色综合站精品国产| 国产在视频线在精品| 中文字幕精品亚洲无线码一区| 亚洲欧美一区二区三区国产| 精品欧美国产一区二区三| 国产在视频线精品| 精品久久久久久久末码| 免费观看性生交大片5| av国产久精品久网站免费入址| 午夜亚洲福利在线播放| 亚洲欧洲国产日韩| 久久99精品国语久久久| 亚洲欧美精品专区久久| 日本午夜av视频| 国产私拍福利视频在线观看| 亚洲精品乱久久久久久| 99在线人妻在线中文字幕| 有码 亚洲区| 18禁在线播放成人免费| 一夜夜www| 久久久久九九精品影院| 成人性生交大片免费视频hd| 春色校园在线视频观看| av在线播放精品| av国产久精品久网站免费入址| 激情 狠狠 欧美| 中文字幕熟女人妻在线| 成年女人永久免费观看视频| 国产高清国产精品国产三级 | 男女边吃奶边做爰视频| 午夜老司机福利剧场| 婷婷色av中文字幕| 国产免费男女视频| 亚洲成人精品中文字幕电影| 乱人视频在线观看| 在现免费观看毛片| 国产久久久一区二区三区| 欧美日本视频| 欧美精品一区二区大全| 色噜噜av男人的天堂激情| 亚洲熟妇中文字幕五十中出| 免费观看的影片在线观看| 99热这里只有精品一区| 免费人成在线观看视频色| 中文亚洲av片在线观看爽| 国产私拍福利视频在线观看| 久久久久九九精品影院| 欧美成人精品欧美一级黄| 99热这里只有精品一区| 尤物成人国产欧美一区二区三区| 国产在视频线精品| 一本一本综合久久| 99在线视频只有这里精品首页| 亚洲国产精品成人久久小说| 最近手机中文字幕大全| 好男人视频免费观看在线| 亚洲高清免费不卡视频| 大香蕉97超碰在线| 久久久久性生活片| 啦啦啦啦在线视频资源| 国产91av在线免费观看| 一级爰片在线观看| 特级一级黄色大片| 黄色日韩在线| 亚洲欧美日韩无卡精品| 日韩强制内射视频| 亚洲欧洲国产日韩| 欧美高清性xxxxhd video| 亚洲丝袜综合中文字幕| 久久久久免费精品人妻一区二区| 久久久久国产网址| 亚洲欧美精品自产自拍| 3wmmmm亚洲av在线观看| 长腿黑丝高跟| 久久久久久大精品| 欧美激情在线99| 日韩av在线大香蕉| 日本爱情动作片www.在线观看| 淫秽高清视频在线观看| a级毛片免费高清观看在线播放| 成人漫画全彩无遮挡| av卡一久久| 久久国内精品自在自线图片| 最近2019中文字幕mv第一页| 老司机福利观看| 一个人免费在线观看电影| 少妇猛男粗大的猛烈进出视频 | 大香蕉久久网| 国产真实伦视频高清在线观看| 婷婷六月久久综合丁香| 日本爱情动作片www.在线观看| 最近2019中文字幕mv第一页| 国产v大片淫在线免费观看| 乱码一卡2卡4卡精品| 日日摸夜夜添夜夜添av毛片| 国产精品1区2区在线观看.| h日本视频在线播放| 国产免费一级a男人的天堂| 2021天堂中文幕一二区在线观| 日韩精品有码人妻一区| 少妇的逼好多水| 麻豆国产97在线/欧美| 18禁在线播放成人免费| 成人无遮挡网站| 成年免费大片在线观看| 中文在线观看免费www的网站| 免费一级毛片在线播放高清视频| 国产伦理片在线播放av一区| 国内精品一区二区在线观看| 欧美三级亚洲精品| 中文欧美无线码| 亚洲自偷自拍三级| 亚洲国产精品成人综合色| 日韩亚洲欧美综合| 91精品国产九色| 免费观看精品视频网站| 国产伦在线观看视频一区| 亚洲中文字幕日韩| 少妇人妻一区二区三区视频| 美女黄网站色视频| 久久99热这里只频精品6学生 | 亚洲av不卡在线观看| 亚洲色图av天堂| 国产精品无大码| 嫩草影院精品99| 国产精品一区www在线观看| 成人午夜精彩视频在线观看| 国产三级中文精品| 国产国拍精品亚洲av在线观看| 国产亚洲精品久久久com| 亚洲人成网站高清观看| 中文亚洲av片在线观看爽| 有码 亚洲区| 免费黄色在线免费观看| 赤兔流量卡办理| 亚洲经典国产精华液单| 国产精品一及| 老司机影院成人| 亚洲成av人片在线播放无| 亚洲真实伦在线观看| 国产亚洲最大av| 非洲黑人性xxxx精品又粗又长| 尾随美女入室| 亚洲精品,欧美精品| 视频中文字幕在线观看| 午夜免费激情av| 亚洲av电影不卡..在线观看| 亚洲欧美成人精品一区二区| 欧美日韩一区二区视频在线观看视频在线 | 亚洲精品456在线播放app| 免费一级毛片在线播放高清视频| 人人妻人人澡人人爽人人夜夜 | 国产大屁股一区二区在线视频| 久久久久免费精品人妻一区二区| 亚洲av日韩在线播放| 99热网站在线观看| 美女xxoo啪啪120秒动态图| 99久久精品国产国产毛片| 久久久久久久国产电影| 麻豆成人av视频| 国产三级在线视频| 一级毛片电影观看 | 成年av动漫网址| 国内精品一区二区在线观看| 亚洲欧美精品专区久久| 亚洲精品国产成人久久av| 国产成人aa在线观看| 看片在线看免费视频| 欧美精品一区二区大全| 91久久精品国产一区二区三区| 桃色一区二区三区在线观看| 大又大粗又爽又黄少妇毛片口| 国产免费男女视频| 久久久久久久久久久免费av| 日韩中字成人| 日本爱情动作片www.在线观看| 麻豆精品久久久久久蜜桃| 搞女人的毛片| 少妇被粗大猛烈的视频| 天美传媒精品一区二区| 超碰97精品在线观看| 欧美日韩国产亚洲二区| 国产真实伦视频高清在线观看| 欧美性感艳星| 成年免费大片在线观看| 伊人久久精品亚洲午夜| 人人妻人人澡欧美一区二区| 九九热线精品视视频播放| 熟女人妻精品中文字幕| 久久久久久久久中文| 一区二区三区高清视频在线| 国产亚洲午夜精品一区二区久久 | 亚洲精品亚洲一区二区| 国产国拍精品亚洲av在线观看| av女优亚洲男人天堂| 亚洲国产精品合色在线| 一级毛片电影观看 | 国产极品精品免费视频能看的| 一夜夜www| 国产淫语在线视频| 亚洲怡红院男人天堂| 久热久热在线精品观看| 亚洲精品日韩在线中文字幕| 免费观看的影片在线观看| 亚洲,欧美,日韩| 热99在线观看视频| 91精品伊人久久大香线蕉| 看片在线看免费视频| 国产精品av视频在线免费观看| 99热6这里只有精品| av.在线天堂| 亚洲av熟女| 成年av动漫网址| 日本免费一区二区三区高清不卡| 91久久精品国产一区二区成人| 日韩 亚洲 欧美在线| 干丝袜人妻中文字幕| 人体艺术视频欧美日本| 黑人高潮一二区| 亚洲美女搞黄在线观看| 国产免费一级a男人的天堂| 亚洲av一区综合| 啦啦啦观看免费观看视频高清| 少妇人妻一区二区三区视频| 18禁在线无遮挡免费观看视频| 亚洲欧美精品综合久久99| 级片在线观看| a级一级毛片免费在线观看| 身体一侧抽搐| 久久久午夜欧美精品| 美女大奶头视频| 亚洲av福利一区| 亚洲av一区综合| 午夜福利视频1000在线观看| 精品久久久久久成人av| 精品国内亚洲2022精品成人| 国产高清不卡午夜福利| 久久精品人妻少妇| 久久这里只有精品中国| 精品久久久久久久久av| 人妻少妇偷人精品九色| a级一级毛片免费在线观看| 亚洲精品久久久久久婷婷小说 | 欧美bdsm另类| 好男人视频免费观看在线| 一卡2卡三卡四卡精品乱码亚洲| 波野结衣二区三区在线| 美女黄网站色视频| 女人十人毛片免费观看3o分钟| 在线免费观看不下载黄p国产| 国产av一区在线观看免费| 91精品一卡2卡3卡4卡| 九九在线视频观看精品| 国产 一区精品| 久久久久久久午夜电影| 啦啦啦啦在线视频资源| 国产成人freesex在线| 久久精品夜夜夜夜夜久久蜜豆| 18禁在线无遮挡免费观看视频| 国产又黄又爽又无遮挡在线| 精品人妻偷拍中文字幕| 久久久久久大精品| 欧美成人午夜免费资源| 国产免费一级a男人的天堂| 噜噜噜噜噜久久久久久91| 干丝袜人妻中文字幕| 久久久久久久久久久免费av| 99热网站在线观看| 精品人妻熟女av久视频| 日本熟妇午夜| 亚洲av.av天堂| 成人性生交大片免费视频hd| 亚洲,欧美,日韩| 日日撸夜夜添| 午夜爱爱视频在线播放| 高清视频免费观看一区二区 | 少妇的逼水好多| 国产日韩欧美在线精品| 午夜福利在线在线| 好男人视频免费观看在线| 不卡视频在线观看欧美| 日韩强制内射视频| 国产精品野战在线观看| 久久婷婷人人爽人人干人人爱| 精品欧美国产一区二区三| 中文字幕av在线有码专区| 男人舔女人下体高潮全视频| 1024手机看黄色片| 内射极品少妇av片p| 久久久久久久久久黄片| 少妇人妻精品综合一区二区| 大香蕉久久网| 欧美精品国产亚洲| 免费看日本二区| 99视频精品全部免费 在线| 成人特级av手机在线观看| 午夜福利成人在线免费观看| 欧美最新免费一区二区三区| 久久精品影院6| 观看免费一级毛片| 日韩制服骚丝袜av| 又爽又黄无遮挡网站| 国产极品精品免费视频能看的| 69人妻影院| 久久久久免费精品人妻一区二区| 亚洲精品色激情综合| 久久久国产成人精品二区| 男女啪啪激烈高潮av片| 亚洲欧美日韩高清专用| 全区人妻精品视频| 一边摸一边抽搐一进一小说| 人妻制服诱惑在线中文字幕| 舔av片在线| 国产成人a区在线观看| 亚洲伊人久久精品综合 | 秋霞伦理黄片| 美女大奶头视频| 国产精品久久视频播放| 成年av动漫网址| 99久久九九国产精品国产免费| 最近最新中文字幕免费大全7| 夫妻性生交免费视频一级片| 91av网一区二区| 一区二区三区乱码不卡18| 小蜜桃在线观看免费完整版高清| 只有这里有精品99| 免费搜索国产男女视频| av在线老鸭窝| 纵有疾风起免费观看全集完整版 | 亚洲图色成人| 欧美区成人在线视频| 在线免费观看的www视频| 男女视频在线观看网站免费| 欧美高清性xxxxhd video| 色噜噜av男人的天堂激情| 国产精品一二三区在线看| 亚洲av日韩在线播放| 亚洲国产最新在线播放| 国产成人freesex在线| 久久综合国产亚洲精品| 看片在线看免费视频| 欧美日韩在线观看h| 免费在线观看成人毛片| 寂寞人妻少妇视频99o| 亚洲精品影视一区二区三区av| 97超碰精品成人国产| 国产成人精品婷婷| 又爽又黄a免费视频| 哪个播放器可以免费观看大片| 亚洲欧美一区二区三区国产| 小说图片视频综合网站| 最近视频中文字幕2019在线8| 又黄又爽又刺激的免费视频.| 精品人妻视频免费看| 午夜爱爱视频在线播放| 美女脱内裤让男人舔精品视频| av在线播放精品| 久久久久久久久久成人| 久久久久久国产a免费观看| 十八禁国产超污无遮挡网站| 日韩精品有码人妻一区| 久久久国产成人精品二区| 午夜福利高清视频| 一区二区三区乱码不卡18| 人妻制服诱惑在线中文字幕| 免费观看在线日韩| 一个人看视频在线观看www免费| 国产精品爽爽va在线观看网站| 日韩视频在线欧美| 在线免费观看的www视频| 高清在线视频一区二区三区 | 草草在线视频免费看| 桃色一区二区三区在线观看| 精品人妻偷拍中文字幕| 大香蕉97超碰在线| av天堂中文字幕网| 日日摸夜夜添夜夜爱| 免费av观看视频| 国产精品三级大全| 欧美一级a爱片免费观看看| 亚洲怡红院男人天堂| 特大巨黑吊av在线直播| 亚洲欧美日韩无卡精品| 国产成人aa在线观看| 大话2 男鬼变身卡| 欧美成人一区二区免费高清观看| 亚洲国产精品合色在线| 午夜福利在线观看吧| 寂寞人妻少妇视频99o| 国产美女午夜福利| 亚洲精品久久久久久婷婷小说 | 欧美日本视频| 天美传媒精品一区二区| 成人高潮视频无遮挡免费网站| 一夜夜www| h日本视频在线播放| 两个人视频免费观看高清| 亚洲,欧美,日韩| 一级毛片久久久久久久久女| 久久精品国产亚洲av涩爱| 一区二区三区四区激情视频| 一个人看视频在线观看www免费| 日本免费a在线| 九草在线视频观看| 日韩成人av中文字幕在线观看| 国产伦理片在线播放av一区| 我要看日韩黄色一级片| 99热精品在线国产| 久久99热这里只有精品18| av黄色大香蕉| 日韩大片免费观看网站 | 亚洲国产精品合色在线| 91在线精品国自产拍蜜月| 久久久久国产网址| 一本久久精品| 欧美三级亚洲精品| 日本欧美国产在线视频| 亚洲精品aⅴ在线观看| 水蜜桃什么品种好| 日本五十路高清| 亚洲一区高清亚洲精品| 3wmmmm亚洲av在线观看| 成年av动漫网址| 日本wwww免费看| 99热精品在线国产| 淫秽高清视频在线观看| 91狼人影院| 国产免费视频播放在线视频 | kizo精华| 99九九线精品视频在线观看视频| 久久亚洲国产成人精品v| 免费看美女性在线毛片视频| 岛国毛片在线播放| 97热精品久久久久久| 寂寞人妻少妇视频99o| 国产老妇伦熟女老妇高清| 在线观看av片永久免费下载| 麻豆av噜噜一区二区三区| 最新中文字幕久久久久| 久久久a久久爽久久v久久| 国产av不卡久久|