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

    改進(jìn)的鍵基正交各向異性近場動力學(xué)模型

    2020-04-09 08:16鄭國君陳瑞申國哲夏陽
    計(jì)算機(jī)輔助工程 2020年1期
    關(guān)鍵詞:單向剪切裂紋

    鄭國君 陳瑞 申國哲 夏陽

    摘要:針對經(jīng)典的鍵基近場動力學(xué)(bond?based peridynamic, BPD)模型受固定泊松比限制的問題,提出一種改進(jìn)BPD模型。該模型可解除泊松比限制,并可用于分析正交各向異性單向板的變形和裂紋擴(kuò)展問題。在改進(jìn)BPD模型中,每根鍵受到軸向和橫向成對力的作用,額外增加的節(jié)點(diǎn)轉(zhuǎn)動可消除由橫向力引起的附加彎矩,從而確保該模型滿足角動量守恒條件。仿真結(jié)果驗(yàn)證所提出的改進(jìn)BPD模型的精度,并展示其預(yù)測碳纖維復(fù)合材料變形和裂紋擴(kuò)展的能力。

    關(guān)鍵詞:

    近場動力學(xué); 鍵基; 剪切影響系數(shù); 泊松比; 顆粒旋轉(zhuǎn); 角動量守恒; 正交各向異性; 碳纖維復(fù)合材料

    中圖分類號:TB334;TP391.99

    文獻(xiàn)標(biāo)志碼:B

    An improved bond?based peridynamic orthotropic model

    ZHENG Guojun, CHEN Rui, SHEN Guozhe, XIA Yang

    (School of Vehicle Engineering, Dalian University of Technology, Dalian 116024, Liaoning, China)

    Abstract:

    As to the problem that the classical bond?based peridynamic(BPD) model is limited by the fixed Poisson′s ratio, an improved BPD model is proposed, which removes the limitation of Poisson′s ratio and can be used to analyze the deformation and crack propagation of orthotropic anisotropic plates. In the improved BPD model, each bond is subjected to axial and transverse forces in pairs, and the additional rotation of particles can eliminate the additional bending momentum caused by transverse forces, so that the balance of angular momentum can be ensured. The simulation results can verify the accuracy of the improved BPD model, and show its ability to predict the deformation and crack propagation in carbon fiber composite material.

    Key words:

    peridynamic; bond?based; shear influence coefficient; Poisson′s ratio; particle rotation; angular momentum balance; orthotropic; carbon fiber composite material

    引?言

    復(fù)合材料是指由2種或者2種以上不同性質(zhì)的材料通過物理或者化學(xué)的方法結(jié)合而成的具有新性能的材料。復(fù)合材料具有質(zhì)量輕、強(qiáng)度高、抗疲勞性好等優(yōu)點(diǎn)[1],在航空航天、汽車、風(fēng)力發(fā)電、醫(yī)療器械等領(lǐng)域的應(yīng)用越來越廣泛,如空中客車A380、波音787等客機(jī)的主要結(jié)構(gòu)(整體機(jī)身、機(jī)翼和艙門等)都采用復(fù)合材料[2]。隨著復(fù)合材料應(yīng)用領(lǐng)域越來越廣泛,其載荷環(huán)境日益復(fù)雜,損傷失效問題日益突出,對其失效模式的分析也成為工程中的重點(diǎn)和難點(diǎn)。

    傳統(tǒng)的連續(xù)介質(zhì)理論具有微分形式的運(yùn)動方程,在面對損傷破壞等不連續(xù)問題時,由于位移場的不連續(xù)性,在不連續(xù)的區(qū)域難以得到位移場的偏微分方程,因此傳統(tǒng)連續(xù)介質(zhì)理論在求解不連續(xù)問題時會遇到困難[3?4]。目前,針對復(fù)合材料結(jié)構(gòu)損傷等問題主要采用有限元法,但是有限元法也是以連續(xù)介質(zhì)理論為基礎(chǔ)的,因此使用有限元法求解損傷失效問題也極具挑戰(zhàn)性,需要借助一些附加的失效準(zhǔn)則重新劃分網(wǎng)格并進(jìn)行求解[5?7],而且損傷只能沿某些特定的方向傳播。為解決傳統(tǒng)有限元在求解損傷裂紋時的缺陷,有學(xué)者提出擴(kuò)展有限元和非連續(xù)有限元等思想[8?14],該思想主要使用獨(dú)立于網(wǎng)格劃分的思想解決裂紋擴(kuò)展問題,不需要對結(jié)構(gòu)內(nèi)部存在的裂紋等缺陷進(jìn)行網(wǎng)格劃分。擴(kuò)展有限元法可提高模型描述復(fù)雜位移場的能力,避免網(wǎng)格的重新劃分,解決大量的斷裂問題。但是,在解決不連續(xù)問題時,擴(kuò)展有限元法仍然需要額外的斷裂準(zhǔn)則。COX等[15]提出分子動力學(xué)理論解決傳統(tǒng)有限元存在的缺陷,但卻增大計(jì)算的消耗。KRNER等[16]根據(jù)原子間的作用力,考慮長程力的影響,提出非局部連續(xù)理論,該理論可不區(qū)分不連續(xù)性,但由于裂紋的存在,其位移導(dǎo)數(shù)仍然是不存在的。

    為解決傳統(tǒng)數(shù)值方法在求解不連續(xù)問題時存在的困難,SILLING[17]和SILLING等[18]提出近場動力學(xué)理論(peridynamic,PD),基于非局部思想,通過求解空間積分方程描述物質(zhì)的力學(xué)行為。

    鍵基近場動力學(xué)(bond?based PD,BPD)模型使用成對力函數(shù)描述節(jié)點(diǎn)間的相互作用,因此對于各向同性材料來說,鍵力僅與兩個節(jié)點(diǎn)間的相對伸長率有關(guān),在平面應(yīng)力狀態(tài)下泊松比被限制為1/3,在平面應(yīng)變狀態(tài)下泊松比被限制為1/4[19]。GERSTLE等[20]通過引入桿單元提出微極模型,通過增加成對力矩模擬具有不同泊松比的線彈性材料,這樣可解除泊松比的限制,但是未考慮鍵的轉(zhuǎn)動,因此不滿足角動量守恒。REN等[21]提出考慮剪切變形的線彈性固體PD模型,從總變形中減去剛體的轉(zhuǎn)動部分,但是未具體處理泊松比的限制。GHAJARI等[22]引入勒讓德多項(xiàng)式,提出正交各向異性材料的連續(xù)微模量函數(shù)模型,能夠預(yù)測復(fù)雜的斷裂現(xiàn)象,但是其泊松比仍然是受限制的。ZHU等[23]考慮鍵轉(zhuǎn)動的影響,重構(gòu)近場動力學(xué)鍵基模型,消除泊松比的限制。ZHOU等[24]提出共軛鍵線彈性模型,鍵能不僅與其法向的伸長有關(guān),也與一對共軛鍵的旋轉(zhuǎn)角度有關(guān),從而克服泊松比的限制。SILLING等[25]提出基于狀態(tài)的PD模型,通過引入節(jié)點(diǎn)的變形狀態(tài),徹底解除泊松比的限制,但其計(jì)算過程比較復(fù)雜。

    對于正交各向異性材料,雖然有些PD模型可解除泊松比情況限制,但未考慮鍵轉(zhuǎn)動的影響,不滿足角動量守恒,在剪切模擬時誤差較大?;跔顟B(tài)的PD模型雖然可有效解除泊松比的限制,但是基于狀態(tài)的PD模型計(jì)算過程比基于鍵的PD模型更復(fù)雜,應(yīng)用更不便。

    本文基于Timoshenko梁理論提出一種改進(jìn)的鍵基梁模型,考慮節(jié)點(diǎn)的旋轉(zhuǎn)效應(yīng)和剪切變形,不僅可以解除泊松比的限制,也可以使橫向剛度計(jì)算更加精確。

    1?BPD基礎(chǔ)理論

    在經(jīng)典的BPD模型中,PD模型的參考構(gòu)型和當(dāng)前構(gòu)型[17]見圖1。參考構(gòu)型中某節(jié)點(diǎn)在位置xi的動力學(xué)方程可以寫為

    式中:ρ為質(zhì)量密度;u為位移矢量場;f為成對力函數(shù),表示節(jié)點(diǎn)xj施加到節(jié)點(diǎn)xi上的每體積的力;V為節(jié)點(diǎn)xi所占據(jù)的體積;b為體力矢量場;Ηxi為節(jié)點(diǎn)xi的作用域。

    為滿足線動量和角動量守恒,對參考構(gòu)型中的任意相對位置矢量ξ和相對位移矢量η有

    f(η,ξ)=-f(-η,-ξ),η,ξ

    (η+ξ)×f(η,ξ)=0,η,ξ

    (2)

    其中:

    ξ=xj-xi

    η=u(xj,t)-u(xi,t)

    (3)

    根據(jù)式(2)和(3)可知,兩個節(jié)點(diǎn)彼此施加的力大小相等、方向相反,且與當(dāng)前構(gòu)型中的相對位置矢量平行。

    在微彈性材料模型中,成對力函數(shù)可以由標(biāo)量微勢能函數(shù)w(η,ξ)推導(dǎo)[17],即

    f(η,ξ,t)=w(η,ξ)η

    (4)

    單個鍵力與相對伸長率之間的線性關(guān)系可以根據(jù)線微彈性理論假設(shè)和微勢能函數(shù)推導(dǎo),即

    w(η,ξ)=c(η,ξ)s2ξ/2(5)

    式中:c(η,ξ)為鍵剛度常數(shù);s為鍵的伸長;ξ為ξ的模。s的定義為

    s=(ξ+η-ξ)/ξ

    (6)

    對于微彈性材料來說,給定節(jié)點(diǎn)xi處的應(yīng)變能密度可通過對作用域內(nèi)的微勢能函數(shù)w(η,ξ)積分得到,即

    WPD(xi)=12∫Ηxiw(η,ξ)dVxj

    (7)

    式中的1/2是因?yàn)槊總€節(jié)點(diǎn)的鍵能只有總鍵能的1/2。相同載荷作用下PD得到的應(yīng)變能密度應(yīng)等于經(jīng)典彈性力學(xué)理論得到的應(yīng)變能密度,即

    WPD(xi)=WCL(xi)

    (8)

    從而可以得到鍵剛度系數(shù)與彈性模量之間的關(guān)系。

    2?改進(jìn)的BPD模型

    2.1?求解鍵常數(shù)

    對于改進(jìn)的BPD模型,由節(jié)點(diǎn)xi和在其近場區(qū)域Ηxi內(nèi)的相鄰點(diǎn)xj組成PD鍵,見圖2。

    與經(jīng)典的BPD模型相同,假設(shè)鍵的軸向力只與鍵的拉伸變形和軸向力密度有關(guān),即

    f^xij=-f^xji=cNs

    (9)

    式中:f^xij為在近場區(qū)域Ηxi內(nèi)的節(jié)點(diǎn)xj施加在節(jié)點(diǎn)xi上的軸向力密度;f^xji為在近場區(qū)域Ηxj內(nèi)的節(jié)點(diǎn)xi施加在節(jié)點(diǎn)xj上的軸向力密度;cN為要求解的PD參數(shù)。

    考慮用鍵的橫向力消除泊松比的限制。鍵的橫向力附著在兩個節(jié)點(diǎn)上,會出現(xiàn)附加的彎矩,必須對其進(jìn)行消除以滿足角動量守恒。類似于有限元模型中的梁單元,可通過增加兩個節(jié)點(diǎn)的旋轉(zhuǎn)自由度抵消由橫向力引起的彎矩。鍵的橫向位移、力和節(jié)點(diǎn)的旋轉(zhuǎn)角度及彎矩關(guān)系見圖2。與軸向力相同,在變形構(gòu)型上建立橫向力和彎矩。假設(shè)PD鍵是具有長度ζ和高度Δ的梁模型,高度相對于長度的比值Δ/ζ沒有小到足以忽略剪切變形的影響。Timoshenko梁理論可描述剪切變形和轉(zhuǎn)動慣量對各種細(xì)長梁的影響。[26]剪切影響系數(shù)可解釋剪切應(yīng)力的變化,適于描述短梁的剪切變形性能。[27]在局部坐標(biāo)系中,基于Timoshenko梁理論的有限元方程為

    p=Kd

    (10)

    式中:p為附著在材料節(jié)點(diǎn)xi和xj之間的鍵的橫向力和彎矩分量;d為鍵的位移分量;K為Timoshenko梁單元的局部剛度矩陣[28?29]。p、d和K分別定義為

    p=[fixfiymizfjxfjymjz]T(11)

    d=[uiviθiujvjθj]T

    (12)

    K=EAξ00-EAξ00012EI(1+b)ξ36EI(1+b)ξ20-12EI(1+b)ξ36EI(1+b)ξ206EI(1+b)ξ2(4+b)EI(1+b)ξ0-6EI(1+b)ξ2(2-b)EI(1+b)ξ-EAξ00EAξ000-12EI(1+b)ξ3-6EI(1+b)ξ2012EI(1+b)ξ3-6EI(1+b)ξ206EI(1+b)ξ2(2-b)EI(1+b)ξ0-6EI(1+b)ξ2(4+b)EI(1+b)ξ

    (13)

    式中:E為彈性模量;G為剪切模量;I為梁的轉(zhuǎn)動慣量;b為剪切影響系數(shù)。b可使梁的橫向剛度更加準(zhǔn)確,其定義為

    b=6EΔ2/5Gξ2(14)

    引入節(jié)點(diǎn)的附加旋轉(zhuǎn)角度和彎矩以計(jì)算鍵的應(yīng)變能,基于Timoshenko梁理論的PD模型見圖3。

    由鍵連接的兩個節(jié)點(diǎn)之間的位移是相互獨(dú)立的,初始的相對位置矢量為ξ,相對位移矢量為η,x′軸與全局坐標(biāo)系的x軸夾角為。

    在單個鍵中,PD的軸向力密度、橫向力密度和彎矩密度可以表示為

    f^y=12d(1+b)ξ2ξ-θ^1+θ^22

    M^z=M^1z+M^2z2=6c(1+b)ξξ-θ^1+θ^22

    (15)

    式中:c和d為要求解的PD參數(shù);和分別為軸向和橫向位移;θ^1和θ^2為局部坐標(biāo)系中節(jié)點(diǎn)的旋轉(zhuǎn)角度;b為剪切影響系數(shù)。因此,

    =2-1

    =2-1

    (16)

    那么,存儲在由相鄰兩個節(jié)點(diǎn)xi和xj形成的PD鍵內(nèi)的應(yīng)變能密度為

    w(η,ξ)=dTKd/2=

    wf^x(η,ξ)+

    wf^y(η,ξ)+wM^z(η,ξ)

    (17)

    式中:wf^x(η,ξ)為由軸向變形引起的應(yīng)變能密度;wf^y(η,ξ)為由剪切變形引起的應(yīng)變能密度;wM^z(η,ξ)為由節(jié)點(diǎn)旋轉(zhuǎn)力矩引起的應(yīng)變能密度。

    wf^x(η,ξ)=c2ξ2

    wf^y(η,ξ)=12d2(1+b)ξ22ξ-(θ^1+θ^2)2

    wM^z(η,ξ)=d2(1+b)ξ-6ξ(θ^1+θ^2)+

    (4+b)(θ^21+θ^22)+(4-2b)θ^1θ^2

    (18)

    節(jié)點(diǎn)的旋轉(zhuǎn)角度與坐標(biāo)系無關(guān),因此定義

    θ=(θ^1+θ^2)/2 (19)

    那么應(yīng)變能密度方程可以化簡為

    w(η,ξ)=

    12cξ2+12d(1+b)ξξ-θ2+dξ(θ^1-θ^2)2

    (20)

    在局部小變形過程中,可以假設(shè)兩個節(jié)點(diǎn)的旋轉(zhuǎn)角度是相等的,即

    θ^1=θ^2(21)

    因此,式(20)可化簡為

    w(η,ξ)=12cξ2+12d(1+b)ξξ-θ2

    (22)

    在x′和y′方向上的應(yīng)變分別為ε1和ε2,那么在局部坐標(biāo)系下的應(yīng)變定義為

    ε1ε2γ12=cos2 sin2 sin cos sin2 cos2 -sin cos -2sin cos 2sin cos cos2 -sin2 εxεyγxy

    (23)

    由局部坐標(biāo)系與全局坐標(biāo)系關(guān)系可知

    12=ξε1=ξ(εxcos2 θ+εysin2 θ+γxysin θcos θ)

    12=ξ12γ12+θ=ξ-εxsin θcos θ+εysin θcos θ+γxy(cos2 θ-sin2 θ)2+θ

    (24)

    能量密度

    WPD=12∫H12cξε21+12d(1+b)ξ12γ122tdH

    (25)

    正交各向異性單向板材料可以作為二維問題處理,其PD模型見圖4,厚度方向僅考慮單層物質(zhì)點(diǎn),面內(nèi)分為纖維方向和基體方向,纖維鍵與基體鍵在拉伸方向的鍵剛度分別設(shè)為cf和cm,纖維鍵與基體鍵在剪切方向的鍵剛度分別設(shè)為df和dm,任意角度的鍵剛度滿足

    c=cf+cm,θ=

    c=cm,θ≠

    d=df+dm,θ=

    d=dm,θ≠

    (26)

    式中:θ為纖維方向與x軸的夾角。

    因此,式(25)中的應(yīng)變能密度可以近似表示為

    WPD=12∑Qq=112cfξε21+12df(1+b2)ξ12γ122tVq+

    12∫H12cmξε21+12dm(1+b1)ξ12γ122tdH

    (27)

    式中:b1和b2分別為基體方向和纖維方向的剪切影響因子。

    在經(jīng)典連續(xù)介質(zhì)力學(xué)中,正交各向異性材料的平面應(yīng)力狀態(tài)主方向的應(yīng)力應(yīng)變關(guān)系為

    σ1σ2τ12=C11C120C21C22000C66ε1ε2γ12

    (28)

    式中:Cij為剛度矩陣。因此,經(jīng)典連續(xù)介質(zhì)理論的應(yīng)變能密度為

    WCL=12Cijεiεj=12C11ε21+C12ε1ε2+

    12C22ε22+12C66γ212

    (29)

    式中:i和j均取1、2、6,ε6=γ12。

    PD得到的應(yīng)變能密度應(yīng)等于經(jīng)典連續(xù)介質(zhì)理論的應(yīng)變能密度,可求得cm、cf、dm和df分別為

    cf=2(E1-E2)(1-ν12ν21)Qq=1ξVq

    df=2G12(1-ν12ν21)-(1-ν12)E2Qq=15G12ξ25G12ξ2+6E1Δ2×3(1-ν12ν21)Qq=1Vqξ

    cm=6(1+ν12)E2(1-ν12ν21)tπδ3

    dm=λ(1-3ν12)E26πδt(1-ν12ν21)

    (30)

    式中:ν為泊松比;λ為調(diào)整系數(shù),

    λ=ΦΦ-arctan Φ,?Φ=δΔ5G126E2

    (31)

    2.2?失效準(zhǔn)則

    對于二維平面問題,破壞每單位斷裂面積的所有鍵所需要的功等于材料的臨界應(yīng)變能釋放率,由此可以計(jì)算臨界拉伸。[30]對于正交各向異性單向板PD模型,鍵失效定義為纖維鍵和基體鍵失效。在纖維方向和基體方向的鍵的臨界伸長示意見圖5。對于纖維方向的鍵,如果鍵伸長超過臨界伸長sft,則纖維鍵發(fā)生斷裂;對于其他方向的所有鍵,如果鍵伸長超過臨界伸長smt,則基體鍵發(fā)生斷裂。令均勻介質(zhì)中單位斷裂區(qū)域的所有鍵斷裂所需要的功等于應(yīng)變能釋放率GI,Cr,1(即垂直于纖維方向的Ι型斷裂應(yīng)變能釋放率[31]),可計(jì)算臨界伸長sft;等效相應(yīng)的功和應(yīng)變能釋放率GI,Cr,2(即沿著纖維方向的Ι型斷裂應(yīng)變能釋放率[31]),可計(jì)算臨界伸長smt。臨界伸長sft和smt可以表示為

    sft=4πGI,Cr,19E1δ,?smt=4πGI,Cr,29E2δ

    (32)

    節(jié)點(diǎn)xi與其近場區(qū)域內(nèi)的任意點(diǎn)xj之間有鍵的相互作用,當(dāng)鍵的伸長超過臨界伸長s0時,鍵會發(fā)生斷裂。定義標(biāo)量函數(shù)μ描述鍵是否發(fā)生斷裂,見式(33)。當(dāng)鍵伸長小于臨界伸長時μ=1,意味著鍵未發(fā)生斷裂,否則鍵斷裂并且相應(yīng)的鍵力變?yōu)?。

    μ(ξ,t)=1,s(ξ,t′)

    (33)

    將節(jié)點(diǎn)xi的損傷定義為函數(shù)φ,見式(34)。若在其近場區(qū)域內(nèi)沒有鍵的斷裂,則損傷值為0;當(dāng)與節(jié)點(diǎn)xi連接的所有鍵都斷開時,節(jié)點(diǎn)的損傷值為1;當(dāng)損傷函數(shù)達(dá)到0.5時產(chǎn)生裂縫。

    φ(x,t)=1-∫H μ(ξ,t)dV′∫HdV′

    (34)

    3?數(shù)值結(jié)果

    3.1?碳纖維單向板拉伸

    為驗(yàn)證模型的有效性,對不同纖維方向的單向板進(jìn)行數(shù)值模擬,分別考慮纖維方向θ為0°、15°、30°、45°和90°的單向板,板長為100 mm、寬為20 mm、厚為1 mm。初始坐標(biāo)設(shè)置在薄板左邊的中心點(diǎn),在左側(cè)建立寬度為3Δ的區(qū)域作為約束區(qū)域,右端施加靜力拉伸載荷,力的總和為F=42 000 N。纖維方向的彈性模量E1=41.0 GPa,基體方向的彈性模量E2=10.4 GPa,泊松比ν12=0.28,剪切模量G12=4.3 GPa,密度ρ=1 970 kg/m3。近場區(qū)域均勻離散為1 mm×1 mm正方形物質(zhì)點(diǎn)的集合,近場區(qū)域設(shè)置為δ=3Δ。使用OpitiStruct進(jìn)行數(shù)值仿真模擬,有限元計(jì)算同樣也采用1 mm×1 mm正方形網(wǎng)格,將均勻分布在矩形板上的7個作用點(diǎn)作為測量點(diǎn),見圖6。

    為驗(yàn)證模型的有效性,繪制具有不同纖維方向的正交各向異性單向板沿x軸的位移云圖,與有限元分析進(jìn)行對比,結(jié)果見圖7~11。

    PD模型位移云圖與有限元模型位移云圖結(jié)果吻合較好,說明該模型準(zhǔn)確。

    3.2?碳纖維單向板裂紋擴(kuò)展

    基于本文提出的BPD模型,引入含有中心裂紋缺陷的碳纖維單向板,模擬含缺陷碳纖維單向板的漸進(jìn)損傷,通過研究單向板的損傷擴(kuò)展路徑和最終破壞模式,驗(yàn)證該方法的可行性和有效性。

    考慮長為70.0 mm、寬為40.0 mm、厚為0.6 mm的單向板,初始坐標(biāo)設(shè)置在薄板的中心點(diǎn),纖維方向的彈性模量E1=105 GPa,基體方向的彈性模量E2=8.4 GPa,泊松比ν12=0.32,剪切模量G12=4 GPa,密度ρ=1 800 kg/m3。近場區(qū)域設(shè)置為δ=3.015Δ,均勻離散為0.5 mm×0.5 mm正方形物質(zhì)點(diǎn)的集合。缺陷在試件的幾何中心,裂紋長度2a=8.0 mm,在兩側(cè)施加2.0 mm的位移載荷。含線裂紋碳纖維單向板受拉伸載荷作用示意見圖12。纖維鍵的臨界伸長sft=0.03,基體鍵的臨界伸長smt=0.02。纖維角度為45°時,其漸進(jìn)損傷結(jié)果對比見圖13。

    當(dāng)纖維方向?yàn)?5°時,基體首先損傷且損傷開始于初始裂紋的兩端;隨著位移的增大,裂紋沿著纖維方向向兩端擴(kuò)展,當(dāng)其達(dá)到纖維的臨界伸長時,纖維發(fā)生斷裂。PD仿真結(jié)果與試驗(yàn)結(jié)果[32]較吻合,可驗(yàn)證該模型的有效性和可行性,同時也表明PD可以模擬碳纖維單向板裂紋損傷,凸顯PD理論在

    模擬裂紋損傷方面的優(yōu)越性。

    4?結(jié)?論

    提出一種改進(jìn)的BPD模型用于克服泊松比受到限制的問題,假設(shè)節(jié)點(diǎn)受到軸向和橫向力的作用,考慮節(jié)點(diǎn)的旋轉(zhuǎn)以消除由橫向力引起的額外彎矩,從而確保該模型滿足角動量守恒條件。由于鍵的高度與長度的比值沒有小到足以忽略的程度,因此在該模型中考慮鍵的剪切影響系數(shù),提高其精度。

    通過模擬不同纖維方向的碳纖維單向板的單向拉伸,驗(yàn)證改進(jìn)的碳纖維單向板BPD模型的準(zhǔn)確性和解除泊松比限制的特點(diǎn)。通過模擬不同纖維角度矩形板的裂紋傳播路經(jīng),驗(yàn)證模型在裂紋損傷模擬中的可行性和準(zhǔn)確性。對比結(jié)果表明,數(shù)值模擬結(jié)果和有限元仿真結(jié)果均與試驗(yàn)結(jié)果一致,可得出以下結(jié)論:

    (1)改進(jìn)的鍵基PD模型通過在鍵的兩個節(jié)點(diǎn)上加入成對的橫向力,可消除泊松比的限制。這是因?yàn)橛?個獨(dú)立的PD參數(shù)對應(yīng)于復(fù)合材料4個宏觀材料常數(shù),分別為彈性模量、泊松比、剪切模量和密度。

    (2)通過考慮節(jié)點(diǎn)的旋轉(zhuǎn),確保該模型滿足角動量守恒,是該模型能夠精確模擬變形的關(guān)鍵。

    (3)鍵的剪切影響系數(shù)可以提高單軸拉伸下板件的橫向變形精度,對模擬裂紋損傷和擴(kuò)展至關(guān)重要。

    參考文獻(xiàn):

    [1]?杜善義. 先進(jìn)復(fù)合材料與航空航天[J]. 復(fù)合材料學(xué)報, 2007, 24(1): 1?12. DOI: 10.3321/j.issn:1000?3851.2007.01.001.

    [2]?任鵬. 復(fù)合材料層合板沖擊損傷及沖擊后疲勞壽命研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 2014.

    [3]?LIEW K M, ZHAO X, FERREIRA A J M. A review of meshless methods for laminated and functionally graded plates and shells[J]. Composite Structures, 2011, 93(8): 2031?2041. DOI: 10.1016/j.compstruct.2011.02.018.

    [4]?NGUYEN V P, RABCZUK T, BORDAS S, et al. Meshless methods: A review and computer implementation aspects[J]. Mathematics and Computers in Simulation, 2008, 79(3): 763?813. DOI: 10.1016/j.matcom.2008.01.003.

    [5]?MUROTANI K, YAGAWA G, CHOI J B. Adaptive finite elements using hierarchical mesh and its application to crack propagation analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 1?14. DOI: 10.1016/j.cma.2012.07.024.

    [6]?SCHLLMANN M, FULLAND M, RICHARD H A. Development of a new software for adaptive crack growth simulations in 3D structures[J]. Engineering Fracture Mechanics, 2003, 70(2): 249?268. DOI: 10.1016/S0013?7944(02)00028?0.

    [7]?EL KHAOULANI R, BOUCHARD P O. An anisotropic mesh adaptation strategy for damage and failure in ductile materials[J]. Finite Elements in Analysis and Design, 2012, 59: 1?10. DOI: 10.1016/j.finel.2012.04.006.

    [8]?BELYTSCHKO T, BLACK T. Elastic crack growth in finite elements with minimal remeshing[J]. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601?620.

    [9]?COMBESCURE A, GRAVOUIL A, GRéGOIRE D, et al. X?FEM a good candidate for energy conservation in simulation of brittle dynamic crack propagation[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(5): 309?318. DOI: 10.1016/j.cma.2007.04.011.

    [10]?SHI J X, CHOPP D, LUA J, et al. Abaqus implementation of extended finite element method using a level set representation for three?dimensional fatigue crack growth and life predictions[J]. Engineering Fracture Mechanics, 2010, 77(14): 2840?2863. DOI: 10.1016/j.engfracmech.2010.06.009.

    [11]?GENIAUT S, GALENNE E. A simple method for crack growth in mixed mode with X?FEM[J]. International Journal of Solids and Structures, 2012, 49(15/16): 2094?2106. DOI: 10.1016/j.ijsolstr.2012.04.015.

    [12]?GRGOIRE D, MAIGRE H, RTHOR J, et al. Dynamic crack propagation under mixed?mode loading: Comparison between experiments and X?FEM simulations[J]. International Journal of Solids and Structures, 2007, 44(20): 6517?6534. DOI: 10.1016/j.ijsolstr.2007.02.044.

    [13]?PRABEL B, MARIE S, COMBESCURE A. Using X?FEM method to model dynamic propagation and arrest of cleavage cracks in ferritic steel[J]. Engineering Fracture Mechanics, 2008, 75 (10): 2984?3009. DOI: 10.1016/j.engfracmech.2008.01.008.

    [14]?MENK A, BORDAS S P A. Crack growth calculations in solder joints based on microstructural phenomena with X?FEM[J]. Computational Materials Science, 2011, 50(3): 1145?1156. DOI: 10.1016/j.commatsci.2010.11.014.

    [15]?COX B N, GAO H J, GROSS D, et al. Modern topics and challenges in dynamic fracture[J]. Journal of Mechanics and Physics of Solids, 2005, 53(3): 565?596. DOI: 10.1016/j.jmps.2004.09.002.

    [16]?KRNER E. Elasticity theory of materials with long range cohesive forces[J]. International Journal of Solids and Structures, 1967, 3(5): 731?742. DOI: 10.1016/0020?7683(67)90049?2.

    [17]?SILLING S A. Reformulation of elasticity theory for discontinuities and long?range forces[J]. Journal of Mechanics and Physics of Solids, 2000, 48(1): 175?209. DOI: 10.1016/S0022?5096(99)00029?0.

    [18]?SILLING S A, ASKARI E. A meshfree method based on peridynamic model of solid mechanics[J]. Computers & Structures, 2005, 83(17/18): 1526?1535. DOI: 10.1016/j.compstruc.2004.11.026.

    [19]?WALTER G, NICOLAS S, STEWART S. Peridynamic modeling of plain and reinforced concrete structures[C]// Proceedings of 18th International Conference on Structural Mechanics in Reactor Technology(SMiRT 18). Bejing: IASMiRT, 2005.

    [20]?GERSTLE W, SAU N, SILLING S. Peridynamic modeling of concrete structures[J]. Nuclear Engineering and Design, 2007, 237(12/13): 1250?1258. DOI: 10.1016/j.nucengdes.2006.10.002.

    [21]?REN H L, ZHUANG X Y, RABCZUK T. A new peridynamic formulation with shear deformation for elastic solid[J]. Journal of Micromechanics and Molecular Physics, 2016, 1(2): 1?28. DOI: 10.1142/S2424913016500090.

    [22]?GHAJARI M, IANNUCCI L, CURTIS P. A peridynamic material model for analysis of dynamic crack propagation in orthotropic media[J]. Computer Methods in Applied Mechanics and Engineering, 2014, 276: 431?452. DOI: 10.1016/j.cma.2014.04.002.

    [23]?ZHU Q Z, NI T. Peridynamic formulations enriched with bond rotation effects[J]. International Journal of Engineering Science, 2017, 121: 118?129. DOI: 10.1016/j.ijengsci.2017.09.004.

    [24]?ZHOU X P, WANG Y T, SHOU Y D, et al. A novel conjugated bond linear elastic model in bond?based peridynamics for fracture problems under dynamic loads[J]. Engineering Fracture Mechanics, 2018, 188: 151?183. DOI: 10.1016/j.engfracmech.2017.07.031.

    [25]?SILLING S A, EPTON M, WECKNER O, et al. Peridynamic states and constitutive modeling[J]. Journal of Elasticity, 2007, 88(2): 151?184. DOI: 10.1007/s10659?007?9125?1.

    [26]?TIMOSHENKO S P. LXVI: On correction for shear of differential equation for transverse vibrations of prismatic bars[J]. Philosophical Magazine, 1921, 41(245): 744?746. DOI: 10.1080/14786442108636264.

    [27]?KIM J U, RENARDY Y. Boundary control of Timoshenko beam[J]. SIAM Journal on Control and Optimization, 1987, 25(6): 1417?1429. DOI: 10.1137/0325078.

    [28]?REDDY J N. On dynamic behaviour of Timoshenko beam finite elements[J]. Sadhana, 1999, 24(3): 175?198. DOI: 10.1007/BF02745800.

    [29]?CHANDRUPATLA T R. Introduction to finite elements in engineering[M]. Gewerbestrasse: Pearson Schweiz Ag, 2010: 69.

    [30]?HA Y D, BOBARU F. Studies of dynamic crack propagation and crack branching with peridynamics[J]. International Journal of Fracture, 2010, 162(1/2): 229?244. DOI: 10.1007/s10704?010?9442?4.

    [31]?PINEDA E, WAAS A. A Thermodynamically?based mesh objective work potential theory for predicting intralaminar progressive damage and failure in fiber?reinforced laminates [C]// Proceedings of 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference. Honolulu: AIAA, 2012. DOI: 10.2514/6.2012?1612.

    [32]?劉貞谷. 基于近場動力學(xué)的復(fù)合材料單向板漸進(jìn)損傷模擬與實(shí)驗(yàn)研究[D]. 哈爾濱: 哈爾濱工業(yè)大學(xué), 2017.

    (編輯?武曉英)

    猜你喜歡
    單向剪切裂紋
    有了裂紋的玻璃
    連退飛剪剪切定位控制研究與改進(jìn)
    T應(yīng)力對脆性材料初始裂紋起裂角影響研究
    單向空間
    心生裂紋
    單向街書店:智力、思想和文化生活的公共空間
    機(jī)械法制備石墨烯球磨過程的能量層分析
    知世界
    剪切與沖切的辨析及工程算例
    一種解決軋機(jī)和剪切線共享一臺WINCC服務(wù)器所產(chǎn)生的弊端的方法
    成人特级av手机在线观看| 国产综合懂色| 国产欧美日韩精品一区二区| 九九在线视频观看精品| 哪个播放器可以免费观看大片| 一边亲一边摸免费视频| 午夜精品在线福利| 午夜精品国产一区二区电影 | 亚洲精品亚洲一区二区| 国产成人福利小说| 男人和女人高潮做爰伦理| 亚洲第一电影网av| 高清午夜精品一区二区三区 | 欧美一级a爱片免费观看看| .国产精品久久| 黄片wwwwww| 日韩强制内射视频| 非洲黑人性xxxx精品又粗又长| 天堂网av新在线| 欧美日韩综合久久久久久| 亚洲欧美精品专区久久| 久久婷婷人人爽人人干人人爱| 国产麻豆成人av免费视频| 亚洲熟妇中文字幕五十中出| 亚洲熟妇中文字幕五十中出| 国产亚洲91精品色在线| 一边亲一边摸免费视频| 国国产精品蜜臀av免费| 青青草视频在线视频观看| 欧美成人a在线观看| 国产亚洲精品久久久com| 麻豆精品久久久久久蜜桃| 99久久人妻综合| 精品国内亚洲2022精品成人| 亚洲性久久影院| 久久久欧美国产精品| 亚洲经典国产精华液单| 亚洲一区二区三区色噜噜| 人人妻人人澡欧美一区二区| 免费av毛片视频| 丰满乱子伦码专区| 免费搜索国产男女视频| 精品人妻一区二区三区麻豆| 日本-黄色视频高清免费观看| 亚洲综合色惰| 精品一区二区三区视频在线| 看黄色毛片网站| 中文字幕av成人在线电影| 床上黄色一级片| 国产精品久久久久久精品电影小说 | 内地一区二区视频在线| 久久久久久久久久久免费av| 久久久久久久久久黄片| 熟女人妻精品中文字幕| 成人永久免费在线观看视频| 国产一区二区在线观看日韩| 午夜老司机福利剧场| av黄色大香蕉| 精品熟女少妇av免费看| av又黄又爽大尺度在线免费看 | 麻豆乱淫一区二区| 国产亚洲精品av在线| 久久久午夜欧美精品| 久久久久久大精品| 久久精品国产99精品国产亚洲性色| 国产色爽女视频免费观看| 久久久久国产网址| 亚洲欧美精品专区久久| 变态另类丝袜制服| 男女视频在线观看网站免费| 最近手机中文字幕大全| 久久人人爽人人片av| 久久久久久久久久久丰满| 99热这里只有是精品50| 麻豆成人av视频| 国产亚洲精品久久久久久毛片| 在线免费十八禁| 99热这里只有是精品在线观看| 欧美高清性xxxxhd video| 赤兔流量卡办理| www.av在线官网国产| 色哟哟哟哟哟哟| 波多野结衣高清作品| 一区二区三区免费毛片| 国内精品久久久久精免费| 日韩欧美 国产精品| 国产一区二区激情短视频| 亚洲人与动物交配视频| 天天一区二区日本电影三级| 爱豆传媒免费全集在线观看| 一级黄片播放器| 国产伦理片在线播放av一区 | 国产亚洲av嫩草精品影院| 国产高清三级在线| av在线蜜桃| 床上黄色一级片| 美女 人体艺术 gogo| 免费av观看视频| 日本与韩国留学比较| 天美传媒精品一区二区| 狂野欧美白嫩少妇大欣赏| 国产精品久久电影中文字幕| 观看美女的网站| 亚洲av.av天堂| 国产高清激情床上av| 在线观看美女被高潮喷水网站| 久久国内精品自在自线图片| 国产精品伦人一区二区| 亚洲av中文字字幕乱码综合| .国产精品久久| 亚洲精品自拍成人| 亚洲国产欧洲综合997久久,| 能在线免费看毛片的网站| av女优亚洲男人天堂| 亚洲av不卡在线观看| 黄色欧美视频在线观看| 欧美人与善性xxx| 亚洲av成人精品一区久久| 欧美人与善性xxx| 两个人的视频大全免费| 在线观看一区二区三区| 国产一区二区亚洲精品在线观看| 天美传媒精品一区二区| 秋霞在线观看毛片| 久久久久免费精品人妻一区二区| 国产成人91sexporn| 亚洲av免费高清在线观看| 久久久久久久亚洲中文字幕| 少妇熟女aⅴ在线视频| 国产精品麻豆人妻色哟哟久久 | 国产亚洲精品av在线| 亚洲av二区三区四区| 2022亚洲国产成人精品| 狠狠狠狠99中文字幕| 99热这里只有精品一区| 亚洲中文字幕日韩| 午夜亚洲福利在线播放| 免费观看精品视频网站| eeuss影院久久| 日韩高清综合在线| 尾随美女入室| 黑人高潮一二区| 最近2019中文字幕mv第一页| 亚洲,欧美,日韩| 在线观看一区二区三区| 亚洲人与动物交配视频| 变态另类成人亚洲欧美熟女| 淫秽高清视频在线观看| 又爽又黄无遮挡网站| 国产黄片美女视频| 国产国拍精品亚洲av在线观看| 一级黄片播放器| 国产精品国产高清国产av| 免费av毛片视频| 亚洲成人久久爱视频| 日韩国内少妇激情av| 日韩制服骚丝袜av| 亚洲国产色片| 免费无遮挡裸体视频| 国产精品人妻久久久影院| 午夜激情欧美在线| 在线a可以看的网站| 国产黄a三级三级三级人| 色视频www国产| 亚洲av二区三区四区| АⅤ资源中文在线天堂| 欧美高清性xxxxhd video| 我的老师免费观看完整版| 听说在线观看完整版免费高清| 亚洲av不卡在线观看| 国产午夜福利久久久久久| 淫秽高清视频在线观看| 日本一本二区三区精品| 亚洲成人中文字幕在线播放| 99在线视频只有这里精品首页| 亚洲国产日韩欧美精品在线观看| 免费不卡的大黄色大毛片视频在线观看 | 国产精品一及| 美女高潮的动态| 国语自产精品视频在线第100页| 高清日韩中文字幕在线| 性欧美人与动物交配| 熟女电影av网| 嫩草影院精品99| 高清毛片免费看| 一进一出抽搐gif免费好疼| 偷拍熟女少妇极品色| 亚洲第一电影网av| 日韩av不卡免费在线播放| 少妇猛男粗大的猛烈进出视频 | 国产黄片视频在线免费观看| 国产精品野战在线观看| 免费电影在线观看免费观看| 国产一级毛片在线| 久久精品人妻少妇| 少妇猛男粗大的猛烈进出视频 | 欧美在线一区亚洲| 看黄色毛片网站| 菩萨蛮人人尽说江南好唐韦庄 | 我的老师免费观看完整版| 综合色丁香网| 成人一区二区视频在线观看| 久久婷婷人人爽人人干人人爱| 国产美女午夜福利| 久久久久久久久久久免费av| 美女xxoo啪啪120秒动态图| 国产久久久一区二区三区| 国产爱豆传媒在线观看| 一本久久中文字幕| 国内精品一区二区在线观看| 亚洲欧洲日产国产| 天堂影院成人在线观看| 村上凉子中文字幕在线| 久久国内精品自在自线图片| 欧美激情国产日韩精品一区| 国产精品三级大全| 国产大屁股一区二区在线视频| 男人舔奶头视频| 亚洲精品影视一区二区三区av| 精品久久久久久久人妻蜜臀av| 亚洲欧洲国产日韩| 亚洲精品久久国产高清桃花| av在线老鸭窝| 国产老妇伦熟女老妇高清| 蜜桃亚洲精品一区二区三区| 欧美高清成人免费视频www| 日本熟妇午夜| 在线观看av片永久免费下载| 长腿黑丝高跟| 能在线免费观看的黄片| 国产亚洲精品久久久久久毛片| 国产精品精品国产色婷婷| 午夜久久久久精精品| 亚洲真实伦在线观看| 给我免费播放毛片高清在线观看| 国产精品.久久久| 亚洲国产欧美人成| 十八禁国产超污无遮挡网站| 日韩精品有码人妻一区| 美女 人体艺术 gogo| 国产熟女欧美一区二区| 欧美在线一区亚洲| 两性午夜刺激爽爽歪歪视频在线观看| 欧美一区二区国产精品久久精品| 久久久精品大字幕| 日韩欧美一区二区三区在线观看| 欧美成人精品欧美一级黄| 精品一区二区三区人妻视频| 国产成人影院久久av| 国产成人一区二区在线| 国产一区二区亚洲精品在线观看| 亚洲欧美日韩高清在线视频| 丰满的人妻完整版| 欧美高清性xxxxhd video| 午夜福利在线在线| 日韩 亚洲 欧美在线| 波野结衣二区三区在线| 狠狠狠狠99中文字幕| 成人午夜精彩视频在线观看| 深夜精品福利| 最近最新中文字幕大全电影3| 亚洲成人精品中文字幕电影| 日本撒尿小便嘘嘘汇集6| 日本黄大片高清| 一个人观看的视频www高清免费观看| 免费av观看视频| av又黄又爽大尺度在线免费看 | 免费看日本二区| 久久久久久九九精品二区国产| 美女内射精品一级片tv| 看免费成人av毛片| 亚洲av一区综合| 免费人成在线观看视频色| www.色视频.com| 久久久国产成人精品二区| 成人毛片a级毛片在线播放| 亚洲精品久久国产高清桃花| 三级经典国产精品| 日日摸夜夜添夜夜爱| 午夜免费男女啪啪视频观看| 国产成人福利小说| 青春草视频在线免费观看| 丰满的人妻完整版| 国内久久婷婷六月综合欲色啪| 少妇人妻一区二区三区视频| 亚洲成a人片在线一区二区| 边亲边吃奶的免费视频| 久久久久网色| 欧美区成人在线视频| 岛国在线免费视频观看| 亚洲,欧美,日韩| a级毛片免费高清观看在线播放| 久久人妻av系列| 亚洲欧美成人综合另类久久久 | 亚洲在久久综合| 国产视频首页在线观看| 99久久精品一区二区三区| 久久久久久伊人网av| 成年免费大片在线观看| 2022亚洲国产成人精品| 久久人人爽人人片av| 国产高潮美女av| 春色校园在线视频观看| 精品不卡国产一区二区三区| 97热精品久久久久久| 黄色配什么色好看| 狠狠狠狠99中文字幕| av在线亚洲专区| 亚洲欧洲日产国产| 国产激情偷乱视频一区二区| 久久久久久国产a免费观看| 国产高清三级在线| 亚洲精品色激情综合| 99久久人妻综合| 看非洲黑人一级黄片| 久久精品久久久久久久性| 亚洲精品亚洲一区二区| 99热6这里只有精品| 欧美精品一区二区大全| 日韩亚洲欧美综合| 亚洲欧美精品自产自拍| 99久国产av精品国产电影| 国产高清不卡午夜福利| 亚洲aⅴ乱码一区二区在线播放| 伦精品一区二区三区| 噜噜噜噜噜久久久久久91| 久久久久久久亚洲中文字幕| 国产精品一区二区三区四区久久| 亚洲人成网站高清观看| 国产精品人妻久久久影院| 中文在线观看免费www的网站| 成人国产麻豆网| 成人二区视频| 3wmmmm亚洲av在线观看| 一卡2卡三卡四卡精品乱码亚洲| 久久人人精品亚洲av| 久久久久性生活片| 国产精品福利在线免费观看| 亚洲成a人片在线一区二区| 亚洲精品日韩在线中文字幕 | 日韩中字成人| 中文字幕制服av| 成人三级黄色视频| 欧美成人精品欧美一级黄| 91aial.com中文字幕在线观看| 看黄色毛片网站| 久久久久久久午夜电影| 亚洲成av人片在线播放无| 色哟哟·www| 欧美一级a爱片免费观看看| 日韩三级伦理在线观看| 大又大粗又爽又黄少妇毛片口| 爱豆传媒免费全集在线观看| 精品久久久噜噜| 免费电影在线观看免费观看| 国产日韩欧美在线精品| 长腿黑丝高跟| 欧美变态另类bdsm刘玥| 舔av片在线| 亚洲人成网站在线观看播放| 国产一区二区三区在线臀色熟女| 国产伦一二天堂av在线观看| 亚洲七黄色美女视频| 少妇熟女aⅴ在线视频| 69人妻影院| 精品少妇黑人巨大在线播放 | 中文字幕av在线有码专区| 麻豆国产97在线/欧美| 成人综合一区亚洲| 97在线视频观看| 小蜜桃在线观看免费完整版高清| 99九九线精品视频在线观看视频| 国产探花极品一区二区| 一个人观看的视频www高清免费观看| 伦理电影大哥的女人| 69人妻影院| 岛国毛片在线播放| 97热精品久久久久久| 九九在线视频观看精品| 亚洲四区av| 身体一侧抽搐| 国产精品美女特级片免费视频播放器| 日本五十路高清| 美女高潮的动态| 97人妻精品一区二区三区麻豆| 欧美成人精品欧美一级黄| 少妇丰满av| 蜜桃久久精品国产亚洲av| 久久这里只有精品中国| 国产在线精品亚洲第一网站| 亚洲国产欧美人成| 色5月婷婷丁香| 搡老妇女老女人老熟妇| 国产精品.久久久| 国产精品久久视频播放| 夜夜看夜夜爽夜夜摸| 一本一本综合久久| 国内少妇人妻偷人精品xxx网站| 热99re8久久精品国产| 久久99精品国语久久久| 精品日产1卡2卡| 欧美bdsm另类| 最近最新中文字幕大全电影3| 国产精品一区二区三区四区久久| 欧美变态另类bdsm刘玥| 一个人免费在线观看电影| 久久精品国产亚洲av涩爱 | 亚洲电影在线观看av| 亚洲欧美日韩高清专用| 午夜福利在线在线| 人妻久久中文字幕网| 草草在线视频免费看| 久久久久久久久久久丰满| 丰满人妻一区二区三区视频av| 国产成人a区在线观看| 久久这里只有精品中国| 婷婷色av中文字幕| av在线蜜桃| 夜夜爽天天搞| 亚洲在久久综合| 久久精品影院6| 精品久久久久久久久久久久久| 亚洲美女视频黄频| 国产单亲对白刺激| av在线蜜桃| 毛片一级片免费看久久久久| 精品久久久噜噜| 99热网站在线观看| 免费看av在线观看网站| 亚洲成人久久性| 亚洲国产精品成人综合色| 小蜜桃在线观看免费完整版高清| 日韩欧美三级三区| 中文字幕免费在线视频6| 亚洲精品久久久久久婷婷小说 | 久久久久久大精品| 成人综合一区亚洲| 免费av毛片视频| 特大巨黑吊av在线直播| 亚洲精品亚洲一区二区| 给我免费播放毛片高清在线观看| 白带黄色成豆腐渣| 淫秽高清视频在线观看| 丝袜喷水一区| 最新中文字幕久久久久| 欧美丝袜亚洲另类| av在线观看视频网站免费| 性欧美人与动物交配| 精品国产三级普通话版| 国产色爽女视频免费观看| 成年女人永久免费观看视频| 免费看光身美女| 国模一区二区三区四区视频| 日本免费一区二区三区高清不卡| 成人鲁丝片一二三区免费| 搡女人真爽免费视频火全软件| 国产在线精品亚洲第一网站| 久久久精品欧美日韩精品| 九草在线视频观看| 国内精品宾馆在线| 人人妻人人澡欧美一区二区| 日本爱情动作片www.在线观看| 大香蕉久久网| 一本一本综合久久| 中文字幕av成人在线电影| 又爽又黄无遮挡网站| 久久99热这里只有精品18| 亚洲三级黄色毛片| 亚洲最大成人av| 日本熟妇午夜| 久久国内精品自在自线图片| 国产亚洲5aaaaa淫片| 成年av动漫网址| 亚洲18禁久久av| 国产午夜精品一二区理论片| 日韩av在线大香蕉| 成人毛片60女人毛片免费| 桃色一区二区三区在线观看| 久久精品夜夜夜夜夜久久蜜豆| 久久人妻av系列| 综合色av麻豆| 色综合色国产| av女优亚洲男人天堂| 日韩欧美精品v在线| 国产在视频线在精品| 哪里可以看免费的av片| 天天一区二区日本电影三级| 给我免费播放毛片高清在线观看| 国产三级在线视频| 久久久久久久久久久丰满| 床上黄色一级片| 伊人久久精品亚洲午夜| 亚洲四区av| 久久99精品国语久久久| 精品人妻一区二区三区麻豆| 淫秽高清视频在线观看| 国产精品1区2区在线观看.| 国产精品久久久久久精品电影小说 | 中文欧美无线码| 又爽又黄a免费视频| 久久精品综合一区二区三区| 看十八女毛片水多多多| 国产真实乱freesex| av国产免费在线观看| 日本一本二区三区精品| 我要搜黄色片| 久久精品91蜜桃| 国产亚洲欧美98| 亚洲av.av天堂| 18+在线观看网站| 六月丁香七月| 日韩av不卡免费在线播放| 午夜精品国产一区二区电影 | 国产精品麻豆人妻色哟哟久久 | 欧美一区二区国产精品久久精品| 成年版毛片免费区| 亚洲成人久久爱视频| 色噜噜av男人的天堂激情| 啦啦啦观看免费观看视频高清| 日本免费a在线| 亚洲电影在线观看av| 啦啦啦啦在线视频资源| 国产精品国产三级国产av玫瑰| 91精品一卡2卡3卡4卡| 欧美成人一区二区免费高清观看| 神马国产精品三级电影在线观看| 国产精品永久免费网站| 九九热线精品视视频播放| 伦理电影大哥的女人| 成人性生交大片免费视频hd| 国产不卡一卡二| 免费人成在线观看视频色| 亚洲av熟女| 麻豆久久精品国产亚洲av| 又爽又黄无遮挡网站| 亚州av有码| 精品少妇黑人巨大在线播放 | 国产欧美日韩精品一区二区| 少妇猛男粗大的猛烈进出视频 | 亚洲图色成人| 亚洲无线在线观看| 夫妻性生交免费视频一级片| 久久人人爽人人片av| 久久久久久久久大av| 日日摸夜夜添夜夜添av毛片| 美女cb高潮喷水在线观看| 国产真实伦视频高清在线观看| 丝袜喷水一区| 日本爱情动作片www.在线观看| 三级国产精品欧美在线观看| 午夜精品国产一区二区电影 | 欧美最黄视频在线播放免费| 欧美在线一区亚洲| 亚洲精品乱码久久久久久按摩| 直男gayav资源| 18禁黄网站禁片免费观看直播| 大又大粗又爽又黄少妇毛片口| 欧美zozozo另类| 欧美精品一区二区大全| 国产精品久久电影中文字幕| 淫秽高清视频在线观看| 插逼视频在线观看| 免费看日本二区| 少妇的逼好多水| 三级经典国产精品| 99久国产av精品国产电影| 天天躁夜夜躁狠狠久久av| 夜夜夜夜夜久久久久| 在线免费观看的www视频| 国产伦在线观看视频一区| 国产高清有码在线观看视频| 在线观看免费视频日本深夜| 在线免费观看不下载黄p国产| 国产高清激情床上av| 在线播放无遮挡| 午夜免费激情av| av免费观看日本| 联通29元200g的流量卡| 国产成人a∨麻豆精品| 99国产极品粉嫩在线观看| 久久久a久久爽久久v久久| 狠狠狠狠99中文字幕| 国产乱人视频| 最近最新中文字幕大全电影3| 一夜夜www| 国产乱人视频| 国产亚洲av片在线观看秒播厂 | 国产不卡一卡二| 久久久久久久久中文| 性色avwww在线观看| 亚洲成人中文字幕在线播放| 看片在线看免费视频| 久久午夜福利片| 自拍偷自拍亚洲精品老妇| 亚洲在线自拍视频| 一个人观看的视频www高清免费观看| 99国产极品粉嫩在线观看| 天堂网av新在线| 天天躁日日操中文字幕| 日本与韩国留学比较| 99热精品在线国产| 欧美xxxx黑人xx丫x性爽| 日韩欧美精品v在线| 99久久精品一区二区三区| 高清毛片免费观看视频网站| 日韩中字成人| 日韩在线高清观看一区二区三区| 一本一本综合久久| 草草在线视频免费看| 欧美人与善性xxx| 国产精品美女特级片免费视频播放器| 我的女老师完整版在线观看| 国产探花极品一区二区| 三级毛片av免费| 国产精品一区二区三区四区久久|