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

含氫原子缺陷晶界的剪切行為1)

2017-07-03 14:59:42趙東偉郁汶山申勝平
力學(xué)學(xué)報(bào) 2017年3期
關(guān)鍵詞:變形

趙東偉 郁汶山 申勝平

(西安交通大學(xué)航天航空學(xué)院機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室/陜西省航天結(jié)構(gòu)振動(dòng)控制工程實(shí)驗(yàn)室,西安710049)

含氫原子缺陷晶界的剪切行為1)

趙東偉 郁汶山2)申勝平3)

(西安交通大學(xué)航天航空學(xué)院機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室/陜西省航天結(jié)構(gòu)振動(dòng)控制工程實(shí)驗(yàn)室,西安710049)

針對(duì)4個(gè)α-Fe對(duì)稱傾斜晶界,采用分子靜力學(xué)考察了4個(gè)晶界中H原子偏析能的分布特征,并采用分子動(dòng)力學(xué)方法研究了晶界內(nèi)植入不同數(shù)量H原子對(duì)其在室溫條件下剪切行為的影響.H原子通過隨機(jī)方式植入界面內(nèi),利用植入H原子數(shù)量與晶界面積的比值來定義H原子面密度ρ.在含H原子晶界剪切行為分析過程中,重點(diǎn)考察了在不同H原子密度ρ下,4個(gè)晶界的初始塑性臨界應(yīng)力和晶界遷移位移的變化趨勢(shì)以及4個(gè)晶界在加載過程中的微觀變形機(jī)理.研究表明:晶界內(nèi)的H原子偏析能明顯偏低,4個(gè)晶界附近的H原子會(huì)自發(fā)向晶界內(nèi)偏析;隨著植入H原子數(shù)量的逐漸增多,晶界的初始塑性臨界應(yīng)力和后續(xù)變形階段應(yīng)力均會(huì)降低.晶界內(nèi)植入H原子會(huì)從本質(zhì)上改變晶界的微觀變形機(jī)理,進(jìn)而影響晶界在外載荷條件下的遷移屬性.與不含H原子晶界的變形機(jī)理對(duì)比發(fā)現(xiàn),加載過程中晶界的微結(jié)構(gòu)會(huì)發(fā)生劇烈的演化,H原子的擴(kuò)散和團(tuán)簇化效應(yīng)會(huì)導(dǎo)致晶界內(nèi)出現(xiàn)納米孔缺陷.

晶界,氫原子,偏析能,剪切行為

引言

金屬材料的氫脆(hydrogen embrittlement)問題最早由Johnson于1874年提出[1],此后,氫脆問題在學(xué)術(shù)界和工程領(lǐng)域得到證實(shí)并受到廣泛的關(guān)注[23].這是因?yàn)榇罅抗こ虒?shí)踐和研究發(fā)現(xiàn),金屬材料中氫(H)的出現(xiàn)不僅會(huì)降低金屬材料韌性導(dǎo)致材料出現(xiàn)脆斷現(xiàn)象,而且會(huì)極大地影響到材料塑性變形行為等其他力學(xué)屬性[49].這些影響直接關(guān)系到金屬構(gòu)件的安全性和使用壽命.鑒于此,目前已有大量借助實(shí)驗(yàn)[5-15]、理論分析[12,16]和微尺度模擬[1720]等手段對(duì)金屬氫脆問題開展的研究工作.通過這些研究結(jié)果可以獲得對(duì)工程材料性能設(shè)計(jì)和結(jié)構(gòu)安全預(yù)測(cè)有益的結(jié)論.然而基于宏觀層面的實(shí)驗(yàn)和理論研究工作,無法從微觀層面上揭示和分析材料氫脆問題所涉及到的與點(diǎn)缺陷相關(guān)微觀機(jī)理.事實(shí)證明,這個(gè)局限性可以借助密度泛函理論方法和分子動(dòng)力學(xué)方法的微尺度模擬來突破[2123].這方面也已有大量的研究報(bào)道,如H原子的擴(kuò)散[2426]、H原子與其他缺陷體的相互作用[24,2728]以及H原子對(duì)材料斷裂性質(zhì)的影響[17,19-20,29-33]等.

H原子在面心或者體心立方結(jié)構(gòu)的金屬單晶體中,會(huì)在不同的Fe原子間隙停留,形成如四面體和八面體等多種不同構(gòu)型.這些H原子在溫度場(chǎng)和應(yīng)力場(chǎng)等其他驅(qū)動(dòng)力作用下,會(huì)從晶體內(nèi)的一個(gè)間隙跳躍至另一個(gè)間隙位置.通過這種運(yùn)動(dòng)機(jī)制,H原子實(shí)現(xiàn)了晶體內(nèi)的擴(kuò)散.大量研究證實(shí),當(dāng)晶體材料內(nèi)存在晶界(grain boundaries,GBs)時(shí),處于晶界附近的H原子會(huì)傾向于向晶界內(nèi)部運(yùn)動(dòng)[5,2526].大量H原子在晶界處的偏析會(huì)提高晶界內(nèi)H原子點(diǎn)缺陷的濃度,影響晶界的力學(xué)行為,改變材料宏觀變形屬性[14,3438].晶界作為金屬材料內(nèi)一類十分常見的面缺陷,其在不含雜質(zhì)原子條件下的力學(xué)行為已有大量的研究工作,這些工作不僅涵蓋了多種金屬材料,還從不同加載模式下對(duì)晶界的微觀變形機(jī)理做了充分的研究[3944].但目前關(guān)于H原子缺陷對(duì)晶界力學(xué)行為影響方面的微尺度模擬研究只有關(guān)于晶界斷裂行為的報(bào)道[17],而關(guān)于H原子缺陷對(duì)晶界的剪切行為鮮見報(bào)告.為此本文以體心立方結(jié)構(gòu)α-Fe材料中的4個(gè)常見晶界作為研究對(duì)象,借助分子動(dòng)力學(xué)模擬方法研究H原子缺陷對(duì)晶界在剪切載荷作用下的微觀變形等力學(xué)響應(yīng).

本文首先借助分子靜力學(xué)方法考察了 4個(gè)以〈110〉為傾斜軸的α-Fe對(duì)稱傾斜晶界(symmetric tilt grain boundaries)結(jié)構(gòu)和其中H原子偏析能分布規(guī)律;然后采用分子動(dòng)力學(xué)方法研究了剪切載荷作用下不同數(shù)量H原子對(duì)晶界剪切行為的影響,同時(shí)還詳細(xì)分析了加載過程中的微觀變形機(jī)理,并與不含H原子晶界剪切變形機(jī)理進(jìn)行了對(duì)比.通過本研究,為進(jìn)一步分析和揭示H原子缺陷對(duì)多晶α-Fe材料力學(xué)行為的影響提供一定理論參考.

