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

Lode相關斷裂準則在6061-T6511H鋁合金Taylor桿斷裂預報中的應用

2018-02-10 02:55:04楊慶年陳孝珍肖新科
振動與沖擊 2018年2期
關鍵詞:有限元模型

楊慶年, 陳孝珍, 肖新科, 李 凡, 張 偉

(1. 南陽理工學院 土木工程學院,河南 南陽 473004; 2. 重慶大學 土木工程學院,重慶 400045; 3. 哈爾濱工業大學 航天學院,哈爾濱 150008)

由于價格低廉、比強度高、承載能力良好和優異的加工性,金屬材料常用于制造戰斗部或構筑軍民防護結構[1]。在與作用對象碰撞過程中,金屬彈或靶往往發生大變形和多種形式和機理的斷裂破壞[2]。近年來,以有限元方法為代表的數值計算在沖擊工程領域獲得了諸多成功應用。

但是,金屬材料的延性斷裂是一種非常局部化的現象,裂紋萌生位置塑性變形大、應力高、應變梯度大。另外,金屬彈或靶還將經歷豐富和多變的應力狀態[3-7]。因此,預測金屬的延性斷裂就要求采用合理精度的本構模型和斷裂準則。正因為如此,Zukas等[8]指出,數值計算的精度和有效性受材料模型的限制。

近期,越來越多的試驗研究及基于細觀力學的數值計算表明,除應力三軸度外,Lode參數對金屬材料延性斷裂有影響[9-13],即斷裂應變同時與應力狀態的兩個參數有關[14]。但目前,一方面沖擊工程界大量采用僅應力三軸相關的斷裂準則(如Johnson-CooK[15])并取得不少成功案例,采用Lode參數相關斷裂準則的數值計算鮮見報道;另一方面,文獻中仍可見到數值計算無法合理預報試驗結果的案例,如文獻[16-19]。

因此,有必要研究應力三軸度和Lode參數同時相關的斷裂準則在沖擊斷裂問題數值預報中的效用。

金屬圓柱桿在撞擊剛性靶(即Taylor撞擊)的過程中可出現大變形和多種斷裂行為。本研究首先在一級輕氣炮上開展Taylor撞擊試驗,獲取寬泛撞擊速度范圍內彈體的變形和斷裂行為;然后,開展材料性能測試,標定本構模型和斷裂準則參數;最后,建立三維有限元模型開展數值Taylor撞擊試驗,得到Lode參數相關和不相關兩種斷裂準則的數值預報效果,通過與試驗結果的對比分析不同斷裂準則的預報效果,說明Lode參數相關斷裂準則的效用。

1 Taylor撞擊試驗及結果

1.1 試驗概況

試驗原材料為直徑35 mm的6061-T6511H鋁合金棒材,生產廠家為Sapa。該材料其他化學成分的質量分數分別為:w(Si)=0.76%,w(Fe)=0.45%,w(Cu)=0.30%,w(Mn)=0.06%,w(Mg)=0.88%,w(Cr)=0.07%,w(Zn)=0.03%,w(Ti)=0.02%。 該材料的基本材料參數為:ρ=2 700 kg/m3, 比熱Cp=890 J/(kg·K), 熔點Tm=925 K, 泊松比v=0.33。

彈體為平頭圓柱狀,名義直徑和長度分別為5.90 mm和29.48 mm。靶板為直徑250 mm、厚度25 mm的高強裝甲鋼,通過三個螺栓牢固的固定在靶架上。

試驗在一級輕氣炮上完成,撞擊過程采用Photron公司的FASTCAM SA5高速相機監控,初始速度由激光測速儀測得。試驗中,通過改變高壓氣室的初始壓力控制子彈的撞擊速度??偣查_展了6發試驗,得到的撞擊速度范圍是163.4~327.7 m/s。試驗中,靶板無明顯變形。

1.2 試驗結果

表1列出了試驗結果。其中,D0、L0和m依次表示子彈體初始的直徑、長度和質量;Df和Lf表示撞擊結束后彈體頭部的最終直徑和彈體最終的長度;V0表示彈體的初始撞擊速度。

表1 6061-T6511H鋁合金Taylor撞擊結果

從表1可見:

(1) 在撞擊速度在163.4~221.2 m/s時,6061-T6511H鋁合金彈體發生鐓粗變形,即彈體頭部直徑變大,長度變短,如圖1所示。

(2) 在撞擊速度在221.2~327.7 m/s時,彈體發生剪切開裂,即裂紋與軸線成45°夾角,如圖2所示。

