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

復材非線性及漸進損傷的態型近場動力學模擬

2016-06-17 06:42:18劉肅肅
浙江大學學報(工學版) 2016年5期
關鍵詞:復合材料

劉肅肅,余 音

(上海交通大學 航空航天學院,上海 200240)

?

復材非線性及漸進損傷的態型近場動力學模擬

劉肅肅,余音

(上海交通大學 航空航天學院,上海 200240)

摘要:為了實現復合材料材料非線性行為在態型近場動力學的模擬,引入單參量非線性本構模型來確定態型近場動力學中力狀態與變形狀態之間的非線性關系.為了實現復合材料漸進損傷行為的態型近場動力學模擬,在力狀態表達式中引入一個標量函數,以實現利用強度準則判斷物質點之間作用力是否失效.采用動力松弛法作為數值求解方法.結果顯示,模型計算得到的單向復合材料板(AS4/PEEK)的非線性應力-應變曲線與試驗數據吻合良好;對中心含孔復合材料板的漸進損傷模擬則表明,該模型模擬的破壞形式與已有研究的模擬結果及實驗觀察結果一致.探索和實現了復合材料非線性本構在態型近場動力學中的引入以及復合材料漸進損傷在態型近場動力學中的模擬.

關鍵詞:近場動力學;復合材料;非線性;漸進損傷;動力松弛法

近場動力學(peridynamic theory, PD)[1-2]是一種新興的基于非局部作用思想的理論.不同于傳統理論的是,該理論采用積分形式構建運動方程,在解決諸多不連續問題方面具有獨特的優勢.近年來,PD理論受到了廣泛的關注并得到了迅速的發展.黃丹等[3]對PD的基本方法及其相關應用進行了介紹,并將PD應用于模擬混凝土的沖擊破壞過程[4];胡祎樂等[5]則將PD應用于分析復合材料層壓板的漸進損傷;孫璐妍[6]在研究中建立了金屬平板和曲板穩定性分析的PD模型.

最初的鍵型近場動力學理論(bond-based peridynamics)[1]在采用對點力函數(pairwise force function)來描述物質點之間的相互作用力時,存在一定的不足[3,7].這種對點力函數假設認為兩物質點之間的相互作用力是等大反向的,并且未考慮鄰域內其他因素(如其他物質點及物質點間作用力)對該作用力的影響,在描述各向同性、線彈性材料時,其泊松比只能局限于0.25.此外,在應用該理論描述金屬材料的塑性時,也會出現與實驗現象相矛盾的結果[3,7].為了解決上述問題,Silling等[7-8]提出態型近場動力學理論(state-based peridynamics),該理論不僅可以應用于不同泊松比的材料,同時也可以在PD理論框架下實現一些經典力學模型的應用:Foster等[9]在態型PD理論中利用金屬塑性理論分析了6061-T6鋁合金的粘塑性響應;Tupek等[10]則在態型PD理論中整合了經典力學破壞模型,模擬了6061-T6鋁合金在沖擊載荷下的破壞過程及單軸拉伸載荷下的力學響應.

已有的復合材料PD模型多采用鍵型PD理論進行建模:Askari等[11]在用鍵型PD理論對復合材料進行建模時,分別采用纖維鍵和基體鍵來考慮纖維和基體的不同性能;Killic等[12]將單向復合材料離散成纖維域和基體域,兩者之間的比例可以反映纖維體積含量;Hu等[13]則采用宏觀偏軸模量來等效確定復合材料層內鍵的大小.這些模型給出了與實驗現象相符的破壞模式,對宏觀等效性能的精確定量分析還比較缺乏,相關的定量分析也多在線彈性范圍內,既能準確預測復材材料非線性行為、又能描述其損傷行為的PD模型有待進一步研究.

本文采用態型PD理論對復合材料進行建模,采用經典單參量非線性本構模型確定力狀態與變形狀態之間的關系,并在損傷描述中引入傳統的強度準則.為驗證模型的有效性,編寫了數值計算程序,并對AS4/PEEK復合材料在單軸拉伸載荷作用下的材料非線性行為及漸進損傷過程進行模擬,與文獻中的實驗結果及有限元的計算結果進行對比.

1態型近場動力學理論簡介