1 晶界模型

考慮以〈110〉傾斜軸的4個(gè)在金屬材料中常見的 α-Fe對(duì)稱傾斜晶界作為研究對(duì)象,它們分別是∑9(114),∑9(221),∑11(113)和∑11(332).這4個(gè)晶界的取向差角度(misorientation angle)分別為38.94?,141.06?,50.48?和129.52?.它們是∑9和∑11這兩大類晶界取向差下的4個(gè)最基本的晶界類型.圖1為兩個(gè)晶粒A和B構(gòu)成的晶界原子模型示意圖,其中晶界處于模型高度方向的中間位置.晶界的幾何和物理參數(shù)見表1.

模型中 x和 y方向分別與晶界的傾斜軸和晶界面法向方向一致.模型在 3個(gè)方向的尺寸記為L(zhǎng)x×Ly×Lz,它們必須是每個(gè)晶粒在相應(yīng)方向?qū)?yīng)晶向指數(shù)長(zhǎng)度的整數(shù)倍[4546].

表1 晶界的幾何和物理參數(shù)Table 1 The geometrical and physical properties of four GBs

圖1 晶界模型Fig.1 Schematic of grain boundary model

本文所有模擬均采用LAMMPS軟件[47]實(shí)現(xiàn),原子結(jié)構(gòu)圖均在Ovito軟件[48]中進(jìn)行后處理.此外,借助嵌入原子勢(shì)函數(shù) (embedded atom method,EAM)[4950]模擬Fe原子之間以及Fe原子與H原子之間的力場(chǎng)作用.采用由Ackland等[51]和Ramasubramaniam等[52]發(fā)展的可模擬Fe-H系統(tǒng)的EAM勢(shì)函數(shù),不僅能較為準(zhǔn)確地描述α-Fe晶體和晶界的力學(xué)性質(zhì),還能描述H原子缺陷在α-Fe晶界中的諸多性質(zhì).此外,為了能獲得處于能量最低狀態(tài)的晶界結(jié)構(gòu),在晶界模型進(jìn)行能量?jī)?yōu)化之前,必須要對(duì)模型在晶界核附近的原子數(shù)量進(jìn)行調(diào)整[45].采用Tschopp等[5354]發(fā)展的方法來調(diào)控晶界中的原子數(shù)量.該方法首先通過兩個(gè)晶粒間的相互剛性滑移來枚舉不同結(jié)合方式的晶界,接著刪除晶界內(nèi)相鄰原子間距小于給定準(zhǔn)則中的其中一個(gè)原子.通過對(duì)所有枚舉出的晶界結(jié)構(gòu)在0K下進(jìn)行能量?jī)?yōu)化,最終篩選出晶界能最低的晶界結(jié)構(gòu).晶界能γGB的計(jì)算公式為[45]

其中,每個(gè) α-Fe原子的內(nèi)聚能量 (cohesive energy)E為-4.013eV;AGB為晶界在xoz面內(nèi)的面積.為了排除自由表面對(duì)晶界能計(jì)算的影響,僅考慮能量?jī)?yōu)化后晶界模型中的包含晶界且厚度為40?的塊體.因此,式(1)中ETot和N分別為此塊體的總勢(shì)能和總Fe原子數(shù)量.所獲得的平均晶界能見表1,與文獻(xiàn)結(jié)果對(duì)比十分吻合.圖2所示為4個(gè)晶界的原子結(jié)構(gòu),可以發(fā)現(xiàn)由于晶界類型的不同,晶界的結(jié)構(gòu)不同.

圖2 能量?jī)?yōu)化后的晶界結(jié)構(gòu);晶界結(jié)構(gòu)的原子著色采用了共緊鄰原子分析方法[55]Fig.2 Relaxed structures of four GBs;Atoms in GB structure are colored using common neighbor analysis(CNA)method[55]

在模擬4個(gè)晶界在不含H和含H原子情況下的z方向剪切行為過程中,僅考慮室溫條件(300K).全部系綜均為NVT,并采用Nos′e-Hoover熱浴[5657]實(shí)現(xiàn)系統(tǒng)控溫.剪切加載通過給模型上下表面附近厚度為10?的塊體施加z方向位移來實(shí)現(xiàn).每一步位移加載量控制為0.2?,積分時(shí)間步長(zhǎng)為1fs,應(yīng)變率控制為分子動(dòng)力學(xué)典型的應(yīng)變率108s-1.

2 結(jié)果及分析

2.1 晶界內(nèi)H原子偏析能

為了便于討論不同晶界內(nèi)不同位置處H原子偏析能力的差異性,定義H原子在晶界內(nèi)的偏析能E為[25,58]

相對(duì)于單晶體而言,晶界結(jié)構(gòu)的復(fù)雜性使得H原子在晶界內(nèi)部的偏析位置十分復(fù)雜.為了獲得更多H原子在晶界內(nèi)的偏析位置,借助Voronoi胞體劃分法實(shí)現(xiàn).首先在模型x,y和z三個(gè)方向施加周期邊界的條件下,借助voro++代碼[59]獲得4個(gè)晶界原子模型中的所有原子的Voronoi胞體.進(jìn)一步將每個(gè)Voronoi胞體的節(jié)點(diǎn)視作H原子的插入位置.在獲得這些位置信息的基礎(chǔ)上,針對(duì)每個(gè)位置插入H原子并借助式(2)來計(jì)算H原子偏析能.

圖 3所示為 ∑9(114)和 ∑11(332)晶界內(nèi)部氫原子偏析能分布趨勢(shì).從晶界結(jié)構(gòu)圖(圖2)中可以看出,晶界核位置Fe原子結(jié)構(gòu)不同于遠(yuǎn)離晶界處的BCC結(jié)構(gòu).從CNA方法分析的結(jié)果來看,晶界核中的Fe原子具有HCP、BCC和其他復(fù)雜結(jié)構(gòu)特征,使得H原子偏析能在這兩個(gè)晶界內(nèi)部呈顯著變化.此外,晶界結(jié)構(gòu)復(fù)雜性導(dǎo)致了H原子偏析能分布的不均勻性.

