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

CoCrFeMnNi 高熵合金沖擊波響應與層裂強度的分子動力學研究1)

2022-08-30 02:41:20袁福平熊啟林闞前華
力學學報 2022年8期
關鍵詞:方向

杜 欣 袁福平 熊啟林 張 波 闞前華 張 旭 ,2)

* (西南交通大學力學與航空航天學院,成都 610031)

? (中國科學院力學研究所非線性力學國家重點實驗室,北京 100190)

** (華中科技大學航空航天學院,武漢 430074)

引言

高熵合金作為新型結構材料,其擁有高硬度、高強度、良好的塑性、抗輻照性能以及出色的熱穩定性和抗腐蝕能力[1-5],這為高熵合金未來在航空航天、太空探索以及深海探測等領域的應用提供了巨大的應用前景.作為一種經典的高熵合金[6],CoCrFeMnNi 高熵合金在極端環境下的力學行為研究備受關注.然而,現有的表征手段無法觀測到動態連續的微結構演變,使得高熵合金在納米尺度的層裂行為缺少深入的理解.分子動力學(molecular dynamics,MD) 方法具有時間尺度短,空間尺寸小,應變率大的特點[7],可以在極短時間尺度下研究材料沖擊動態力學響應,且容易捕捉受沖擊過程中的微結構演化,例如孔洞形核過程[8].因此,利用納米空間尺度和皮秒時間尺度的MD 方法分析材料的動態力學行為具有天然優勢.

目前,MD 模擬方法已成為研究材料沖擊動態力學行為的重要手段之一,并已有大量的研究報道.Li 等[9]利用MD 研究了缺陷可逆性對碳化硅動態拉伸強度的影響.結果表明,在低于臨界壓縮應變的情況下,SiC 在沿[001]晶向加載時產生可逆的變形孿晶,這種可逆性使得晶體擁有高強度.在超過臨界應變時,該過程變得不穩定,并觸發孿晶界運動,產生不可逆的纏結缺陷,導致強度顯著降低.Jian等[10]通過MD 模擬研究了沖擊引起CoCrNi 中熵合金非晶化的閾值和原因,并與采用平均原子勢的模擬結果進行對比,發現中熵合金由于晶格畸變效應而產生了更多的無序結構,造成其層裂強度降低.Cekil 等[11]研究了碳化硼在動態載荷作用下的力學性能,在沿[0001]和兩個方向分別以2.0 和3.0 km/s 的速度進行沖擊時碳化硼發生了非晶化,C-B-C 鏈的屈曲以及晶格旋轉是發生非晶化的主要原因,其中方向不發生晶格旋轉并且其C-B-C鏈發生屈曲的臨界沖擊速度大于[0001]方向.Zhu等[12]通過MD 研究了Cu/Ni 納米多層體系在不同沖擊脈沖持續時間下的動態響應和層裂,模擬結果發現層裂只出現在銅層區域,而不會出現在鎳層區域;并且隨著沖擊持續時間的增加,層裂機制從空隙在銅層內均勻成核轉變為在Cu/Ni 界面處成核,導致整體層裂強度下降.Xie 等[13]使用蒙特卡洛以及MD 方法研究了單晶和多晶CoCrNi 在沖擊加載下的力學行為、缺陷演化和層裂破壞特征.結果表明單晶中熵合金的Hugoniot 彈性極限表現出異常的各向異性,并且孔洞成核于層錯相交處或者沖擊壓縮產生的非晶態區域.在多晶中,孔洞形核則位于晶界處的Ni 聚集區域,一般為三叉晶界處,并隨著晶界分離逐漸擴大.Thürmera 等[14]利用實驗和MD 模擬研究了CoCrFeMnNi 高熵合金沖擊波作用下層裂現象,發現層裂面產生了納米孿晶,MD 模擬中證實CoCrFeMnNi 高熵合金擁有較Cu,Fe 和Ta 更高的層裂強度.

以上研究雖已對傳統合金,中熵合金和高熵合金的沖擊動態損傷進行了報道,但現有文獻中尚未考慮彈塑性雙波分離對層裂強度的影響.此外,在CoCrFeMnNi 高熵合金中沿不同方向沖擊時的微結構演化以及微結構對層裂強度的影響仍需進一步分析.利用MD 模擬軟件Lammps[15],研究沖擊方向以及速度對單晶CoCrFeMnNi 高熵合金沖擊波響應以及層裂行為的影響,并利用后處理軟件OVITO 分析二者對沖擊時微觀結構演化的影響[16-17].相關研究成果有助于從原子尺度上理解CoCrFeMnNi 高熵合金的動態力學行為,促進CoCrFeMnNi 高熵合金未來作為抗沖擊材料在航空航天等領域內的應用.

