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

激光沖擊下CoCrFeMnNi 高熵合金微觀塑性變形的分子動力學模擬1)

2021-10-12 08:55:06熊啟林周留成闞前華蔣雖合
力學學報 2021年12期
關鍵詞:方向結構

杜 欣 熊啟林 周留成 闞前華 蔣雖合 張 旭 ,2)

* (西南交通大學力學與工程學院,成都 610031)? (華中科技大學航空航天學院,武漢 430074)

** (空軍工程大學等離子體動力學重點實驗室,西安 710038)

?? (北京科技大學新金屬材料國家重點實驗室,北京 100083)

引言

激光沖擊強化技術利用激光穿過約束層,作用于燒蝕層產生大量等離子體,等離子體繼續吸收激光能量發生爆炸從而產生沖擊波作用于金屬表面.當沖擊波峰值壓力大于材料的屈服強度時,激光沖擊使材料發生超高應變率(>106s-1)的塑性變形,造成材料表面出現殘余壓應力分布并且發生表面晶粒細化[1-3].與傳統的機械處理技術相比,激光沖擊強化后材料的粗糙度得到很好的改善,并且磨損壽命和疲勞強度得到了提高[4-5].因此,目前激光沖擊強化技術已經在航空工業中得到廣泛應用[3].

高熵合金是具有多個化學主元(通常大于4 種)的新型結構材料[6],多主元并未使高熵合金變脆,反而表現出更好的性能,例如高硬度、高強度、良好的塑性、出色的熱穩定性和耐腐蝕性[7-11].這些優異的性能使得高熵合金在工程中的潛在適用性受到了廣泛關注.CoCrFeMnNi 作為一種典型的高熵合金,最早由Cantor 等[12]在2004 年提出,所以也稱為“Cantor 合金”,其晶體結構為面心立方,在室溫(300 K),低溫(77 K)和超低溫(15 K)下具有出色的強度和塑性變形能力[13-15].這些優異的性能使得CoCrFeMnNi 高熵合金未來有望在天然氣輸運、航空航天等領域運用.由于激光沖擊強化技術已廣泛應用于航空航天等領域,所以研究CoCrFeMnNi 高熵合金的激光沖擊強化有重要的意義[16].

分子動力學模擬是基于經典牛頓力學對原子或者分子的運動進行模擬的方法,有助于揭示材料在原子尺度上的力學行為和變形機理[17].利用分子動力學來模擬激光沖擊,可以直觀觀察到沖擊過程中以及沖擊之后材料微觀組織的演化過程.許多學者針對激光沖擊強化問題,開展了分子動力學相關研究.例如,Meng 等[18]研究了激光沖擊波在Al-Cu 合金中的傳播過程,發現溫度會影響擴展位錯的形成和演化,并且彈、塑性波速度隨溫度的升高而減小.Xiong 等[19]研究了銅單晶受沖擊壓縮的彈塑性雙波結構,指出由沖擊波引起的缺陷形態表現出明顯的晶體取向依賴性.陳亞洲等[20]發現在激光沖擊下純鈦內部產生了孿生變形,變形孿晶生長經歷了沿垂直加載方向生長、無序生長、形成孿晶柵3 個過程,得到了彈塑性波分離的雙波結構.徐高峰等[21]研究了激光沖擊純鈦的溫度效應,指出在深冷條件下(77 K),沖擊波的速度高于常溫條件,能夠產生穩定的彈、塑性雙波結構,產生的高密度堆垛層錯釘扎了位錯,從而實現材料強化.Germann 等[22]發現在某些沖擊方向下,在波陣面附近由于原子面間的彈性振動產生了獨立的波列.Bringa 等[23]發現激光沖擊產生的高應力會阻礙晶界運動,限制了材料的軟化.上述研究都采用分子動力學模擬方法,分析了沖擊波傳播過程以及沖擊波作用于材料的微觀機理.然而,目前針對CoCrFeMnNi 高熵合金激光沖擊強化的分子動力學模擬鮮見報道,這限制了從原子尺度上理解CoCrFeMnNi 高熵合金激光沖擊過程中的微觀結構演化、沖擊波響應.