圖3 ∑9(114)和∑11(332)內(nèi)部氫原子偏析能分布趨勢(shì)(全部結(jié)構(gòu)圖為[110]方向投影圖)Fig.3 Distributions of hydrogen segregation energies in GBs∑9(114)and∑11(332)

圖4所示為4個(gè)晶界中的H原子偏析能隨晶界法向距離變化的趨勢(shì).從圖4可以看出,H原子偏析能分布特征在4個(gè)晶界中明顯不同.由于4個(gè)晶界均為傾斜對(duì)稱晶界,因而偏析能在晶界位置(GB location)兩側(cè)的分布具有一定的對(duì)稱性.針對(duì)∑9(114),∑11(221),∑11(113)和∑11(332)四個(gè)晶界,最小H原子偏析能分別為:-0.527 eV,-0.491 eV,-0.539 eV和-0.457 eV.顯然,它們彼此間差別不顯著.此外,H偏析能僅僅在晶界核附近的-5?~5?區(qū)域內(nèi)顯著變化,而在遠(yuǎn)離晶界核的位置處偏析能趨向于零.這說明當(dāng)H原子處于晶界核附近一定距離范圍內(nèi)時(shí),極易被晶界吸收.文獻(xiàn)中常常將此距離范圍定義為晶界附近點(diǎn)缺陷的特征尺度,將其與最小偏析能和最小缺陷形成結(jié)合來衡量不同晶界對(duì)各類點(diǎn)缺陷的吸收效率的差異[25,4546].關(guān)于H原子在晶界內(nèi)的其他偏析性質(zhì),詳見文獻(xiàn)[25].從以上H偏析能的分析來看,本文所選取的4個(gè)晶界對(duì)其附近的H原子具有自發(fā)吸收能力.

圖4 氫原子偏析能隨晶界面法向方向距離變化趨勢(shì)Fig.4 Variations of hydrogen segregation energies in four GBs vs.distance from GBs

2.2 晶界的耦合剪切變形

圖5給出了在300K條件下4個(gè)晶界沿z方向施加剪切載荷的應(yīng)力--應(yīng)變曲線.從圖5可以看出,加載過程中雙晶粒晶界首先發(fā)生彈性變形,表現(xiàn)為應(yīng)力呈線性增大.當(dāng)外載荷達(dá)到晶界發(fā)生塑性變形的臨界值時(shí),應(yīng)力突然降低.繼續(xù)加載應(yīng)力再次逐漸線性增大直至應(yīng)力再突然降低.這種現(xiàn)象常被稱為黏滑過程(stick-slip)[44,60].特定的晶界變形機(jī)理誘發(fā)了加載過程中的每次應(yīng)力釋放.可以明顯看出,4個(gè)晶界在剪切載荷作用下的應(yīng)力--應(yīng)變響應(yīng)明顯不同.這些不同點(diǎn)主要表現(xiàn)在塑性變形發(fā)生時(shí)對(duì)應(yīng)的臨界應(yīng)力和后續(xù)變形階段的應(yīng)力響應(yīng).在加載過程中,與其他3個(gè)晶界相比,∑11(332)的應(yīng)力較小.

圖5 剪切載荷作用下不含H原子晶界的應(yīng)力--應(yīng)變曲線(4個(gè)晶界的初始塑性加載點(diǎn)分別用A,B,C和D標(biāo)注)Fig.5 Stress vs.strain curves of four GBs without hydrogen during the shear

實(shí)際上,本文選取的4個(gè)晶界在加載過程中均發(fā)生了耦合剪切變形(shear-coupling).這種特殊的晶界變形機(jī)理最早由Cahn等[44,60]發(fā)現(xiàn)于Cu晶界中.簡(jiǎn)言之,在這種變形模式下,晶界遷移(migration)和晶界附近的材料剪切變形同時(shí)發(fā)生.圖6所示為4個(gè)晶界在加載過程中的遷移位移與加載應(yīng)變之間的變化關(guān)系.從圖6可以看出,除了∑9(114),其他3個(gè)晶界均向上遷移(正位移).此外,∑11(332)的遷移是一個(gè)持續(xù)不斷的過程,且從加載起始階段便發(fā)生.而其他3個(gè)晶界,只有當(dāng)加載應(yīng)變達(dá)到一定值才發(fā)生遷移,且晶界遷移位移變化為跳躍增大模式,這意味著它們的遷移必須經(jīng)歷一定的應(yīng)變積累階段.

圖6 剪切載荷作用下的不含H原子晶界遷移位移曲線;晶界的遷移發(fā)生在晶界的初始塑性加載點(diǎn)A,B,C和DFig.6 GB displacements of four GBs without Hydrogen during the shear

以∑11(332)晶界為例,對(duì)耦合剪切變形機(jī)理進(jìn)行詳細(xì)分析.圖7所示為∑11(332)晶界在4個(gè)加載應(yīng)變下的原子結(jié)構(gòu)圖.首先,在初始晶界模型中的一個(gè)帶狀區(qū)域里選取部分原子,然后在加載過程中跟蹤這些原子的位移,以跟蹤晶界的滑移變形.從圖7可以看出,加載過程晶界不斷向上遷移.被跟蹤的帶狀區(qū)域除了晶界遷移所經(jīng)歷的區(qū)域呈現(xiàn)出錯(cuò)位現(xiàn)象之外,其他帶狀區(qū)域始終保持垂直狀態(tài).說明此晶界在加載過程中幾乎未發(fā)生彈性變形.加載過程中晶界中出現(xiàn)的塑形應(yīng)變可以由圖7(c)中us與模型高度的比值獲得.

圖7 不同剪切載應(yīng)變下∑11(332)晶界的變形構(gòu)型圖(全部結(jié)構(gòu)圖為[110]方向投影圖)Fig.7 Snapshots of GB∑11(332)at the di ff erent shear strains(all structures are viewed along[110]direction)

為了便于闡述∑11(332)晶界的遷移機(jī)理,圖8所示為晶界遷移過程中的兩個(gè)相鄰狀態(tài)的原子結(jié)構(gòu)圖.首先,[110]方向投影的BCC單晶體可以被視作由很多B結(jié)構(gòu)單元(structure units,SUs)周期復(fù)制生成.而∑11(332)晶界的結(jié)構(gòu),可以視作由很多A和E結(jié)構(gòu)單元構(gòu)成.在加載過程中,原始的結(jié)構(gòu)單元E會(huì)逐漸演化為單元A和B,而晶界中原始的單元A會(huì)與上半部分晶粒中的B結(jié)合生成E單元.在加載過程中,這樣的演化重復(fù)不斷發(fā)生,使得晶界持續(xù)不斷地向上遷移.在這個(gè)晶界局部結(jié)構(gòu)的演化過程中,晶界結(jié)構(gòu)微結(jié)構(gòu)演化會(huì)導(dǎo)致上半部分晶粒中的原子具有向右的位移,這正是晶界模型塑性應(yīng)變產(chǎn)生的根源.通過類似方法,還可以討論其他3個(gè)晶界在塑性變形階段的微觀變形機(jī)理.