近場動力學將物體離散成一系列物質點,在物體參考構型R內,物質點x存在半徑為δ的鄰域(Horizon)H,在該鄰域內物質點x與其他物質點x′存在相互作用.記物質點x與物質點x′的初始坐標分別為x和x′,位移分別為u和u′,則兩物質點的相對位置和相對位移可分別表示為ξ=x′-x,η=u′-u,如圖1所示.

圖1 參考構型與鄰域Fig.1 Reference configuration and horizon

在近場動力學中,物質點x上建立運動方程為

(1)

(2)

(3)

(4)

Silling等[7]通過與經典理論中應變能密度等效的方式得到了變形梯度的表達式:

(5)

(6)

2單參量模型

2.1單參量非線性本構

在表征復合材料物理非線性的宏觀力學理論模型中,由Sun等[14]提出的單參量模型能夠較為準確地描述單向復合材料在偏軸載荷作用下的非線性響應.在平面應力狀態下,該模型在假設纖維方向為線彈性的基礎上,獲得了僅含單參量的塑性勢能函數表達式:

(7)

(8)

式中:A和n為擬合參數.將總的應變增量分為彈性部分增量dεe和塑性部分增量dεp,dε=dεe+dεp.結合流動法則,塑性應變增量為

(9)

(10)

該塑性功增量又可表示為等效應力與等效塑性應變增量的乘積形式:

(11)

(12)

塑性應變增量的表達式可聯立式(7)、(8)、(9)和(12)得到

(13)

dε=Sdσ,

式中:S為柔度陣.根據式(13)結合復合材料彈性理論可得到柔度陣為

(14)

式中:E1、E2、E12為彈性模量,ν12為泊松比.

定義塑性泊松比為

(15)

(16)

2.2強度準則

Sun等[15]在對AS4/PEEK單向復合材料在各偏軸拉伸載荷作用下的強度進行測試后,認為其強度準則為

(17)

(18)

式中:Y和S分別為橫向拉伸強度及剪切強度.

3含非線性本構及損傷描述的態型PD模型

3.1本構方程

在態型近場動力學理論中,力狀態與應力張量σ之間存在以下關系[7]:

(19)

將式(19)表示為增量形式,應力增量可用應變增量表示,相應的力狀態增量為

(20)

位移梯度可以用變形梯度式(5)表示

▽u=F-I.

(21)

式中:I為單位陣.則應變張量為

(22)

將式(5)代入式(22),再代入式(20),得到

(23)

3.2損傷描述

一些復合材料PD模型[11-13]在采用鍵理論分析損傷問題時,一般采用鍵的臨界伸長率作為鍵是否失效的判斷依據.然而材料的損傷特性也有可能取決于其他狀態參量[10,15],如應力狀態等.本文嘗試了一種利用強度準則判斷損傷是否發生的方法,為損傷在態型PD理論中的描述提供參考.

(24)

式中:D和D′分別為鍵ξ兩端物質點x與物質點x′的損傷參量.首先判斷鍵ξ是纖維鍵還是基體鍵(本文中鍵的方向與纖維鋪層方向一致則為纖維鍵,否則為基體鍵):若為纖維鍵,則采用式(17)中的“纖維斷裂”判據判斷鍵兩端物質點處的應力狀態,達到判據的物質點對應的D(或D′)值取為1,否則取值為0;若為基體鍵,則采用式(17)中的“基體破壞”判據對鍵兩端物質點的D(D′)值進行判斷.

將式(24)引入式(20)得到

(25)

如此,即可根據式(24)、 (25),依據式(17)的強度準則對兩物質點之間作用力是否失效作出判斷.式(25)的力狀態表達式中包含了損傷和斷裂的描述,可用于描述裂紋萌生和擴展過程.

進一步,物質點x處損傷程度可用損傷指數φ進行表征:

(26)

4數值方法

4.1離散化方法

對模型進行均勻離散化,并采用黎曼和對積分項進行求解.根據式(5)和 (6),得到物質點xj處的變形梯度和形狀張量分別為

(27)

(28)

式中:xn和xj分別為物質點n和物質點j的位置矢量,m為物質點j鄰域H范圍內的物質點總數,Vn為物質點n的體積.將式(3)代入式(27)、再代入式(22)得到物質點xj處的非局部應變張量為

(29)

式中:un和uj分別為物質點n和物質點j的位移矢量.式(29)即為物質點j處的非局部應變,由該物質點鄰域H范圍內的位移場計算得到.

4.2數值求解

