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

平面冪律塑性I型裂紋尖端應力場全解1)

2023-02-25 02:24:36王志強蔡力勛黃茂波
力學學報 2023年1期
關(guān)鍵詞:裂紋效應

王志強 蔡力勛 黃茂波

(西南交通大學力學與航空航天學院,應用力學與結(jié)構(gòu)安全四川省重點實驗室,成都 610031)

引言

在航空航天、船舶、核電和化工等工程中,裂紋構(gòu)件的裂紋尖端應力場與斷裂強度是結(jié)構(gòu)安全評價的重要基礎(chǔ).

Cherepanov[1]和Rice[2]分別在1967 年和1968年針對含I 型裂紋無限體各自獨立提出與裂紋尖端環(huán)向積分路徑無關(guān)的J積分.隨后,Hutchinson[3]、Rice 等[4](HRR) 針對冪律硬化材料,基于J積分和塑性解析求解在理論上給出描述裂尖應力場的HRR 奇異解,該解關(guān)于角度分布的解析描述一直未得到解決.1983 年,Shih[5]考慮冪律塑性應力?應變關(guān)系,通過有限元分析計算了HRR 解中的積分常數(shù)IN和無量綱角度函數(shù),相關(guān)結(jié)果均由數(shù)值表格形式給出.1991 年,Betegon 等[6]利用邊界層方法研究小范圍屈服條件下不同T應力(T11,作用于裂紋面且垂直于裂紋前沿的應力分量)對彈塑性應力場的影響,提出J-T11方法,該方法以T11作為表征裂尖約束程度的彈性參量用于描述I 型裂紋有限尺寸試樣的彈塑性裂尖場.后繼研究[7-9]表明,具有彈性屬性的T應力不適用于大范圍或全范圍屈服條件.1991 年,對于中心裂紋板和單邊裂紋彎曲兩個有限幾何尺寸試樣,O’Dowd等[10-11]為反映有限尺寸對裂尖場的幾何約束效應,考慮在HRR 場基礎(chǔ)上進行約束補償,提出J-Q方法,該方法對不同材料和試樣的不同幾何尺寸,完全依賴有限元分析,所獲得Q值有不同結(jié)果,無普適性Q方程,應用不便.1986 年Li 等[12]提出含高次展開項的冪律塑性漸進解作為HRR 場約束補償項的新思路,1991 年Sharma 等[7]及1993 年Yang 等[13]采用一次高階漸近項作為HRR 場補償項提出了裂尖場描述模型,該模型通過復雜數(shù)值分析僅對特定材料和幾何工況給出漸進解參數(shù),相較HRR 場描述精度上有了提升,不過,對于不同材料和不同幾何工況需給出不同漸進解參數(shù),普適性仍然不足.Yang等[8,14]基于裂紋尖端應力場的3 項漸近展開式,給出了第二和第三展開項的系數(shù)之間的近似關(guān)系,提出了J-A2理論.Nikishkov 等[9,15]采用類似方法,提出J-A理論,本質(zhì)上和J-A2理論一樣[16-17].隨后,Ding 等[18-19]進行廣泛的有限元分析,獲得經(jīng)驗方程,預測約束參數(shù)A,Q,A2等參數(shù).上述應力漸進解均無關(guān)于角度的描述函數(shù).事實上,任何復雜函數(shù)關(guān)系都可以由不同函數(shù),如高階級數(shù)函數(shù)、特殊函數(shù)等函數(shù)同時表征,描述形式并非唯一.Ji 等[20-21]通過數(shù)值分析提出特殊的指數(shù)函數(shù)來表征特定材料具有奇異性的彈塑性應力場,雖然該修正的奇異性指數(shù)并沒有適用于不同材料和幾何條件的統(tǒng)一表達式,但這是針對塑性裂尖奇異場研究提出的一個重要概念和新方法.

