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

基于晶體塑性理論的疲勞裂紋起始數(shù)值模擬

2016-09-09 01:52:10劉俊卿曹書文
航空材料學(xué)報(bào) 2016年2期
關(guān)鍵詞:裂紋變形模型

劉俊卿, 李 蒙, 左 帆, 劉 紅, 曹書文

(西安建筑科技大學(xué) 理學(xué)院,西安 710055)

?

基于晶體塑性理論的疲勞裂紋起始數(shù)值模擬

劉俊卿,李蒙,左帆,劉紅,曹書文

(西安建筑科技大學(xué) 理學(xué)院,西安 710055)

在疲勞荷載作用早期,材料在晶體尺度出現(xiàn)裂紋的萌生,為了研究疲勞短裂紋的發(fā)展,利用MonteCarlo法建立了多晶體晶粒集合的Voronoi有限元模型,并在AnsysUsermat子程序接口下編寫了晶體塑性本構(gòu)方程子程序,修正了以拉伸硬化為主的疲勞裂紋起裂的計(jì)算機(jī)模擬方法,并結(jié)合TANAKA和MURA的位錯(cuò)偶極子模型,模擬了晶粒集合在疲勞荷載作用下的裂紋萌生,最后與文獻(xiàn)實(shí)驗(yàn)數(shù)據(jù)對比,這一修正的方法更加符合實(shí)驗(yàn)結(jié)果和宏觀現(xiàn)象。

疲勞;晶體塑性本構(gòu);裂紋起始;有限元模擬

材料在疲勞載荷作用下,裂紋在微觀小尺度上萌生,然后擴(kuò)展為宏觀裂紋,并最終在疲勞壽命的最后循環(huán)中導(dǎo)致試樣斷裂[1-3],微觀研究表明[4-14],在大多數(shù)疲勞問題中,裂紋的微觀活動占據(jù)疲勞壽命的大部分,因此對疲勞裂紋形核的研究具有重要意義。材料疲勞裂紋形核往往是從材料多處起裂[6-8],由材料應(yīng)力集中點(diǎn)、滑移帶或者微觀缺陷處起始,并在疲勞荷載反復(fù)作用下,短裂紋群擴(kuò)展并交匯形成宏觀裂紋。以前的研究往往僅局限于單一短裂紋的擴(kuò)展分析[2-5,9-13],因而對裂紋群體作用關(guān)注不夠,或提出以拉伸硬化為主的多裂紋萌生模擬[6],但精確度不高。本工作考慮晶體塑性本構(gòu)模型[15-18],改進(jìn)以往的模擬方法,模型使用Usermat編寫,結(jié)合有限元方法模擬了晶粒集合在疲勞荷載作用下的裂紋萌生,與已有的實(shí)驗(yàn)數(shù)據(jù)對比,其模擬結(jié)果更加接近真實(shí)實(shí)驗(yàn)數(shù)據(jù)。

1 多晶體晶粒集合模型的建立

1.1晶體模型的模擬

疲勞裂紋的萌生環(huán)境即多晶體模型,多采用MonteCarlo直晶界組織仿真模型[15],該模型模擬金屬晶粒組織成核、結(jié)晶、晶粒生長及晶界成型過程,具有金屬結(jié)晶的主要特征,圖1所示為一個(gè)金屬多晶體模型實(shí)例,本研究采用這一模型。

圖1 多晶體晶體晶界Voronoi圖Fig.1 Voronoi chart of multi crystal

1.2單個(gè)晶粒開裂模型

假定該模型中單個(gè)晶粒為立方晶體系,因此每個(gè)晶粒都是彈性各向異性的。TANAKA和MURA[16]提出了疲勞裂紋的位錯(cuò)偶極子模型,基于位錯(cuò)自由能的增加,揭示了裂紋萌生及短裂紋擴(kuò)展規(guī)律。TANAKA和MURA認(rèn)為,在材料局部滑移帶施加剪切應(yīng)力循環(huán)下,位錯(cuò)偶極子不斷向滑移層兩邊堆積,引起自由能的增加,首次達(dá)到最大應(yīng)力幅時(shí),位錯(cuò)自由能表示為

(1)

式中:μ是剪切模量,κ是位錯(cuò)摩擦力,τ為剪切應(yīng)力峰值,l為滑移帶中心距離位錯(cuò)運(yùn)動受阻處即晶界距離。