1 原子模型與加載方式

為研究沖擊方向和沖擊速度對層裂強度以及沖擊過程中微觀結構演化的影響,利用Atomsk 軟件[18]構建了取向分別為的三個單晶CoCrFeMnNi 高熵合金模型,三個模型的尺寸均為74.7 nm × 13.3 nm ×7.5 nm,晶格常數為a=3.595 ?.其中x方向為[100]的單晶CoCrFeMnNi 高熵合金模型如圖1 所示.在沖擊之前對模型進行弛豫,使系統達到平衡狀態.在弛豫過程中,模型的x,y和z方向均采用周期性邊界條件,并在NPT 系綜下弛豫50 ps (1 ps=1×10-12s).在誘導沖擊波產生和沖擊波傳播階段采用NVE 系綜,并且沖擊方向(x方向)為自由邊界條件.采用活塞法誘導沖擊波的產生[19],即在模型的一端設置虛擬墻作為活塞并給予一定的速度進行沖擊.在本文中,為探究沖擊速度對層裂強度的影響,施加的沖擊速度Up為0.6,0.9,1.2,1.5 和1.8 km/s.虛擬墻作為活塞進行沖擊的過程持續5 ps.之后,在NVE 系綜下沖擊波在模型內傳播并在自由面發生反射形成拉應力區,該過程持續時間為35 ps.在整個模擬過程中,時間步長為1 fs (1 fs=1×10-15s),采用的勢函數為Choi等[20]提出的第二近鄰修正嵌入原子(2 NN MEAM) 勢.該勢函數已被用來模擬CoCrFeMnNi 高熵合金的沖擊動態力學行為[14,19].此外,該勢函數也被用來模擬CoCrFeMnNi 高熵合金的拉伸、循環、納米壓痕以及輻照行為[21-24].因此,該勢函數具有一定可靠性.

圖1 x 方向為[100]取向的單晶CoCrFeMnNi 高熵合金原子模型.黃色和紅色平面分別代表虛擬墻和自由面Fig.1 Atomistic model of single-crystal CoCrFeMnNi high-entropy alloy with [100] orientation in x direction.The yellow and red planes represent virtual wall and free surface,respectively

2 結果與分析

2.1 自由面速度

自由面速度可以反映層裂強度和動態損傷的發生.當沖擊波到達樣品自由面時,自由面速度增大到最大值.隨后,沖擊波在自由面發生反射,反射后的稀疏波在與衰減沖擊波尾部的稀疏波相遇后,自由面發生卸載,造成自由面速度下降并且形成拉應力區域.當拉應力達到閾值后發生層裂,而反射后的稀疏波作用于層裂面時會發生再次反射,使得自由面速度增加,產生速度回跳[10,25-26].

沿[100],[110]和[111]方向進行沖擊加載時的自由面速度演化如圖2 所示.可以看出,在以0.6 km/s 的速度沿[100]方向進行沖擊時自由面速度未發生回跳,表明在該速度下沿[100]方向沖擊未出現動態損傷[27].但是,沿[110] 和[111] 方向以0.6 km/s 的速度沖擊時自由面速度發生了多次回跳,表明沿[110]和[111]方向沖擊時動態損傷發生的臨界沖擊速度小于[100]方向.

圖2 自由面速度演化Fig.2 Evolutions of free surface velocity

在虛擬墻以0.6 km/s 的速度沿三個方向沖擊5 ps 后,樣品內部的塑性變形未隨著沖擊波的傳播而增加,表明此時只產生了彈性波,未產生塑性波.但是,由于在CoCrFeMnNi 高熵合金中拉伸強度低于壓縮強度,表現出了明顯的拉-壓非對稱性[22],使得在[110]和[111]方向下沖擊波反射后形成的拉應力區產生了塑性變形,如圖3 所示.其中,綠色原子為面心立方(face-centered cubic,FCC)結構;紅色原子為密排六方(hexagonal close-packed,HCP)結構,并且主要為層錯;藍色原子為體心立方(bodycentered cubic,BCC)結構;白色原子為無序結構[28].采用了共鄰分析(common neighbor analysis,CNA)方法[29]分析.可以看出,在虛擬墻沖擊5 ps 時,由于模型沖擊表面原子受到虛擬墻的作用發生畸變,使得沖擊表面附近產生了少量層錯.在沖擊波正向傳播的過程中(10,15 ps),塑性變形并未隨著沖擊波的傳播而增加.但是,當沖擊波在自由面發生反射后(20 ps),由于拉-壓非對稱性,造成了在拉應力區發生塑性變形,產生了層錯.此時,[110]方向在當前的二維視角中未產生貫穿模型的孔洞.但隨著沖擊波的進一步作用,在25 ps 時[110]方向出現了貫穿模型的孔洞.