圖8 ∑11(332)晶界耦合剪切變形機(jī)理(黑色和白色原子處于兩個(gè)相鄰的(110)面)Fig.8 Mechanism of shear-coupling deformation in GB∑11(332)(black and white atoms are on the two successive(110)planes)

2.3 含氫原子晶界的力學(xué)行為

為了研究室溫條件下4個(gè)晶界中含有不同數(shù)量H原子的剪切變形行為,定義晶界內(nèi)的H原子密度ρ為

其中,NH為晶界內(nèi)植入的H原子數(shù)量.在給定晶界范圍內(nèi),每個(gè)H原子的x和z坐標(biāo)借助隨機(jī)函數(shù)生成.需要指出的是,H原子y坐標(biāo)需要晶界位置附近給予一定的微擾量.其目的是為了避免在分子模擬運(yùn)行過程中出現(xiàn)Fe和H原子過分接近,使得勢(shì)函數(shù)無法捕捉兩者之間相互作用關(guān)系.此外,為了促使植入的H原子均勻地分布于晶界內(nèi),在剪切加載前需讓植入H原子的晶界模型在NVT系綜下演化10 ps.針對(duì)每個(gè)晶界,考慮了ρ從0.02 ?-2到0.18 ?-2變化的9種情況,在不同H原子密度下4個(gè)晶界中植入H原子數(shù)量見表2.

表2 每個(gè)晶界中植入的H原子密度對(duì)應(yīng)的H原子數(shù)量NHTable 2 The number of hydrogen atoms placed in each GB for a specifie hydrogen density

圖 9中給出了 4個(gè)晶界在 ρ = 0.0?-2,0.04?-2,0.08?-2,0.12?-2和0.16?-2下的應(yīng)力--應(yīng)變曲線.顯然,不同H原子密度下的應(yīng)力--應(yīng)變曲線在彈性階段的斜率幾乎不變,這說明H原子的植入對(duì)晶界彈性變形沒有顯著影響.然而,H原子的植入對(duì)晶界塑性變形產(chǎn)生了很大影響.這些影響表現(xiàn)在初始塑性變形所對(duì)應(yīng)的臨界應(yīng)力值變化和后續(xù)流動(dòng)階段曲線構(gòu)型的變化.

圖9 含有不同面密度H原子晶界在剪切載荷作用下的應(yīng)力--應(yīng)變曲線Fig.9 Stress vs.strain curves of four GBs for di ff erent hydrogen densities in them

圖9 含有不同面密度H原子晶界在剪切載荷作用下的應(yīng)力--應(yīng)變曲線(續(xù))Fig.9 Stress vs.strain curves of four GBs for di ff erent hydrogen densities in them(continued)

圖10所示為4個(gè)晶界的初始塑性變形臨界應(yīng)力與不同H原子面密度之間的關(guān)系曲線.從圖10可以看出,隨著H原子密度的增大,∑9(114),∑9(221)和∑11(113)晶界的初始塑性變形臨界應(yīng)力呈現(xiàn)整體減小趨勢(shì).∑9(114)和∑11(113)晶界的臨界應(yīng)力雖然在H原子面密度增大過程中,在某些數(shù)據(jù)點(diǎn)出現(xiàn)了震蕩現(xiàn)象,但震蕩幅度相對(duì)較小.∑9(221)晶界的臨界應(yīng)力數(shù)據(jù)點(diǎn)震蕩則比較強(qiáng)烈,且數(shù)據(jù)震蕩的幅度隨著H原子面密度增大而逐漸變?nèi)?這說明當(dāng)H原子面密度相對(duì)較小時(shí),H原子在晶界分布的不均勻性對(duì)∑9(221)晶界初始塑性變形臨界應(yīng)力影響比較顯著,而這種分布的不均勻性對(duì)∑9(114)和∑11(113)兩個(gè)晶界的數(shù)據(jù)點(diǎn)震蕩性影響則很弱.

圖10 晶界初始塑性臨界應(yīng)力與不同H原子面密度之間的關(guān)系曲線Fig.10 The variations of the critical stress corresponding to the incipient plasticity of four GBs vs.the hydrogen densities

由于晶界內(nèi)的H原子是通過隨機(jī)方式植入,且植入H原子的晶界在加載前首先經(jīng)歷了長(zhǎng)達(dá)10 ps的構(gòu)型演化過程;而這個(gè)過程可以在一定程度上削弱由于隨機(jī)方式植入H原子對(duì)模擬結(jié)果所帶來的隨機(jī)性,這一點(diǎn)可以由∑9(114),∑11(113)和∑11(332)晶界所表現(xiàn)出的微弱震蕩證實(shí).因而,∑9(221)晶界結(jié)果所呈現(xiàn)的具有一定規(guī)律性的減幅震蕩特點(diǎn),與該晶界變形機(jī)理的特殊性有很大關(guān)系.同時(shí),通過觀察圖9(a)~圖9(c)還可以發(fā)現(xiàn),植入H原子會(huì)導(dǎo)致晶界初始塑性變形在較小的剪切應(yīng)變下發(fā)生.由此得出結(jié)論:晶界中H原子的植入會(huì)導(dǎo)致這些晶界易于進(jìn)入塑性變形階段.從圖10還可以看出,在后續(xù)塑性變形階段,植入H原子會(huì)降低降低晶界的流動(dòng)應(yīng)力.

∑11(332)晶界在不含H原子情況(ρ=0.0 ?-2)下,其塑性變形在非常小的載荷下會(huì)持續(xù)不斷地發(fā)生.然而,H原子的植入會(huì)提高其初始塑性變形臨界應(yīng)力,并且隨著H原子密度的增大,∑11(332)晶界的初始塑性變形臨界應(yīng)力會(huì)逐漸增大.這主要是因?yàn)閳D8所示的∑11(332)晶界結(jié)構(gòu)演化的關(guān)鍵在于晶界結(jié)構(gòu)體單元之間的順利演化;而當(dāng)H原子植入該晶界時(shí),這些結(jié)構(gòu)單元的構(gòu)型會(huì)遭到破壞,致使晶界結(jié)構(gòu)無法順利地按照?qǐng)D8所示機(jī)理進(jìn)行演化.