針對以上問題,利用分子動力學模擬軟件Lammps[24],對單晶CoCrFeMnNi 高熵合金進行了激光沖擊的分子動力學模擬.通過位錯提取算法(dislocation analysis,DXA) 以及共鄰分析方法(common neighbor analysis,CNA)來分析不同沖擊方向下微結構演變過程[25-27],并且通過粒子速度以及沖擊波波速來分析不同沖擊方向下的沖擊波特性,此外分析了沖擊后殘余應力與位錯密度沿沖擊深度的分布.

1 原子模型與加載方式

為研究沖擊時雙波結構以及微結構演化的取向相關性,利用Atomsk 軟件[28]分別建立尺寸為74.7 nm × 13.3 nm × 7.5 nm,取向為x=[100],y=[010],z=[001],x=[110],y=,z=[001]和x=[111],的3 個單晶CoCrFeMnNi 高熵合金模型(晶格常數a=3.595 ?),其中x方向為[100]取向的單晶模型如圖1 所示.此外,作為對比分析,建立了與CoCrFeMnNi 高熵合金具有相同尺寸及取向的純Ni 模型.在沖擊之前,采用NPT 系綜使體系在300 K 下弛豫50 ps 來達到熱力學平衡,弛豫時,x,y和z方向采用周期性邊界條件,時間步長為1 fs.沖擊加載時,沖擊方向(x方向)為自由邊界條件并采用NVE 系綜.

圖1 x 方向為[100]取向的單晶CoCrFeMnNi 高熵合金原子模型Fig.1 Atomistic model of single crystal CoCrFeMnNi high-entropy alloy with [100] orientation in x direction

CoCrFeMnNi 高熵合金的分子動力學模擬采用Choi 等[29]提出的第二近鄰修正嵌入原子(2 NN MEAM)勢函數,該勢函數已被用于模擬CoCrFeMnNi高熵合金的循環變形、壓痕、溫度以及應變率效應[30-33].基于該勢函數構建的高熵合金模型徑向分布函數如圖2(a)所示,可以看出其徑向分布呈現平滑分布的特征,這與傳統的面心立方結構(facecentered cubic,FCC)合金徑向分布函數中表現出的單峰特征有所不同,體現了高熵合金的的晶格畸變效應[34-35].此外,在該勢函數下構建的CoCrFeMnNi高熵合金原子模型的所有第一近鄰鍵對的鍵長滿足高斯分布,如圖2(b)所示.

圖2 徑向分布函數與第一近鄰所有鍵對的鍵長分布Fig.2 The radial distribution function and the bond length distribution of all bond pairs in the first nearest neighbor

在分子動力學模擬中,有3 種方法產生沖擊波:(1)收縮性邊界條件法,通過收縮邊界產生由兩側向內部傳播的沖擊波,常用在流體沖擊波的模擬;(2)對稱碰撞法,通過將材料與飛片發生碰撞來產生沖擊波;(3)活塞法,通過將邊界處一定層數的原子作為活塞來誘導沖擊波的產生[1].在本文中采用活塞法來進行激光沖擊的分子動力學模擬,取x方向邊界處三層原子作為活塞,其中活塞原子速度為Up,產生的沖擊波速度為Us.只有沖擊波引起的應力大于其Hugoniot 彈性極限,才會使材料產生動態塑性變形,從而有可能產生彈性波和塑性波分離的雙波結構.材料的Hugoniot 彈性極限 σHEL可表示為[19]

式中,υ 為泊松比,σY為屈服應力.由于彈、塑性雙波結構的產生,沖擊強化材料分為5 個區域,即活塞區、塑性區、彈塑性區、彈性區和未受沖擊區[19],如圖3 所示.