之后每半周期增加的自由能為

(2)

即n次循環(huán)載荷下總自由能為

(3)

結(jié)合荷載譜,即可以給定自由能的總值,本研究討論單向恒應(yīng)力比為-1的情況,此時(shí)忽略U1,公式化為

(4)

取Ws為沿滑移帶單位面積斷裂能,當(dāng)令

U=4lWs

(5)

則可得到恒幅荷載下單個(gè)裂紋萌生的疲勞壽命ns

(6)

2 晶體塑性本構(gòu)方程

2.1晶體的變形分解

采用單晶體塑性本構(gòu)方程來描述晶體變形及其塑性滑移[17-18],晶體的變形包括晶格畸變和位錯(cuò)的滑移,其中晶格畸變(包括其剛性轉(zhuǎn)動)可看作連續(xù)介質(zhì)力學(xué)中的彈性變形,采用彈性力學(xué)處理;其不可逆部分,即位錯(cuò)的滑移,在理想狀態(tài)下,可假設(shè)位錯(cuò)滑移在晶粒內(nèi)部是均勻的,并用變形梯度場變量來描述。

在小變形的情況下,晶體變形可分解為彈性變形和塑性變形的疊加,其中彈性變形反映了晶體的伸縮和旋轉(zhuǎn),是晶格畸變的表示,塑性變形是由晶體沿著特定的滑移系滑移產(chǎn)生的,一個(gè)滑移系由固定的滑移面和滑移方向組成,其具體的乘法分解圖如圖2所示。

晶體變形梯度張量按照晶體塑性理論乘法分解為隨晶體伸縮和旋轉(zhuǎn)的彈性梯度張量及隨滑移系滑移的塑性梯度張量,其具體變形梯度張量為

圖2 晶體的變形分解Fig.2 Decomposition of deformation of crystal

(7)

式中:Fe為彈性變形梯度張量,F(xiàn)p為塑性變形梯度張量。

(8)

(9)

同理依照(9)式,變形率張量和旋轉(zhuǎn)率張量也可以分解為彈性部分和塑性部分,將變形率張量D分解為

(10)

其中塑性部分

(11)

式中:

(12)

旋轉(zhuǎn)率張量W分解為

(13)

其中塑性部分

(14)

式中:

(15)

2.2晶體塑性本構(gòu)方程

假設(shè)單晶體的彈性性質(zhì)不受滑移系變形的影響,因此可采用經(jīng)典的彈性模量張量和變形率張量二點(diǎn)乘積形式的彈性本構(gòu)方程:

(16)

(17)

考慮到

(18)

以初始構(gòu)形為基準(zhǔn)的剛體導(dǎo)數(shù)是

(19)

定義:

βα=Wασ-σWα

(20)

即有

(21)

代入(19)式,即可得到晶體塑性本構(gòu)方程

(22)

2.3硬化規(guī)律

文獻(xiàn)[17]采用Schmid定律描述流變應(yīng)力,對于臨界的滑移系,具體硬化規(guī)律為

(23)

對于非臨界的滑移系,恒有

(24)

同時(shí)采用經(jīng)典的線性臨界剪應(yīng)力、滑移剪切率關(guān)系,有

(25)

式中:σc為屈服臨界剪應(yīng)力,其初始值記為σ0,hαβ是滑移硬化系數(shù),下標(biāo)表示取自不同滑移系組合,其中當(dāng)α=β時(shí),hαβ稱為自硬化系數(shù),當(dāng)α≠β時(shí),hαβ為潛在硬化系數(shù),對于不同滑移面,其值一般不同,這里采用Peirce提出的硬化系數(shù)簡化模型,即

(26)

式中:h是自硬化系數(shù),q為潛在硬化率,均為常量,q的取值一般建議為1到1.5。

3 裂紋萌生的有限元模擬及分析

模擬金屬材料在室溫下的疲勞裂紋萌生,在疲勞早期其晶體滑移方向主要沿{110}族晶面,按照這個(gè)開動的滑移面分析,文中采用Ansys有限元軟件及其接口Fortran的子程序。

3.1參數(shù)的選取

為了與文獻(xiàn)[6]實(shí)驗(yàn)結(jié)果對照,同樣選取F82H金屬材料,其彈性模量取為單晶純鐵測定值,見表1。