此外,也有學者在研究斷裂韌性時嘗試將裂尖塑性區(qū)與裂尖約束聯(lián)系起來.Anderson 等[22]提出A-D 法,即用裂尖前最大正應力σ1與屈服應力σy的比值來表征約束,比值越大,約束越大.Mostafavi等[23-26]對A-D 法進行了改進,用試樣或結(jié)構(gòu)斷裂時的裂尖塑性區(qū)面積Ac與標準高約束平面應變斷裂韌性試樣斷裂時的裂尖塑性區(qū)面積Assy(參考塑性區(qū)面積)的比值 φ來表征約束水平來表征約束水平,用φ可以將任何約束(不同的面內(nèi)和面外約束)水平下的斷裂韌性JIC試驗數(shù)據(jù)關(guān)聯(lián)起來,建立一致性的斷裂韌性與約束的關(guān)聯(lián)線.

本文考慮冪律塑性本構(gòu)關(guān)系,基于能量密度等效原理[27-29]和量綱分析提出具有中值能量密度的代表性體積單元 (representative volume element,RVE)的等效應力(σM)解析方程,以能量密度等效體積單元的應力因子為應力特征量,提出平面冪律塑性I 型裂紋尖端應力場半解析模型.針對有限平面應變和平面應力緊湊拉伸(compact tension,CT)試樣和單邊裂紋彎曲(single edge bend,SEB)試樣,給出應力場半解析模型參數(shù)的確定方法,并進行數(shù)值分析的模型驗證.

1 理論模型

1.1 等效應力解析方程

假設(shè)延性材料均勻連續(xù)、各向同性,對于受單向加載的任意構(gòu)元(見圖1),Chen 等[27-29]基于復雜應力與單軸應力下 RVE 應變能密度等效(圖1中任意A點RVE)及應變能密度中值等效(圖1 中值M點RVE),提出彈塑性應變能U控制方程為

圖1 能量密度等效示意圖Fig.1 Schematic diagram of energy density equivalence

式中,Veff為有效變形域體積,εeq?M為RVE 等效應變.

冪律塑性條件下,設(shè)材料應力應變關(guān)系符合冪律

式中,εp?eq為等效塑性應變,σeq為等效應力,K為材料的應變硬化系數(shù),N為應力硬化指數(shù).進一步假定Veff,εp?eq與加載線位移hp之間符合冪律關(guān)系

式中,特征體積V*為A*h*,A*為特征面積;h*為特征位移;k1和k2分別為有效體積系數(shù)和有效體積指數(shù),k3和k4分別為等效塑性應變系數(shù)和等效塑性應變指數(shù).

由式(1)、式(2)和式(3)得到冪律塑性應變能Up控制方程

根據(jù)文獻[27-28],受單向加載構(gòu)元的載荷?位移關(guān)系表示為

式中,P為載荷,P*為特征載荷,ζp和mp分別為冪律塑性條件下的無量綱載荷?位移關(guān)系系數(shù)和指數(shù),m為體積折減系數(shù).

能量密度中值點RVE 的等效應力σM與載荷P及加載線位移hp與構(gòu)元特征尺寸分別構(gòu)成相似判據(jù)π1,π2,設(shè)為

假定π1和π2之間關(guān)系符合冪律

式中,k5和k6分別為等效應力系數(shù)和等效應力指數(shù).由式(6)和式(7)可得

聯(lián)立式(4)和式(8),整理可得

在準靜態(tài)加載和冪律塑性條件下,根據(jù)功-能原理,載荷對構(gòu)元試樣所做的功等于構(gòu)元試樣的變形能(U),進而由上式可得

將式(5)的載荷P代入式(10)中整理可得

上式等號左右項積分整理可得

考慮到等號右端項為常數(shù),則等號左端項的無量綱變量hp/h*的指數(shù)項必為0,則1?(k2+k4+k6)=0,此時右端項必為1,進而可得

結(jié)合冪律應變式(4)、關(guān)于P與hp的無量綱式(8)及k5,k6表達式(13),可得σM解析方程,定義該應力為應力因子