圖11為在后續(xù)加載過程中含H原子晶界的結(jié)構(gòu)演化模式,其中,H原子密度為0.18 ?-2且剪應(yīng)變分別為0.0和0.15.從圖11可以看出,∑9(114)晶界在剪應(yīng)變?yōu)?.15時(shí)僅僅發(fā)生了滑移變形.然而,從圖6的結(jié)果來看,不含H原子∑9(114)晶界在剪應(yīng)變?yōu)?.15時(shí)晶界會(huì)向下遷移.由此可見,H原子的植入會(huì)導(dǎo)致晶界的變形機(jī)理發(fā)生本質(zhì)變化.而其他3個(gè)晶界(∑9(221),∑11(113)和∑11(332))在發(fā)生剪切滑移變形的同時(shí),還發(fā)生了向上遷移.此外,觀察初始狀態(tài)時(shí)含有H的4個(gè)晶界的局部結(jié)構(gòu),并與圖2中不含H晶界結(jié)構(gòu)進(jìn)行對(duì)比可以發(fā)現(xiàn),H原子的植入并未導(dǎo)致晶界局部結(jié)構(gòu)發(fā)生很大變化.然而,晶界結(jié)構(gòu)會(huì)隨著加載而發(fā)生劇烈變化.同時(shí)H原子不再像初始狀態(tài)時(shí)均勻分布于晶界內(nèi),而是發(fā)生劇烈擴(kuò)散現(xiàn)象,進(jìn)而逐漸匯聚成團(tuán)簇結(jié)構(gòu).這不僅使得原來平整的晶界不再平整,而且導(dǎo)致晶界內(nèi)出現(xiàn)類似納米孔的微結(jié)構(gòu).

從圖6可以看出,晶界遷移位移隨剪應(yīng)變的變化率因晶界不同而呈現(xiàn)出不同的趨勢(shì).由于不同的剪應(yīng)變與不同的加載時(shí)間相互對(duì)應(yīng),晶界遷移位移隨移隨加載時(shí)間的變化率即為晶界在加載過程中的遷移速率.

為了進(jìn)一步探究植入H原子對(duì)晶界遷移速率,圖12所示為∑9(114)和∑11(113)兩個(gè)晶界遷移位移隨剪應(yīng)變的變化曲線.從圖12可以看出,在ρ=0.0時(shí),∑9(114)晶界以恒定的速率向下遷移,當(dāng)ρ>0.0時(shí),該晶界無遷移現(xiàn)象發(fā)生,遷移速率為零.顯然,植入H原子對(duì)∑9(114)晶界的遷移行為影響巨大.然而在∑11(113)晶界內(nèi)植入H原子,其遷移速率在不同H原子面密度下或增大、或減小,且基本保持勻速遷移.可以發(fā)現(xiàn),植入H原子對(duì)∑11(113)晶界的耦合剪切變形并不徹底,僅僅影響了晶界的遷移速率;對(duì)∑9(114)晶界則徹底改變了變形機(jī)制.∑9(114)晶界在含H原子條件下的遷移速度為零,晶界塑性應(yīng)變僅僅由晶粒間滑移承擔(dān),因而其變形機(jī)理由耦合剪切變形改變?yōu)榫Ы玳g滑移.

圖11 H原子密度為0.18 ?-2且剪應(yīng)變分別為0.0和0.15時(shí),4個(gè)晶界的原子結(jié)構(gòu)圖(灰色和紅色原子分別對(duì)應(yīng)Fe和H,綠色帶狀區(qū)域用來跟蹤晶界變形過程)Fig.11 Structures of four GBs at shear strains 0.0 and 0.15 for Hydrogen density of 0.18 ?-2(Fe and hydrogen atoms are colored in grey and red,the green strip is used to track the GB deformation)

圖12 植入H原子對(duì)∑9(114)和∑11(113)晶界遷移位移的影響Fig.12 Migration displacements of GBs∑9(114)and∑11(113)for di ff erent hydrogen densities

3 結(jié)論

本文針對(duì)4個(gè)α-Fe傾斜對(duì)稱晶界建立了對(duì)應(yīng)的原子模型.首先借助分子靜力學(xué)方法,通過調(diào)整晶界內(nèi)原子數(shù)量獲得晶界能量最低狀態(tài);其次考察了氫原子在4個(gè)晶界中偏析能的分布特點(diǎn)和變化趨勢(shì):最后借助分子動(dòng)力學(xué)方法,研究了室溫條件下含不同數(shù)量H原子和不含H原子晶界的剪切變形,并詳細(xì)地分析了加載過程中晶界的微觀變形機(jī)理.通過本文研究,獲得以下結(jié)論:

(1)從4個(gè)晶界中的H原子偏析能分布規(guī)律以及較小的H原子偏析能來看,晶界附近的H原子會(huì)向晶界內(nèi)部自發(fā)偏析.

(2)H原子的植入不會(huì)改變晶界的彈性變形行為,但會(huì)改變晶界的塑形變形行為.除了特例∑11(332)晶界之外,在其他3個(gè)晶界中植入H原子均會(huì)降低初始塑性變形和后續(xù)流動(dòng)階段的應(yīng)力值.

(3)H原子的植入會(huì)從本質(zhì)上改變晶界遷移和晶界微結(jié)構(gòu)演化機(jī)理.此外,植入H原子會(huì)改變晶界的遷移速率,使得原來不含H原子晶界的耦合剪切變形模式會(huì)變?yōu)槊黠@的晶界滑移變形模式.

(4)晶界中植入的H原子在加載起始階段基本上均勻分布于晶界內(nèi)部,隨著加載會(huì)發(fā)生劇烈遷移和匯聚現(xiàn)象,這會(huì)導(dǎo)致晶界結(jié)構(gòu)中出現(xiàn)原子混亂區(qū)域和類似納米孔.

文中僅選取了4個(gè)α-Fe傾斜對(duì)稱晶界作為研究對(duì)象.最后需要指出的是,若考慮不同傾斜軸、不同取向差角度和晶界對(duì)稱與否等因素,晶界空間實(shí)際上十分龐大[25,58].因此,要獲得較為普適的關(guān)于H原子植入對(duì)晶界剪切行為影響的結(jié)論,尚需進(jìn)一步擴(kuò)大研究對(duì)象并考慮其他類型的晶界.此外,本文針對(duì)性地研究了晶界在室溫下的剪切行為,至于它們?cè)谄渌麥囟葪l件下以及拉伸載荷作用下的力學(xué)行為尚需進(jìn)一步的研究.