圖1 6061-T6511H鋁合金Taylor桿的鐓粗Fig.1 Mushrooming of 6061-T6511H aluminum alloy Taylor rods

(a) V0=232.8 m/s

(b) V0=327.7 m/s圖2 6061-T6511H鋁合金Taylor桿的剪切開裂Fig.2 Shear cracking of 6061-T6511H aluminum alloy Taylor rods

圖3給出了其中一發試驗的高速相機圖像,其中(a)和(b)顯示的是著靶前,圖(c)顯示的是著靶中,(d)~(e)顯示的是著靶后的反彈過程??梢娮訌椩谧矒羟昂妥矒艉缶哂辛己玫娘w行姿態。

圖3 T3試驗T3的高速攝像圖像Fig.3 Pictures recorded by high speed camera for Test T3

2 材料性能測試和表征

2.1 材料模型

本構模型擬采用肖新科提出的修正形式的Johnson-Cook模型[20],寫為

(1)

與原始Johnson-Cook模型相比, 公式(1)所示的MJC模型多一個參數p,該參數的增加使得該模型可以表征材料屈服強度隨溫度增加而發生非線性軟化的現象。同時使得原始J-C模型變為MJC模型的一種特殊情況,即p=1。

斷裂準則擬采用文獻[21]提出的修正形式的Johnson-Cook準則(MJC準則)以及文獻[22]提出的L-Y-H準則。后者相對后者,不但包含了應力三軸度的影響,也包含了Lode參數的影響。

MJC準則寫為

(2)

式中:D1~D6為模型參數。η為應力三軸度,定義為

η=(σ1+σ2+σ3)/(3σeq)

(3)

式中:σ1~σ3依次為第一,中間和第三主應力。

原始L-Y-H準則[22]不包含溫度和應變率項。為了應用在沖擊問題中,采用式(2)所示的公式描述溫度和應變率的影響,最終L-Y-H準則寫為

(4)

式中:l1~l4及D4~D6是模型參數,L為Lode參數,定義為

(5)

另外,采用絕熱假定,溫升ΔT寫為

(6)

式中:χ為塑性功轉熱系數, 取為0.9;ρ為材料密度;CP為比熱。

2.2 模型參數標定

為了標定上述模型的相關參數,開展了如下四組試驗:

(a) 光滑圓棒試樣的拉伸試驗。該試驗用于標定MJC本構模型中的參數A,B和n。

(b) 不同缺口半徑的拉伸試驗。由于不同缺口半徑的光滑圓棒試樣(包括光滑圓棒試樣,缺口半徑為∞)的Lode參數L=-1,但應力三軸度不同[23]。因此,本組試驗與(a)試驗一起可標定MJC準則的參數D1~D3。

(c) 光滑圓棒試樣25~500 ℃的拉伸試驗。標定MJC模型的參數p和m和MJC準則及L-Y-H準則的參數D4~D6。

(d) 光滑圓棒試樣的扭轉試驗。該試驗中Lode參數為L=0[23]。因此,試驗(d)與試驗(a)和(b)一起可以標定L-Y-H準則的參數l1~l4。

由于篇幅所限,模型參數的標定方法主要參考文獻[24],這里僅列出要點。

光滑圓棒和缺口試樣缺口處的名義直徑為6 mm,缺口半徑R有三組:R=2 mm、3 mm、9 mm。使用20 mm長引伸計記錄含缺口段試樣在軸向的伸長量,試驗得到的載荷位移曲線如圖4所示。

圖4 光滑圓棒和缺口圓棒試樣的載荷位移曲線Fig.4 Load-displacement curves of smooth round and notched round bars

圖5給出了不同溫度下光滑圓棒試樣的載荷位移曲線,除常溫外,其他溫度下的位移為十字頭位移。試樣直徑6 mm,標距段長度為40 mm。圖6給出了不同溫度下材料的屈服強度和失效應變。失效應變定義為εf=ln(A0/Af), 其中A0和Af分別為試樣的初始橫截面面積和拉斷后的端口面積。需要說明的是,試驗中除500 ℃外,其他拉伸試驗均為兩組平行試驗,為了便于顯示,圖5中僅畫出其中一組試驗的結果。

圖5 不同溫度下光滑圓棒試樣的載荷位移曲線Fig.5 Load-displacement curves of smooth round bars at various temperatures

圖6 屈服強度和斷裂應變隨溫度的變化Fig.6 Variations of yield stress and fracture strain versus temperature