圖3 采用CNA 方法分析以0.6 km/s 的沖擊速度沿[110]和[111]方向沖擊時的微觀組織演化Fig.3 Microstructure evolutions by CNA at the 0.6 km/s shock velocity along the [110] and [111] directions

此外,從圖2(b)和圖2(c)中可以看出,隨著沖擊速度的增加,在沿[110]和[111]方向進行沖擊時,自由面速度達到峰值前經歷了兩次增加的過程.其中,第一次增加是由于彈性波作為先驅波首先到達自由面,第二次增加是由于塑性波到達自由面.由于沿[100]方向未發生彈塑性雙波分離現象,所以其自由面速度在達到峰值前只經歷了一次增加的過程.[100] 方向未出現雙波分離現象與該方向擁有較[110]和[111]方向更多可開動滑移系有關,多個滑移系的開動造成彈性波波速較低,甚至低于塑性波波速[30].當沖擊速度為0.9 km/s 時,在沿[110]和[111]方向沖擊時樣品內部產生了塑性波,但此時在[111]方向的自由面速度中未體現出彈塑性雙波分離,并且[110]方向彈塑性雙波分離現象不明顯,只對自由面速度造成微小擾動.此外,在以1.2 km/s 的沖擊速度沿[111]方向沖擊時,自由面速度曲線反映出的彈塑性雙波分離現象同樣不明顯.隨著沖擊速度的進一步增加,在以1.2 km/s 和1.5 km/s 分別沿[110]和[111]方向沖擊時,自由面速度在達到峰值前明顯經歷了兩次增加的過程,表明了彈塑性雙波分離現象的出現.由于沿[111]方向加載時Schmid因子小于[110]方向(分別為0.272 和0.408),所以沿[111]方向進行沖擊時塑性變形所需的外加應力更大[19],導致了[111]較[110]方向擁有更大的彈塑性雙波分離現象臨界沖擊速度.之后,彈塑性雙波分離現象隨著沖擊速度的增加而減弱.

為進一步分析彈塑性雙波結構,得到沿[100],[110]和[111]方向以0.6~1.8 km/s 沖擊5 ps 時的粒子速度剖面圖,如圖4 所示.可以看出,在[110]和[111]方向以0.6 km/s 的速度沖擊時曲線并未出現臺階段,表明此時未出現彈塑性雙波分離現象.在[110]和[111]方向以0.9 km/s 的速度沖擊時,曲線中出現了臺階段,但是臺階段的粒子速度與沖擊表面粒子速度相差不大,使得在圖3 中[110]方向的自由面速度未表現出明顯的彈塑性雙波分離現象,而[111]方向未出現彈塑性雙波分離現象.在1.2 km/s的沖擊速度下,[110]和[111]方向出現了明顯的臺階段,但是由于彈性波的單獨擾動,使原子面在波前附近處彈性振動[31],導致[111]方向彈性波前附近粒子速度上升并與沖擊表面附近粒子速度相差仍然不大,造成圖3 中[111]方向的自由面速度未表現出明顯的彈塑性雙波分離現象.由于臺階段粒子速度與沖擊表面速度相差不大,而此時已發生雙波分離,使得自由面受到了分離的彈性波和塑性波的兩次作用,造成了在圖3 中自由面速度第一次回跳點滯后.隨著沖擊速度的進一步增加,塑性變形加劇,使塑性波傳播速度增加,造成塑性波與彈性先驅波之間的速度差減小,從而使得臺階段的寬度逐漸減小,表明在發生明顯的彈塑性雙波分離現象后,彈塑性雙波分離現象隨著沖擊速度的增加而減弱.在[111]方向發生雙波分離現象后(Up≥1.2 km/s),[111]方向先驅波的波前附近粒子速度明顯大于[110]方向.此外,在[100]方向下以0.6~1.8 km/s 速度沖擊時,曲線中均未產生臺階段,表明[100]方向未出現雙波分離現象.

圖4 沖擊5 ps 時粒子速度剖面圖Fig.4 Stress profile at 5 ps