1 Johnson WH.On some remarkable changes produced in iron and steel by the action of hydrogen and acids.Proceedings of the Royal Society of London,1874,23:168-179

2南雲(yún)道彥.鋼的氫脆的新研究方向.熱處理,2010,25(3):1-6(Michihiko Nagumo.Turning of the research direction on hydrogen embrittlement of steel.Heat Treatment,2010,25(3):1-6(in Chinese))

3羅潔,郭正洪,戎詠華.先進(jìn)高強(qiáng)度鋼氫脆的研究進(jìn)展.機(jī)械工程材料,2015,39(8):1-9(Luo Jie,Guo Zhenghong,Rong Yonghua.Research progess on hydrogen embrittlement in advanced high strength steels.Materials for Mechanical Engineering,2015,39(8):1-9(in Chinese))

4 Nagumo M.Characteristic features of deformation and fracture in hydrogen embrittlement//Fundamentals of Hydrogen Embrittlement,Springer Singapore,Singapore,2016,137-165

5宋仁國(guó),張寶金,曾梅光.高強(qiáng)鋁合金晶界偏析與氫致斷裂機(jī)理的研究.航空材料學(xué)報(bào),1997,17(1):31-38(Song Renguo,Zhang Baojin,Zeng Meiguang.The research on the hydrogen segregation andhydrogeninducedfracturemechanismsinhighstrengthAlalloy.Journal of Aeronautical Materials,1997,17(1):31-38(in Chinese))

6張立文,劉寶璋,哈寬富.α-Fe中氫致裂紋的研究.材料科學(xué)進(jìn)展,1988,2(5):44-49(Zhang Liwen,Liu Baozhang,Ha Kuanfu.Hydrogen cracking in α-Fe.Progess in Material Science,1988,2(5):44-49(in Chinese))

7 Neeraj T,Srinivasan R,Li J.Hydrogen embrittlement of ferritic steels:Observations on deformation microstructure,nanoscale dimples and failure by nanovoiding.Acta Materialia,2012,60(13-14):5160-5171

8 Sanchez J,Lee SF,Martin-Rengel MA,et al.Measurement of hydrogen and embrittlement of high strength steels.Engineering Failure Analysis,2016,59:467-477

9 Demetriou V,Robson JD,Preuss M,et al. Study of the e ff ectofhydrogen charging on thetensilepropertiesand microstructure of four variant heat treatments of nickel alloy 718. International Journal of Hydrogen Energy,2017(http://doi.org/10.1016/j.ijhydene.2017.02.149)

10 Martin ML,Robertson IM,Sofronis P.Interpreting hydrogeninduced fracture surfaces in terms of deformation processes:A new approach.Acta Materialia,2011,59(9):3680-3687

11 Nagao A,Smith CD,Dadfarnia M,et al.The role of hydrogen in hydrogen embrittlement fracture of lath martensitic steel.Acta Materialia,2012,60(13-14):5182-5189

12 Xing X,Chen W,Zhang H.Prediction of crack propagation under cyclic loading based on hydrogen di ff usion.Materials Letters,2015,152:86-89

13 AbrahamDP,AltstetterCJ.Hydrogen-enhancedlocalizationofplasticity in an austenitic stainless steel.Metallurgical and Materials Transactions A,1995,26(11):2859-2871

14張建,李秀艷,趙明久等.晶界相對(duì)Fe-Ni-Cr奧氏體合金氫脆的影響.金屬學(xué)報(bào),2008,44(9):1095-1098(Zhang Jian,Li Xiuyan,Zhao Mingjiu,et al.E ff ects of grain boundary phases on hydrogen embrittlement of Fe-Ni-Cr austenitic alloys.Acta Metallurgica Sinica,2008,44(9):1095-1098(in Chinese))

15陳業(yè)新.金屬間化合物的環(huán)境氫脆.上海大學(xué)學(xué)報(bào)(自然科學(xué)版),2011,17(4):487-502(Chen Yexin.Environmental hydrogen embrittlement of intermetallics.Journal of Shanghai University(Natural Science Edition),2011,17(4):487-502(in Chinese))

16吳友義,卞歡,毛彩云等.納晶金屬材料氫脆的微觀力學(xué)模型.工程力學(xué),2014,31(12):228-233(Wu Youyi,Bian Huan,Mao Caiyun,et al.A micro-mechanical model of hydrogen-induced embrittlement in nanocrystalline metals.Engineering Mechanics,2014,31(12):228-233(in Chinese))

17 Tehranchi A,Curtin WA.Atomistic study of hydrogen embrittlement of grain boundaries in nickel:I.Fracture.Journal of the Mechanics and Physics of Solids,2017,101:150-165

18 Xing X,Chen W,Zhang H.Atomistic study of hydrogen embrittlement during cyclic loading:Quantitative model of hydrogen accumulation e ff ects.International Journal of Hydrogen Energy,2017,42(7):4571-4578

19 Tehranchi A,Yin B,Curtin WA.Softening and hardening of yield stress by hydrogen–solute interactions.Philosophical Magazine,2017,97(6):400-418

20 Song J,Curtin WA.Atomic mechanism and prediction of hydrogen embrittlement in iron.Nature Materials,2013,12(2):145-151

21何昱辰,劉向軍.基于粗粒化水分子模型的Cu-H2O納米流體黏度模擬.力學(xué)學(xué)報(bào),2014,46(6):871-878(He Yuchen,Liu Xiangjun.Simulation studies of viscosities of Cu-H2O nanofluidbased on coarse graining water molecules.Chinese Journal of Theoretical and Applied Mechanics,2014,46(6):871-878(in Chinese))

22劉海,李啟楷,何遠(yuǎn)航.高速?zèng)_擊壓縮梯恩梯的分子動(dòng)力學(xué)模擬.力學(xué)學(xué)報(bào),2015,47(1):174-179(Liu Hai,Li Qikai,He Yuanhang.Molecular dynamics simulations of high velocity shock compressed TNT.Chinese Journal of Theoretical and Applied Mechanics,2015,47(1):174-180(in Chinese))