表1 單晶純鐵的彈性模量/GPaTable 1 Elastic modulus of monocrystalline iron

按照文獻(xiàn)[19-20]取位錯(cuò)摩擦力κ取自純鐵的初始屈服應(yīng)力,即當(dāng)且僅當(dāng)位錯(cuò)驅(qū)動力大于位錯(cuò)摩擦力時(shí)發(fā)生滑移系的切向滑移,材料滑移系長度取為平均晶粒尺寸,材料剪切模量和單位面積斷裂功均為常數(shù),則公式(6)等價(jià)于

(27)

式中Ceq為一常數(shù)。

(28)

具體參數(shù)列于表2。

3.2駐留滑移帶

分析模型在循環(huán)荷載作用下的初始塑性滑移云圖如圖3所示,為單向恒應(yīng)力比為-1的0.6%應(yīng)變范圍作用下的初始滑移云圖。云圖顯示了最大駐留滑移帶出現(xiàn)在應(yīng)力集中的個(gè)別晶粒,滑移沿固定方向,晶粒之間應(yīng)力差別較大,滑移程度亦不相同,且晶界是塑性滑移及應(yīng)力集中的天然屏障,塑性滑移多集中在晶界范圍內(nèi),應(yīng)力集中多出現(xiàn)在多個(gè)晶體的晶界交集處。按照這個(gè)云圖模型,位錯(cuò)將優(yōu)先發(fā)生于滑移程度大的晶粒中,并在滑移過程中集中于晶界處。

圖3 晶粒集合的滑移系云圖Fig.3 Chart of multicrystal slip band

通過體積平均化的方法,提取每個(gè)晶粒的平均滑移剪切應(yīng)力,代入TANAKA和MURA模型公式中,可以求出晶體萌生一個(gè)晶界長度裂紋時(shí)的荷載循環(huán)數(shù)。Ansys在引入裂紋后,更新應(yīng)力,重復(fù)以上步驟,可以模擬疲勞循環(huán)初始階段裂紋萌生以及其裂紋群體之間相互作用、相互影響的情況。

3.3裂紋萌生及其走向

分析循環(huán)荷載進(jìn)程中裂紋群分布情況,如圖4示例為0.6%荷載作用下裂紋群演化云圖。多晶體在沒有裂紋出現(xiàn)的初始情況下,塑性滑移梯度較大,晶粒和晶粒之間有很大的應(yīng)力差,應(yīng)力集中在少數(shù)晶粒中;隨著裂紋首先出現(xiàn)在應(yīng)力集中的晶粒中,晶粒開裂后多晶體模型出現(xiàn)明顯的應(yīng)力重新分布情況,塑性滑移主要集中在短裂紋群周圍晶粒,遠(yuǎn)離短裂紋群的晶粒群則表現(xiàn)出較低滑移變形,且塑性滑移梯度很小,裂紋群的發(fā)展主要集中在初始裂紋周圍,并沿垂直于疲勞荷載方向發(fā)展,這一結(jié)果和宏觀疲勞實(shí)驗(yàn)相一致,印證了模擬的準(zhǔn)確性。

圖4 不同循環(huán)次數(shù)下裂紋群有限元模擬Fig.4 Finite element simulation of crack group at different fatigue cycles(a)N=1193;(b)N=2333;(c)N=3223;(d)N=4506

3.4采用晶體塑性模型構(gòu)建的多晶體裂紋萌生模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對比

圖5 0.5%應(yīng)變荷載作用下模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對比Fig.5 Comparison of simulation results and experimental data under 0.5% strain load

圖6 0.6%應(yīng)變荷載作用下模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對比Fig.6 Comparison of simulation results and experimental data under 0.6% strain load

圖7 0.76%應(yīng)變荷載作用下模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對比Fig.7  Comparison of simulation results and experimental data under 0.76% strain load