圖7給出了扭轉試樣的幾何形狀、尺寸以及在MTS 809試驗機上獲得的扭矩-轉角曲線。

圖7 光滑圓棒試樣的扭矩-轉角曲線Fig.7 Torque-twist curve of smooth round bar

利用試驗(a)標定A,B和n。易知,A為屈服強度。從圖4知,試樣發生了頸縮,且頸縮后仍發生了較大的伸長變形。頸縮后變形將集中在頸縮區域,頸縮處的應力狀態將由單向應力狀態向多軸應力狀態轉變,此后的試件單向真應力應變與試件的等效應力應變將不存在對等關系。對于考慮頸縮后的等效應力應變關系,即參數B和n的獲取,可以利用有限元迭代法。計算時,不斷調整參數B和n的取值,直至有限元計算得到的載荷位移曲線與試驗結果相吻合為止,具體過程見肖新科的研究。得到的模型參數見表2,預測的載荷位移曲線見圖4。

表2 6061-T6511H鋁合金的材料參數

由圖6所示的屈服強度隨溫度的變化可擬合出參數p和m,擬合效果見圖6。

至此,除參數C外,MJC本構模型的全部參數已獲得。參數C采用反向計算的辦法得到,即建立Taylor撞擊試驗的二維軸對稱計算模型,迭代C,直至獲得的Taylor桿的最終幾何尺寸Df和Lf與試驗吻合為止。計算中,以發生鐓粗變形的3發試驗為依據。除彈體單元尺寸不同外,計算的有限元模型同肖新科的研究結果。本研究中彈體的單元尺寸為0.1 mm×0.1 mm。得到的C=0.049,模擬效果如表3所示。

缺口圓棒的斷口直徑不宜測量,因此采用有限元計算的方法獲取斷裂應變:在有限元計算中獲得對稱軸上等效塑性應變最大單元的應力狀態和等效塑性應變歷程。

由于無合適理論公式計算實心圓柱扭轉試樣的等效塑性應變,也用數值計算的方法獲取失效應變。計算思路同缺口試樣,所不同的是扭轉試樣的破壞由試樣表層開始,因此斷裂應變取自試樣表面等效塑性應變最大的單元。

表3 鐓粗Taylor桿的試驗與數值計算結果對比

還需要指出的是,在上述缺口圓棒和扭轉試樣的有限元計算中,首先要適當調整輸入本構模型參數A,B和n,直至獲得的載荷-位移曲線與試驗一致為止。這種方法在獲取材料斷裂應變中經常采用,如肖新科、Gilioli等研究結果。需要調整本構模型參數是因為材料的塑性流動行為受應力狀態及塑性變形大小的影響[26]。

上述圓棒拉伸和扭轉試驗全部建立二維軸對稱模型,標距段內單元尺寸不超過0.15 mm×0.2 mm。有限元網格剖分情況如圖8。獲得的斷裂應變見表4和圖9。可見,除應力三軸度外,Lode參數對材料的斷裂應變也有影響。

(a) 圓棒試樣拉伸

(b) 圓棒試樣扭轉圖8 準靜態加載時各試樣的有限元計算模型Fig.8 Finite element models of various specimens under quasi-static loading

表4 各應力狀態下的斷裂應變

J-C斷裂準則經常通過拉伸試驗標定,因此標定J-C準則時使用光滑圓棒和缺口圓棒的數據,標定結果如圖9。標定L-Y-H準則時,同時還采用扭轉試驗的結果,標定結果如圖9所示。從圖9可知:

(1) MJC和L-Y-H準則都可以對圓棒拉伸試驗的斷裂應變給出合理預報。

(2) 由于沒有考慮Lode參數的影響,MJC準則不能預報材料斷裂應變隨Lode參數的變化,特別是剪切狀態時材料斷裂應變遠低于拉伸試驗的特性無法得到合理預報。

(3) L-Y-H準則可很好地預報所研究試驗狀態下材料的斷裂應變。

由上述結論可知,若材料的塑性變形接近剪切,即L=0,采用MJC準則時必然會高估材料的延性,從而低估材料的損傷。

(a) 二維

(b) 三維及L-Y-H準則預測結果圖9 6061-T6511H鋁合金的斷裂應變及預報效果Fig.9 Fracture data and the predictions of fracture criteria

根據圖6中所列不同溫度下材料的斷裂應變,可擬合出參數D5和D6,擬合效果如圖6所示。需注意的是圖6中的應變為材料失效時的近似主應變,本文假定失效時的等效應變與該主應變相同。