圖3 活塞法沖擊強化材料的結構分區示意圖Fig.3 Schematic diagram of the piston method

2 結果與分析

2.1 粒子速度

以1.2 km/s(高于Hugoniot 彈性極限)對單晶沿3 個方向進行沖擊加載,沖擊總時間為5 ps,并沿沖擊方向進行分層處理,取每層質心位置作為沖擊距離[36].在高應變率下,局部區域的粒子速度由原子的平移和熱運動這兩部分構成,而大多情況下,原子的熱運動速度是其運動速度[37].所以,將除去質心速度后的平均粒子速度作為局部區域內的粒子速度,得到1 ps 和5 ps 時的粒子速度剖面圖,如圖4 所示.在沖擊時間為1 ps 的沖擊初期,由于受沖擊區域較淺,導致沿各沖擊方向沖擊都未有明顯的彈、塑性雙波分離現象.在沖擊時間為5 ps 時,沿[100]方向沖擊時曲線并未出現臺階段,即并未出現彈塑性雙波分離現象,而沿[110]和[111]方向沖擊時產生了彈性波與塑性波分離的雙波結構.這是由于[100]方向有更多的滑移系(8 個滑移系)激活,導致滑移更容易發生,以至于有更顯著的塑性變形,造成塑性波的傳播加快,從而引起彈、塑性波的分離不明顯.而在[110]密排方向(4 個滑移系激活)和[111](6 個滑移系激活)沖擊時塑性波較當前沖擊方向下的彈性波波速慢,表現出了明顯的彈、塑性雙波分離現象.但是,當沖擊速度提升到一定程度之后,沿[110]以及[111]沖擊時的塑性波將追趕上彈性波,從而也不再表現出彈、塑性波分離的雙波結構[38].

圖4 粒子速度剖面圖Fig.4 Particle velocity profile

由圖4 可以明顯觀察到沿不同方向進行沖擊時的沖擊波波速的差異.在沿[100]方向進行沖擊時其塑性波波速明顯快于其余方向.沿[110]方向進行沖擊時其塑性波波速略快于沿[111] 方向進行沖擊時的塑性波波速,這與[110]方向為其密排方向有關.但是,沿[111]方向進行沖擊時,由于彈性波的單獨擾動,使得原子面在波前附近處彈性振動,導致附近粒子速度上升,從而產生了獨立的波列[23].CoCrFeMnNi 高熵合金表現出的彈塑性雙波分離現象和原子面彈性振動的取向效應與目前所研究的Cu 和Al 等面心立方結構金屬結果一致[2,39].此外,可以看出沿[111]方向進行沖擊時產生塑性波的臨界沖擊速度(約0.9 km/s) 大于[110] 方向(約0.8 km/s).臨界沖擊速度的差異與可開動滑移系的Schmid 因子大小有關.沿[111]方向加載時Schmid因子為0.272,而[110]方向的可開動滑移系Schmid因子為0.408,所以沿[111]方向進行沖擊時所需要的外加應力更大,導致了其臨界沖擊速度大于[110]方向.但是,由于沿[100]方向沖擊時的粒子速度剖面圖未表現出彈、塑性雙波分離現象,所以不能從中判斷出[100]沖擊方向下塑性變形發生的臨界沖擊速度.

2.2 局部應力

沖擊加載時的應力狀態可以反映沖擊波特性以及微觀組織變化,在這里計算了沖擊1 ps,5 ps 時各區域的正應力 σxx、最大分切應力 σshear、等效應力σvm,如圖5 所示.其中,最大分切應力 σshear和馮·米塞斯等效應力 σvm可表示為[1]

圖5 局部應力Fig.5 Local stress