在圖5~7中所示三種應(yīng)力水平下,當(dāng)循環(huán)數(shù)不超過500~700次時(shí),短裂紋數(shù)密度均與實(shí)驗(yàn)結(jié)果符合較好,但當(dāng)循環(huán)數(shù)較大時(shí),實(shí)驗(yàn)數(shù)據(jù)出現(xiàn)了拐點(diǎn),與模擬結(jié)果出現(xiàn)較大偏差,根據(jù)文獻(xiàn)[21]等的研究,這是由于裂紋在實(shí)際循環(huán)荷載作用下出現(xiàn)了穿晶擴(kuò)展,即當(dāng)短裂紋萌生階段結(jié)束,短裂紋群密度達(dá)到臨界值時(shí),在疲勞荷載持續(xù)作用下將導(dǎo)致部分裂紋穿過晶界,擴(kuò)展至相鄰晶粒中,這一部分裂紋中的絕大部分將在達(dá)到二倍于晶粒尺寸時(shí)停止擴(kuò)展,出現(xiàn)宏觀裂紋。

4 結(jié)論

(1)多晶體金屬在疲勞荷載作用早期,其裂紋是從個(gè)別應(yīng)力集中較高的晶粒中萌生,并擴(kuò)展至晶界處停止擴(kuò)展。

(2)晶界作為晶粒的屏障,將位錯(cuò)的滑移限制在晶粒內(nèi)部,導(dǎo)致晶體變形不均勻。

(3)當(dāng)晶粒開裂后,應(yīng)力及滑移變形重新分布,集中在短裂紋周圍晶粒處,加速了其裂紋的出現(xiàn)和擴(kuò)展,導(dǎo)致短裂紋出現(xiàn)集中,并沿垂直于荷載方向分布。

(4)裂紋數(shù)密度達(dá)到臨界時(shí),裂紋起始階段結(jié)束,疲勞過程進(jìn)入宏觀裂紋擴(kuò)展階段。

[1]SCHIJVEJ.Fatigueofstructureandmaterials[M].[S.l.]:Springer,2009:11-43.

[2] 吳楠,張顯程,王正東,等.GH4169合金在650℃下疲勞小裂紋萌生和擴(kuò)展行為[J].航空材料學(xué)報(bào),2015,35(5):71-76.

(WUN,ZHANGXC,WANGZD,et al.InitiationandpropagationofsmallfatiguecrackofGH4169alloyat650℃[J].JournalofAeronauticalMaterials,2015,35(5):71-76.)

[3] 王凱,閆志峰,王文先,等.循環(huán)荷載作用下鎂合金溫度演化及高周疲勞性能預(yù)測[J].材料工程,2014(1):85-89.

(WANGK,YANZF,WANGWX,et al.Temperatureevolutionandfatiguepropertiespredictionforhighcyclefatigueofmagnesiumalloy[J].JournalofMaterialsEngineering,2014(1):85-89.)

[4]BJERKéNC,MELINS.Atooltomodelshortcrackfatiguegrowthusingadiscretedislocationformulation[J].InternationJournalofFatigue,2003,25(6):559-566.

[5]CHRISTHJ,F(xiàn)RITZENCP,K?STERP.Micromechanicalmodelingofshortfatiguecracks[J].CurrentOpinioninSolidStateandMaterialsScience,2014,18(4):205-211.

[6] 黃新躍,BRUECKNER-FOITA.多點(diǎn)裂紋起始過程的計(jì)算機(jī)模擬[J].航空材料學(xué)報(bào),2008,28(2):30-33.

(HUANGXY,BRUECKNER-FOITA.Computersimulationonprocessofmulti-crackinitiation[J].JournalofAeronauticalMaterials,2008,28(2):30-33.)

[7]AZARAS,SVENSSONLE,NYHUSB.Effectofcrystalorientationandtextureonfatiguecrackevolutioninhighstrengthsteelwelds[J].InternationalJournalofFatigue,2015(77):95-104.

[8]KRUPPU,ALVAREZ-ARMASBI.Shortfatiguecrackpropagationduringlow-cycle,highcycleandvery-high-cyclefatigueofduplexsteel-anunifiedapproach[J].InternationalJournalofFatigue,2014(65):78-85.

[9]STRUBBIAR,HERES,ALVAREZ-ARMASAI,et al.ShortfatiguecracksnucleationandgrowthinleanduplexstainlesssteelLDX2101[J].MaterialsScienceandEngineeringA,2014,615:169-174.

[10]譚曉明,張丹峰,陳躍良,等.基于疲勞裂紋萌生機(jī)理的鋁合金疲勞壽命可靠性評估方法[J].航空材料學(xué)報(bào),2014,34(2):84-89.