1.2 平面I 型裂紋尖端應力場半解析模型

對于I 型裂紋構(gòu)元,有效變形體積在加載過程中保持不變,則式(4)中k2為0;由文獻[30]可知冪律塑性條件下中值點RVE 等效應變與試樣加載線位移之間呈線性關(guān)系,即,式(4)中k4為1.則k5=(k1k3)?1,k6=0.進一步由式(4)、式(5) RVE 等效應力的解析方程可簡化為

對于有限尺寸試樣,裂紋尖端應力場不再符合HRR 解.考慮HRR 解的奇異性,對裂尖應力,以應力因子σM作為特征應力,則基于特征應力的無量綱裂尖應力表為

式中,當sub 為eq 時,σeq是等效應力,當sub 為ij時,σij是應力張量(ij=11 時為rr,ij=22 時為θθ,ij=12 時為rθ),σM是特征應力,L是特征長度(通常取為1 mm),f和λ為關(guān)于N與θ的函數(shù),m為試樣尺寸a/W的綜合指數(shù).

考慮到HRR 場奇異性指數(shù)為1/(N+1),對分子、分母各疊加一個關(guān)于θ增量β’sub,即(1+β’sub)/(N+1+β’sub),可進一步假設(shè)fsub=βsub/(N+βsub),βsub=1+β’sub,βsub為與三角函數(shù)相關(guān)的特殊函數(shù),即

式中,dsubp(p=0,1,2,3)為待定參數(shù).

有限元計算表明冪律塑性條件下裂尖等效應力等值線在平面應變和平面應力下分別具有蝶翅輪廓和扇貝輪廓的特征,為描述這兩類特殊的等應力輪廓線,設(shè)應力分布系數(shù)λsub為三角特殊函數(shù)

式中,psubk(k=0,1,2)與qsubm(m=0,1,2,3)為待定參數(shù),psubk(k=0,1,2)與N相關(guān),csubk(k=0,1,2),bsubk(k=0,1,2) 為psubk的待定參數(shù).

修正綜合指數(shù)msub與N相關(guān),假設(shè)

式中,msub0和msub1為待定參數(shù).

聯(lián)立式(16)~ 式(19)可得平面冪律塑性I 型裂紋尖端應力場半解析模型

若式(16) 中 λeq·(L/r)feq取為不同程度常數(shù)C,則由式(17)、式(18)有

選擇合適的級數(shù)項數(shù)和相應參數(shù),可在數(shù)學上使上述函數(shù)能恰當描述蝶翅輪廓式或扇貝輪廓式等應力線,進而由式(20)可推導出平面I 型裂紋有限尺寸試樣等效應力等值線表達式

2 有限元分析

2.1 分析條件

本文采用有限元分析(finite element analysis,FEA)軟件ANSYS19.0 對CT 試樣和SEB 試樣進行研究,本文所針對CT 試樣和SEB 試樣的建模參數(shù)標定及模型驗證都針對有限尺寸試樣.幾何構(gòu)型如圖2 所示,a為裂紋長度,b為裂紋剩余韌帶長度.

圖2 試樣示意圖Fig.2 Schematic diagram of specimens

由于對稱性,圖3 給出了試樣的有限元分析的1/2模型,CT 試樣的寬度W為50 mm,SEB 試樣的寬度W為20 mm,其跨度H和寬度W之比H/W為2.裂紋尖端采用1/4 節(jié)點奇異單元(KSCON),其余部分均采用4 節(jié)點單元(Plane182),裂紋剩余韌帶部分施加對稱約束邊界條件.CT 試樣采用銷釘孔上半圓內(nèi)邊緣施加法向等效均布壓載荷(合力F即為試樣的載荷,F=0.25qW),考慮銷釘加載時加載孔X方向固定,對加載孔上方接觸點施加X向約束;SEB 試樣采用壓頭位移加載,壓頭(或加載輥)均采用Target169單元,構(gòu)元試樣表面均采用contact 172單元.