從圖5 可以看出,沿[100]方向進行沖擊時,首先出現彈性先驅波,正應力、最大分切應力、馮·米塞斯等效應力在波前附近急劇下降,未表現出彈、塑性雙波分離現象.其余方向下正應力隨著深度的增加有明顯的臺階段出現,表明彈塑性雙波已分離,此時臺階段只存在彈性波的擾動.而最大分切應力和馮·米塞斯等效應力在靠近沖擊表面處較小,這是由于受到沖擊載荷的作用,材料未表現出明顯的泊松效應,導致材料處于靜水應力狀態[19].正應力、最大分切應力和馮·米塞斯等效應力在波前的急劇下降表明了模型中受沖擊區域向未受沖擊區域的轉變.沿[111]方向進行沖擊時,波前附近的局部應力增加同樣表示著由于原子面之間的彈性振動導致出現的獨立的波列.

沿[100],[110],[111] 3 個方向進行沖擊時壓力分布情況如圖6 所示.可以看出,在沿[111]方向進行沖擊時其Hugoniot 彈性極限(Hugoniot elastic limit,HEL)大于[110]方向,這同樣與[111]方向可開動滑移系的Schmid 因子大小有關.但是,由于沿[100]方向進行沖擊時的壓力分布同樣未表現出彈、塑性雙波分離現象,所以不能直觀獲得其Hugoniot彈性極限[19].

圖6 壓力分布Fig.6 Pressure distribution

2.3 微觀結構演變

采用位錯提取算法(DXA)分析位錯密度沿深度方向的分布,如圖7 所示.CoCrFeMnNi 高熵合金沖擊過程中沿深度方向位錯密度先增后減,表現出明顯的位錯密度梯度分布.在沖擊表面附近處位錯密度先增加,是由于在在沖擊表面附近處最大分切應力較小.較小的最大分切應力造成位錯誘導的塑性變形受限,從而影響沖擊表面附近的位錯成核以及增殖,使得位錯密度呈現沿深度方向先增加后減小的梯度分布.分析位錯密度的大小以及峰值出現的位置可以看出,在[110]沖擊方向下位錯密度峰值所對應的深度最深,沿[100]方向進行沖擊會有更多的滑移系開動,但是其位錯密度最小.其中,由于[110]沖擊方向為密排方向,所以其塑性波傳播速度略快于[111]沖擊方向,導致其位錯密度峰值所對應的深度最深.而[111]沖擊方向較[110]沖擊方向擁有更多的可開動滑移系,所以其塑性區域的位錯密度大于[110]沖擊方向的位錯密度.此外,CoCrFeMnNi高熵合金由于晶格畸變降低了位錯成核的能壘[40],導致了高熵合金沖擊過程中的位錯密度大于純Ni.但是,在純Ni 中以1.2 km/s 沿[111]方向進行沖擊時只出現了彈性波,未發生塑性變形.在純Ni 中沿[111]方向以1.2 km/s 進行沖擊未產生塑性變形反映出了在此沖擊方向下的臨界沖擊速度大于[100]和[110] 方向,這與在CoCrFeMnNi 高熵合金在[111]方向下擁有更高的臨界沖擊速度相同.同樣,在純Ni 的沖擊加載過程中也產生了沿沖擊深度方向先增后減的梯度位錯密度.

圖7 位錯密度分布Fig.7 Dislocation density distribution

沿[100]方向進行沖擊時擁有最小的位錯密度是由于其產生了體心立方(body-centered cubic,BCC)中間相,導致了更少的不全位錯形核,如圖8所示.圖中3 個及以上藍色原子層表示為BCC 中間相,綠色原子為FCC 結構,紅色原子為密排六方(hexagonal close-packed,HCP)結構.圍繞兩個或多個FCC 原子層的一層HCP 原子層是一個孿晶界,兩個相鄰的HCP 原子層構成內層錯,兩個HCP 原子層夾一層FCC 原子層構成外層錯,白色原子部分表示無序結構[28].在沖擊過程中材料內塑性波傳播的區域產生了大量的層錯以及無序結構[27].此外,在沖擊過程中產生了大量的擴展位錯,進一步提高位錯滑移的阻力,從而促進位錯的進一步增殖.高密度位錯使周邊晶格發生畸變,導致原子排列呈無序化,產生了無序結構.從缺陷在材料內的分布深度可以看出在沿[110]方向沖擊時塑性波傳播快于[111]方向.