依據文獻[27],應變率對6061-T6斷裂應變的影響可以忽略。另外,Gilioli在研究6061-T6鋁合金的抗侵徹性能時也忽略了應變率對斷裂應變的影響。在無試驗數據的情況下,本文認為D4=0。

至此,獲得了MJC本構模型和MJC斷裂準則及L-Y-H準則的全部模型參數,列于表2中。

2.3 靶板的材料模型及參數

由于試驗中靶板的變形不明顯,借用文獻[28]中的雙線性硬化模型來描述其本構關系,并忽略應變率對其強度的影響,不考慮其在撞擊過程中的溫升,即:

(7)

式中:E和Et分別為材料的彈性模型和切線模量;σ0為材料的屈服強度;ε0為材料發生初始屈服時所對應的應變。模型參數除密度外同文獻[28],如表5所示。

表5 裝甲鋼的材料參數

3 有限元計算及結果

3.1 有限元模型

有限元計算在Abaqus/Explicit中進行,Taylor桿和靶板均采用三維實體,如圖10所示。初始時刻Taylor桿與靶板有0.5 mm間隙,同時為了與實驗條件更為接近,對靶板周圍20 mm范圍進行完全固定。

Taylor桿與靶板均采用C3D8R單元。為提高計算效率,Taylor桿前端區域單元劃分較密,大小約為0.13 mm×0.13 mm×0.15 mm;遠離靶板部分的單元稍微稀疏一些,其大小變至約0.13 mm×0.13 mm×0.4 mm。靶板的網格尺寸相對較大,與Xiao的結論相同。

圖10 Taylor桿撞擊試驗有限元模型Fig.10 FE calculation model of the Taylor impact test

3.2 計算結果

依據試驗的速度范圍,分別采用MJC模型和MJC斷裂準則或L-Y-H準則,通過改變子彈的初始撞擊速度開展了數值Taylor撞擊試驗,記錄了彈體打靶后的變形和斷裂模式。圖11和12給出了對應于試驗中最高撞擊速度的三發試驗的計算結果,圖中SDV1和SDV5分別為等效塑性應變和損傷。從圖中可以看出:

(1) 采用MJC斷裂準則時,彈體全部發生鐓粗變形,無斷裂發生。而試驗中僅V0=221.2 m/s時,彈體發生鐓粗變形,更高速度下彈體均發生剪切破壞。

(2) 采用L-Y-H斷裂準則時,彈體都發生了剪切開裂。與試驗結果相比,預測的裂紋形式合理(見圖2),即主要裂紋與彈體軸線方向呈45°角,但稍微過高預報了彈體的損傷程度,如V0=221.2 m/s時預報到了彈體頭部邊緣發生開裂,V0=327.7 m/s時預報了較試驗更長的裂紋長度。

兩者預測結果的差別可由Taylor桿頭部材料的應力狀態來進行解釋。圖13給出了V0=232.8 m/s時彈體頭部邊緣一個失效單元的損傷累積及應力狀態歷程,采用的斷裂準則為L-Y-H準則預測??梢姡瑩p傷主要在2.5~7.5 μs內造成,該時間內應力三軸度的水平介于0和0.18之間,而羅德參數處于0與-0.5之間(如圖所示)。依據圖9,在該應力狀態下,L-Y-H斷裂準則預報的斷裂應變值遠小于MJC的預測結果。因此采用L-Y-H斷裂準則預報的損傷值一定高于MJC的預測結果,預報的破壞程度也更強。

圖11 采用MJC斷裂準則時數值模擬預報的彈體變形與斷裂Fig.11 Numerical predicted deformation and fracture pattern by using MJC fracture criterion

圖12 采用L-Y-H斷裂準則時數值模擬預報的彈體變形與斷裂Fig.12 Numerical predicted deformation and fracture pattern by using MJC fracture criterion

圖13 彈體頭部邊緣剪切裂紋上失效單元的損傷及應力狀態歷程Fig.13 Damage and stress state history of an element located in the shear crack of projectile nose rim

4 結 論

在一級輕氣炮上開展了直徑5.9 mm的6061-T6511H鋁合金圓柱桿撞擊剛性靶的Taylor試驗,得到了鐓粗和剪切開裂兩種變形和斷裂模式;通過開展材料性能測試標定了本構模型和斷裂準則;最后通過有限元計算開展了數值打靶試驗,得到了Lode參數無關和相關斷裂準則的預報效果,得到如下結論:

(1) 6061-T6511H鋁合金的斷裂行為與Lode參數相關。

(2) 采用修改形式的Johnson-Cook準則時,由于沒考慮Lode參數對材料斷裂應變的影響,有限元計算過高估計了材料的延性,低估了Taylor桿的斷裂。

(3) Lode參數相關的L-Y-H斷裂準則預報到了與試驗一致的裂紋形式,但略微過高估計了Taylor桿的斷裂。

本研究沒能通過試驗手段得到應變率對材料斷裂應變的影響,也忽略了應力狀態對材料本構關系的影響。后期將開展這兩方面的研究,以期更加準確地評估Lode參數相關斷裂準則的有效性。

[ 1 ] DEY S. High-strength steel plates subjected to projectile impact-An experimental and numerical study [D]. Trondheim: Norwegian University of Science and Technology, 2004.

[ 2 ] 肖新科. 雙層金屬靶的抗侵徹性能和Taylor桿的變形與斷裂[D]. 哈爾濱:哈爾濱工業大學,2010.

[ 3 ] TENG X, WIERZBICKI T, HIERMZIER S, et al. Numerical Prediction of Fracture in the Taylor Test[J]. International Journal of Solids and Structures, 2005, 42: 2929-2948.

[ 4 ] TENG X, WIERZBICKI T. Evaluation of six fracture models in high velocity perforation [J]. Engineering Fracture Mechanics, 2006, 73: 1653-1678.

[ 5 ] XIAO X, ZHANG W, WEI G, et al. Experimental and numerical investigation on the deformation and failure behavior in the Taylor test[J]. Materials and Design, 2011, 32: 2663-2674.

[ 6 ] GILIOLI A, MANES A, GIGLIO M, et al. Predicting ballistic impact failure of aluminium 6061-T6 with the rate-independent Bao-Wierzbicki fracture mode[J]. International Journal of Impact Engineering, 2015, 76: 207-220.

[ 8 ] ZUKAS J A, NICHOLAS T, SWIFT H F, et al. Impact dynamics [M]. Wiley, 1982.

[ 9 ] WIERZBICKI T, BAO Y, LEE Y W, et al. Calibration and Evaluation of Seven Fracture Models[J]. International Journal of Mechanical Sciences, 2005, 47: 719-743.

[10] BARSOUM I, FALESKOG J. Rupture mechanisms in combined tension and shear-experiments [J]. International Journal of Solids and Structures, 2007, 44: 1768-1786.

[11] GAO X, ZHANG T, HAYDEN M, et al. Effects of the stress state on plasticity and ductile failure of an aluminum 5083 alloy[J]. International Journal of Plasticity, 2009, 25: 2366-2382.

[12] LIAN J, SHARAF M, ARCHIE F, et al. A hybrid approach for modelling of plasticity and failure behaviour of advanced high-strength steel sheets [J]. International Journal of Damage Mechanics, 2012, 22: 188-218.

[13] ERICE B, PéREZ-MARTN M J, GLVEZ F. An experimental and numerical study of ductile failure under quasi-static and impact loadings of Inconel 718 nickel-base superalloy[J]. International Journal of Impact Engineering, 2014, 69: 11-24.

[14] BAI Y, WIERZBICKI T. A comparative study of three groups of ductile fracture loci in the 3D space[J]. Engineering Fracture Mechanics, 2015, 135: 147-167.

[15] JOHNSON G R, COOK W H. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures [J]. Engineering Fracture Mechanics, 1985, 21: 31-48.

[16] ANDERSON Jr C E, CHOCRON I S, NICHOLLS A E. Damage modeling for taylor impact simulations [J]. Journal De Physique IV France, 2006, 134: 331-337.

[17] 陳剛. 半穿甲戰斗部彈體穿甲效應數值模擬與實驗研究[D]. 中國工程物理研究院, 2006.

[18] DEY S, B?RVIK T, HOPPERSTAD O S, et al. The effect of target strength on the perforation of steel plates using three different projectile nose shapes [J]. International Journal of Impact Engineering, 2004, 30: 1005-1038.

[19] KANE A, B?RVIK T, HOPPERSTAD O S, et al. Finite element analysis of plugging failure in steel plates struck by blunt projectiles [J]. Journal of Applied Mechanics, 2009, 76: 051302-1-11.

[20] JOHNSON G R, COOK W H. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures [C]//Proceeding of the Seventh International Symposium on Ballistic. The Netherlands: The Hague, 1983: 541-547.