圖3 FEA 模型圖Fig.3 FEA model of specimens

2.2 FEA 模型驗證

為在網(wǎng)格模型設(shè)計中選用合適網(wǎng)格密度,需考察裂尖附近不同網(wǎng)格密度對 FEA 結(jié)果的影響,針對有限平面應變(plane strain)、平面應力(plane stress)條件的CT 試樣,有限元分析預設(shè)材料本構(gòu)關(guān)系為冪律,即式(2),預設(shè)K:1000 MPa,N:10,a/W: 0.5 完成冪律塑性分析,考察裂尖附近不同網(wǎng)格密度下的應力分布.定義圖3 中所示距裂尖5 mm扇形范圍內(nèi)單元尺寸為 1 倍網(wǎng)格密度(扇形單元徑向尺寸S為0.08 mm),2 倍、4 倍、8 倍網(wǎng)格密度定義為對裂尖扇形區(qū)網(wǎng)格加密為1 倍網(wǎng)格密度的相應倍數(shù)(扇形單元徑向尺寸S相應為0.04 mm,0.02 mm,0.01 mm).圖4 和圖5 中圖例的標記符號NE和NN分別表示不同密度下扇形網(wǎng)格區(qū)沿徑向單元總數(shù)和節(jié)點數(shù).圖4 和圖5 分別給出了有限平面應變和平面應力條件下不同網(wǎng)格密度CT 試樣 的FEA 獲取的載荷位移曲線和裂尖等效應力分布,結(jié)果表明,采用 1 倍網(wǎng)格密度時已經(jīng)滿足計算要求.

圖4 平面應變條件下裂尖附近網(wǎng)格密度對FEA 結(jié)果的影響Fig.4 Effects of element meshing size around the crack tip on the FEA results under plane strain conditions

圖5 平面應力條件下裂尖附近網(wǎng)格密度對FEA 結(jié)果的影響Fig.5 Effects of element meshing size around the crack tip on the FEA results under plane stress conditions

2.3 模型參數(shù)的確定

為了更加了解式(20)的特征,以下部分以平面應變條件下CT 試樣的應力分布為例,給出了參數(shù)確定方法.

2.3.1 綜合指數(shù)m和修正綜合指數(shù)msub

選取無量綱裂紋長度a/W為0.45~ 0.75(間隔0.05),應變硬化系數(shù)K為1000 MPa,應力硬化指數(shù)N分別為12,10,6 和4,在冪律塑性條件下進行計算.文獻[30]給出了CT 和SEB 試樣在平面應變條件下的體積折減系數(shù)m與參數(shù)k1~k4,進而通過式(15)得到k5,用同樣方法可標定出平面應力條件下體積折減系數(shù)m與等效應力系數(shù)k5,方程(15)參數(shù)如表1 所示.

表1 方程(15)模型參數(shù)Table 1 Parameters of Eq.(15)

進一步以應力因子σM為特征量,對不同無量綱裂紋長度a/W與應力硬化指數(shù)N下無量綱應力σsub/σM~r/L冪律擬合得到系數(shù)λsub,在給定應力硬化指數(shù)N時,將系數(shù)λsub與無量綱裂紋長度(1?a/W)冪律擬合得到msub,接著在不同N條件擬合出式(19)的待定參數(shù)msub0和msub1.圖6 為修正綜合指數(shù)參數(shù)標定過程,表2 給出了CT 和SEB 試樣分別在平面應變和平面應力條件下的待定參數(shù)msub0和msub1.在以下研究中,以σM作為特征應力,獲得裂尖的無量綱應力分布分析參數(shù)βsub和λsub.

圖6 參數(shù)標定Fig.6 Calibration parameter

表2 方程(19)參數(shù)Table 2 Parameters of Eq.(19)

2.3.2 奇異性指數(shù)fsub及其參數(shù)βsub