圖8 采用共鄰分析方法分析微觀組織變化Fig.8 Microstructure evolutions by CNA

將CoCrFeMnNi 高熵合金與純Ni 在沖擊時間為5 ps 時受沖擊區域內的微結構進行對比,如圖9所示.可以看出,由于高熵合金中存在的晶格畸變效應,使得其在沖擊過程中更容易發生晶格失配[40].導致了在高熵合金中產生的無序結構含量高于純Ni,但在[100]沖擊方向下高熵合金產生的BCC 中間相含量低于純Ni.

圖9 沖擊時間為5 ps 時CoCrFeMnNi 高熵合金與純Ni 受沖擊區域內微結構含量Fig.9 The microstructure content in the impact area of CoCrFeMnNi HEA and pure Ni at 5 ps

FCC 晶胞中的面心原子和與之相鄰的FCC 晶胞中的面心原子以及共邊頂角原子和共面的面心原子可以構成一個體心四方 (body-centered tetragon,BCT)結構,如圖10 所示.圖中綠色直線構成了兩個相鄰的FCC 晶胞,藍色直線構成了一個BCT 晶胞.

圖10 FCC 結構中的BCT 結構示意圖Fig.10 Schematic diagram of BCT structure in FCC structure

在沿[100]方向進行沖擊時,由于沖擊波高壓作用導致FCC結構中BCT晶胞的高度從a壓縮至與其橫向長度相同的高熵合金的晶格常數),使得BCT 結構轉變為BCC結構,進而表現出FCC 結構向BCC 中間相轉變的現象[39].在沿FCC 結構的〈100〉晶向方向進行加載時,FCC 結構在沖擊波作用下是否穩定可通過修正的Born 穩定性準則來判斷,即[41]

式中,M1為自旋穩定準則,M2為剪切穩定準則,M3為Born 穩定準則,P為壓力,Cij為完美晶體的彈性常數.局部壓力沿沖擊深度方向分布如圖11 所示,可以看出,在受沖擊區域的壓力超過了Born 穩定準則的臨界壓力,但是未達到剪切穩定準則的臨界壓力.因此,在沿[100]方向進行沖擊加載時FCC 結構中的BCT 結構不穩定并轉變為BCC 中間相.由于裂紋尖端的應力集中,使得這種壓力相關的相轉變機制也出現在裂紋擴展的分子動力學模擬中[42-43].

圖11 沿[100]方向沖擊時壓力分布Fig.11 Pressure distribution at [100] shocking direction

為了探究高熵合金沿[100]方向沖擊時BCC 中間相的演化情況,提取沖擊為1 ps 時受沖擊區域的微結構在隨后沖擊過程中的演化情況,如圖12 所示.可以看出,BCC 中間相在隨后的沖擊加載過程中會逐漸轉變為FCC 結構,并且有部分FCC 結構會轉變為HCP 結構(層錯).BCC 相在結構上不穩定并且兩相間能壘隨著壓縮應力的減小而減小,導致了在隨后的沖擊加載過程中部分BCC 相轉變回FCC 結構[41].高位錯密度以及位錯芯區域原子嚴重錯排,導致了在沖擊過程中產生了大量無序原子.

圖12 高熵合金中沿[100]沖擊方向下微結構演化Fig.12 Microstructure evolution at [100] direction in HEA

2.4 沖擊后殘余應力與位錯密度