23華軍,候燕,段志榮等.石墨烯輻照損傷及力學(xué)性能研究.力學(xué)學(xué)報(bào),2016,48(5):1080-1087(HuaJun,HouYan,DuanZhirong,etal.Study on irradiation damage and mechanical property of graphene.ChineseJournalofTheoreticalandAppliedMechanics,2016,48(5):1080-1087(in Chinese))

24 Fernandez N,Ferro Y,Kato D.Hydrogen di ff usion and vacancies formation in tungsten:Density functional theory calculations and statistical models.Acta Materialia,2015,94:307-318

25 SolankiKN,TschoppMA,BhatiaMA,etal.Atomisticinvestigation of the role of grain boundary structure on hydrogen segregation and embrittlement in α-Fe.Metallurgical and Materials Transactions A,2013,44(3):1365-1375

26 Bhatia MA,Groh S,Solanki KN.Atomic-scale investigation of point defects and hydrogen-solute atmospheres on the edge dislocation mobility in alpha iron.Journal of Applied Physics,2014,116(6):064302

27 Ismer L,Park MS,Janotti A,et al.Interactions between hydrogen impurities and vacancies in Mg and Al:A comparative analysis based on density functional theory.Physical Review B,2009,80(18):184110

28 Jiang B,Wan FR,Geng WT.Strong hydrogen trapping at helium in tungsten:Density functional theory calculations.Physical Review B,2010,81(13):134112

29 Fu CL,Painter GS.First principles investigation of hydrogen embrittlement in FeAl.Journal of Materials Research,1991,6(4):719-723

30 JiangDE,CarterEA.Firstprinciplesassessmentofidealfractureenergies of materials with mobile impurities:implications for hydrogen embrittlement of metals.Acta Materialia,2004,52(16):4801-4807

31 Serebrinsky S,Carter EA,Ortiz M.A quantum-mechanically informed continuum model of hydrogen embrittlement.Journal of the Mechanics and Physics of Solids,2004,52(10):2403-2430

32 Song J,Curtin WA.A nanoscale mechanism of hydrogen embrittlement in metals.Acta Materialia,2011,59(4):1557-1569

33 Zhou X,Ouyang B,Curtin WA,et al.Atomistic investigation of the influenc of hydrogen on dislocation nucleation during nanoindentation in Ni and Pd.Acta Materialia,2016,116:364-369

34 Tabata T,Birnbaum HK.Direct observations of the e ff ect of hydrogen on the behavior of dislocations in iron.Scripta Metallurgica,1983,17(7):947-950

35 Eberhart ME,Johnson KH,Latanision RM.A molecular orbital model of intergranular embrittlement.Acta Metallurgica,1984,32(6):955-959

36 Robertson IM,Tabata T,Wei W,et al.Hydrogen embrittlement and grainboundaryfracture.ScriptaMetallurgica,1984,18(8):841-846

37 RiceJR,WangJS.Embrittlementofinterfacesbysolutesegregation.Materials Science and Engineering:A,1989,107(89):23-40

38 Martin ML,Somerday BP,Ritchie RO,et al.Hydrogen-induced intergranular failure in nickel revisited.Acta Materialia,2012,60(6-7):2739-2745

39 Spearot DE,Tschopp MA,Jacob KI,et al.Tensile strength of〈100〉and 〈110〉tilt bicrystal copper interfaces.Acta Materialia,2007,55(2):705-714

40 Tschopp MA,McDowell DL.Tension-compression asymmetry in homogeneous dislocation nucleation in single crystal copper.Applied Physics Letters,2007,90(12):121916

41 Tschopp MA,McDowell DL.Dislocation nucleation in sigma 3 asymmetrictiltgrainboundaries.InternationalJournalofPlasticity,2008,24:191-217

42盧磊,陳先華,黃曉旭等.納米孿晶純銅的極值強(qiáng)度及納米孿晶提高金屬材料綜合強(qiáng)韌性.中國(guó)基礎(chǔ)科學(xué),2010,12(1):16-18(Lu Lei,Chen Xianhua,Huang Xiaoxu,et al.Revealing the maximum strengthinnanotwinnedcopper.ChinaBasisScience,2010,12(1):16-18(in Chinese))

43魏宇杰.納米金屬材料的界面力學(xué)行為研究.金屬學(xué)報(bào),2014,50(2):183-190(Wei Yujie.Investigation of mechanical behavior of interfaces in nanostructured metals.Acta Metallurgica Sinica,2014,50(2):183-190(in Chinese))

44 Cahn JW,Mishin Y,Suzuki A.Coupling grain boundary motion to shear deformation.Acta Materialia,2006,54(19):4953-4975.

45 Yu WS,Demkowicz MJ.Non-coherent Cu grain boundaries driven by continuous vacancy loading.Journal of Materials Science,2015,50(11):4047-4065

46 Yu WS,Shen SP,Liu QF.Energetics of point defect interacting with bi-crystal Σ3 copper grain boundaries.Computational Materials Science,2016,118:47-55

47 Plimpton S.Fast parallel algorithms for short-range molecular dynamics.Journal of Computational Physics,1995,117(1):1-19

48 Alexander S.Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool.Modelling and Simulation in Materials Science and Engineering,2010,18(1):015012

49 Daw MS,Baskes MI.Embedded-atom method derivation and application to impurities,surfaces,and other defects in metals.Physical Review B,1984,29:6443-6453

50 Daw MS,Foiles SM,Baskes MI.The embedded-atom method:a review of theory and applications.Materials Science Reports,1993,9:251-310

51 Ackland GJ,Mendelev MI,Srolovitz DJ,et al.Development of an interatomic potential for phosphorus impurities in α-iron.Journal of Physics:Condensed Matter,2004,16:S2629

52 Ramasubramaniam A,Itakura M,Carter EA.Interatomic potentials for hydrogen in α iron based on density functional theory.Physical Review B,2009,79:174101

53 Tschopp MA,McDowell DL.Structures and energies of∑3 asymmetric tilt grain boundaries in copper and aluminium.Philosophical Magazine,2007,87:3147-3173

54 Tschopp MA,McDowell DL.Asymmetric tilt grain boundary structure and energy in copper and aluminium.Philosophical Magazine,2007,87:3871-3892

55 Faken D,J′onsson H.Systematic analysis of local atomic structure combined with 3D computer graphics.Computational Materials Science,1994,2:279-286

56 Nose S.A unifie formulation of the constant temperature molecular dynamicsmethods.JournalofChemicalPhysics,1984,81:511-519

57 Hoover WG.Canonical dynamic equilibrium phase-space distributions.Physical Review A,1985,31:1695-1697