參考常用材料,選取無量綱裂紋長度a/W為0.6,應變硬化系數(shù)K為1000 MPa,應力硬化指數(shù)N分別為12,10,6 和4,在冪律塑性條件下進行計算.在距離裂紋尖端0.06b范圍內(nèi)的應力分布可用關(guān)于L/r的冪函數(shù)進行描述,其奇異性分布指數(shù)fsub可由回歸確定;考慮到應力的奇異性主要體現(xiàn)在裂紋尖端附近的右扇形區(qū),故主要關(guān)注應力在θ=0 到θ=π/2 的等效應力分布,在平面應變條件下得到各個角度的fsub和對應的βsub.參數(shù)βsub與θ關(guān)系如圖7 所示,由(17)可得三角函數(shù)關(guān)系如下

圖7 角度θ 對β 影響Fig.7 Effect of θ on β

式中,dsubp(p=0,1,2,3)為待定參數(shù).

類似可得到平面條件下βsub與θ的關(guān)系,表3給出了CT 和SEB 試樣分別在平面應變和平面應力條件下的參數(shù)dsubp(p=1,2,3,4)

表3 方程(23)模型參數(shù)Table 3 Parameters of Eq.(23)

2.3.3 無量綱應力分布系數(shù)λsub

為考慮各類因素對無量綱應力分布系數(shù)λsub的影響,以平面應變條件下試樣裂紋面(θ=0)的應力分布為例,將應力硬化指數(shù)N固定為10,以σM作為特征應力,獲得裂尖的無量綱應力分布,以此分析λsub.圖8 給出了載荷P和應力硬化系數(shù)K對λsub的影響,表明無量綱應力分布系數(shù)λsub與載荷P和應力硬化系數(shù)K無關(guān).

圖8 載荷P 和應力硬化系數(shù)K 對λsub 的影響Fig.8 Effects of loading P and strain hardening coefficient K on the λsub

將無量綱裂紋長度a/W和應變硬化系數(shù)K分別固定為0.6 和1000 MPa,以σM作為特征應力,獲得無量綱應力分布系數(shù)λsub如圖9 所示,在不同硬化指數(shù)下,λsub隨θ呈規(guī)律性變化;psub0與psub1均可以近似表示為與硬化指數(shù)N冪律相關(guān)的函數(shù),如圖10所示.類比三角特殊函數(shù)表征線彈性裂尖應力場[31],式(18)整理可得不同應力的λsub如下

圖9 角度θ 與N 對λsub 的影響Fig.9 Effects of θ and N on λsub

圖10 psubk(k=0,1)與N 的關(guān)系Fig.10 Relationships between psubk(k=0,1) and N

式中,psubk(k=0,1)與qsubm(m=0,1,2,3)為待定參數(shù),psubk(k=0,1)與N相關(guān),csubk(k=0,1),bsubk(k=0,1)為psubk(k=0,1)的待定參數(shù).

綜上,平面冪律塑性I 型裂紋尖端應力場半解析模型可顯式表達為

式中,參數(shù)csubk,bsubk(k=0,1) 與參數(shù)qsubm(k=0,1,2,3)列于表4.

表4 模型參數(shù)Table 4 Parameters of model

3 模型驗證

3.1 應力分布的驗證

圖8 表明應力分布系數(shù)λsub與應力硬化系數(shù)K無關(guān),而本文標定參數(shù)選用硬化指數(shù)N=4,6,10,12,故本文模型材料參數(shù)使用范圍為N=4~ 12,表1~表4 給出了相應的公式參數(shù)及其適用范圍,根據(jù)式(25)可預測CT 試樣和SEB 試樣的裂尖應力分布.圖11~ 圖14 分別給出了CT 試樣和SEB 試樣在平面應變、平面應力條件下的不同裂紋長度和不同硬化指數(shù)下,FEA 和本文公式預測的徑向應力分布對比;圖15~ 圖18 分別給出了CT 試樣和SEB 試樣在平面應變、平面應力條件下的不同裂紋長度和不同硬化指數(shù)下,裂紋尖端徑向距離不同情況下的FEA 和本文公式預測的環(huán)向應力分布對比.結(jié)果表明本文模型預測的結(jié)果與FEA 結(jié)果均吻合良好.