[21] 郭子濤. 彈體入水特性及不同介質中金屬靶的抗侵徹性能研究[D]. 哈爾濱: 哈爾濱工業大學, 2012.

[22] LOU Y, YOON J W, HUH H. Modeling of shear ductile fracture considering a changeable cut-off value for stress triaxiality[J]. International Journal of Plasticity, 2014, 54: 56-80.

[23] BAI Y, WIERZBICKI T. A new model of metal plasticity and fracture with pressure and Lode dependence[J]. International Journal of Plasticity, 2008, 24: 1071-1096.

[24] BAO Y, WIERZBICKI T. A comparative study on various ductile crack formation criteria [J]. Journal of Engineering Materials and Technology, 2004, 126: 314-324.

[25] BAO Y, WIERZBICKI T. On fracture locus in the equivalent strain and stress triaxiality space[J]. International Journal of Mechanical Sciences, 2004, 46(1): 81-98.

[26] ALGARNI M, BAI Y, CHOI Y. A study of inconel 718 dependency on stress triaxiality and lode angle in plastic deformation and ductile fracture [J]. Engineering Fracture Mechanics, 2015, 147: 140-157.

[27] LESUER D R, KAY G J, LEBLANC M M. Modeling large-strain, high rate deformation in metals [C]. Third Biennial Tri-Laboratory Engineering Conference on Modeling and Simulation, Pleasanton, CA. 2001.

[28] B?RVIK T, HOPPERSTAD O S, BERSTAD T, et al. A computational model of viscoplasticity and ductile damage for impact and penetration [J]. European Journal of Mechanics-A/Solids, 2001, 20: 685-712.

猜你喜歡
有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
新型有機玻璃在站臺門的應用及有限元分析
上海節能(2020年3期)2020-04-13 13:16:16
基于有限元的深孔鏜削仿真及分析
基于有限元模型對踝模擬扭傷機制的探討
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 99热亚洲精品6码| 欧美午夜在线播放| 视频在线观看一区二区| 国产麻豆永久视频| 色综合综合网| 激情网址在线观看| 欧洲高清无码在线| 91在线精品麻豆欧美在线| 国产91丝袜在线播放动漫 | 国产精品成人一区二区不卡| 亚洲精品成人福利在线电影| 狠狠ⅴ日韩v欧美v天堂| 国产激情无码一区二区免费| 91丝袜美腿高跟国产极品老师| 在线观看av永久| 国产偷国产偷在线高清| 欧美日韩一区二区在线播放| 国内精品伊人久久久久7777人| 色婷婷成人网| 久久久久青草大香线综合精品 | 91久久偷偷做嫩草影院免费看 | 九一九色国产| 亚洲天堂网在线播放| 欧美不卡在线视频| 亚洲国产精品VA在线看黑人| h网址在线观看| 日韩人妻无码制服丝袜视频| 欧美成人a∨视频免费观看 | 好紧太爽了视频免费无码| 久久 午夜福利 张柏芝| 国产乱肥老妇精品视频| 国产欧美中文字幕| 国产91小视频在线观看| 欧美在线免费| 久久久久免费看成人影片| 奇米影视狠狠精品7777| 免费观看欧美性一级| 亚洲性一区| 精品欧美视频| yy6080理论大片一级久久| 有专无码视频| 亚洲AV无码乱码在线观看裸奔| 国产一二三区在线| 欧美色视频网站| 在线观看精品自拍视频| 日韩黄色精品| 成人在线亚洲| 综合网天天| 多人乱p欧美在线观看| 在线精品亚洲一区二区古装| 亚洲bt欧美bt精品| 欧美日韩中文字幕在线| 国产区成人精品视频| 欧洲一区二区三区无码| 欧美成人区| 欧美日韩专区| 97久久人人超碰国产精品| 国产一级视频久久| 久久香蕉国产线看观看式| 婷婷综合色| 天堂中文在线资源| 国产精品入口麻豆| 黄网站欧美内射| 内射人妻无套中出无码| 亚洲av日韩av制服丝袜| 91毛片网| 91av成人日本不卡三区| 91丝袜在线观看| 热伊人99re久久精品最新地| 日本91视频| 久久九九热视频| 草草影院国产第一页| 日本91视频| 黄色在线不卡| 国产成人调教在线视频| 伊人久综合| 美女免费精品高清毛片在线视| 国产成人综合在线视频| 成人字幕网视频在线观看| 久久semm亚洲国产| 老司国产精品视频| 婷婷色中文|