將模型兩端固定,中心區域采用NVE 系綜使沖擊波在模型內充分傳播,以達到保載的目的,保載過程中沖擊波經歷了多次反射的過程.沿[100]方向沖擊后保載過程中沖擊波的傳播造成微結構變化的情況如圖13 所示.BCC 中間相的產生與沖擊波的繼續傳播以及反射密切相關,在塑性波進一步傳播時的波前附近會產生BCC 中間相,而在沖擊波已傳播過的區域中,BCC 中間相會轉變為FCC 結構以及層錯.并且,在沖擊波到達模型另一端的固定端時會發生反射,造成該區域內的壓力增加,從而產生大量的BCC 中間相.但是,隨著沖擊波的衰減,后續產生的BCC 中間相含量減少,甚至在保載完成后不存在BCC 中間相,只存在層錯以及無序結構.

圖13 沿[100]方向沖擊后保載過程的微結構演化Fig.13 The microstructure evolution of the holding process after shocking in the direction of [100]

當全局應力達到平穩水平時認為保載完成,隨后取消模型兩端約束,在NVE 系綜下弛豫以達到卸載的目的,當全局應力為零時認為卸載完成[44].卸載后獲得沖擊方向殘余應力分量 σxx隨沖擊深度的變化情況如圖14 所示.可以看出,在沖擊表面為殘余壓應力,芯部為拉應力.由于沖擊波的反射,造成了殘余應力表現出了雙向沖擊的分布情況[45].沿單晶[100]方向沖擊時芯部的殘余拉應力最大.沿[111]方向沖擊時較[110]方向擁有更多的可開動滑移系,導致[111]方向受沖擊區域的塑性變形更加劇烈,造成在[111]方向沖擊表面的殘余壓應力大于[110]方向.雖然沿[100]方向進行沖擊時可開動滑移系最多,但中間相的產生導致了沖擊后沖擊表面附近的殘余壓應力小于沿[111]方向沖擊后的殘余壓應力.值得注意的是,沖擊以及保載過程中的約束層原子對殘余應力的大小產生了一定的影響[46].

圖14 殘余應力 σxx 沿深度分布曲線Fig.14 The enolution curves of residual stress σxx with the depth

分析卸載之后位錯密度沿沖擊深度分布情況如圖15 所示.可以看出,卸載之后位錯密度同樣表現出沿深度先增加后減小的梯度分布情況.但是,由于沖擊波的反復反射導致塑性變形中缺陷的動態回復,造成了卸載后的位錯密度大小較沖擊過程中的有所減小.此外,由于[110]方向下產生了較多的無序結構,造成位錯在無序結構處釘扎,使沖擊波反射后位錯的反向運動受阻,以至于位錯動態回復更難發生,導致卸載之后[110]方向下的位錯密度略大于[111]方向下的位錯密度.

圖15 卸載后位錯密度分布Fig.15 Dislocation density distribution after unloading

激光沖擊強化會產生梯度位錯密度、梯度殘余應力以及梯度晶粒尺寸3 種梯度結構.梯度位錯密度可以使材料強度得到有效提高[47-49];梯度晶粒尺寸可以保持材料的應變強化能力[50-53]梯度殘余應力可以降低疲勞裂紋的發生幾率和抑制裂紋的擴展,從而延長零件的服役壽命[45,54-55].但限于分子動力學模擬的尺寸限制,在這里無法獲得晶粒細化的模擬結果.以上研究成果可對分析實際工程中的沖擊波傳播特性和沖擊波對材料的加工應用進行理論指導.

3 結論

基于Atomsk 和Lammps 軟件對CoCrFeMnNi高熵合金進行了激光沖擊的模擬,研究了不同沖擊方向下沖擊過程中的沖擊波的傳播特性、局部應力的分布情況、微觀結構的變化以及沖擊后的殘余應力和位錯密度分布情況,并將CoCrFeMnNi 高熵合金沖擊過程微結構演化與純Ni 進行對比,獲得以下主要結論.

(1) 激光沖擊強化誘導的彈塑性雙波分離現象表現出了明顯的取向效應.[100]方向擁有更多的可動滑移系,在沿[100]方向進行沖擊時,滑移更容易發生,塑性變形更容易.這導致了在[100]沖擊方向下塑性波傳播速度變快,從而未出現雙波結構.而沿[110]和[111]方向進行沖擊加載時,可動滑移系較少,其塑性變形也更困難,進而導致了彈、塑性雙波分離的現象.