(TANXM,ZHANGDF,CHENYL,et al.Probabilisticmethodtopredictfatiguelifebasedoncrackinitiatingmicro-mechanismofaluminumalloy[J].JournalofAeronauticalMaterials,2014,34(2):84-89.)

[11]童弟華,吳學(xué)仁,劉建中,等.基于小裂紋理論的鑄造鈦合金ZTC4疲勞壽命預(yù)測[J].材料工程,2015,43(6):60-65.

(TONGDH,WUXR,LIUJZ,et al.FatiguelifepredictionofcasttitaniumalloyZTC4basedonthesmallcracktheory[J].JournalofMaterialsEngineering,2015,43(6):60-65.)

[12]ZHANGTT,JIANGJ,SHOLLOCKBA,et al.Sliplocalizationandfatiguecracknucleationnearanon-metallicinclusioninpolycrystallinenickel-basedsuperalloy[J].MaterialsScienceandEngineeringA,2015,641:328-339.

[13]張麗,吳學(xué)仁.基于小裂紋理論的GH4169高溫合金的疲勞全壽命預(yù)測[J].航空材料學(xué)報(bào),2014,34(6):75-83.

(ZHANGL,WUXR.Fatigue-lifepredictionmethodbasedonsmall-cracktheoryinGH4169superalloy[J].JournalofAeronauticalMaterials,2014,34(6):75-83.)

[14]ZHANGKS,JUJW,LIZH,et al.Micromechanicsbasedfatiguelifepredictionofapolycrystallinemetalapplyingcrystalplasticity[J].MechanicsofMaterials,2015,85:16-37.

[15]司良英,鄧關(guān)宇,呂程,等.基于voronoi圖的晶體塑性有限元多晶幾何建模[J].材料與冶金學(xué)報(bào),2009,8(3):193-197.

(SILY,DENGGY,LVC,et al.Polycrystalgeometrymodelingofcrystalplasticityfiniteelementmethodwithvoronoidiagram[J].JournalofMaterialsandMetallurgy,2009,8(3):193-197.)

[16]TANAKAK,MURAT.Atheoryoffatiguecrackinitiationatinclusions[J].MetallurgicalandMaterialsTransactionsA.1982,13(1):117-123.

[17]MANONUKULA,DUNNEFPE.High-andlow-cyclefatiguecrackinitiationusingpolycrystalplasticity[J].TheRoyalSociety,2003,460:1881-1903.

[18]王自強(qiáng),段祝平.塑性細(xì)觀力學(xué)[M].北京:科學(xué)出版社,1995:109-115.

[19]孔金星,陳輝,何寧,等.純鐵材料動態(tài)力學(xué)性能測試及本構(gòu)模型[J].航空學(xué)報(bào),2013,35(7):2063-2070.

(KONGJX,CHENH,HEN,et al.Dynamicmechanicalpropertytestsandconstitutivemodelofpureironmaterial[J].ActaAeronauticaetAstronauticasinica,2013,35(7):2063-2070.)

[20]孫賓,李兆霞.描述微裂紋成核與擴(kuò)展的疲勞損傷多尺度模型及其應(yīng)用[J].東南大學(xué)學(xué)報(bào),2014,44(2):333-338.

(SUNB,LIZX.Multi-scalefatiguedamagemodelformicro-cracksnucleationandgrowthanditsapplication[J].JournalofSoutheastUniversity,2014,44(2):333-338.)

[21]白以龍,柯孚久,夏蒙棼.固體中微裂紋系統(tǒng)統(tǒng)計(jì)演化的基本描述[J].力學(xué)學(xué)報(bào),1991,23(3):290-298.

(BAIYL,KEFJ,XIAMF.Formulationofstatisticalevolutionofmicro-cracksinsolids[J].JournalofMechanics,1991,23(3):290-298.)

(SchoolofScience,Xi′anUniversityofArchitecture&Technology,Xi′an710055,China)

(責(zé)任編輯:徐永祥)

Numerical Simulation of Fatigue Crack Initiation Using Crystallographic Constitutive Equation

LIU Junqing,LI Meng,ZUO Fan,LIU Hong,CAO Shuwen