整體而言,隨著沖擊速度的增加,[110] 和[111]方向的彈塑性雙波分離現象先增強后減弱,自由面速度多次回跳現象減弱,自由面峰值速度增加.此外,沿[100]方向進行沖擊時產生動態損傷的臨界沖擊速度大于[110]和[111]方向,表明了[100]方向擁有更強的抗沖擊能力.

2.2 層裂強度

式中,ρ0和cB分別為初始密度和體波聲速,ΔU=umax-upull-back為自由面速度演化曲線中峰值速度與第一個回跳點速度之差.在這里,初始密度ρ0=7.856 g/cm3[14],體波聲速cB=4.5 km/s[33].

然而,在沖擊速度超過Hugoniot 彈性極限時,在沖擊壓縮過程中已發生塑性變形,此時需要考慮彈塑性響應和層裂片厚度的影響,修正后的層裂強度可表示為[26,34-35]

此外,層裂強度還可通過沖擊方向最大拉應力計算,即[35-36]

為了進一步分析層裂強度與沖擊方向和沖擊速度的關系,采用經典和修正聲學近似方法以及最大拉應力方法,得到[100],[110]和[111]方向下層裂強度與沖擊速度的關系,如圖5 所示.

圖5 采用三種方法得到沿[100],[110]和[111]方向沖擊時層裂強度與沖擊速度的關系Fig.5 Three methods are used to obtain the relationship between spall strength and shock velocity when the shocking along the [100],[110] and[111] directions

從圖5 中可以看出,在沿[100]方向進行沖擊時以三種方法獲得的層裂強度都表現出隨沖擊速度的增加而減小的變化趨勢.在[110]和[111]方向下,當處于彈塑性雙波分離的臨界速度時([110] 和[111]方向下分別為0.9 和1.2 km/s),三種方法獲得的層裂強度都發生了明顯減小.其中,在雙波分離現象出現時,由于彈性波前粒子速度與塑性波速度相差不大,造成自由面速度曲線中的峰值速度下降,從而使得通過經典和修正聲學近似方法計算層裂強度時自由面速度演化曲線中峰值速度與第一個回跳點速度之差減小,造成層裂強度大幅減小,使得在該速度下通過經典和修正聲學近似方法計算的層裂強度小于更快沖擊速度下的層裂強度.此外,由于雙波分離現象的出現,使得在塑性波前出現了熱尖峰,如圖6所示.塑性波前的熱尖峰會在塑性波反射后再次對模型的溫度產生影響,從而使得通過最大拉應力得到的層裂強度在雙波分離現象出現時也表現出了明顯的減小趨勢.

圖6 以0.6 和0.9 km/s 的速度沿[110]沖擊時的溫度剖面圖Fig.6 Stress profile in the [110] direction at the shock velocities of 0.6 and 0.9 km/s

在采用經典和修正聲學近似方法計算層裂強度時所需要的一些關鍵物理量,如體波聲速、初始密度、縱波聲速均不是動態變化的,因此只能近似地得到動態力學行為中的層裂強度.然而,MD 可以實時提取沖擊方向的應力分布.因此,采用在沖擊加載時模型沖擊方向出現的最大局部拉伸應力作為層裂強度更準確.在這里,采用最大拉伸應力來描述層裂強度與沖擊方向和沖擊速度之間的關系,如圖7 所示.

圖7 采用最大拉伸應力方法得到沿[100],[110]和[111]方向沖擊時層裂強度與沖擊速度的關系Fig.7 The maximum tensile stress method is used to obtain the relationship between spall strength and shock velocity when shocked along the [100],[110] and [111] directions

可以看出,沿三個方向沖擊時層裂強度都表現出了隨著沖擊速度增加而減小的變化趨勢.彈塑性雙波分離的出現使得[110]和[111]方向的層裂強度大幅減小.整體而言,[100]方向擁有更大的層裂強度.在塑性變形較小時(Up≤0.9 km/s),沿[111]方向沖擊時層裂強度大于[110]方向.但是,隨著沖擊速度的增加產生較大的塑性變形后(Up≥1.2 km/s),沿[110]方向沖擊時層裂強度大于[111]方向.