圖11 CT 試樣平面應變條件下裂尖無量綱應力徑向分布對比Fig.11 Comparison with FEA and predicted results by Eq.(25) for normalized stress radial distributions near the crack tip of CT specimens under plane strain conditions

圖12 CT 試樣平面應力條件下裂尖無量綱應力徑向分布對比Fig.12 Comparison with FEA and predicted results by Eq.(25) for normalized stress radial distributions near the crack tip of CT specimens under plane stress conditions

圖13 SEB 試樣平面應變條件下裂尖無量綱應力徑向分布對比Fig.13 Comparison with FEA and predicted results by Eq.(25) for normalized stress radial distributions near the crack tip of SEB specimens under plane strain conditions

圖14 SEB 試樣平面應力條件下裂尖無量綱應力徑向分布對比Fig.14 Comparison with FEA and predicted results by Eq.(25) for normalized stress radial distributions near the crack tip of SEB specimens under plane stress conditions

圖15 CT 試樣平面應變條件下裂尖無量綱應力環(huán)向分布對比Fig.15 Comparison with FEA and predicted results by Eq.(25) for normalized stress angular distributions near the crack tip of CT specimens under plane strain conditions

圖16 CT 試樣平面應力條件下裂尖無量綱應力環(huán)向分布對比Fig.16 Comparison with FEA and predicted results by Eq.(25) for normalized stress angular distributions near the crack tip of CT specimens under plane stress conditions

圖17 SEB 試樣平面應變條件下裂尖無量綱應力環(huán)向分布對比Fig.17 Comparison with FEA and predicted results by Eq.(25) for normalized stress angular distributions near the crack tip of SEB specimens under plane strain conditions

圖18 SEB 試樣平面應力條件下裂尖無量綱應力環(huán)向分布對比Fig.18 Comparison with FEA and predicted results by Eq.(25) for normalized stress angular distributions near the crack tip of SEB specimens under plane stress conditions

3.2 等效應力等值線驗證

根據(jù)式(25)可推導出平面I 型裂紋有限尺寸試樣等效應力等值線表達式

考慮極限載荷PL表達式為[32]

圖19、圖20 和圖21、圖22 分別給出了CT試樣和SEB 試樣加載到極限載荷PL時,不同裂紋長度和硬化指數(shù)下,有限元分析得到的等效應力等值線與方程(26)預測的等效應力等值線對比,當?shù)刃Ζ襡q達到屈服應力σ0時,等效應力等值線表示塑性區(qū)邊界.結(jié)果顯示,預測的結(jié)果與有限元分析結(jié)果吻合良好.

圖19 CT 試樣平面應變條件下等效應力σeq 等效應力等值線對比Fig.19 Comparison with FEA results predicted by Eq.(26) for contour lines of the equivalent stress σeq near the crack tip of CT specimen under plane strain conditions

圖20 SEB 試樣平面應變條件下等效應力σeq 等效應力等值線對比Fig.20 Comparison with FEA results predicted by Eq.(26) for contour lines of the equivalent stress σeq near the crack tip of SEB specimen under plane strain conditions

圖21 CT 試樣平面應力條件下等效應力σeq 等效應力等值線對比Fig.21 Comparison with FEA results predicted by Eq.(26) for contour lines of the equivalent stress σeq near the crack tip of CT specimen under plane stress conditions

圖22 SEB 試樣平面應力條件下等效應力σeq 等效應力等值線對比Fig.22 Comparison with FEA results predicted by Eq.(26) for contour lines of the equivalent stress σeq near the crack tip of SEB specimen under plane stress conditions