動力松弛法是一種應用廣泛的數值求解方法.在應用于靜力問題時,動力松弛法可采用虛擬質量和虛擬阻尼,進而將靜力問題轉化成動力問題來進行求解.在施加初始增量步之后,系統首先在不平衡力的作用下發生運動,而后在迭代過程中不平衡力逐漸減小,系統動能也在阻尼作用下漸趨于零,最后達到靜力平衡狀態.依據達朗貝爾原理,第n次迭代時物質點的運動方程為

(30)

(31)

整理得

(32)

(33)

在求解靜力問題時,虛擬密度、虛擬阻尼與計算結果無關,但會影響計算穩定性和收斂速度[16].一般取Δt=1,將虛擬密度陣設為如下對角形式:

(34)

則對角元素與Δt需滿足以下關系式[17]:

(35)

式中:Q為平衡方程中的剛度陣.第n次迭代的虛擬阻尼可按以下表達式確定[17]:

(36)

i=1,2,3.

(37)

4.3收斂準則

采用相對位移收斂準則判斷物質點是否達到平衡狀態:

(38)

5算例與結果

在第3、4章基礎上編寫了復合材料近場動力學模型的數值計算程序.在對模型進行均勻離散化之后,先給定離散點的位移邊界條件及載荷條件,而后按照動力松弛法進行迭代求解.

為驗證模型及程序的有效性,分別對AS4/PEEK復合材料中心無孔層壓板和中心含孔層壓板受單軸拉伸載荷作用的情況進行建模計算,并與試驗及有限元結果進行對比.試驗數據及材料參數來自文獻[15,18],材料參數如表1所示.其中參數A1和A2分別為式(8)中算例5.1和算例5.2所取的A值.

表1 室溫下AS4/PEEK復合材料的材料參數

5.1中心無孔層壓板受偏軸拉伸

圖2 中心無孔層壓板受偏軸拉伸載荷示意圖Fig.2 Laminate without a hole under off-axis tension load

中心無孔層壓板受偏軸拉伸載荷作用的示意圖如圖2所示,單向復合材料板的纖維方向與載荷作用方向的夾角為θ,取層壓板有效尺寸為190.5 mm×19.0 mm.在離散化建模時,取物質點間距Δx=0.5 mm,鄰域半徑按照文獻[19]的建議取δ=3Δx,模型最終包含物質點總數為14 478個.邊界條件為一邊固支,另一邊施加縱向固定位移,邊界條件施加在兩端的物質點上.

該模型的離散物質點數雖然不多,但是考慮到非線性迭代計算需要保證一定的增量步數,故而通過將各物質點的計算映射到計算機的單個線程的方式來實施并行運算,提高計算效率.

最終得到的偏軸拉伸應力-應變曲線與試驗數據的對比結果如圖3所示.圖中εL為根據式(29)計算得到的層壓板中心處物質點在拉伸方向的非局部應變,σL為相應物質點在拉伸方向上的應力.由圖3可見,本模型的計算結果與試驗數據吻合良好,誤差小于9.2%.

圖3 不同偏軸拉伸角度下應力-應變對比結果Fig.3 Stress-strain comparison under different off-axis tension angles

如圖4所示給出了不同偏軸角度下層壓板中心處物質點的塑性泊松比的對比結果.可見,各偏軸角度下由本模型計算所得的塑性泊松比相對式(16)的理論值更接近于試驗值.

圖4 不同偏軸角度下塑性泊松比的對比結果Fig.4 Plastic Poisson’s ratios comparison under different off-axis tension angles

5.2中心含孔層壓板受軸向拉伸

中心含孔層壓板受軸向拉伸載荷作用的示意圖如圖5所示,復合材料層壓板的纖維方向與載荷作用方向一致,取層壓板有效尺寸177.5 mm×45.5 mm.在離散化建模時,本算例仍然取物質點間距Δx=0.5 mm,鄰域半徑δ=3Δx,最終模型包含物質點總數為30 853個.中心圓孔半徑為20.5 mm,圓孔周邊p1處沿環向貼有應變片,p1與圓孔中心的位置關系見圖6.按模型離散點位置,本模型實際所取角φ=32.27°.根據本模型得到的遠端應力σR-孔邊應變εH曲線與試驗數據的對比結果如圖6所示,兩者在應變較小處吻合良好,但在應變較大處出現了差異.這種差異可能是由于單參量塑性模型的參數是在比例加載情況下得到的,只能近似適用于孔邊的受力狀態[18].