為探究沖擊速度以及沖擊方向帶來的層裂強度的差異,采用CNA 方法[29]分析以0.9 km/s 和1.5 km/s 的沖擊速度沿[100],[110]和[111]方向沖擊時的微觀組織演化,如圖8 所示.從圖8(a) 和8(b) 中可以看出,在沿[100] 方向的沖擊階段(5 ps)產生了BCC 中間相,并且在高速沖擊下會產生更多的BCC 中間相.BCC 中間相的產生與沖擊方向和壓力相關,其過程為沿[100] 方向加載時,FCC 結構中的體心四方(body-centered tetragon,BCT)結構轉變為BCC 結構,進而表現出FCC 結構向BCC 中間相的轉變,該過程遵循Born 準則[19,37].BCC 中間相會抑制層錯和無序結構的產生[19].此外,由于BCC 相的結構不穩定,以及FCC 和BCC 兩相間能壘隨著壓縮應力的減小而減小[37],使得在沖擊波到達自由面并發生反射后(15 ps),在稀疏波的卸載作用下BCC 中間相轉變為FCC 結構.

由圖8(c)~圖8(f) 可以看出,在沿[110] 和[111]方向的沖擊階段主要產生了層錯以及無序結構,并且產生的無序結構隨著沖擊速度的增加而增加.在層裂發生初期階段(30 ps),微孔洞形核于大量無序結構處[38],表明該區域發生了非晶化.此外,隨著沖擊速度的增加,產生了更多的微孔洞,并且微孔洞的分布更廣.

圖8 采用CNA 方法分析以0.9 km/s 和1.5 km/s 的沖擊速度沿[100],[110]和[111]方向沖擊時的微觀組織演化Fig.8 Microstructure evolutions by CNA at the shock velocities of 0.9 km/s and 1.5 km/s along the [100],[110] and [111] directions

非晶化是由拉伸應力和溫升共同導致的[13],為了探究非晶化與拉伸應力和溫度的關系,得到沿[100]方向以0.9 km/s 的速度沖擊時的應力剖面圖與溫度演化過程,分別如圖9 和圖10 所示.可以看出,在沖擊過程中(5,10 ps)模型處于壓應力狀態,壓應力和溫升隨著波前沿沖擊方向的傳播而出現在更深的區域.隨后,沖擊波到達自由表面并發生反射,反射后的稀疏波在與衰減沖擊波尾部的稀疏波相遇后(20 ps),形成拉應力區域,并且在拉應力區中產生微孔洞以及發生非晶化.此外,在非晶化區域中擁有更高的溫度.

圖9 沿[100]方向以0.9 km/s 的速度沖擊時的應力剖面圖Fig.9 Stress profile in the [100] direction at the shock velocities of 0.9 km/s

圖10 沿[100]方向以0.9 km/s 的速度沖擊時的微結構和溫度演化Fig.10 Microstructure and temperature evolutions in the [100]direction at the shock velocities of 0.9 km/s

層裂強度與非晶化程度密切相關[10].在這里,提取自由面速度達到峰值時的無序結構含量,如圖11所示.可以看出,[100]方向產生的無序結構含量最少,這使得微孔洞較[110]和[111]方向更難以形核,所以[100] 方向擁有較大的層裂強度.[110] 和[111]方向產生了較多的無序結構,并且[111]方向產生的無序結構含量小于[110]方向.但是,只有在較低速度沖擊(Up≤0.9 km/s)時[111]方向的層裂強度大于[110]方向.此外,隨著沖擊速度的增加,模型內溫度增加(圖6),甚至局部溫度超過熔點,這使得在模型內發生了非晶化,導致產生的無序結構含量隨著沖擊速度的增加而增加.

圖11 自由面速度達到峰值時無序結構含量Fig.11 The content of disordered structure at the peak free surface velocity stage

[110]和[111]方向層裂強度大小關系的轉變與微孔洞形核區域的無序結構含量有關[10].提取了沿沖擊方向寬度為10 nm 的微孔洞形核區域中無序結構含量,如圖12 所示.可以看出,在Up≤0.9 km/s 時,[110] 方向在層裂發生初期的微孔洞形核區域產生了更多的無序結構.但是,在Up≥1.2 km/s 時,[111]方向產生了更多的無序結構.由于在[111]方向雙波分離現象出現后,彈性先驅波的波前附近粒子速度明顯大于[110] 方向(圖4),導致[110] 和[111]方向下層裂初期微孔洞形核區域中無序結構含量大小關系的轉變.此外,多主元使得高熵合金擁有晶格畸變效應并產生更多的無序結構,進而降低了層裂強度[10,13].

圖12 沿[110]和[111]方向沖擊時微孔洞形核區域的無序結構含量Fig.12 The content of disordered structure in the microvoid nucleation region along the [110] and [111] directions

3 結論

