• <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在线直播| 麻豆乱淫一区二区| 免费黄色在线免费观看| 欧美3d第一页| 久久人人爽人人片av| 久久精品夜色国产| 水蜜桃什么品种好| 国产精品久久久久成人av| 色婷婷久久久亚洲欧美| 男人和女人高潮做爰伦理| 国产美女午夜福利| 狂野欧美激情性bbbbbb| 狂野欧美激情性bbbbbb| 男人狂女人下面高潮的视频| 亚洲精品中文字幕在线视频 | 免费在线观看成人毛片| 午夜福利视频精品| 午夜激情久久久久久久| 国产伦精品一区二区三区视频9| 亚洲av男天堂| 国产精品人妻久久久影院| 有码 亚洲区| 亚洲精品国产成人久久av| 丰满乱子伦码专区| 亚洲欧洲精品一区二区精品久久久 | 亚洲内射少妇av| 一本—道久久a久久精品蜜桃钙片| 日韩三级伦理在线观看| 国产精品人妻久久久影院| 国产亚洲最大av| 视频中文字幕在线观看| 国产免费视频播放在线视频| 91在线精品国自产拍蜜月| 美女福利国产在线| 最近手机中文字幕大全| 乱系列少妇在线播放| 少妇 在线观看| 成人二区视频| 国产91av在线免费观看| 91成人精品电影| 国产在线男女| 国产精品99久久久久久久久| 色94色欧美一区二区| 熟女电影av网| 午夜福利网站1000一区二区三区| 精品人妻偷拍中文字幕| 日本与韩国留学比较| 街头女战士在线观看网站| 免费黄频网站在线观看国产| 简卡轻食公司| 欧美日韩av久久| 丰满饥渴人妻一区二区三| 国产欧美日韩一区二区三区在线 | 亚洲精品国产av蜜桃| 亚洲精品乱码久久久v下载方式| 精品亚洲成国产av| 男人添女人高潮全过程视频| 亚洲国产精品999| h日本视频在线播放| 欧美日韩一区二区视频在线观看视频在线| 精品熟女少妇av免费看| 99精国产麻豆久久婷婷| 欧美xxⅹ黑人| 日韩av不卡免费在线播放| av播播在线观看一区| 国产成人91sexporn| 在线观看三级黄色| 男女国产视频网站| 成人18禁高潮啪啪吃奶动态图 | av又黄又爽大尺度在线免费看| 日本vs欧美在线观看视频 | 欧美日韩亚洲高清精品| 国产精品麻豆人妻色哟哟久久| 另类亚洲欧美激情| 日韩熟女老妇一区二区性免费视频| 国产精品一区二区在线观看99| 我的老师免费观看完整版| 亚洲美女黄色视频免费看| 亚洲美女视频黄频| 免费久久久久久久精品成人欧美视频 | 久久精品久久精品一区二区三区| 综合色丁香网| 亚洲婷婷狠狠爱综合网| 狂野欧美激情性bbbbbb| 极品少妇高潮喷水抽搐| h视频一区二区三区| 亚洲无线观看免费| 黑丝袜美女国产一区| 亚洲av中文av极速乱| 国产 一区精品| 日韩中文字幕视频在线看片| av在线播放精品| av女优亚洲男人天堂| 免费看av在线观看网站| 日韩三级伦理在线观看| 国产高清国产精品国产三级| 国精品久久久久久国模美| 午夜老司机福利剧场| 中文字幕精品免费在线观看视频 | 亚洲色图综合在线观看| 一级毛片我不卡| 啦啦啦啦在线视频资源| 亚洲第一av免费看| 成年人免费黄色播放视频 | 黑人猛操日本美女一级片| 午夜老司机福利剧场| 国产一级毛片在线| 大又大粗又爽又黄少妇毛片口| 大香蕉久久网| 在线观看免费日韩欧美大片 | 丰满迷人的少妇在线观看| 国模一区二区三区四区视频| 最近2019中文字幕mv第一页| 亚洲图色成人| 秋霞伦理黄片| 日日摸夜夜添夜夜添av毛片| 欧美日韩综合久久久久久| 亚洲欧美一区二区三区黑人 | h视频一区二区三区| 免费少妇av软件| 亚洲人成网站在线观看播放| 一级毛片aaaaaa免费看小| 成年女人在线观看亚洲视频| 亚洲三级黄色毛片| 嫩草影院入口| 久久精品久久精品一区二区三区| a级毛片免费高清观看在线播放| 久久久久久久精品精品| 99九九线精品视频在线观看视频| 激情五月婷婷亚洲| 亚洲av福利一区| 一级二级三级毛片免费看| 国产综合精华液| 国产精品免费大片| 夫妻午夜视频| 少妇的逼好多水| 国产高清有码在线观看视频| 色哟哟·www| 精品国产一区二区久久| av免费在线看不卡| 国产淫语在线视频| 纯流量卡能插随身wifi吗| 欧美 日韩 精品 国产| 亚洲精品乱码久久久v下载方式| 成年女人在线观看亚洲视频| 久久久精品免费免费高清| 国产伦精品一区二区三区四那| 国产色婷婷99| av黄色大香蕉| 久久亚洲国产成人精品v| 色5月婷婷丁香| 不卡视频在线观看欧美| .国产精品久久| 亚洲精品视频女| 久久久久久久久久久久大奶| 免费观看av网站的网址| av.在线天堂| 欧美日韩一区二区视频在线观看视频在线| 国产亚洲5aaaaa淫片| 日本午夜av视频| 国产精品伦人一区二区| 特大巨黑吊av在线直播| 丝瓜视频免费看黄片| 99久久人妻综合| 亚洲欧美日韩另类电影网站| 丁香六月天网| 少妇精品久久久久久久| 丰满饥渴人妻一区二区三| 国产一级毛片在线| 国产91av在线免费观看| 亚洲精品一二三| 人妻一区二区av| 免费看光身美女| av有码第一页| 久久久亚洲精品成人影院| 午夜福利网站1000一区二区三区| 蜜臀久久99精品久久宅男| 免费人妻精品一区二区三区视频| 亚洲自偷自拍三级| 五月玫瑰六月丁香| 一级片'在线观看视频| 成人特级av手机在线观看| 日本猛色少妇xxxxx猛交久久| 一区二区三区四区激情视频| 六月丁香七月| 高清欧美精品videossex| 精品国产国语对白av| 欧美日韩国产mv在线观看视频| 一区二区三区免费毛片| 丰满饥渴人妻一区二区三| 国产精品伦人一区二区| 蜜桃久久精品国产亚洲av| 中文字幕免费在线视频6| av视频免费观看在线观看| av视频免费观看在线观看| 久热这里只有精品99| 国产色婷婷99| 熟女av电影| 午夜福利,免费看| 永久网站在线| 亚洲精品亚洲一区二区| 中文欧美无线码| 美女主播在线视频| 欧美少妇被猛烈插入视频| 有码 亚洲区| 偷拍熟女少妇极品色| 一级毛片久久久久久久久女| 色视频在线一区二区三区| 在线 av 中文字幕| 欧美日韩av久久| 国产无遮挡羞羞视频在线观看| 又爽又黄a免费视频| www.av在线官网国产| 久久久国产精品麻豆| 亚洲经典国产精华液单| 国产在线免费精品| 伊人久久精品亚洲午夜| 女性生殖器流出的白浆| 九九在线视频观看精品| 国产精品一区二区在线观看99| av.在线天堂| 99久久精品热视频| 精品一区二区免费观看| 多毛熟女@视频| 自拍偷自拍亚洲精品老妇| 日本午夜av视频| 9色porny在线观看| 国产精品一区二区在线不卡| 国产精品久久久久久av不卡| 国产淫片久久久久久久久| 在现免费观看毛片| 下体分泌物呈黄色| 免费人成在线观看视频色| 久久国产精品男人的天堂亚洲 | 青春草视频在线免费观看| 在线观看人妻少妇| 九草在线视频观看| 亚洲精品,欧美精品| 在线观看免费日韩欧美大片 | av免费观看日本| 免费看光身美女| 色吧在线观看| 97超视频在线观看视频| 多毛熟女@视频| 久久6这里有精品| 美女福利国产在线| 免费黄网站久久成人精品| 午夜精品国产一区二区电影| 日韩电影二区| 美女中出高潮动态图| 久久久国产精品麻豆| 在线天堂最新版资源| av福利片在线| 亚洲欧美成人综合另类久久久| 亚洲欧洲日产国产| 久久 成人 亚洲| 久久韩国三级中文字幕| 如何舔出高潮| 久久久国产欧美日韩av| 亚洲欧美成人精品一区二区| 美女视频免费永久观看网站| 夜夜爽夜夜爽视频| 免费人成在线观看视频色| 亚洲中文av在线| 又大又黄又爽视频免费| 免费av中文字幕在线| 一二三四中文在线观看免费高清| 国产精品熟女久久久久浪| 国产成人精品一,二区| kizo精华| 亚洲成色77777| 22中文网久久字幕| 2018国产大陆天天弄谢| 午夜免费鲁丝| 国产精品熟女久久久久浪| 午夜激情福利司机影院| 99热6这里只有精品| 人妻 亚洲 视频| 午夜视频国产福利| 久久国产精品大桥未久av | 久久婷婷青草| 成人毛片60女人毛片免费| 精品国产乱码久久久久久小说| 一边亲一边摸免费视频| 久久人人爽av亚洲精品天堂| 狂野欧美激情性bbbbbb| 中文字幕免费在线视频6| 午夜免费观看性视频| 国产伦在线观看视频一区| 亚洲精品日韩在线中文字幕| 伦理电影大哥的女人| 日韩,欧美,国产一区二区三区| 国产成人精品无人区| av播播在线观看一区| www.色视频.com| 97超视频在线观看视频| 夜夜骑夜夜射夜夜干| 热re99久久精品国产66热6| 国产成人精品一,二区| 欧美老熟妇乱子伦牲交| 2022亚洲国产成人精品| 日日撸夜夜添| 日日撸夜夜添| 观看美女的网站| 下体分泌物呈黄色| 中文字幕制服av| 日本欧美视频一区| 国产白丝娇喘喷水9色精品| 久久热精品热| 伦理电影免费视频| 国内精品宾馆在线| 国产精品国产av在线观看| 国产精品一区二区在线观看99| 日韩一本色道免费dvd| 汤姆久久久久久久影院中文字幕| 18禁在线播放成人免费| av国产久精品久网站免费入址| 欧美xxxx性猛交bbbb| 成人美女网站在线观看视频| 国产永久视频网站| 99热国产这里只有精品6| a级毛片免费高清观看在线播放| 黄色怎么调成土黄色| 欧美xxxx性猛交bbbb| 伊人久久精品亚洲午夜| 永久免费av网站大全| 成人免费观看视频高清| 久久国产乱子免费精品| 丰满少妇做爰视频| 日本欧美视频一区| 亚洲经典国产精华液单| 亚洲av.av天堂| 久久青草综合色| 免费黄色在线免费观看| 成年人午夜在线观看视频| 亚洲四区av| 亚洲情色 制服丝袜| 色婷婷久久久亚洲欧美| 久久国产亚洲av麻豆专区| 热re99久久国产66热| 免费看av在线观看网站| 国产av一区二区精品久久| 欧美三级亚洲精品| 久久精品国产自在天天线| 在线观看免费高清a一片| 麻豆乱淫一区二区| 日韩欧美一区视频在线观看 | 一本色道久久久久久精品综合| 国产一区亚洲一区在线观看| 美女福利国产在线| 人妻系列 视频| 亚洲第一av免费看| 午夜日本视频在线| 午夜激情久久久久久久| 亚洲成人手机| 国产精品蜜桃在线观看| 少妇高潮的动态图| 日本av免费视频播放| 免费久久久久久久精品成人欧美视频 | 久久久国产欧美日韩av| 看十八女毛片水多多多| 黄片无遮挡物在线观看| 成人无遮挡网站| 狂野欧美激情性bbbbbb| 久久狼人影院| 久久久久久久久久久免费av| a级毛片在线看网站| 国产高清有码在线观看视频| 国产乱来视频区| tube8黄色片| 五月开心婷婷网| 亚洲精品久久久久久婷婷小说| 91久久精品国产一区二区三区| 黑人猛操日本美女一级片| 亚洲国产最新在线播放| 亚洲欧美日韩另类电影网站| 欧美一级a爱片免费观看看| 日韩欧美 国产精品| 亚洲国产av新网站| 一级二级三级毛片免费看| 丝袜脚勾引网站| 新久久久久国产一级毛片| 久久久久久久久久久久大奶| 日日撸夜夜添| 男男h啪啪无遮挡| 亚洲精品国产av蜜桃| 亚洲精华国产精华液的使用体验| 欧美日韩视频精品一区| 亚洲精品国产色婷婷电影| 日本欧美国产在线视频| 久久久国产欧美日韩av| 亚洲av欧美aⅴ国产| 色吧在线观看| 亚洲国产最新在线播放| 熟妇人妻不卡中文字幕| 午夜免费男女啪啪视频观看| 精品熟女少妇av免费看| 91久久精品国产一区二区三区| av不卡在线播放| 伊人久久国产一区二区| 欧美激情国产日韩精品一区| 精品久久久久久电影网| 在线观看免费高清a一片| 国产无遮挡羞羞视频在线观看| 日韩在线高清观看一区二区三区| 日韩制服骚丝袜av| 毛片一级片免费看久久久久| 国产深夜福利视频在线观看| 三级经典国产精品| 一级毛片电影观看| 国产淫片久久久久久久久| 欧美97在线视频| 在线看a的网站| 久久综合国产亚洲精品| 国产美女午夜福利| 精品久久久噜噜| 夫妻性生交免费视频一级片| 国产精品免费大片| 久久久久精品性色| 亚洲精品一区蜜桃| 成人免费观看视频高清| 这个男人来自地球电影免费观看 | av在线老鸭窝| 黄色日韩在线| 亚洲美女黄色视频免费看| 国产亚洲精品久久久com| 成人综合一区亚洲| 我的老师免费观看完整版| 卡戴珊不雅视频在线播放| 欧美3d第一页| 一个人看视频在线观看www免费| 亚洲欧美精品专区久久| 一本一本综合久久| 久久精品熟女亚洲av麻豆精品| av天堂中文字幕网| 国产一区二区三区av在线| 在线观看国产h片| 另类精品久久| 男人添女人高潮全过程视频| 在线看a的网站| 免费av中文字幕在线| 夫妻性生交免费视频一级片| 最新的欧美精品一区二区| 欧美国产精品一级二级三级 | 欧美精品高潮呻吟av久久| 黄色配什么色好看| 69精品国产乱码久久久| 久久久久久久精品精品| 久久综合国产亚洲精品| 国产乱来视频区| 黄色日韩在线| 久久99蜜桃精品久久| 国产日韩一区二区三区精品不卡 | 全区人妻精品视频| 久久久久视频综合| 欧美bdsm另类| 亚洲美女黄色视频免费看| 久久久午夜欧美精品| .国产精品久久| 国产一区二区三区综合在线观看 | a级毛色黄片| 亚洲国产精品国产精品| 国产欧美亚洲国产| 久久久精品94久久精品| 色婷婷久久久亚洲欧美| 精品一区二区三卡| 亚洲精品自拍成人| 大香蕉97超碰在线| 国产免费又黄又爽又色| 国产免费一级a男人的天堂| 少妇的逼好多水| 99久国产av精品国产电影| 日韩欧美一区视频在线观看 | 亚州av有码| 中国国产av一级| 成人漫画全彩无遮挡| 国产精品伦人一区二区| 五月开心婷婷网| 亚洲av在线观看美女高潮| 国产日韩欧美在线精品| 午夜福利网站1000一区二区三区| 久久人人爽人人爽人人片va| 中文字幕av电影在线播放| 亚洲av在线观看美女高潮| 亚洲av日韩在线播放| av一本久久久久| 在线免费观看不下载黄p国产| 99re6热这里在线精品视频| 欧美日韩精品成人综合77777| 久久 成人 亚洲| 自线自在国产av| 精品国产一区二区三区久久久樱花| 精品酒店卫生间| 国产黄频视频在线观看| 男女国产视频网站| 性色av一级| 在线精品无人区一区二区三| 女人久久www免费人成看片| 国产欧美日韩综合在线一区二区 | 最黄视频免费看| 国产精品久久久久久av不卡| 国产高清三级在线| 免费大片黄手机在线观看| 日本与韩国留学比较| 欧美日韩视频高清一区二区三区二| 99久久精品国产国产毛片| 国产日韩欧美亚洲二区| 美女大奶头黄色视频| 中国国产av一级| 精品99又大又爽又粗少妇毛片| 成年人免费黄色播放视频 | 成人亚洲精品一区在线观看| 亚洲精品色激情综合| a级毛色黄片| 看十八女毛片水多多多| 久久久久久伊人网av| 高清毛片免费看| 国内揄拍国产精品人妻在线| 妹子高潮喷水视频| 日产精品乱码卡一卡2卡三| 在线 av 中文字幕| 久久99一区二区三区| 国产亚洲精品久久久com| 99视频精品全部免费 在线| 好男人视频免费观看在线| .国产精品久久| 好男人视频免费观看在线| 成人亚洲精品一区在线观看| 免费观看无遮挡的男女| 亚洲精品一区蜜桃| 午夜福利网站1000一区二区三区| 我要看日韩黄色一级片| 久久精品国产亚洲av天美| 精品人妻熟女毛片av久久网站| 一区在线观看完整版| 最后的刺客免费高清国语| 美女中出高潮动态图| 日韩av在线免费看完整版不卡| 噜噜噜噜噜久久久久久91| 99热这里只有是精品在线观看| 亚洲精品一二三| 黄片无遮挡物在线观看| 日韩免费高清中文字幕av| 尾随美女入室| 看非洲黑人一级黄片| 亚洲色图综合在线观看| 国产精品蜜桃在线观看| www.av在线官网国产| 日韩人妻高清精品专区| 精品一区在线观看国产| 久久精品国产鲁丝片午夜精品| 色5月婷婷丁香| 观看av在线不卡| 一级毛片我不卡| 久久狼人影院| 一区二区三区精品91| 啦啦啦视频在线资源免费观看| 亚洲国产成人一精品久久久| 亚洲在久久综合| 草草在线视频免费看| 亚洲色图综合在线观看| 亚洲不卡免费看| 我的老师免费观看完整版| 久久午夜综合久久蜜桃| 水蜜桃什么品种好| 97超碰精品成人国产| 国产伦精品一区二区三区视频9| 51国产日韩欧美| 国产白丝娇喘喷水9色精品| 国产一区亚洲一区在线观看| 国产中年淑女户外野战色| a级一级毛片免费在线观看| 免费播放大片免费观看视频在线观看| 亚洲综合精品二区| 亚洲人成网站在线观看播放| 久久精品国产自在天天线| 国产在视频线精品| 亚洲精品456在线播放app| 人人妻人人澡人人爽人人夜夜| 极品人妻少妇av视频| 人妻少妇偷人精品九色| av专区在线播放| 成人毛片a级毛片在线播放| 精品久久久久久久久亚洲| 精品一区在线观看国产| 久久久久网色| 欧美老熟妇乱子伦牲交| 一级a做视频免费观看| 久久久国产欧美日韩av| 只有这里有精品99| 日韩强制内射视频| 久久国产精品男人的天堂亚洲 | 亚洲国产精品专区欧美| 国精品久久久久久国模美| 国产色婷婷99| 自拍偷自拍亚洲精品老妇| 久久97久久精品| 在线观看免费视频网站a站| 国产精品秋霞免费鲁丝片| 夜夜骑夜夜射夜夜干| 久久鲁丝午夜福利片| 黄色一级大片看看| 一级片'在线观看视频| 51国产日韩欧美| 日韩熟女老妇一区二区性免费视频| 成人特级av手机在线观看| 亚洲欧美一区二区三区国产| 国产精品人妻久久久久久| 欧美精品一区二区免费开放| 久久久精品94久久精品| 中文字幕免费在线视频6| 一级,二级,三级黄色视频| 人人澡人人妻人| 乱码一卡2卡4卡精品|