4 結(jié)論

(1) 基于能量密度等效和量綱分析方法,推導了受單向加載構(gòu)元試樣中值能量密度RVE 的等效應力解析方程,并定義其為σM應力因子;

(2) 以應力因子σM作為特征應力,采用可用于裂尖等效應力等值線表征的蝶翅輪廓式(平面應變)和扇貝輪廓式(平面應力)的三角特殊函數(shù),提出描述平面應變和平面應力冪律塑性I 型裂紋尖端應力場的半解析模型;

(3) 針對平面應變和平面應力CT 和SEB 試樣,本文模型預測的應力場與有限元分析結(jié)果均吻合良好.

猜你喜歡
裂紋效應
裂紋長度對焊接接頭裂紋擴展驅(qū)動力的影響
鈾對大型溞的急性毒性效應
一種基于微帶天線的金屬表面裂紋的檢測
懶馬效應
場景效應
Epidermal growth factor receptor rs17337023 polymorphism in hypertensive gestational diabetic women: A pilot study
微裂紋區(qū)對主裂紋擴展的影響
應變效應及其應用
偶像效應
預裂紋混凝土拉壓疲勞荷載下裂紋擴展速率
主站蜘蛛池模板: 国产在线观看91精品| 久久综合成人| 久久综合伊人 六十路| 国产中文在线亚洲精品官网| 女人av社区男人的天堂| 999精品在线视频| 亚洲二三区| 欧洲亚洲欧美国产日本高清| 欧美精品xx| 色综合久久综合网| 乱人伦中文视频在线观看免费| 亚洲精品视频免费| 成人国产小视频| 久久国产精品麻豆系列| 国产鲁鲁视频在线观看| 四虎综合网| 亚洲免费毛片| 伊人久久婷婷| a色毛片免费视频| 久久精品这里只有国产中文精品| 热思思久久免费视频| 精品久久香蕉国产线看观看gif| 日韩免费中文字幕| 午夜啪啪网| AⅤ色综合久久天堂AV色综合| 日韩久草视频| 制服丝袜亚洲| 人妻夜夜爽天天爽| 免费jjzz在在线播放国产| 人人爽人人爽人人片| 正在播放久久| 男人天堂亚洲天堂| av无码一区二区三区在线| 国产在线观看成人91| 亚洲一区二区精品无码久久久| 亚洲A∨无码精品午夜在线观看| 国产毛片久久国产| 国产精品理论片| 欧美亚洲日韩中文| 免费毛片全部不收费的| 71pao成人国产永久免费视频 | 91福利免费视频| 夜夜高潮夜夜爽国产伦精品| 欧美一级高清片久久99| 国产成人精品2021欧美日韩| 亚洲精品无码av中文字幕| 在线免费观看a视频| 无码精品国产VA在线观看DVD| 国产精品视频a| 真人免费一级毛片一区二区| 亚洲妓女综合网995久久| 国产波多野结衣中文在线播放| 成人小视频在线观看免费| 亚洲最大福利视频网| 国产女人在线| 亚洲永久色| 一区二区偷拍美女撒尿视频| 欧美 国产 人人视频| 一区二区三区精品视频在线观看| 久久天天躁夜夜躁狠狠| 国产精品内射视频| 成人年鲁鲁在线观看视频| 欧美啪啪网| 一级爱做片免费观看久久| 亚洲精品中文字幕无乱码| 久久久精品无码一区二区三区| 114级毛片免费观看| 2021精品国产自在现线看| 国产精品入口麻豆| 一区二区欧美日韩高清免费| 色精品视频| 青青草综合网| 国产一区二区三区在线观看视频| 国产成人精品免费av| 人妻中文字幕无码久久一区| 欧美国产在线看| 亚洲天堂伊人| 国产欧美综合在线观看第七页| 欧美α片免费观看| 国产精品毛片一区| 国产香蕉97碰碰视频VA碰碰看| 18禁色诱爆乳网站|