本文對單晶CoCrFeMnNi 高熵合金進行了沖擊加載的MD 模擬,研究了沖擊方向與沖擊速度對沖擊波響應、層裂強度以及微結構演化的影響,并從微觀結構演化的角度分析了造成層裂強度改變的原因,獲得以下主要結論.

(1) 沿[100]方向進行沖擊時自由面速度未表現出彈塑性雙波分離現象,而以一定速度沿[110]和[111]方向沖擊時產生了雙波結構,并且雙波分離現象隨著沖擊速度的增加表現出先增強后減弱的變化趨勢.

(2) 在層裂發生的初期,微孔洞的形核產生了大量無序結構.隨著沖擊速度的增加,模型內溫度增加,在高沖擊強度和高溫下產生的無序結構更多,從而使得沿[100],[110]和[111]方向沖擊時層裂強度隨著沖擊速度的增加而降低.此外,雙波分離使得[110]和[111]方向的層裂強度下降明顯.

(3) 由于BCC 中間相的存在阻礙了無序結構的產生,使得[100]方向擁有較高的層裂強度.在沖擊速度較小時(Up≤0.9 km/s),沿[111]方向沖擊的層裂強度大于[111] 方向.但是,在沖擊速度較大時(Up≥1.2 km/s),沿[110]方向沖擊的層裂強度略大于[111]方向.層裂初期微孔洞形核區域無序結構含量大小關系的轉變,造成了[110]和[111]方向下層裂強度大小關系的轉變.

猜你喜歡
方向
2023年組稿方向
計算機應用(2023年1期)2023-02-03 03:09:28
方向
青年運動的方向(節選)
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2022年組稿方向
計算機應用(2022年1期)2022-02-26 06:57:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
如何確定位置與方向
2021年組稿方向
計算機應用(2021年3期)2021-03-18 13:44:48
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
大自然中的方向
主站蜘蛛池模板: 免费国产不卡午夜福在线观看| 国产原创第一页在线观看| 日韩AV手机在线观看蜜芽| 国产地址二永久伊甸园| 99免费在线观看视频| 人妻无码一区二区视频| 99视频国产精品| 国产一区二区三区日韩精品| 欧美综合激情| 亚洲毛片在线看| 99视频精品在线观看| 国产自在线播放| 久久伊伊香蕉综合精品| 无码在线激情片| 99re66精品视频在线观看| 免费一级毛片在线观看| 四虎永久在线精品影院| 自拍中文字幕| 福利小视频在线播放| 成人伊人色一区二区三区| 青青青国产免费线在| 欧美亚洲香蕉| 国产在线八区| 日韩亚洲高清一区二区| 色综合中文| 91精品国产综合久久香蕉922| 欧美影院久久| 欧美在线黄| 红杏AV在线无码| 欧美69视频在线| 国产午夜在线观看视频| 天天综合亚洲| 亚洲精品成人片在线观看| 国产精品第一区| 中文国产成人精品久久| 色噜噜中文网| 亚洲美女高潮久久久久久久| 国产精品第三页在线看| 国产成人乱码一区二区三区在线| 成人午夜天| 国产国产人在线成免费视频狼人色| 日韩精品一区二区三区免费| 污污网站在线观看| 亚洲第一在线播放| 成人福利在线视频免费观看| 中文字幕 日韩 欧美| 亚洲h视频在线| 久爱午夜精品免费视频| 一本久道久综合久久鬼色| 国产精品成人免费视频99| 亚洲VA中文字幕| 日韩少妇激情一区二区| 久久精品视频亚洲| 国产成人高清精品免费5388| 亚欧成人无码AV在线播放| 久久天天躁狠狠躁夜夜躁| 精品久久久久成人码免费动漫 | 国产视频资源在线观看| 亚洲欧美日韩另类在线一| 真人高潮娇喘嗯啊在线观看| 又粗又硬又大又爽免费视频播放| 国产黑丝一区| 99福利视频导航| 国产精品视频a| 91精品网站| 久久96热在精品国产高清| 久久久久88色偷偷| av在线人妻熟妇| 国产精品午夜电影| 波多野结衣一区二区三区四区| 亚洲欧美另类色图| 亚洲一区精品视频在线| 五月婷婷综合网| 欧美中文字幕在线视频| 日a本亚洲中文在线观看| 99久久国产精品无码| 日韩黄色在线| 欧洲免费精品视频在线| 全部无卡免费的毛片在线看| 日本人妻一区二区三区不卡影院| 精品国产美女福到在线不卡f| 园内精品自拍视频在线播放|