圖5 中心含孔層壓板受軸向拉伸載荷示意圖Fig.5 Laminate with a centre-hole under axis tension load

圖6 中心含孔層壓板遠端應力—孔邊應變曲線對比Fig.6 Remote stress-hoop strain curve comparison for laminate with a centre-hole

如圖7和8所示分別給出了端部產生總長度1%的位移時,縱向(載荷作用方向)與橫向的位移場分布云圖.損傷指數φ采用式(26)進行表征.(a)和(b)分別為有限元(FEM)的計算結果和本模型的計算結果.由圖7、8可見兩者的縱向位移場分布幾乎沒有差別.在橫向的位移場分布上,由有限元計算所得的正向與反向的最大位移分別為0.320和-0.324 mm,PD計算得到的則對應為0.387和-0.391 mm,但在[-0.324,0.320]范圍外的位移僅分布在靠近圓孔周邊的少量物質點上.總體而言, PD計算所得的位移場分布與有限元的計算結果相一致.

圖7 縱向位移場分布云圖Fig.7 Displacement filed distribution in longitudinal direction

圖8 橫向位移場分布云圖Fig.8 Displacement filed distribution in transverse direction

5.3中心含孔層壓板的漸進損傷模擬

由于近場動力學的優勢主要體現在分析不連續問題,本文采用第3節所述的損傷描述方法對中心含孔層壓板受拉伸載荷作用下的漸進損傷進行了模擬.損傷程度采用式(26)進行表征,沒有損傷時損傷指數為0,完全損傷后為1.選用5.2節所述的層壓板尺寸及離散化參數,纖維方向與拉伸載荷作用方向分別取0°、45°、90°.如圖9所示給出了本模型各角度鋪層層壓板的漸進損傷模擬結果,可見在單軸拉伸載荷作用下各角度鋪層層壓板的損傷均起始于圓孔周圍,而后損傷沿纖維方向發生擴展.對比已有研究[20]中對纖維增強樹脂基復合材料的模擬結果以及實驗觀察結果(見圖10的(a)、(b)),本模型模擬的各鋪層角度下的破壞形式與文獻[20]基本相符合.

圖9 各角度下的漸進損傷模擬結果Fig.9 Progress damage results under different angles

圖10 已有研究中的有限元模擬及實驗結果[20]Fig.10 FEM and experimental results of a relevant research[20]

6結論

(1)本研究采用了態型近場動力學對復合材料進行建模,通過引入單參量非線性本構來確定力狀態與變形狀態之間的關系.采用該模型計算了中心無孔AS4/PEEK單向復合材料板在偏軸拉伸載荷下的力學響應.結果顯示,由本模型計算得到各偏軸角度下的應力-應變曲線與試驗數據吻合良好.

(2)對中心含孔層壓板受軸向拉伸情況的計算表明,本模型得到的遠端應力-孔邊應變曲線在應變較小處吻合良好.計算得到的位移場分布與有限元的計算結果相一致.

(3)在力狀態表達式中引入一個標量函數,以實現利用強度準則判斷物質點之間作用力是否失效,從而模擬復合材料的漸進損傷過程.對中心含孔層壓板的漸進損傷模擬表明,損傷最初起始于孔邊,而后沿纖維方向發生擴展.由本模型模擬的各角度下的破壞形式與有關研究中對纖維增強樹脂基復合材料的模擬結果以及實驗觀察結果相一致.

(4)通過以上算例,模型的有效性得以驗證,實現了復合材料非線性及漸進損傷在態型近場動力學中的模擬,為復合材料在態型近場動力學中的建模提供新的思路與方法.