58 Tschopp MA,Solanki KN,Baskes MI,et al.Generalized framework for interatomic potential design:Application to Fe–He system.Journal of Nuclear Materials,2012,425:22-32

59 RycroftCH,GrestGS,LandryJW,etal.Analysisofgranularfl win a pebble-bed nuclear reactor.Physical Review E,2006,74:021306

60 Ivanov VA,Mishin Y.Dynamics of grain boundary motion coupled to shear deformation:An analytical model and its verificatio by molecular dynamics.Physical Review B,2008,78:064106

SHEAR RESPONSE OF GRAIN BOUNDARIES WITH HYDROGEN DEFECTS1)

Zhao Dongwei Yu Wenshan2)Shen Shengping3)
(State Key Laboratory for Strength and Vibration of Mechanical Structures,Shanxi Engineering Laboratory for Vibration Control of Aerospace Structures,School of Aerospace,Xi’an Jiaotong University,Xi’an 710049,China)

The segregation energy distributions of hydrogen in four α-Fe symmetric tilt grain boundaries(GBs)are analyzed by using molecular statics(MS),and then the shear responses of four GBs embedded with di ff erent number of hydrogen atoms at the room temperature are investigated by using molecular dynamics(MD)methods.To facilitate our quantitative analysis,the hydrogen density ρ is define as the ratio between the number of hydrogen atoms and the GB area.At di ff erent hydrogen densities,the variations of initial critical stress of GB plasticity and GB migration displacements are analyzed,and the micro-deformation mechanisms in each GB with the presence of hydrogen atoms are analyzed as well.It is found that the hydrogen segregation energies are generally lower in GB than those inside grain,which lead four GBs to absorb hydrogen atoms in the vicinity of GBs.With the increase of hydrogen density ρ the critical stress of incipient plasticity as well as the fl w stress could be reduced.Moreover,the micro-deformation mechanisms of fours GBs with hydrogen atoms are di ff erent from those of GBs without hydrogen atoms.In particular,presence of hydrogen atoms remarkably a ff ects GB migration velocity.Thus,GB with hydrogen atoms may undergo a pure slidingdeformation instead of the shear-coupling deformation for GB without hydrogen atoms.Meanwhile,in contrast to GBs without hydrogen atoms,the micro-structures of GB with hydrogen atoms drastically evolve upon loading.In addition,the di ff usion and agglomeration of hydrogen atoms may lead to the formation of nanovoid in GBs.

grain boundary,hydrogen,segregation energy,shear behavior

O341

:A

10.6052/0459-1879-17-132

2017–04–21 收稿,2017–05–23 錄用,2017–05–23 網(wǎng)絡(luò)版發(fā)表.

1)國(guó)家自然科學(xué)基金資助項(xiàng)目(11502191,11372238,11632014).

2)郁汶山,副教授,主要研究方向:納米力學(xué)、輻照損傷材料力學(xué).E-mail:wenshan@mail.xjtu.edu.cn

3)申勝平,教授,主要研究方向:力化學(xué)耦合、撓曲電理論及應(yīng)用.E-mail:sshen@mail.xjtu.edu.cn

趙東偉,郁汶山,申勝平.含氫原子缺陷晶界的剪切行為.力學(xué)學(xué)報(bào),2017,49(3):605-615

Zhao Dongwei,Yu Wenshan,Shen Shengping.Shear response of grain boundaries with hydrogen defects.Chinese Journal of Theoret-ical and Applied Mechanics,2017,49(3):605-615

猜你喜歡
變形
變形記
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
柯西不等式的變形及應(yīng)用
“變形記”教你變形
不會(huì)變形的云
“我”的變形計(jì)
會(huì)變形的折紙
童話世界(2018年14期)2018-05-29 00:48:08
變形巧算
例談拼圖與整式變形
會(huì)變形的餅
主站蜘蛛池模板: 美女扒开下面流白浆在线试听| 色婷婷在线影院| 国产福利小视频在线播放观看| 国模私拍一区二区| 九九精品在线观看| 毛片在线播放a| 无码AV日韩一二三区| 亚洲欧美综合精品久久成人网| 人妻免费无码不卡视频| 久久亚洲高清国产| av在线人妻熟妇| 国产成人高清在线精品| 国产农村1级毛片| 嫩草在线视频| 爱爱影院18禁免费| 国产正在播放| 色综合中文字幕| 午夜日b视频| 国产一级做美女做受视频| 久久精品国产免费观看频道| 色偷偷一区二区三区| 国产精选小视频在线观看| 中文字幕乱码中文乱码51精品| 538精品在线观看| 日韩国产黄色网站| 国产一级毛片在线| 国产一区二区三区在线无码| 国产制服丝袜无码视频| 亚洲va在线∨a天堂va欧美va| 91精品网站| 精品视频一区在线观看| 色窝窝免费一区二区三区| 亚洲一区二区成人| 一级高清毛片免费a级高清毛片| 久久成人国产精品免费软件| 亚洲成网777777国产精品| 国产网站免费看| 国产大片黄在线观看| 麻豆精品在线视频| 好吊色妇女免费视频免费| 中文字幕天无码久久精品视频免费 | 国产在线观看精品| 国产午夜福利在线小视频| 午夜在线不卡| 老司机午夜精品视频你懂的| 香蕉蕉亚亚洲aav综合| 亚洲精品国产自在现线最新| 色综合综合网| 婷婷色中文| 日本国产精品| 无码久看视频| 国产女人在线| 精品国产自| 亚洲第一极品精品无码| 毛片免费试看| 国产成人精品午夜视频'| 久久精品女人天堂aaa| 97视频免费在线观看| 午夜国产精品视频黄| 日韩色图在线观看| 久久大香香蕉国产免费网站| hezyo加勒比一区二区三区| 久久精品国产免费观看频道| 国产精品成人啪精品视频| 91无码人妻精品一区| 97精品伊人久久大香线蕉| 亚洲系列无码专区偷窥无码| a毛片在线免费观看| 人妻少妇乱子伦精品无码专区毛片| 欧洲精品视频在线观看| 亚洲精品不卡午夜精品| 青青操视频在线| 亚洲一区毛片| 亚洲人成色在线观看| 91在线一9|永久视频在线| 亚洲av无码牛牛影视在线二区| 伊伊人成亚洲综合人网7777| 香蕉久人久人青草青草| 亚洲国产精品日韩专区AV| 亚洲国产精品一区二区第一页免 | 欧美日韩精品综合在线一区| 亚洲无码视频喷水|