(2) 沿[100]方向進行沖擊時產生了BCC 中間相,中間相在后續的塑性變形過程中演變為FCC 結構、層錯與少量的無序結構.在沿[110]以及[111]方向進行沖擊時會產生大量的無序結構,這是由于高密度位錯以及位錯芯區域原子錯排,使得原子排列呈無序化.

(3) CoCrFeMnNi 高熵合金由于晶格畸變效應使得位錯更容易形核,導致在沖擊過程中產生了較純Ni 更高的位錯密度以及更多的無序結構.

(4) 由于沖擊波產生的塑性變形不均勻,導致在卸載后模型兩端產生殘余壓應力,芯部產生殘余拉應力.并且,卸載之后的位錯密度同樣呈現沿深度先增加后減小的梯度分布情況.

猜你喜歡
方向結構
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2022年組稿方向
計算機應用(2022年1期)2022-02-26 06:57:42
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
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
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
主站蜘蛛池模板: 国产高潮流白浆视频| 在线观看91精品国产剧情免费| 狠狠亚洲婷婷综合色香| 欧美a在线视频| 日本一本在线视频| 国产日韩精品欧美一区灰| av在线手机播放| 红杏AV在线无码| 欧美在线免费| 亚洲精品国产成人7777| 热热久久狠狠偷偷色男同| 亚洲成人77777| 宅男噜噜噜66国产在线观看| 中日韩一区二区三区中文免费视频| 人妻丰满熟妇AV无码区| 欧美成人手机在线观看网址| 亚洲日韩AV无码精品| 国产农村1级毛片| 91尤物国产尤物福利在线| 无码中文字幕乱码免费2| 亚洲精品黄| 国产视频欧美| 影音先锋亚洲无码| 干中文字幕| 精品一区二区三区无码视频无码| 亚洲成a人片7777| 2021无码专区人妻系列日韩| 婷婷激情五月网| 四虎在线高清无码| 亚洲二区视频| 中文字幕精品一区二区三区视频| 日韩在线欧美在线| 亚洲综合色婷婷| 亚洲精品欧美日本中文字幕| 国产一区二区网站| 日韩av高清无码一区二区三区| 本亚洲精品网站| 国内精自视频品线一二区| 日韩a级毛片| 成人av专区精品无码国产 | 国产大片喷水在线在线视频| 亚洲综合九九| 尤物国产在线| 国产亚洲男人的天堂在线观看| 久久精品最新免费国产成人| 国产精品亚洲一区二区在线观看| 中文字幕在线不卡视频| 亚洲午夜综合网| V一区无码内射国产| 91小视频在线观看免费版高清| 在线免费观看AV| 精品福利视频导航| 亚洲国产精品久久久久秋霞影院 | 波多野结衣一区二区三区四区| 在线亚洲精品自拍| 91无码网站| 毛片a级毛片免费观看免下载| 国产精品私拍99pans大尺度| 欧洲在线免费视频| lhav亚洲精品| 久久久噜噜噜久久中文字幕色伊伊| 国产一区二区人大臿蕉香蕉| 久久6免费视频| 538国产在线| 亚洲第一成年人网站| 中文字幕免费在线视频| 色综合中文字幕| 精品一区二区三区四区五区| 一级全黄毛片| 国产迷奸在线看| 国产中文在线亚洲精品官网| 热久久这里是精品6免费观看| 精品第一国产综合精品Aⅴ| 成人伊人色一区二区三区| 九九这里只有精品视频| 久草网视频在线| 国产熟女一级毛片| 欧美a在线视频| 再看日本中文字幕在线观看| 亚洲视频三级| 激情六月丁香婷婷四房播| 国产精品密蕾丝视频|