Attheearlystageoffatigueloading,shortcrackinitiatedandthenmacrocracknucleationappearedoncrystalscaleinmaterials.Inordertostudythedevelopmentoffatiguecrack,theVoronoifiniteelementmodelofmulticrystalgrainswasestablishedbyMonteCarlomethod;crystallographicconstitutiveequationwasprogrammedwiththefiniteelementmodelbasedontheusersubroutineinANSYSUsermat,acorrectionofthecomputersimulationmethodwhichmainlyconsideredtensionstiffeningeffectforfatiguecrackinitiationwasdeveloped.WiththeconsiderationofthedislocationdipolemodelsuggestedbyTANAKAandMURA,fatiguecrackinitiationofmulticrystalgrainswassimulated.Acomparisonwithcurrentlycollecteddatashowsthatthecorrectedmethodwasmoreconsistedwithtestresultsandmacrophenomenon.

fatigue;crystallographicconstitutiveequation;crackinitiation;finiteelementsimulation

2015-08-10;

2015-10-08

陜西省教育廳專項(xiàng)基金資助項(xiàng)目(15JK1382)

劉俊卿(1957—),男,工學(xué)博士,教授,主要從事工程力學(xué)的教學(xué)與研究,(E-mail)37825pf@sina.com。

10.11868/j.issn.1005-5053.2016.2.012

TG113.25+5

A

1005-5053(2016)02-0074-06

猜你喜歡
裂紋變形模型
一半模型
裂紋長度對焊接接頭裂紋擴(kuò)展驅(qū)動力的影響
重要模型『一線三等角』
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
“我”的變形計(jì)
例談拼圖與整式變形
會變形的餅
3D打印中的模型分割與打包
主站蜘蛛池模板: 成年人国产视频| 日韩免费毛片视频| 国产精品2| 国产一区二区福利| 国产日韩欧美视频| 国产极品粉嫩小泬免费看| 久久精品免费看一| 久久精品这里只有国产中文精品| 欧洲熟妇精品视频| 午夜精品久久久久久久无码软件 | 日本一本正道综合久久dvd| 97久久精品人人| 美女内射视频WWW网站午夜| 国产在线小视频| 玖玖精品在线| 亚洲精品无码高潮喷水A| 青青草原国产一区二区| 亚洲大尺码专区影院| 国产人碰人摸人爱免费视频| 国产亚洲精品资源在线26u| 久久99蜜桃精品久久久久小说| 国产一区在线观看无码| 久久黄色影院| 欧美日在线观看| 欧美一级高清片欧美国产欧美| 国产精品自拍露脸视频| 亚洲AV成人一区二区三区AV| 亚洲国产在一区二区三区| 欧美日韩高清在线| 日韩av电影一区二区三区四区| 国产在线一区视频| 四虎国产精品永久在线网址| 污网站在线观看视频| 婷婷综合在线观看丁香| 四虎国产精品永久一区| 国产一区亚洲一区| 午夜激情婷婷| 久久精品波多野结衣| 2021精品国产自在现线看| 丰满人妻一区二区三区视频| 国产精品制服| 夜夜操天天摸| 国产剧情伊人| 国产亚洲欧美在线视频| 国产主播福利在线观看| 亚洲区欧美区| 青青草原国产| 国产尤物在线播放| 国产成人禁片在线观看| 无码精品福利一区二区三区| 日韩国产综合精选| 暴力调教一区二区三区| 日韩精品少妇无码受不了| 国产日韩欧美成人| 欧美一区二区精品久久久| 91网站国产| 伊人天堂网| 久久 午夜福利 张柏芝| 亚洲香蕉在线| 一区二区三区在线不卡免费| 国产成人综合亚洲欧洲色就色| 国产产在线精品亚洲aavv| 欧美成人一区午夜福利在线| 亚洲婷婷丁香| 亚洲综合极品香蕉久久网| www.99在线观看| 久久国产精品影院| 精品黑人一区二区三区| 久久这里只有精品国产99| 人妻丰满熟妇啪啪| 婷婷六月综合网| 99re在线观看视频| 毛片免费视频| 午夜视频免费试看| 国产丝袜第一页| 欧美日韩高清| 日韩无码黄色网站| 国产亚洲精品97在线观看| 国产大全韩国亚洲一区二区三区| AV在线天堂进入| 欧美日韩久久综合| 国产va在线|