參考文獻(References):[1] SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces [J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175-209.

[2] EMMRICH E, LEHOUCQ R B, PUHST D. Meshfree Methods for Partial Differential Equations VI[M]. Berlin Heidelberg:Springer, 2013: 45-65.

[3] 黃丹, 章青, 喬丕忠,等. 近場動力學方法及其應用[J]. 力學進展, 2010, 40 (4):448-459.

HUANG Dan,ZHANG Qing,QIAO Pei-zhong,et al. A review on peridynamics(PD) method and its application [J].Advances in Mechanics, 2010, 40 (4):448-459.

[4] 沈峰, 章青, 黃丹,等. 沖擊荷載作用下混凝土結構破壞過程的近場動力學模擬[J]. 工程力學, 2012, 29(A01):12-15.

SHEN Feng,ZHANG Qing,HUANG Dan,et al. Peridynamic Modeling of failure process of concrete structures subjected to impacting load [J].Engineering Mechanics, 2012, 29(A01):12-15.

[5] 胡祎樂, 余音, 汪海. 基于近場動力學理論的層壓板損傷分析方法[J]. 力學學報, 2013, 45(4):624-628.

HU Yi-le,YU Yin,WANG Hai. Damage analysis method for laminates based on peridynamic theory [J].Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(4):624-628.

[6] 孫璐妍. 基于近場動力學的金屬板穩定性分析方法研究[D]. 上海: 上海交通大學, 2014.

SUN Lu-yan. Stability analysis method study of metallic plates based on peridynamics[D]. Shanghai: Shanghai Jiao Tong University, 2014.

[7] SILLING S A, EPTON M, WECKNER O, et al. Peridynamic states and constitutive modeling [J]. Journal of Elasticity, 2007, 88(2):151-184.

[8] SILLING S A. Linearized theory of peridynamic states[J]. Journal of Elasticity, 2010, 99(1):85-111.

[9] FOSTER J T, SILLING S A, CHEN W W. Viscoplasticity using peridynamics[J]. International Journal for Numerical Methods in Engineering, 2010, 81(10):1242-1258.

[10] TUPEKA M R, RIMOLIA J J ,RADOVITZKY R. An approach for incorporating classical continuum damage models in state-based peridynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2013, 263:20-26.

[11] ASKARI E, XU J F. Peridynamic analysis of damage and failure in composites [C]∥44th AIAA Aerospace Sciences Meeting and Exhibition. Reno, Nevada,Reston, VA:AIAA, 2006: 1-12.

[12] KILLIC B, AGWAI A, MADENCI E. Peridynamic theory for progressive damage prediction in center-cracked composite laminates [J]. Composite Structures, 2009, 90(2):141-151.

[13] HU Y L, YU Y, WANG H. Peridynamic analytical method for progressive damage in notched composite laminates [J]. Composite Structures, 2014, 108: 801-810.

[14] SUN C T, CHEN J L. A simple flow rule for characterizing nonlinear behavior of fiber composites[J]. Journal of Composite Materials, 1989, 23(10): 1009-1020.

[15] SUN C T, YOON K J. Characterization of elastic-plastic properties of AS4/APC-2 thermoplastic composite[R]. NASA-CR-183145, West Lafayette: School of Aeronautics and Astronautics, Purdue University, 1988.

[16] 單建, 藍天. 動力松弛法在張拉結構靜力分析中的應用[J]. 東南大學學報, 1994, 24(3):94-98.

SHAN Jian,LAN Tian. On dynamic relaxation and its application to static analysis of tension structures [J].Journal of Southeast University, 1994, 24(3):94-98.

[17] KILLIC B. Peridynamic theory for progressive failure prediction in homogeneous and heterogeneous materials[D]. Tucson: The University of Arizona, 2008.

[18] SUN C T, YOON K J. Elastic-plastic analysis of AS4/PEEK composite laminate using a one-parameter plasticity model [J]. Journal of Composite Materials, 1992, 26(2):293-308.

[19] SILLING S A, ASKARI E. A meshfree method based on the peridynamic model of solid mechanics [J]. Computers and Structures, 2005, 83(17):1526-1535.

[20] LAS V, ZEMCIK R. Progressive damage of unidirectional composite panels [J]. Journal of Composite Materials, 2008, 42(1):25-44.

State-based peridynamic modeling of nonlinear behavior and progressive damage of composites

LIU Su-su, YU Yin

(SchoolofAeronauticsandAstronautics,ShanghaiJiaoTongUniversity,Shanghai200240,China)

Abstract:One-parameter nonlinear constitutive model was utilized to get the relationship between deformation vector state and force vector state for the sake of nonlinear behavior modeling of composites in state-based peridynamics. A scalar function was used so that strength criterion could be utilized for the validity judgment of the interaction between particles in order to model the progressive damage of composites in state-based peridynamics. Dynamic relaxation method was utilized as the solution method. The calculated nonlinear stress-strain relations of unidirectional composite panels (AS4/PEEK) accorded well with the test data. The progressive damage analysis results agreed with the experimental and simulating results of a relevant research. The method for nonlinear behavior and progressive damage modeling of composites in state-based peidynamics was explored and verified.

Key words:peridynamics; composite materials; nonlinear behavior; progressive damage; dynamic relaxation method

收稿日期:2015-10-26.浙江大學學報(工學版)網址: www.journals.zju.edu.cn/eng

基金項目:國家自然科學基金資助項目(11372192);上海市自然科學基金資助項目(13ZR1422200);上海市科委重大項目科技攻關資助項目(12DZ1100302).

作者簡介:劉肅肅(1991-),男,碩士生,從事復合材料結構強度等研究.ORCID: 0000-0002-2111-2074. E-mail: fowyaring@sjtu.edu.cn通信聯系人:余音,女,副教授.ORCID: 0000-0003-3587-0260. E-mail: yuyin@sjtu.edu.cn

DOI:10.3785/j.issn.1008-973X.2016.05.025

中圖分類號:TB 332

文獻標志碼:A

文章編號:1008-973X(2016)05-0993-08

猜你喜歡
復合材料
淺談現代建筑中新型復合材料的應用
金屬復合材料在機械制造中的應用研究
敢為人先 持續創新:先進復合材料支撐我國國防裝備升級換代
民機復合材料的適航鑒定
復合材料無損檢測探討
電子測試(2017年11期)2017-12-15 08:57:13
復合材料性能與應用分析
PET/nano-MgO復合材料的性能研究
中國塑料(2015年6期)2015-11-13 03:02:54
ABS/改性高嶺土復合材料的制備與表征
中國塑料(2015年11期)2015-10-14 01:14:14
聚乳酸/植物纖維全生物降解復合材料的研究進展
中國塑料(2015年8期)2015-10-14 01:10:41
TiO2/ACF復合材料的制備及表征
應用化工(2014年10期)2014-08-16 13:11:29
主站蜘蛛池模板: 国产精品区视频中文字幕| 日韩在线网址| 亚洲中文制服丝袜欧美精品| 亚洲aaa视频| 一级毛片不卡片免费观看| 亚洲精品在线影院| 国产91色| 中文字幕色在线| 中文字幕伦视频| 一级全免费视频播放| 97免费在线观看视频| 亚洲第一成人在线| 波多野结衣中文字幕一区二区| 亚洲中文字幕日产无码2021 | 国产欧美亚洲精品第3页在线| 影音先锋丝袜制服| 国产成人综合日韩精品无码首页 | 在线观看国产精美视频| 亚洲精品国产自在现线最新| 久久九九热视频| 亚洲日产2021三区在线| 日韩欧美国产精品| 青青久视频| 专干老肥熟女视频网站| 久久国产精品麻豆系列| 在线观看精品国产入口| 三区在线视频| 国产精品jizz在线观看软件| 日本一区二区三区精品视频| 中文成人在线| 国产不卡在线看| 国产一级二级在线观看| 国产成人艳妇AA视频在线| 国产激情无码一区二区免费| 熟女成人国产精品视频| 久久成人免费| 成人午夜视频网站| 国产精品熟女亚洲AV麻豆| 色综合天天娱乐综合网| 国产sm重味一区二区三区| 无码精油按摩潮喷在线播放 | 美女被操91视频| 久久综合激情网| 久久夜色精品| 欲色天天综合网| 国产激情无码一区二区APP | 亚洲综合香蕉| 国产亚洲日韩av在线| 欧美午夜在线观看| h网址在线观看| 亚洲av色吊丝无码| 毛片网站在线看| 欧美黄网站免费观看| 欧美一级高清片久久99| 专干老肥熟女视频网站| 青青青视频蜜桃一区二区| 亚洲成人网在线播放| 香蕉视频在线观看www| 91黄视频在线观看| 国产成人免费视频精品一区二区| 一级高清毛片免费a级高清毛片| 亚洲乱强伦| 国产在线观看一区精品| 久久永久精品免费视频| 色综合天天综合| 国产成人AV综合久久| 青青草国产在线视频| 新SSS无码手机在线观看| 久久一级电影| 免费在线看黄网址| 国产JIZzJIzz视频全部免费| 国产传媒一区二区三区四区五区| 中日韩欧亚无码视频| 秋霞一区二区三区| 欧美第一页在线| 欧美日韩一区二区在线播放 | 精品国产自在现线看久久| 美女一级毛片无遮挡内谢| 国产精品原创不卡在线| 777午夜精品电影免费看| 亚洲欧洲自拍拍偷午夜色无码| 精品中文字幕一区在线|