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

    基于離散宏單元的砌體結(jié)構(gòu)高效非線性分析方法

    2022-08-01 00:58:16余丁浩
    工程力學(xué) 2022年8期
    關(guān)鍵詞:結(jié)點(diǎn)砌體彈簧

    郭 勇,余丁浩,李 鋼

    (大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,遼寧,大連 116024)

    砌體結(jié)構(gòu)因其造價(jià)低廉、施工方便、性能優(yōu)異而廣泛存在于工業(yè)與民用建筑中,此類結(jié)構(gòu)的主要材料為塊體和砂漿組成的非均質(zhì)復(fù)合材料,其力學(xué)行為由塊體、砂漿和二者之間粘結(jié)界面的力學(xué)性能共同確定。塊體和砂漿為典型的脆性材料,抗拉強(qiáng)度遠(yuǎn)小于抗壓強(qiáng)度。二者之間粘結(jié)界面的力學(xué)性能差,界面受拉行為與材料粘性屬性有關(guān),界面剪切行為同時(shí)與材料粘性和摩擦屬性有關(guān)。在極端環(huán)境荷載作用下,砌體結(jié)構(gòu)將表現(xiàn)出顯著的非線性特征。

    目前,砌體結(jié)構(gòu)非線性分析模型按照建模方式的不同可分為分離式和整體式兩種。分離式模型將砂漿、塊體及二者粘結(jié)界面分開建模,各部分的材料性質(zhì)通過描述各自力學(xué)性能的本構(gòu)關(guān)系表征,模型精細(xì)化程度高,但計(jì)算代價(jià)大。Louren?o和Rots[1]在塊體-砂漿粘結(jié)界面模型中采用復(fù)合屈服準(zhǔn)則考慮砌體墻的拉、壓、剪非線性行為,并在非線性行為中考慮了材料的塑性軟化特征;Sandoval 等[2]利用該模型完成了多孔黏土磚砌體墻和灌漿配筋砌體墻[3]的精細(xì)化數(shù)值模擬,與試驗(yàn)結(jié)果的對(duì)比分析結(jié)果顯示該模型能夠反映砌體的真實(shí)破壞形態(tài); Louren?o[4]和Roca 等[5]對(duì)砌體結(jié)構(gòu)非線性分析的數(shù)值計(jì)算問題進(jìn)行了綜述,認(rèn)為分離式模型最突出的問題在于計(jì)算代價(jià)高昂,嚴(yán)重地限制了其在大型工程結(jié)構(gòu)中的應(yīng)用;Minga等[6]針對(duì)分離式模型應(yīng)用于大型砌體結(jié)構(gòu)計(jì)算代價(jià)高昂的缺點(diǎn),采取將結(jié)構(gòu)分解為若干部分,分配給不同處理器并行計(jì)算的方式進(jìn)行非線性分析。整體式分析模型將塊體和砂漿假定為連續(xù)的勻質(zhì)體,建模過程簡(jiǎn)單,非線性分析所需計(jì)算量小,此類模型屬于宏觀分析模型,可以較好地平衡結(jié)構(gòu)分析的計(jì)算效率和精度,其關(guān)鍵在于建立能夠合理反映結(jié)構(gòu)實(shí)際震害特點(diǎn)和關(guān)鍵失效模式的簡(jiǎn)化建模策略,當(dāng)前已有眾多學(xué)者開展了此方面研究,并提出了多種整體式分析模型,如勻質(zhì)化模型[7-9]、等效框架模型[10-18]等,Grande 等[19]在梁?jiǎn)卧亩瞬吭O(shè)置剪切界面,提出了可考慮砌體墻肢剪切失效的梁?jiǎn)卧P?,將所提模型與節(jié)點(diǎn)區(qū)剛性單元聯(lián)合使用模擬砌體墻的彎曲和剪切非線性行為;Addessi 等以Timoshenko 梁理論為基礎(chǔ),分別在分布式非線性梁?jiǎn)卧P蚚20]和集中非線性梁?jiǎn)卧P蚚21]中設(shè)置塑性剪切鉸,并用該塑性鉸模擬墻肢的剪切失效機(jī)制;Caliò等[22]提出了一種考慮砌體墻面內(nèi)主要失效模式的平面離散宏單元模型,該模型用非線性彈簧構(gòu)造剪切單元和無厚界面單元,通過耦合兩種單元的變形模擬砌體墻的彎曲和剪切非線性行為,目前,該模型經(jīng)不斷改進(jìn)被應(yīng)用于無筋砌體墻在靜力作用下[23]和動(dòng)力作用下[24]的平面外受力特性分析、填充墻鋼筋混凝土框架結(jié)構(gòu)的平面內(nèi)抗震分析[25-26]和平面外抗震分析[27-28]、無筋砌體結(jié)構(gòu)的地震易損性分析[29]等,研究結(jié)果均表明該模型具有良好的分析精度和計(jì)算效率。綜上所述,精細(xì)化的分離式模型適用于結(jié)構(gòu)構(gòu)件的細(xì)部分析和模擬,而整體式模型在大型砌體結(jié)構(gòu)的分析和模擬方面具有更大優(yōu)勢(shì)。

    無論整體模型或分離模型,通常采用傳統(tǒng)變剛度法進(jìn)行非線性求解,不可避免地需對(duì)切線剛度矩陣進(jìn)行實(shí)時(shí)更新和分解,這導(dǎo)致隨問題規(guī)模增大計(jì)算效率急劇降低。為了避免結(jié)構(gòu)大規(guī)模整體切線剛度矩陣的實(shí)時(shí)更新和分解,國(guó)內(nèi)外學(xué)者利用土木工程結(jié)構(gòu)材料非線性的局部化特征,發(fā)展出一系列高效數(shù)值分析方法,例如多尺度方法、子結(jié)構(gòu)方法、重分析方法等。其中,多尺度有限元法通過構(gòu)造多尺度基函數(shù),將局部非線性區(qū)域的精細(xì)化模型與其余部位的宏觀模型進(jìn)行有效耦合,避免了傳統(tǒng)有限元法采用單一小尺度建模,大幅縮減結(jié)構(gòu)自由度,Drosopoulos 和Stavroulakis[30]、 Giambanco 等[31]、 Silva 等[32]和Krejci 等[33]采用多尺度方法對(duì)砌體結(jié)構(gòu)進(jìn)行了非線性分析,分析結(jié)果證明該方法計(jì)算效率明顯優(yōu)于傳統(tǒng)單一尺度建模方法。Clough 和Wilson[34]提出子結(jié)構(gòu)方法,該方法依據(jù)邊界力和變形協(xié)調(diào)條件將整體結(jié)構(gòu)劃分為線彈性子結(jié)構(gòu)和非線性子結(jié)構(gòu),并通過靜力凝聚的方式消去線彈性子結(jié)構(gòu)的內(nèi)部自由度,大幅降低結(jié)構(gòu)自由度數(shù)目,提高了計(jì)算效率,但該方法需要預(yù)判非線性子結(jié)構(gòu)的劃分區(qū)域,且彈、塑性子結(jié)構(gòu)一旦劃分就無法改變,因而適用性有限。為克服其局限性,孫寶印和古泉等[35-36]對(duì)整體結(jié)構(gòu)的單元彈塑性狀態(tài)進(jìn)行實(shí)時(shí)判斷,僅將進(jìn)入非線性狀態(tài)的構(gòu)件作為隔離塑性子結(jié)構(gòu),使大規(guī)模的非線性分析問題轉(zhuǎn)化為整體結(jié)構(gòu)的線彈性分析和規(guī)模較小的局部隔離子結(jié)構(gòu)的非線性分析。重分析方法是一類利用初始結(jié)構(gòu)計(jì)算信息快速求解修改后結(jié)構(gòu)響應(yīng)的計(jì)算方法,能夠避免結(jié)構(gòu)局部修改后進(jìn)行耗時(shí)的完全分析,此類方法最初用于結(jié)構(gòu)設(shè)計(jì),隨后被眾多學(xué)者用來提升局部材料非線性分析的計(jì)算效率。重分析方法主要包含直接法和近似法,直接重分析法[37]將局部結(jié)構(gòu)修改或局部非線性變形產(chǎn)生的剛度矩陣表示為初始彈性剛度矩陣的低秩修正形式,并基于Shearman-Morrison-Woodbury 公式進(jìn)行問題求解,避免了整體切線剛度矩陣耗時(shí)的實(shí)時(shí)更新和分解,組合近似法[38]是基于預(yù)處理共軛梯度法發(fā)展而來的一種高效的近似重分析方法,該方法通過對(duì)變化的剛度矩陣進(jìn)行降階處理來縮減問題的計(jì)算規(guī)模,可快速求解大范圍非線性結(jié)構(gòu)的響應(yīng)。擬力法[39-41]將構(gòu)件的變形分為彈性和塑性兩部分,可在控制方程中保持初始彈性剛度矩陣不變并通過僅更新小規(guī)模塑性剛度矩陣實(shí)現(xiàn)局部非線性問題的高效求解,受其概念啟發(fā),Li等[42-44]提出了隔離非線性法,該方法通過對(duì)非線性材料的應(yīng)變進(jìn)行分解,并與有限元理論結(jié)合,可將切線剛度矩陣轉(zhuǎn)化為初始彈性剛度矩陣的低秩修正形式,結(jié)合Woodbury 公式實(shí)現(xiàn)了一般局部非線性有限元分析問題的高效求解,為克服該方法隨非線性區(qū)域規(guī)模擴(kuò)大時(shí)效率降低的不足,Li等[45-46]通過引入組合近似法的思想對(duì)該方法使用的Woodbury 求解公式進(jìn)行了深度優(yōu)化。上述學(xué)者針對(duì)結(jié)構(gòu)局部非線性特征發(fā)展了多種高效數(shù)值分析方法,極大推動(dòng)了結(jié)構(gòu)分析領(lǐng)域的發(fā)展。

    地震作用下砌體結(jié)構(gòu)的非線性變形狀態(tài)往往集中于樓梯間、縱橫墻交接處等局部區(qū)域,具有典型的局部化特征,因此,利用其局部非線性特征發(fā)展適用的高效分析方法對(duì)提升計(jì)算效率具有積極意義。本文基于整體式空間離散宏單元模型,并充分利用結(jié)構(gòu)局部非線性特征,提出了一種砌體結(jié)構(gòu)高效非線性分析方法,首先建立剪切單元和無厚界面單元,并通過耦合兩者之間的變形模擬砌體失效模式,隨后,結(jié)合隔離非線性思想將剪切單元等效斜向彈簧的軸向變形和無厚界面單元上下表面的相對(duì)變形分解為線彈性和非線性兩部分,并通過附加的塑性自由度描述其中的非線性部分,從而在控制方程中實(shí)現(xiàn)了線性與非線性剛度矩陣的分離,以該控制方程為基礎(chǔ)將結(jié)構(gòu)切線剛度矩陣表示為初始彈性剛度的低秩修正形式,最后,通過引入Woodbury 公式求解控制方程,使非線性分析過程的主要計(jì)算量集中于對(duì)一個(gè)小規(guī)模的非線性矩陣進(jìn)行更新與分解,有效提升了非線性分析效率。數(shù)值模擬與試驗(yàn)結(jié)果的對(duì)比驗(yàn)證了本文方法的準(zhǔn)確性,基于時(shí)間復(fù)雜度理論的計(jì)算效率分析結(jié)果表明該文方法在非線性分析效率方面較傳統(tǒng)變剛度法具有明顯優(yōu)勢(shì)。

    1 離散宏單元模型

    Caliò等[22]基于整體式建模策略提出了砌體結(jié)構(gòu)的空間離散宏單元模型,如圖1 所示,其根據(jù)砌體墻窗洞構(gòu)造將墻體離散成一系列剪切單元和無厚界面單元,剪切單元為由四個(gè)剛性板和兩個(gè)對(duì)角斜向非線性彈簧構(gòu)成的鉸接四邊形結(jié)構(gòu),無厚界面單元由一系列不同功能的非線性彈簧構(gòu)成。不同剪切單元之間或剪切單元與支座之間通過無厚界面單元相互連接,砌體墻的平面內(nèi)、外失效模式通過耦合剪切單元和無厚界面單元的變形進(jìn)行模擬。

    圖1 空間離散宏單元模型Fig. 1 Spatial discrete macro-element model

    砌體墻斜截面剪切破壞是受損砌體結(jié)構(gòu)中最普遍的失效模式,當(dāng)墻體在豎向荷載和水平地震作用下產(chǎn)生的主拉應(yīng)力超過墻體的抗拉強(qiáng)度時(shí),平行于地震作用方向的墻體產(chǎn)生斜向裂縫,如圖2(a)所示。該破壞模式通過剪切單元進(jìn)行模擬,剪切單元將墻肢等效為發(fā)生純剪變形的均質(zhì)板,非線性行為通過兩個(gè)對(duì)角斜向非線性彈簧進(jìn)行描述,如圖3(a)所示。彈簧的初始彈性剛度可基于墻肢與均質(zhì)板線彈性變形范圍內(nèi)的能量等效[22]確定,屈服強(qiáng)度可根據(jù)Turnsek-Cacovic 準(zhǔn)則[22]確定,如圖4(a)所示為斜向彈簧的單軸滯回模型[22,24],模型中考慮了材料屈服軟化和卸載剛度退化等非線性行為,卸載剛度退化系數(shù)可取為β=0.8[24]。

    圖2 砌體墻平面內(nèi)外主要失效模式Fig. 2 Main in-plane and out-of-plane failure mechanisms of a masonry portion

    圖3 離散宏單元模型模擬砌體墻平面內(nèi)外主要失效模式(為簡(jiǎn)化說明,圖(b)、(c)、(d)僅畫出一根斜向彈簧)Fig. 3 Simulation of the main in-plane and out-of-plane failure mechanisms of a masonry portion by means of the macro-element (For simplicity, only one diagonal spring is shown in figures (b), (c) and (d))

    砌體墻正截面彎曲破壞模式表現(xiàn)為墻體受拉區(qū)的開裂和受壓區(qū)的壓碎,如圖2(b)所示。該破壞模式通過在墻肢截面處劃分的若干勻質(zhì)彈塑性纖維模擬,纖維的受力狀態(tài)通過界面單元中m行n列法向非線性彈簧進(jìn)行描述(圖3(b)),在多個(gè)剪切單元與界面單元組合而成的砌體墻模型中,這些界面彈簧僅影響與其相連的剪切單元的一半?yún)^(qū)域(如圖1 所示,墻體左上豎向界面單元僅影響單元左右兩側(cè)陰影區(qū)域),圖4(b)為這些彈簧的單軸滯回模型[22,24],模型中考慮了材料卸載剛度退化行為,卸載剛度退化系數(shù)取為 β=0.8[24]。

    當(dāng)砌體墻豎向荷載較小或砌體與砂漿界面摩擦系數(shù)較低時(shí),墻體易發(fā)生正截面剪切滑移破壞,表現(xiàn)為沿砌體的灰縫出現(xiàn)貫通截面的滑動(dòng)摩擦行為,如圖2(c)所示。為此,沿界面單元縱向設(shè)置一根剪切非線性彈簧描述墻體的正截面剪切滑移破壞,如圖3(c)所示,該剪切彈簧的本構(gòu)關(guān)系采用理想剛塑性模型[23-24](本文數(shù)值計(jì)算時(shí)將圖4(c)所示理想彈塑性模型的彈性剛度取為大數(shù)),屈服強(qiáng)度根據(jù)Mohr-Coulomb 準(zhǔn)則[22]確定,模型中未考慮材料的卸載剛度退化行為,即 β=0。

    圖4 砌體材料的單軸本構(gòu)關(guān)系Fig. 4 Uniaxial constitutive laws of masonry material

    當(dāng)砌體結(jié)構(gòu)縱、橫墻間連接薄弱時(shí),在沿橫墻方向的水平地震作用下結(jié)構(gòu)內(nèi)外墻交接處易產(chǎn)生豎向裂縫,此時(shí)縱墻向外甩出發(fā)生平面外彎曲、剪切和扭轉(zhuǎn)變形,此時(shí)墻體變形如圖2(d)所示。模型中沿界面單元橫向設(shè)置兩根剪切非線性彈簧(圖3(d))模擬墻體的平面外剪切、扭轉(zhuǎn)變形,每根彈簧的單軸滯回模型如圖4(c)所示,屈服強(qiáng)度根據(jù)Coulomb 失效準(zhǔn)則確定[23-24]。

    上述剪切單元通過無厚界面單元連接建立砌體墻片模型,如圖1 所示,進(jìn)一步將墻片模型通過共結(jié)點(diǎn)的方式連接,并根據(jù)結(jié)構(gòu)樓蓋形式及支承條件設(shè)定邊界約束,即可建立整體結(jié)構(gòu)分析模型,進(jìn)而實(shí)現(xiàn)整體結(jié)構(gòu)破壞形式和損傷演化過程的模擬。

    2 非線性隔離宏單元模型

    基于隔離非線性法基本思想,將材料應(yīng)變或單元變形分解為線彈性和非線性兩部分,圖5 以單軸材料的非線性本構(gòu)關(guān)系為例給出了應(yīng)變分解的基本原理,假設(shè)在某個(gè)增量步內(nèi),材料的非線性狀態(tài)由B點(diǎn)到C點(diǎn),根據(jù)初始彈性模量Ee將材料應(yīng)變?cè)隽糠纸鉃閮刹糠郑?/p>

    圖5 非線性材料應(yīng)變分解Fig. 5 Decomposition of nonlinear material strain

    式中: Δε′為材料的線彈性應(yīng)變?cè)隽浚?Δε′′為材料的非線性應(yīng)變?cè)隽?。與應(yīng)變?cè)隽?Δε相對(duì)應(yīng)的應(yīng)力增量 Δσ可表示為:

    式(2)表明,無論何種本構(gòu)模型,應(yīng)變分解后任意時(shí)刻的材料應(yīng)力均可表示為恒常的初始彈性模量與線彈性應(yīng)變的乘積。

    2.1 剪切單元模型

    選擇剪切單元的形心作為基本結(jié)點(diǎn),每個(gè)基本結(jié)點(diǎn)k有7 個(gè)自由度,其中6 個(gè)為剛體運(yùn)動(dòng)自由度,1 個(gè)為剪切變形自由度,如圖6 所示,其中oxyz為單元局部坐標(biāo)系。在迭代過程的任一線性化增量步中,基本結(jié)點(diǎn)k位移增量向量可表示為:

    如圖6 所示,單元每個(gè)剛性板的中心設(shè)置1 個(gè)外部結(jié)點(diǎn),共4 個(gè)外部結(jié)點(diǎn)。其中:①號(hào)、②號(hào)結(jié)點(diǎn)沿局部x軸布置;③號(hào)、④號(hào)結(jié)點(diǎn)沿局部z軸布置。每個(gè)外部結(jié)點(diǎn)r有6 個(gè)自由度,即有6 個(gè)廣義位移和6 個(gè)廣義力,其中僅在xz面內(nèi)的剪力使單元發(fā)生剪切變形。外部結(jié)點(diǎn)r位移增量向量為:

    圖6 剪切單元結(jié)點(diǎn)自由度Fig. 6 Shear element nodal degrees of freedom

    圖7 給出了單元的剪切變形模式(考慮小變形假設(shè)),將③、④號(hào)外部結(jié)點(diǎn)的連線定義為單元剪切變形模式的主線,該主線在xz面內(nèi)的轉(zhuǎn)動(dòng)變形Δγy(基本結(jié)點(diǎn)剪切變形自由度)將引起主線端部外部結(jié)點(diǎn)③、④的平動(dòng)位移dz·Δγy和其他外部結(jié)點(diǎn)①、②的轉(zhuǎn)動(dòng)位移 Δγy。

    圖7 剪切單元剪切變形模式Fig. 7 Shear deformation mode of shear element

    在單元剪切變形模式下,等效斜向彈簧的變形如圖8 所示,根據(jù)圖中幾何關(guān)系,可建立剪應(yīng)變 Δ γy和彈簧變形 Δd之間的關(guān)系:

    圖8 等效斜向彈簧變形圖Fig. 8 Deformation of equivalent diagonal spring

    式中: Δd為剪切單元斜向彈簧變形增量; Φdiag為變形矩陣。

    基于圖5 所示的隔離非線性方法應(yīng)變分解思想,可將等效斜向彈簧變形 Δd分解為線彈性和非線性兩部分:

    式中,線彈性部分 Δd′等于單元等效斜向彈簧初始彈性剛度的倒數(shù)乘以等效斜向彈簧內(nèi)力增量ΔN,即:

    非線性部分 Δd′′等于彈簧變形增量 Δd與線彈性變形增量 Δd′的差值。斜向彈簧內(nèi)力增量可表示為其初始彈性剛度與線彈性變形增量的乘積,即:

    式中,Ddiag,e為剪切單元的等效斜向彈簧的初始彈性剛度。

    2.2 無厚界面單元模型

    基于多點(diǎn)線性約束方法(multi-point constraints,MPC),剪切單元每一外部結(jié)點(diǎn)r的位移可以用基本結(jié)點(diǎn)k剛體運(yùn)動(dòng)引起的位移和剪切變形(圖7 所示剪切變形模式)引起的轉(zhuǎn)動(dòng)位移或平動(dòng)位移表示。在局部坐標(biāo)系中,求解外部結(jié)點(diǎn)r的位移時(shí),定義以基本結(jié)點(diǎn)k為起點(diǎn)、外部結(jié)點(diǎn)r為終點(diǎn)的向量=(dx,dy,dz)(如圖9 給出了k=2,r=④2時(shí)的向量,④2表示基本結(jié)點(diǎn)2 的④號(hào)外部結(jié)點(diǎn)),則單元外部結(jié)點(diǎn)r的位移函數(shù)為:

    圖9 無厚纖維界面單元Fig. 9 Zero-thickness fiber interface element

    式中,等號(hào)右端第一項(xiàng)為基本結(jié)點(diǎn)k剛體運(yùn)動(dòng)引起的位移項(xiàng),第二項(xiàng)和第三項(xiàng)分別為剪切變形引起的轉(zhuǎn)動(dòng)位移項(xiàng)和平動(dòng)位移項(xiàng),其中 Δγrot和Δγdisp分別為剪切變形模式下控制外部結(jié)點(diǎn)r轉(zhuǎn)動(dòng)位移和平動(dòng)位移的轉(zhuǎn)動(dòng)自由度,erot=(,,)和edisp=(,,) 分 別為 Δγrot和 Δγdisp的單位方向向量,在圖7 所示的單元剪切變形模式下,控制外部結(jié)點(diǎn)①、②位移的轉(zhuǎn)動(dòng)自由度 Δγrot=Δγy,Δγdisp=0,erot=(0,1,0),控制外部結(jié)點(diǎn)③、④位移的轉(zhuǎn)動(dòng)自由度 Δγrot=0 , Δγdisp=Δγy,edisp=(0,1,0)。

    圖9 表示出無厚纖維界面單元連接基本結(jié)點(diǎn)1 的③號(hào)外部結(jié)點(diǎn)(即圖中③1結(jié)點(diǎn))和基本結(jié)點(diǎn)2 的④號(hào)外部結(jié)點(diǎn)(即圖中④2結(jié)點(diǎn))。從圖9 可知相鄰③1、④2號(hào)外部結(jié)點(diǎn)的相對(duì)位移產(chǎn)生了界面相對(duì)變形,基于此可將界面單元相對(duì)位移函數(shù)表示為:

    式中: Δuo={ΔuΔvΔwΔθxΔθyΔθz}T為界面單元相對(duì)變形增量;I6為6 階單位矩陣; ΦI為變形矩陣。

    將無厚界面單元的相對(duì)變形 Δuo分解為線彈性和非線性兩部分:

    為方便后文對(duì)圖3(b)、圖3(c)、圖3(d)描述的單元變形分量展開討論,對(duì)界面單元相對(duì)變形增量進(jìn)行初等變換后寫為:

    式中: ΔuCB={ΔwΔθyΔθx}T為單元法向彈簧模擬壓彎受力行為的變量; ΔuinS={Δu}T為局部x向切向彈簧控制面內(nèi)剪切滑移變形的變量;ΔuoutST={ΔvΔθz}T為局部y向切向彈簧模擬面外剪扭受力行為的變量?;诖丝蓪⒎纸夂蟮慕缑鎲卧冃螌憺槿缦滦问剑?/p>

    式(16)中,DCB需通過對(duì)單元中各法向彈簧的彈性剛度集成得到,即:

    式中:ECB和A分別為彈簧的彈性模量和影響面積;xi和yi為第i根彈簧的局部坐標(biāo),彈簧沿局部y向布置m行,沿局部x向布置n列,m×n為彈簧的數(shù)量,DinS根據(jù)單元中沿局部x向切向彈簧的彈性剛度得到:

    式中:EinS代表彈簧的彈性模量;沿單元局部x向設(shè)置一根切向彈簧,影響面積AinS為單元的面積;DoutST需對(duì)單元中沿局部y向兩根彈簧的彈性剛度集成得到:

    式中:EoutST和A分別為彈簧的彈性模量和影響面積;xj表示第j根彈簧的局部坐標(biāo)。

    1) 根據(jù)隔離非線性法中的變形分解原理,單元中第i根法向彈簧的相對(duì)位移增量可分解為:

    實(shí)際分析中采用Newton-Raphson 迭代法求解,理想彈塑性模型會(huì)導(dǎo)致剛度矩陣奇異而出現(xiàn)求解失敗,為了改善數(shù)值計(jì)算的收斂性,本文將本構(gòu)關(guān)系中屈服后剛度與彈性剛度的比值取為0.001。當(dāng)彈簧本構(gòu)關(guān)系采用理想彈塑性模型時(shí),根據(jù)圖5 所示非線性區(qū)段內(nèi)的材料應(yīng)變分解思想,相應(yīng)的彈簧應(yīng)力增量可表示為:

    對(duì)單元中各法向彈簧的應(yīng)力增量進(jìn)行集成,可得與 ΔmCB對(duì)應(yīng)的單元法向內(nèi)力和彎矩表達(dá)式:

    式中: ΔNz為單元法向內(nèi)力; ΔMy和 ΔMx為單元彎矩。將式(22)代入式(16)中可建立與各法向彈簧的相對(duì)位移的關(guān)系:

    對(duì)比式(23)和式(12)可知,單元相對(duì)變形分量 ΔuCB和非線性相對(duì)變形分量可分別表示為:

    式(25)表明,當(dāng)單元中所有法向彈簧的非線性相對(duì)位移均為零值時(shí),單元的相對(duì)非線性變形分量才為零向量(此時(shí)相應(yīng)塑性自由度處于未激活狀態(tài)),若部分彈簧的非線性相對(duì)位移為非零值,則無論這些彈簧數(shù)量的多寡,單元非線性相對(duì)變形分量均為非零向量(此時(shí)相應(yīng)的塑性自由度被激活)。

    2) 沿單元局部x向僅設(shè)置一根切向彈簧,與ΔminS對(duì)應(yīng)的單元剪力為:

    式中, ΔVx為單元剪力。根據(jù)單元線彈性變形表達(dá)式(12)可知當(dāng)單元中該切向彈簧處于彈性狀態(tài)時(shí),單元非線性相對(duì)變形分量為零向量(相應(yīng)塑性自由度未激活),若為非零值,則為非零向量(相應(yīng)塑性自由度被激活)。

    3) 同樣地,對(duì)單元中沿局部y向切向彈簧的剪切相對(duì)位移進(jìn)行分解,將兩根彈簧的應(yīng)力增量進(jìn)行集成,得與 ΔmoutST對(duì)應(yīng)的單元剪力和扭矩為:[ΔVyΔMz]T=

    式中: ΔVy為單元剪力; ΔMz為單元扭矩。進(jìn)一步根據(jù)單元線彈性變形表達(dá)式(12)可知單元相對(duì)變形分量 ΔuoutST和非線性相對(duì)變形分量分別為:

    式(29)表明,當(dāng)單元中兩根切向彈簧的非線性剪切相對(duì)位移均為零值時(shí),單元的非線性相對(duì)變形分量才為零向量(相應(yīng)塑性自由度未激活),若其中一根彈簧的非線性剪切相對(duì)位移為非零值,單元非線性相對(duì)變形分量為非零向量(相應(yīng)塑性自由度被激活)。

    2.3 單元控制方程

    上述剪切單元和界面單元的變形分解后,單元內(nèi)力增量與其變形增量和非線性變形增量相關(guān),式(5)和式(9)、式(10)分別給出了剪切單元和界面單元變形增量的計(jì)算模型,單元非線性變形通過 Δd′′和進(jìn)行描述,其中每個(gè)元素代表一個(gè)塑性自由度。

    圖9 建立的兩基本結(jié)點(diǎn)單元可作為模型的基礎(chǔ)單元,包含兩個(gè)剪切子單元及一個(gè)界面子單元,建立其虛功方程如下:

    式中: δu為單元的虛變形向量;ΔT=[ΔNT]T為單元內(nèi)力; δa為單元虛結(jié)點(diǎn)位移向量; Δf為結(jié)點(diǎn)荷載增量。將單元位移函數(shù)式(5)、式(9) 、式(10)和內(nèi)力增量表達(dá)式(8)、式(16)代入虛功方程式(30),整理可得:

    此外,考慮單元中各彈簧的非線性本構(gòu)關(guān)系可建立補(bǔ)充方程,對(duì)于進(jìn)入非線性狀態(tài)的單元內(nèi)力使用下式求解:

    將式(32)代入式(31),可得基礎(chǔ)單元的隔離非線性控制方程:

    式中:ke為單元的初始彈性剛度矩陣;k′為描述單元結(jié)點(diǎn)力與非線性變形之間聯(lián)系的系數(shù)矩陣;為描述單元非線性行為的系數(shù)矩陣。具體通過下式計(jì)算:

    式中,B矩陣為代表基本結(jié)點(diǎn)自由度與兩剪切子單元斜向彈簧變形 Δd和界面子單元相對(duì)變形 Δuo之間關(guān)系的變形矩陣。對(duì)B矩陣,可先將式(9)代入式(10)得到界面單元的變形矩陣,再將式(5)擴(kuò)充為與其統(tǒng)一階數(shù)的矩陣形式,最后進(jìn)行合并得到。De=diag(DI,e,Ddiag,e),且Ddiag,e=diag(,),和分別為圖9 中相鄰1、2 號(hào)剪切單元的等效斜向彈簧的初始彈性剛度,Dt為相應(yīng)的切向剛度矩陣。

    3 結(jié)構(gòu)控制方程求解

    將基礎(chǔ)單元控制方程進(jìn)行集成,可得整體結(jié)構(gòu)的隔離非線性控制方程:

    式中: ΔX和 ΔF分別為整體結(jié)構(gòu)的結(jié)點(diǎn)位移增量向量和結(jié)點(diǎn)荷載增量向量;為整體結(jié)構(gòu)的非線性變形增量向量,其中僅包含了結(jié)構(gòu)中呈激活狀態(tài)的m個(gè)塑性自由度;Ke為整體結(jié)構(gòu)的初始彈性剛度矩陣,其規(guī)模為n×n階,階數(shù)等于結(jié)構(gòu)位移自由度數(shù)n;K′為描述整體結(jié)構(gòu)結(jié)點(diǎn)力與非線性變形之間聯(lián)系的系數(shù)矩陣,其規(guī)模為n×m階;為描述整體結(jié)構(gòu)材料非線性行為的系數(shù)矩陣,其規(guī)模為m×m階。當(dāng)結(jié)構(gòu)發(fā)生局部材料非線性變形時(shí),激活的塑性自由度數(shù)目m將遠(yuǎn)小于結(jié)構(gòu)自由度數(shù)目n(即m<<n)。

    式(38)表明:當(dāng)m遠(yuǎn)小于n時(shí),任意時(shí)刻結(jié)構(gòu)的切線剛度矩陣可表示為初始彈性剛度矩陣的低秩攝動(dòng)形式。

    采用完全Newton-Raphson 迭代法進(jìn)行結(jié)構(gòu)非線性分析,引入Woodbury 公式(式(39))在每個(gè)迭代步中對(duì)相應(yīng)式(38)進(jìn)行高效求解:

    式中,KSchur為Schur 補(bǔ)矩陣:

    從式(39)可以看出,結(jié)構(gòu)整體剛度項(xiàng)Ke在整體分析過程中始終保持初始彈性狀態(tài),僅需在分析前進(jìn)行一次分解即可,當(dāng)材料非線性集中在結(jié)構(gòu)的局部區(qū)域時(shí),控制方程中塑性自由度數(shù)目遠(yuǎn)小于結(jié)構(gòu)自由度數(shù)目(即m<<n),此時(shí)Woodbury公式在求解結(jié)構(gòu)切向位移反應(yīng)時(shí)的主要計(jì)算開銷將僅集中于小規(guī)模的Schur 補(bǔ)矩陣(m×m階)的更新分解[43],從而可避免傳統(tǒng)變剛度法中對(duì)整體切線剛度矩陣(n×n階)的反復(fù)合成與分解運(yùn)算,實(shí)現(xiàn)計(jì)算效率的大幅提升。

    4 數(shù)值算例

    4.1 試驗(yàn)驗(yàn)證

    為了驗(yàn)證本文提出的砌體結(jié)構(gòu)非線性分析方法的正確性,以Anthoine 等[47]的擬靜力試驗(yàn)數(shù)據(jù)為依據(jù),對(duì)2 片不同高度的磚砌體墻進(jìn)行數(shù)值模擬。墻體的材料屬性和幾何屬性分別如表1 和表2所示,2 片墻體在豎向壓應(yīng)力下的低周反復(fù)荷載試驗(yàn)結(jié)果如圖10 所示,其中圖10(a)為Wall-1 正截面彎曲破壞的試驗(yàn)滯回曲線,可以看出,曲線捏縮較為顯著,墻體耗能能力較弱但強(qiáng)度在變形較大時(shí)沒有明顯退化;圖10(b)為Wall-2 斜截面剪切破壞的試驗(yàn)滯回曲線,曲線的包圍面積較大,墻體耗能能力相對(duì)較強(qiáng),但存在嚴(yán)重的剛度和強(qiáng)度退化,試驗(yàn)詳細(xì)信息見參考文獻(xiàn)[47]。

    表1 墻體材料屬性Table 1 Material properties of masonry walls

    表2 墻體幾何屬性Table 2 Geometry properties of masonry walls

    圖10 砌體墻片試驗(yàn)滯回曲線Fig. 10 Experimental behaviour of simple piers

    采用本文方法模擬2 片砌體墻在豎向荷載和水平荷載共同作用下的面內(nèi)滯回行為,每個(gè)模型使用1 個(gè)剪切單元和1 個(gè)設(shè)置于底端的無厚界面單元。模型中各類非線性彈簧的力學(xué)參數(shù)取值根據(jù)表1 中砌體墻力學(xué)參數(shù)確定。兩個(gè)模型數(shù)值模擬結(jié)果如圖11 所示,分別為Wall-1 和Wall-2 基底剪力與頂點(diǎn)位移之間的滯回曲線,此外,圖12給出了Caliò和Pantò等求解的Wall-1 和Wall-2 兩片砌體墻宏單元模型在低周反復(fù)荷載作用下的滯回曲線(CP 方法求解結(jié)果)[22]。對(duì)比砌體墻數(shù)值模擬與試驗(yàn)結(jié)果:三者剛度和強(qiáng)度整體上吻合良好,其中CP 方法求解Wall-1 滯回曲線捏縮更為顯著,這是由模型中非線性彈簧本構(gòu)關(guān)系中滯回規(guī)則的差異所導(dǎo)致的。為更充分驗(yàn)證本文方法的正確性,圖13 提取了三組滯回曲線的骨架曲線,從圖中知,對(duì)比Wall-1 頂端最大位移處荷載P1,本文方法求解結(jié)果和CP 方法求解結(jié)果與試驗(yàn)結(jié)果的相對(duì)誤差分別為12.38%和16.14%;對(duì)比Wall-2 頂端最大位移處荷載P2,本文方法求解結(jié)果和CP 方法求解結(jié)果與試驗(yàn)結(jié)果的相對(duì)誤差分別為9.51%和25.73%??梢姡罕疚姆椒ㄇ蠼庹`差小于CP 方法,本文所提方法可以有效地模擬砌體墻的失效模式,驗(yàn)證了該分析方法的正確性。

    圖11 砌體墻片本文方法計(jì)算滯回曲線Fig. 11 Numerical simulation results of piers(proposed method)

    圖12 砌體墻片CP 方法計(jì)算滯回曲線Fig. 12 Numerical simulation results of piers (CP method)

    圖13 砌體墻片骨架曲線Fig. 13 Skeleton curves of simple piers

    4.2 動(dòng)力非線性分析

    為進(jìn)一步對(duì)本文方法進(jìn)行驗(yàn)證,對(duì)某4 層約束砌體結(jié)構(gòu)進(jìn)行地震反應(yīng)分析,該結(jié)構(gòu)平面簡(jiǎn)圖如圖14 所示。結(jié)構(gòu)層高為3 m,總高度為12 m,縱墻上窗洞的尺寸為1.5 m ×1.8 m,門洞的尺寸為1.0 m ×1.8 m。在結(jié)構(gòu)縱橫墻交接處設(shè)置構(gòu)造柱(如圖14 所示),結(jié)構(gòu)所有縱橫墻在樓層標(biāo)高處設(shè)置鋼筋混凝土圈梁,構(gòu)造柱和圈梁采用C20 混凝土,截面尺寸為250 mm ×250 mm,縱筋采用4 根HRB335直徑16 mm 的鋼筋;墻體厚度為250 mm,材料為燒結(jié)普通磚,強(qiáng)度等級(jí)為MU10,混合砂漿強(qiáng)度等級(jí)為M10。結(jié)構(gòu)樓面恒載為5.49 kN/m2,活載為2.0 kN/m2,屋面恒載為6.371 kN/m2。設(shè)防烈度為8 度(0.2g),Ⅱ類場(chǎng)地,地震動(dòng)選取汶川地震中臥龍波(幅值為400 cm/s2),加速度時(shí)程如圖15 所示。圖16 為本文方法建立的該結(jié)構(gòu)數(shù)值分析模型,其中,模型假設(shè)在結(jié)構(gòu)首層底端設(shè)置固定約束,圈梁和構(gòu)造柱采用隔離非線性纖維梁?jiǎn)卧M(jìn)行模擬,共劃分1894 個(gè)隔離非線性纖維梁?jiǎn)卧猍48],砌體墻劃分為4198 個(gè)剪切單元和9770 個(gè)無厚界面單元,該空間模型共有39 934 個(gè)位移自由度。圈梁與構(gòu)造柱和砌體墻的相互連接通過在接觸面設(shè)置界面單元的方式實(shí)現(xiàn),另外,在樓板的每一區(qū)格設(shè)置兩根斜向剛性桿傳遞水平荷載和協(xié)調(diào)樓面結(jié)點(diǎn)位移。為簡(jiǎn)化模擬,本例中僅考慮纖維界面單元的軸向和彎曲非線性變形,其剪切和扭轉(zhuǎn)變形假設(shè)為彈性。

    圖14 結(jié)構(gòu)平面簡(jiǎn)圖Fig. 14 Plane of the masonry structure

    圖15 臥龍地震波加速度記錄Fig. 15 Acceleration record of Wolong earthquake

    圖16 約束砌體結(jié)構(gòu)數(shù)值分析模型Fig. 16 Discrete macro-element model of the masonry structure

    結(jié)構(gòu)動(dòng)力響應(yīng)分析之前對(duì)結(jié)構(gòu)進(jìn)行模態(tài)分析,發(fā)現(xiàn)結(jié)構(gòu)第一陣型是一階Y向平動(dòng),一階陣型的周期為T1=0.14 s,第二陣型是一階X向平動(dòng),二階陣型的周期是T2=0.1 s。分析結(jié)構(gòu)在Y向水平地震作用下的動(dòng)力響應(yīng),結(jié)構(gòu)外縱墻與外橫墻交接處各樓層(圖16 示)位移時(shí)程響應(yīng)如圖17所示,結(jié)構(gòu)的基底剪力時(shí)程響應(yīng)如圖18 所示,根據(jù)各樓層位移時(shí)程曲線發(fā)現(xiàn)該結(jié)構(gòu)底層層間位移角最大,峰值達(dá)到0.6%,依據(jù)文獻(xiàn)[49]中約束砌體結(jié)構(gòu)0.4%的倒塌限值判定該砌體結(jié)構(gòu)已達(dá)到倒塌的破壞狀態(tài)。本砌體結(jié)構(gòu)主要震害集中在底層,外橫墻墻肢的損傷分布如圖19 所示,墻體下部1 層、2 層的剪切單元均發(fā)生剪切破壞,下部損傷明顯比上部嚴(yán)重。

    圖17 各樓層位移時(shí)程曲線Fig. 17 Time-history curves of all floors

    圖18 結(jié)構(gòu)基底剪力時(shí)程曲線Fig. 18 Base shear force curve of the masonry structure

    圖19 外橫墻損傷分布圖Fig. 19 Damage distribution of outer transversal wall

    圖20 給出了塑性自由度數(shù)隨增量步的變化曲線,在3.265 s(增量步為653 步)之前結(jié)構(gòu)處于彈性狀態(tài),塑性自由度數(shù)為0。非線性分析過程中產(chǎn)生的最大塑性自由度數(shù)為6545,占總自由度數(shù)的16.4%,而每個(gè)增量步的平均塑性自由度為1997,僅占總自由度數(shù)的5%,說明結(jié)構(gòu)具有顯著的局部非線性特征,采用式(39)進(jìn)行求解時(shí),由于需進(jìn)行實(shí)時(shí)更新和分解的Schur 補(bǔ)矩陣階數(shù)等于塑性自由度數(shù)目,局部非線性引發(fā)的較少的塑性自由度可保證本文方法實(shí)現(xiàn)大幅度的效率提升。

    圖20 塑性自由度時(shí)程曲線Fig. 20 Time-history curves of plastic degrees of freedom

    為了進(jìn)一步論證本文方法的高效性,將傳統(tǒng)非線性分析方法(矩陣分解采用LDLT法)和本文方法的時(shí)間復(fù)雜度[50]進(jìn)行統(tǒng)計(jì)對(duì)比,如圖21 所示,傳統(tǒng)非線性分析方法和本文方法的平均時(shí)間復(fù)雜度分別為1.54 ×1010和1.40 ×109,傳統(tǒng)方法的平均時(shí)間復(fù)雜度是本文方法的11 倍,證明了本文方法的高效性。

    圖21 時(shí)間復(fù)雜度時(shí)程曲線Fig. 21 Time-history curves of time complexity

    5 結(jié)論

    本文通過耦合剪切單元和無厚界面單元來模擬砌體結(jié)構(gòu)失效行為,以隔離非線性法理論為基礎(chǔ),通過將剪切單元中等效斜向彈簧軸向變形和無厚界面單元上下表面的相對(duì)變形分解為線彈性和非線性兩部分,推導(dǎo)出砌體結(jié)構(gòu)離散宏單元的隔離非線性控制方程,采用Woodbury 公式求解控制方程,得到如下結(jié)論:

    (1) 本文利用砌體結(jié)構(gòu)材料非線性的局部化特征建立了結(jié)構(gòu)非線性隔離宏單元模型,該模型基于多點(diǎn)線性約束方法使每個(gè)基本結(jié)點(diǎn)的自由度大幅縮減,相比于已有宏觀模型和精細(xì)模型,其在大型砌體結(jié)構(gòu)非線性分析時(shí)既能保證較高的計(jì)算精度也能使模型規(guī)模維持在相對(duì)較低的程度。同時(shí),模型中激活的塑性自由度數(shù)目等于結(jié)構(gòu)控制方程中小規(guī)模非線性矩陣的維度,其可直觀地衡量砌體結(jié)構(gòu)中產(chǎn)生增量非線性變形的區(qū)域規(guī)模。

    (2) 通過分離結(jié)構(gòu)彈塑性狀態(tài)使砌體結(jié)構(gòu)非線性分析時(shí)僅對(duì)小規(guī)模非線性矩陣進(jìn)行更新和分解,可以避免傳統(tǒng)變剛度法中對(duì)整體切向剛度矩陣的更新與分解,該方法的時(shí)間復(fù)雜度明顯低于傳統(tǒng)變剛度法,可以實(shí)現(xiàn)砌體結(jié)構(gòu)非線性問題的高效求解。

    猜你喜歡
    結(jié)點(diǎn)砌體彈簧
    砌體墻上安裝摩擦型阻尼器施工技術(shù)探討
    豎向開槽砌體墻燃?xì)獗▌?dòng)力響應(yīng)及加固
    析彈簧模型 悟三個(gè)性質(zhì)
    Ladyzhenskaya流體力學(xué)方程組的確定模與確定結(jié)點(diǎn)個(gè)數(shù)估計(jì)
    如何求串聯(lián)彈簧和并聯(lián)彈簧的勁度系數(shù)
    時(shí)間彈簧
    論建筑工程中砌體結(jié)構(gòu)現(xiàn)狀及前景
    采高對(duì)砌體梁關(guān)鍵層位置的影響分析
    基于Raspberry PI為結(jié)點(diǎn)的天氣云測(cè)量網(wǎng)絡(luò)實(shí)現(xiàn)
    基于DHT全分布式P2P-SIP網(wǎng)絡(luò)電話穩(wěn)定性研究與設(shè)計(jì)
    精品免费久久久久久久清纯| 成年女人毛片免费观看观看9| 国产精品久久电影中文字幕| 伦精品一区二区三区| 12—13女人毛片做爰片一| 国内揄拍国产精品人妻在线| 啪啪无遮挡十八禁网站| 亚洲人成网站在线播| 欧美激情国产日韩精品一区| 国产 一区精品| 长腿黑丝高跟| 国产淫片久久久久久久久| 亚洲不卡免费看| 欧美激情国产日韩精品一区| 亚洲国产精品sss在线观看| 欧美性猛交黑人性爽| 国内精品宾馆在线| 人妻丰满熟妇av一区二区三区| 久久香蕉精品热| 男女下面进入的视频免费午夜| 露出奶头的视频| 成人性生交大片免费视频hd| 真人做人爱边吃奶动态| 全区人妻精品视频| 在线免费观看的www视频| videossex国产| 免费黄网站久久成人精品| 国产免费男女视频| 国产一区二区三区av在线 | 国产免费av片在线观看野外av| 成人午夜高清在线视频| 亚洲久久久久久中文字幕| 天天一区二区日本电影三级| 在线观看一区二区三区| 嫩草影院入口| 日本一本二区三区精品| 国产午夜福利久久久久久| 亚洲av.av天堂| 欧美黑人巨大hd| 精品日产1卡2卡| 一个人看视频在线观看www免费| 亚洲真实伦在线观看| 国产一区二区三区在线臀色熟女| 成人鲁丝片一二三区免费| 欧美日本视频| 18禁黄网站禁片免费观看直播| 亚洲性夜色夜夜综合| 非洲黑人性xxxx精品又粗又长| 波多野结衣巨乳人妻| ponron亚洲| 麻豆国产97在线/欧美| 日韩亚洲欧美综合| 中文字幕高清在线视频| 精品人妻1区二区| 我要看日韩黄色一级片| 久久久国产成人精品二区| 欧美+日韩+精品| 2021天堂中文幕一二区在线观| 99九九线精品视频在线观看视频| 麻豆av噜噜一区二区三区| 在线播放无遮挡| 神马国产精品三级电影在线观看| 观看美女的网站| 免费看日本二区| 日日摸夜夜添夜夜添av毛片 | 中文字幕免费在线视频6| 日本熟妇午夜| 色噜噜av男人的天堂激情| 亚洲国产精品成人综合色| 老司机福利观看| 日韩,欧美,国产一区二区三区 | 精品一区二区三区av网在线观看| 老熟妇仑乱视频hdxx| 欧美一区二区国产精品久久精品| 亚洲成a人片在线一区二区| 国产伦在线观看视频一区| 精品人妻视频免费看| 搡老妇女老女人老熟妇| 极品教师在线视频| 婷婷精品国产亚洲av| 国产伦精品一区二区三区视频9| 夜夜爽天天搞| 国产高清三级在线| 日韩精品青青久久久久久| 国产一区二区三区在线臀色熟女| 在线免费观看的www视频| 欧美xxxx黑人xx丫x性爽| 尤物成人国产欧美一区二区三区| 亚洲成人久久性| 女生性感内裤真人,穿戴方法视频| 99精品久久久久人妻精品| 国产色爽女视频免费观看| 日本黄大片高清| 国产国拍精品亚洲av在线观看| 国产真实乱freesex| 可以在线观看毛片的网站| 此物有八面人人有两片| а√天堂www在线а√下载| 校园春色视频在线观看| 动漫黄色视频在线观看| 88av欧美| 天天躁日日操中文字幕| 好男人在线观看高清免费视频| 99热这里只有是精品在线观看| 亚洲在线自拍视频| 国内精品一区二区在线观看| 中文字幕av成人在线电影| 国产蜜桃级精品一区二区三区| 听说在线观看完整版免费高清| 性色avwww在线观看| 免费看av在线观看网站| 99热精品在线国产| .国产精品久久| 久久人人精品亚洲av| 88av欧美| 日本 欧美在线| 欧美中文日本在线观看视频| 在线国产一区二区在线| 三级国产精品欧美在线观看| avwww免费| 精品久久久久久久末码| 淫妇啪啪啪对白视频| ponron亚洲| 亚洲性夜色夜夜综合| 91麻豆av在线| 1000部很黄的大片| 精品免费久久久久久久清纯| 国产高清三级在线| 亚洲自拍偷在线| 日韩av在线大香蕉| 免费观看精品视频网站| 精华霜和精华液先用哪个| 欧美xxxx黑人xx丫x性爽| 成人综合一区亚洲| 人人妻人人澡欧美一区二区| 天天躁日日操中文字幕| 小蜜桃在线观看免费完整版高清| 无遮挡黄片免费观看| 国产精品久久久久久av不卡| 窝窝影院91人妻| 久久人人精品亚洲av| 少妇猛男粗大的猛烈进出视频 | 精品日产1卡2卡| 久久香蕉精品热| 国国产精品蜜臀av免费| 精品99又大又爽又粗少妇毛片 | 成人av在线播放网站| 麻豆精品久久久久久蜜桃| 丰满的人妻完整版| 舔av片在线| 精品免费久久久久久久清纯| 在线免费观看不下载黄p国产 | 成年女人看的毛片在线观看| 欧美bdsm另类| 深爱激情五月婷婷| 国产成人aa在线观看| 真实男女啪啪啪动态图| 婷婷精品国产亚洲av| 美女免费视频网站| 噜噜噜噜噜久久久久久91| 午夜精品久久久久久毛片777| 日本在线视频免费播放| 黄色配什么色好看| 中文亚洲av片在线观看爽| 亚洲自拍偷在线| 国产高清激情床上av| 嫩草影院新地址| 中文在线观看免费www的网站| 老熟妇仑乱视频hdxx| 天美传媒精品一区二区| 99精品在免费线老司机午夜| 成人国产麻豆网| 日韩亚洲欧美综合| 欧美不卡视频在线免费观看| 久久精品综合一区二区三区| 中文字幕免费在线视频6| 成年女人看的毛片在线观看| 亚洲av熟女| 久久久久久久亚洲中文字幕| 搡老岳熟女国产| 久久国内精品自在自线图片| 欧美+日韩+精品| 日日撸夜夜添| 蜜桃亚洲精品一区二区三区| 三级男女做爰猛烈吃奶摸视频| 变态另类丝袜制服| 免费黄网站久久成人精品| 两个人的视频大全免费| 一级毛片久久久久久久久女| 联通29元200g的流量卡| 欧美+亚洲+日韩+国产| 哪里可以看免费的av片| 一区二区三区免费毛片| 成年免费大片在线观看| 麻豆一二三区av精品| 99久久久亚洲精品蜜臀av| 亚洲国产欧洲综合997久久,| 亚洲av免费高清在线观看| 十八禁国产超污无遮挡网站| 在线a可以看的网站| 丰满的人妻完整版| 亚洲va在线va天堂va国产| 午夜久久久久精精品| 女同久久另类99精品国产91| 日日撸夜夜添| av女优亚洲男人天堂| 99久久精品一区二区三区| 免费观看在线日韩| 久久精品影院6| 禁无遮挡网站| 高清在线国产一区| 国产黄片美女视频| 精品国产三级普通话版| 免费人成视频x8x8入口观看| 女的被弄到高潮叫床怎么办 | 美女黄网站色视频| 99精品久久久久人妻精品| 国产色爽女视频免费观看| av.在线天堂| 免费高清视频大片| 在现免费观看毛片| 99热精品在线国产| 久久久久久九九精品二区国产| 国产av不卡久久| 18禁黄网站禁片午夜丰满| 久久久久久久久中文| 欧美日韩综合久久久久久 | 成人特级av手机在线观看| 黄色视频,在线免费观看| 久久久久性生活片| 欧美高清成人免费视频www| 国产男靠女视频免费网站| 亚洲av一区综合| 91麻豆av在线| 亚洲国产高清在线一区二区三| 亚洲精品色激情综合| 国产成人av教育| .国产精品久久| 听说在线观看完整版免费高清| 国产综合懂色| 精品久久久久久,| 精品一区二区三区视频在线| 国产高清不卡午夜福利| 欧美激情在线99| 国产激情偷乱视频一区二区| 国产精品一区二区三区四区久久| 日韩欧美国产在线观看| avwww免费| 国内精品久久久久久久电影| 成年女人看的毛片在线观看| 色精品久久人妻99蜜桃| 欧美激情久久久久久爽电影| 少妇裸体淫交视频免费看高清| 国产色婷婷99| www.色视频.com| 啦啦啦啦在线视频资源| 亚洲精品亚洲一区二区| 麻豆国产av国片精品| 亚洲经典国产精华液单| 亚洲av第一区精品v没综合| 窝窝影院91人妻| 国产精品久久久久久精品电影| 美女 人体艺术 gogo| 一个人看视频在线观看www免费| 亚洲欧美日韩无卡精品| 国产伦人伦偷精品视频| 日本三级黄在线观看| 精品一区二区三区视频在线观看免费| 嫩草影院精品99| 欧美高清性xxxxhd video| av天堂中文字幕网| 美女黄网站色视频| bbb黄色大片| 看黄色毛片网站| 亚洲成人中文字幕在线播放| 一区二区三区免费毛片| 99热这里只有是精品50| 日韩人妻高清精品专区| 欧美性感艳星| 亚洲avbb在线观看| 中文亚洲av片在线观看爽| 亚洲精品粉嫩美女一区| 久久99热6这里只有精品| 免费无遮挡裸体视频| 亚洲18禁久久av| 中文字幕精品亚洲无线码一区| 亚洲最大成人中文| 国产精品免费一区二区三区在线| 毛片女人毛片| 国产精品不卡视频一区二区| 搡女人真爽免费视频火全软件 | 国产精品女同一区二区软件 | 中文字幕av成人在线电影| 深爱激情五月婷婷| 内地一区二区视频在线| 成年女人看的毛片在线观看| 久久久久久国产a免费观看| 亚洲精品456在线播放app | 国产精品久久视频播放| 国产精品一区二区免费欧美| 亚洲av中文av极速乱 | 亚洲av五月六月丁香网| 日韩欧美三级三区| 99热6这里只有精品| 网址你懂的国产日韩在线| 欧美精品国产亚洲| 免费看av在线观看网站| 在线观看免费视频日本深夜| 亚洲精品国产成人久久av| 日本撒尿小便嘘嘘汇集6| 亚洲国产高清在线一区二区三| 99久久无色码亚洲精品果冻| 亚洲色图av天堂| 欧美日韩乱码在线| 内射极品少妇av片p| 九九热线精品视视频播放| 人妻丰满熟妇av一区二区三区| 国产精品一及| 色尼玛亚洲综合影院| 欧美bdsm另类| 亚洲人成网站在线播| 精品人妻视频免费看| 国产美女午夜福利| 日韩欧美 国产精品| 久久国内精品自在自线图片| 国产高清不卡午夜福利| 狂野欧美激情性xxxx在线观看| 久久久国产成人精品二区| 国产精品不卡视频一区二区| 久久九九热精品免费| 国产欧美日韩精品一区二区| av天堂在线播放| 变态另类成人亚洲欧美熟女| 99热精品在线国产| 人人妻,人人澡人人爽秒播| 国产精品av视频在线免费观看| 午夜福利视频1000在线观看| 国产精品不卡视频一区二区| www.色视频.com| 精品无人区乱码1区二区| 黄色配什么色好看| 精品一区二区免费观看| 在现免费观看毛片| 成人高潮视频无遮挡免费网站| 搡女人真爽免费视频火全软件 | 精品久久久噜噜| 日韩中文字幕欧美一区二区| 久久99热这里只有精品18| 国产欧美日韩一区二区精品| 亚洲最大成人中文| 久久精品影院6| 午夜激情欧美在线| 国产午夜精品论理片| 精品国产三级普通话版| 最近在线观看免费完整版| 亚洲性久久影院| 热99在线观看视频| 久久九九热精品免费| 18+在线观看网站| 精品久久久久久久久亚洲 | 国产一级毛片七仙女欲春2| 琪琪午夜伦伦电影理论片6080| 国产日本99.免费观看| 亚洲欧美日韩东京热| 婷婷六月久久综合丁香| 精品无人区乱码1区二区| 白带黄色成豆腐渣| 亚洲18禁久久av| 欧美丝袜亚洲另类 | 色综合站精品国产| 男人狂女人下面高潮的视频| 日韩在线高清观看一区二区三区 | 深爱激情五月婷婷| 老师上课跳d突然被开到最大视频| av天堂在线播放| 国产精华一区二区三区| 亚洲,欧美,日韩| 99在线人妻在线中文字幕| 免费黄网站久久成人精品| 国产高清视频在线播放一区| 精品人妻一区二区三区麻豆 | 国产探花在线观看一区二区| 一区二区三区高清视频在线| 色精品久久人妻99蜜桃| 免费人成在线观看视频色| 亚洲精品日韩av片在线观看| 在现免费观看毛片| 日韩精品中文字幕看吧| 久久精品综合一区二区三区| 亚洲天堂国产精品一区在线| 亚洲国产精品sss在线观看| 简卡轻食公司| 99视频精品全部免费 在线| 我要看日韩黄色一级片| 狠狠狠狠99中文字幕| 欧美xxxx黑人xx丫x性爽| 啦啦啦观看免费观看视频高清| 男人舔奶头视频| 欧美成人性av电影在线观看| 麻豆av噜噜一区二区三区| 少妇的逼水好多| 一进一出好大好爽视频| 国产真实乱freesex| 搞女人的毛片| 一进一出好大好爽视频| 麻豆av噜噜一区二区三区| 日本黄大片高清| 欧美激情在线99| 亚洲中文日韩欧美视频| 97超视频在线观看视频| 国产精品伦人一区二区| 国产乱人视频| 男人的好看免费观看在线视频| 国内精品一区二区在线观看| 91久久精品国产一区二区三区| 国产av一区在线观看免费| 一级a爱片免费观看的视频| 日韩av在线大香蕉| 国产伦一二天堂av在线观看| 女生性感内裤真人,穿戴方法视频| 日本三级黄在线观看| av在线观看视频网站免费| 久久午夜亚洲精品久久| 网址你懂的国产日韩在线| 少妇裸体淫交视频免费看高清| 成人一区二区视频在线观看| 成人国产麻豆网| 蜜桃亚洲精品一区二区三区| 国产免费av片在线观看野外av| 国产 一区 欧美 日韩| 亚洲熟妇中文字幕五十中出| 国产精品不卡视频一区二区| 69av精品久久久久久| 成年版毛片免费区| 99久久精品国产国产毛片| 一进一出抽搐gif免费好疼| 韩国av在线不卡| 精华霜和精华液先用哪个| 国内久久婷婷六月综合欲色啪| 丰满人妻一区二区三区视频av| 亚洲欧美日韩高清在线视频| 乱人视频在线观看| 国产毛片a区久久久久| 性插视频无遮挡在线免费观看| 看十八女毛片水多多多| 在线观看美女被高潮喷水网站| 18禁黄网站禁片免费观看直播| 少妇猛男粗大的猛烈进出视频 | 免费观看在线日韩| 噜噜噜噜噜久久久久久91| 国产在线精品亚洲第一网站| 日日啪夜夜撸| 国产精品嫩草影院av在线观看 | 在线观看午夜福利视频| 国产高潮美女av| 两人在一起打扑克的视频| 少妇丰满av| 国产又黄又爽又无遮挡在线| 欧美日韩中文字幕国产精品一区二区三区| 久久精品国产亚洲av涩爱 | 一区二区三区四区激情视频 | 老熟妇乱子伦视频在线观看| 男人和女人高潮做爰伦理| 日日啪夜夜撸| 欧美高清性xxxxhd video| 美女高潮喷水抽搐中文字幕| netflix在线观看网站| 岛国在线免费视频观看| 免费观看在线日韩| 久久99热6这里只有精品| 日本 av在线| 99热这里只有是精品50| 欧美性感艳星| 国产精品亚洲一级av第二区| 国产成人av教育| 日韩中字成人| 国内精品久久久久精免费| 精品久久久久久久久av| 亚洲最大成人中文| 亚洲av免费高清在线观看| 婷婷精品国产亚洲av在线| 欧美极品一区二区三区四区| 欧美精品国产亚洲| 老司机深夜福利视频在线观看| 亚洲精品久久国产高清桃花| 91午夜精品亚洲一区二区三区 | 亚洲黑人精品在线| 久久人妻av系列| 免费av毛片视频| 国内久久婷婷六月综合欲色啪| 国产黄a三级三级三级人| 内射极品少妇av片p| 国产成人影院久久av| 国产三级中文精品| 蜜桃亚洲精品一区二区三区| 免费不卡的大黄色大毛片视频在线观看 | 国产91精品成人一区二区三区| 亚洲男人的天堂狠狠| 精品人妻视频免费看| 国产视频内射| 嫩草影视91久久| 亚洲成人久久爱视频| 两个人的视频大全免费| 麻豆国产av国片精品| 22中文网久久字幕| 麻豆成人午夜福利视频| 国产高清视频在线播放一区| 国产淫片久久久久久久久| 1000部很黄的大片| 美女cb高潮喷水在线观看| 麻豆成人av在线观看| 欧美成人a在线观看| 内射极品少妇av片p| 日韩欧美免费精品| av在线亚洲专区| 少妇裸体淫交视频免费看高清| 国产伦在线观看视频一区| 日本a在线网址| 欧美三级亚洲精品| 国产69精品久久久久777片| 国产私拍福利视频在线观看| netflix在线观看网站| 欧美xxxx黑人xx丫x性爽| 嫩草影视91久久| 久久精品国产亚洲av天美| 国产精品日韩av在线免费观看| 日韩精品中文字幕看吧| 精品人妻熟女av久视频| 久久这里只有精品中国| 人人妻人人看人人澡| 午夜福利在线观看吧| 全区人妻精品视频| 欧美色视频一区免费| 91久久精品电影网| 国产人妻一区二区三区在| 免费观看在线日韩| 国产精品嫩草影院av在线观看 | 亚洲熟妇熟女久久| 天堂网av新在线| АⅤ资源中文在线天堂| 亚洲三级黄色毛片| 精品午夜福利在线看| 黄色欧美视频在线观看| 99在线人妻在线中文字幕| 真人做人爱边吃奶动态| 国产一区二区三区在线臀色熟女| 18+在线观看网站| 超碰av人人做人人爽久久| 看黄色毛片网站| 琪琪午夜伦伦电影理论片6080| 麻豆一二三区av精品| 久久久久久国产a免费观看| 小蜜桃在线观看免费完整版高清| 国产一区二区激情短视频| 中文字幕av在线有码专区| 黄色视频,在线免费观看| 熟妇人妻久久中文字幕3abv| 日本a在线网址| 黄色日韩在线| 欧美成人一区二区免费高清观看| 成年版毛片免费区| 国产精华一区二区三区| 国产乱人视频| 天堂影院成人在线观看| 制服丝袜大香蕉在线| 在线天堂最新版资源| 欧美丝袜亚洲另类 | 看黄色毛片网站| 俺也久久电影网| 久久国产乱子免费精品| 窝窝影院91人妻| 精品久久久久久久末码| 婷婷精品国产亚洲av在线| 男女那种视频在线观看| 两个人视频免费观看高清| 中文字幕av成人在线电影| 日韩欧美国产在线观看| 在线观看美女被高潮喷水网站| av在线蜜桃| 99热这里只有精品一区| 嫩草影院精品99| 国产午夜福利久久久久久| 欧美日韩亚洲国产一区二区在线观看| 色哟哟·www| 黄色日韩在线| 啦啦啦啦在线视频资源| 亚洲内射少妇av| 国产高清不卡午夜福利| 久久久国产成人精品二区| 九色成人免费人妻av| 国产高清不卡午夜福利| 男人狂女人下面高潮的视频| 麻豆国产97在线/欧美| 欧美日韩综合久久久久久 | 97碰自拍视频| 男女那种视频在线观看| 欧美一区二区亚洲| 国产精品1区2区在线观看.| 精品欧美国产一区二区三| 国产三级中文精品| 国产单亲对白刺激| 他把我摸到了高潮在线观看| 熟妇人妻久久中文字幕3abv| 能在线免费观看的黄片| 我的女老师完整版在线观看| 人妻制服诱惑在线中文字幕| 99热6这里只有精品| 亚洲自偷自拍三级| 亚洲经典国产精华液单| 干丝袜人妻中文字幕| 久久99热这里只有精品18| 免费在线观看日本一区| 国产精品女同一区二区软件 | 欧美xxxx黑人xx丫x性爽| 国内少妇人妻偷人精品xxx网站|