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

多元Al合金的熱力學和熱物理數據庫的建立及凝固過程顯微組織演變的模擬

2010-10-24 09:20:18徐洪輝劉樹紅張利軍趙冬冬李一為王培生金展鵬黃伯云
中國材料進展 2010年6期
關鍵詞:擴散系數界面數據庫

杜 勇,徐洪輝,孔 毅,劉樹紅,張利軍,趙冬冬,李一為,王培生,金展鵬,黃伯云

(中南大學粉末冶金國家重點實驗室,湖南長沙410083)

多元Al合金的熱力學和熱物理數據庫的建立及凝固過程顯微組織演變的模擬

杜 勇,徐洪輝,孔 毅,劉樹紅,張利軍,趙冬冬,李一為,王培生,金展鵬,黃伯云

(中南大學粉末冶金國家重點實驗室,湖南長沙410083)

介紹了相圖熱力學及熱物理性能數據庫的建立方法;綜述了國內外學者構筑多元Al合金相圖熱力學和熱物理性能數據庫的研究工作。并介紹了使用相場方法,使用熱力學和熱物理性能數據庫模擬多元Al合金凝固時顯微組織演變的幾個實例。最后展望了新一代熱力學和熱物理性能數據庫建立及凝固過程顯微組織演變定量模擬未來的可能發展方向。

相圖熱力學;熱物理性能;相場方法;顯微組織;凝固

前 言

大多數工業材料都是多組元、多相合金。這些材料的性能很大程度上取決于其微觀結構,即在凝固和后續熱處理過程中所生成相的種類和含量、枝晶間距、微觀偏析等。弄清楚怎樣通過選擇合金成分和工藝參數來控制多相微觀組織的形成,是材料領域一個富有挑戰性的研究方向。獲得液相和固相精確的熱力學和熱物理性能數據是描述凝固過程微觀組織演變的必備條件。主要的熱物理性能數據包括擴散系數、界面能、熱導率、密度等,它們都是溫度和成分的函數。

不論在宏觀還是微觀尺度,凝固過程是模擬計算領域一個越來越受關注的研究方向。不論從哪個尺度出發,都必須先知道凝固體系關鍵的熱力學和熱物理性能數據。目前,由于缺乏可靠的關鍵數據,特別是液相擴散系數、固液界面能等熱物理性能數據,凝固過程模擬的發展受到了很大限制。這種狀況對于多組元體系來說更是如此。此外,由于缺乏微觀結構的實驗數據,模擬結果往往得不到驗證,這也阻礙了凝固過程模擬的發展。至今為止,對多組元合金微觀結構參數的定量了解相當有限。

近十多年來,相場方法廣泛用于模擬微觀結構演變,該方法已成功應用于凝固、晶粒長大和粗化、薄膜體系的相變等領域[1-4]。相場方法是一種基于經典熱力學和動力學的唯象學模型。其總體思路是用場變量來構造體系的總能量(包括塊體化學自由能、界面能、彈性應變能、磁性能等),而體系的微觀結構演變是總能量的最小化過程。相場方法使用的場變量分為二種:一類是整個演變過程中守恒的物理量(如成分);另一類是非守恒的物理量(如長程有序參數)。該方法通過求解成分場的Cahn-Hilliard非線性擴散方程及長程有序參數的Ginzburg-Landau方程來描述微觀結構的時空演變規律。同其它方法相比,相場方法在描述微觀結構演變時具有不少優點。由于采用擴散相界面(Diffuse Interface),即通過場變量的梯度變化來表征相界面,因而可對各種復雜的微觀結構進行二維和三維模擬。德國科學家Steinbach等[3]開發了世界上第一個基于相場方法的商業程序MICRESS(MICRostructureEvolutionSimulation Sof tware),并主要用于凝固過程時微觀結構演變的模擬。

本文研究內容包括:①多元Al合金熱力學和熱物理性能數據庫的建立;②使用相場方法模擬多元Al合金凝固過程顯微組織演變的幾個實例;③新一代熱力學和熱物理性能數據庫建立及凝固過程顯微組織演變定量模擬未來的可能發展方向。

多元合金相圖熱力學數據庫

. 相圖熱力學研究的主要實驗方法

利用熱力學數據庫和相關熱力學計算軟件可以獲取多組元體系相平衡、亞穩相平衡、仲平衡、各種熱力學性質、熱力學因子等信息。多組元熱力學數據庫涵蓋各相熱力學模型及其熱力學參數。通常,為了優化出這些熱力學參數,需要輸入相圖、化合物生成焓、溶體混合焓和活度等相圖熱力學實驗數據。高精度的熱力學數據庫源于合理的熱力學模型和精確可靠的相圖熱力學實驗數據。下面將簡要介紹相圖測定的主要實驗方法。

相圖的實驗測定方法分為靜態法和動態法兩種,前者有淬冷法、擴散偶法、X射線結構分析法等,后者有熱分析法、膨脹法、電阻法等。相圖的精確測定需要多種方法綜合使用。有關相圖實驗測定方法的書籍和文獻較多[5-7],故在此只簡單介紹最常用的靜態法中的淬冷法和動態法中的熱分析法。另外,還將簡要介紹測定化合物生成焓和溶體混合焓的量熱法。

淬冷法 靜態法是把一系列已知成分的試樣在特定溫度下進行恒溫熱處理,使試樣達到平衡狀態或局部平衡狀態,然后測定試樣的某些性質在該溫度下隨成分的依賴關系。淬冷法是靜態法中應用最廣泛的方法。其主要原理是將選定的不同組成的試樣長時間在一系列預定的溫度下保溫,使它們達到對應溫度下的平衡狀態,然后快速淬冷試樣。由于相變來不及進行,淬冷后的試樣保持了高溫下的平衡狀態。用光學顯微鏡(OM)、掃描電子顯微鏡(SEM)、電子探針微區分析(EPMA)和X射線衍射(XRD)等方法對淬冷試樣進行顯微組織、顯微成分和物相分析,可確定給定溫度下的平衡相成分及相平衡關系,從而繪制出相圖。

從目標溫度淬冷到室溫的過程中,試樣的高溫平衡狀態能否得到保持,是影響測定結果的重要因素。近年來,高溫顯微鏡和高溫X射線衍射技術的發展,使得直接觀察試樣在高溫下的相平衡狀態成為可能。

熱分析法 相圖測定中最常用的熱分析方法是差熱分析和差示掃描量熱分析。差熱分析(Differential ThermalAnalysis,DTA)是在程序控制溫度下測定試樣和參比物(在測量溫度范圍內不發生任何熱效應的物質)之間的溫度差和溫度關系的一種技術[8]。差示掃描量熱分析(Differential Scanning Calorimetry,DSC)是在程序控制溫度下,測量輸入到試樣和參比物的能量差隨溫度或時間變化的一種技術。和差熱分析一樣,DSC主要用于測定液相線溫度、液相面的投影圖、固相線以及固態相變的溫度。此外,DSC可以直接測量熱量,這是DSC與DTA的一個重要區別。

一般情況下,相圖的精確測定需要多種方法配合使用,靜態法和動態法相結合。諸如,首先制備一系列已知成分的合金,均勻化處理后,用DTA、DSC、熱膨脹等動態法確定相變溫度(適合變溫截面、液相面投影圖等的測定),用XRD、顯微鏡、EPMA、金相法等靜態法確定等溫截面上的相界位置。

量熱法 測定熱力學性質的方法包括量熱法、電動勢法、等靜壓法等。所測量的熱力學性質包括形成焓、室溫時的熵、熱容和相變焓等。量熱法是應用最為廣泛的一種測量熱力學性質的方法,幾乎所有熱力學性質都可由量熱法測量。量熱計(Calorimeter)是一種測量反應熱的裝置,可分為恒容式和恒壓式量熱計。形成焓是吉布斯自由能的主要部分,高溫量熱法成為測定形成焓的主要技術手段,該技術所用的測試儀器從溶液量熱法發展到高溫反應量熱法。Kleppa[9]從Calvet型量熱計[10]改進發明了高溫反應量熱計,該儀器與傳統的Setaram測試原理相同,但前者的測試溫度可以超過1 400 K。Kleppa量熱計是利用間接法測量化合物在室溫下的形成焓:首先測定由室溫的元素在某一溫度生成化合物的熱量變化,然后測定化合物由室溫加熱到同一高溫的熱量變化;兩種熱量之差即為化合物的形成焓。

. 相圖計算方法和第一性原理方法

20世紀70年代以來,隨著熱力學基礎理論與計算機技術的發展,由Kaufman,Hillert和Ansara等人奠基的熱化學與相圖計算耦合研究,發展成為一門涵蓋熱化學、相圖熱力學與計算技術的交叉學科CALPHAD(CALculation of PHAse Diagrams)[11-12]。1977年,國際上相圖計算領域唯一的專業刊物CALPHAD雜志創刊。CALPHAD旨在開發熱力學模型和熱力學優化計算軟件,評估實驗數據并建立自洽的數據庫,進而通過熱力學計算預示多組元材料的性質,從而推進計算熱力學的發展。

目前,國際上已有多個成熟的商業化相圖計算軟件,例如Thermo-Calc[13]、Factsage[14]和Pandat[15]等。Thermo-Calc軟件是由瑞典皇家工學院在Hillert,Sundman,Jansson等人的工作基礎上,于1981年推出的相圖熱力學計算軟件。經過近30年的發展,Thermo-Calc軟件現已成為數據齊全、功能強大、結構較為完整的計算系統,目前是世界上享有相當聲譽的熱力學計算軟件。

使用相圖計算方法即CALPHAD技術可以高效構筑多元相圖。但CALPHAD技術依賴于有限的相平衡及熱力學數據(特別是二元及三元系的數據)優化熱力學模型中的參數。第一性原理(First-Principles Calculation)計算方法僅僅需要原子序數以及晶體結構信息作為輸入,便可獲得穩定相及亞穩相的熱力學數據(如生成焓、熱容等)。很多在CALPHAD方法中需要的參數無法由實驗直接獲得,但是可以直接從第一性原理計算獲得。

近幾年來,國內外學者在結合第一性原理計算和CALPHAD技術研究領域做了很多工作。將第一性原理計算的晶體結構、熱力學和相圖信息輸入到熱力學數據庫可以擴大CALPHAD方法預測多元合金的熱力學性質和相圖的應用范圍[16-18]。

. 合金相圖熱力學研究介紹

鋁合金數據庫COST507是歐共體上世紀90年代開發的商用數據庫。由于其使用的很多實驗數據過時、欠準確,因此用COST507數據庫計算出的相平衡經常不準確。近6年來,中南大學粉末冶金國家重點實驗室通過實驗和計算相結合的方法建立了多組元Al基體系熱力學數據庫,并圍繞多元Al基合金凝固及隨后時效強化過程的微結構演變模擬開展了系統研究。

張利軍等[18]利用量熱法、第一性原理計算和CALPHAD方法對Al-Fe-Ni體系在整個成分和溫度范圍內相圖熱力學性質進行了系統研究。他們通過反應量熱法和第一性原理計算獲得三元化合物的生成焓。借助第一性原理計算獲得了亞點陣模型中的三元端際化合物和L12有序相的能量。所計算的能量有助于隨后的熱力學建模和熱力學優化。他們還修正了多元系fcc-Al/L12有序-無序轉變的熱力學模型,計算預測出的三元bcc-B2有序相的溶解度隙與實驗值吻合得很好。圖1是用第一性原理計算的Al-Fe-Ni系三元化合物τ1相和τ2相的生成焓和熵與溫度的關系。圖2是計算的垂直截面相圖與實測數據的比較,圖a,b,c和d分別為:80%Al,60%Al,50%Al和75%Ni(均原子分數)。

圖1 第一性原理計算的Al-Fe-Ni系三元化合物τ1相和τ2相生成焓(a)及τ1相和τ2相的熵(b)與溫度的關系Fig.1 Effect of finite temperature on enthalpy of formation(a)and entropy(b)calculated via fist-principle for theτ1andτ2 phases in the Al-Fe-Ni system

多元合金熱物理性能數據庫

. 多元合金的擴散系數矩陣

微觀結構演變過程中的很多現象,如凝固過程樹枝晶的生長、平衡相或亞穩相在基體上的析出、相的長大和消失等,都同擴散現象密切相關。

擴散理論表明:為描述含n個組元相的擴散行為,需要一個由(n-1)2個擴散系數Dij(i=1…n-1,j=1…n-1)組成的矩陣。為減少描述多元相的擴散所需的參數個數,Anderson等[27]提出了由實測擴散系數獲得各元素在相中的原子移動性參數,然后再用原子移動性參數計算擴散系數矩陣的方法。根據這一思想,他們開發了D ICTRA軟件。D ICTRA軟件是通過原子移動性和由熱力學數據庫獲得的熱力學因子來描述各種擴散系數隨溫度和成分的變化規律。最近,Zhang等[28]利用D ICTRA軟件描述了Al-Ni體系中有序和無序相的擴散系數。通過考慮L12(γ′)有序相的成分范圍及缺陷濃度,他們建立了能顯著減少描述γ′相原子移動性參數個數的擴散動力學模型。他們所獲得的原子移動性參數可以很好地描述各種擴散系數的實驗數據(圖3),并首次準確描述了氣/固擴散偶中與時間相關的濃度曲線(圖4)。

圖2 計算的Al-Fe-Ni系垂直截面相圖與實驗值的比較:(a)80%Al,(b)60%Al,(c)50%Al和(d)75%Ni(均原子分數)Fig.2 Comparison between calculated vertical section phase diagram and experimental values from the literature forAl-Fe-Ni system:(a)80%Al,(b)60%Al,(c)50%Al,and(d)75%Ni(all atomic percentage)

圖3 計算的Ni在L12相中的雜質擴散系數同實驗數據的比較(M代表一個常數,加入此常數使得不同成分下的曲線分離,以便比較。)Fig.3 Comparison between calculated tracer diffusion coefficients ofNi in L12 phase and the experimental data

由于原子移動性參數只有通過對實測擴散系數進行優化獲得,因此測定擴散系數是很有意義的工作。但測定擴散系數的各種實驗方法都有局限性,因此利用第一性原理計算等方法預測擴散系數越來越受到國內外學者的關注。第一性原理計算自擴散系數是基于空位機制(包括單空位和雙空位機制)來計算自擴散系數。對于fcc結構的純金屬自擴散系數一般采用公式(1)來描述。

圖4 模型預測的Al/Ni氣/固擴散偶的濃度曲線與實驗數據的比較Fig.4 Comparison between model-predicted concentrationcurves of vapor/solid Al/Ni and experimental data

其中,f為關聯因子(對于fcc體系為0.781 5),α為晶格常數,C為空位濃度,ω為振動頻率。使用第一性原理可以準確地計算出晶格常數,空位濃度以及振動頻率,從而得到自擴散系數。

空位濃度C一般采用公式(2)來描述

振動頻率ω的一般表述形式為:

其中ΔSm為空位遷移熵,ΔHm為空位遷移能,v為有效振動頻率,結合TST(Transition State Theory過渡態理論),第一性原理能夠準確地預測振動頻率。

雜質擴散系數(fcc體系)的第一性原理計算是基于含有稀薄雜質原子的fcc體系的空位機制進行的。對于fcc體系來講,5振動頻率模型(Five Jump Frequency)是一個很好描述雜質擴散系數的模型,其表示式為:

式中,D2為體系的雜質擴散系數,D0為純fcc體系的自擴散系數,f2為雜質擴散系數的關聯因子,f0為自擴散系數的關聯因子,ω0-4為體系的5個振動頻率,其中ω0為自擴散系數的振動頻率。結合TST,第一性原理可以準確地預測fcc體系的雜質擴散系數。

互擴散系數的第一性原理計算也是基于空位機制。根據擴散理論,擴散系數可表示為:

其中,L為動力學因子,Θ為熱力學因子。由Kubo-Green機理可知,

由D的矩陣對角化可以得到互擴散系數。通過式(6)可以得到熱力學和動力學因子從而預測互擴散系數。

圖5 Al自擴散系數的第一性原理計算與實驗結果比較(GGA,LDA分別是Generalized Gradient Approximation,LocalDensity Approximation的簡寫)Fig.5 Comparison between first-principle calculation for selfdiffusion coefficients of pure-Al and experimental results

最近,Mantina等[38]使用第一性原理十分準確地預測了純Al的自擴散系數,得到的計算結果可以很好地預測實驗數據(圖5)。此前,他們還報道了用第一性原理預測雜質擴散系數的方法[44]。應用該方法,他們計算了Mg,Si及Cu在稀薄fcc相Al合金中的雜質擴散系數(圖6)。Ven和Ceder[56]報道了用第一性原理計算空位控制的二元合金擴散系數的方法,并成功應用于Al-Li二元系,獲得了Al-Li二元系的互擴散系數(圖7)。

. 界面能

除相圖熱力學、擴散系數外,界面能是描述多元多相材料微觀結構演變必不可少的重要信息。例如時效強化過程材料強度的增量同析出相的形貌密切相關,而析出相同基體間的界面能和彈性應變能決定了析出相的形貌。用實驗方法測定界面能的難度很大,目前測定金屬固-液界面能的常用方法是基于均質形核理論的形核過冷度法和基于Gibbs-Thomson關系式的晶界法。形核過冷度法的準確性取決于所測的過冷度是不是均質形核的過冷度。但長期以來人們一直沒有解決如何判定一個形核過程是均質形核過程,還是異質形核過程的問題。晶界法測定固-液界面能時,因其包含的環節很多,故很容易造成系統誤差。由于以上兩個方面的原因,目前金屬固-液界面能實驗測定數據相差較大。

圖6 計算的Mg(a)和Si(b)在fcc Al中的雜質擴散系數與實驗數據的比較Fig.6 Comparison between calculated diffusion coefficient of magnesium(a)and silicon(b)in fcc Al and experimental data

圖7 Al-Li互擴散系數的第一性原理計算(600 K)Fig.7 First-principles calculated interdiffusion coefficient of the Al-Li system at 600 K

近10年來,從界面結構出發,采用第一性原理或嵌入原子方法[57]計算界面能得到了國內外學者的廣泛關注[58-59]。最近,魯曉剛等根據實測的析出物的體積和平均尺寸隨時間的變化曲線,通過在一維框架內描述析出物的形核、長大及粗化過程得到了鋼鐵材料中較準確的界面能數據[60]。劉梓葵等[61]應用第一性原理超胞方法和密度泛函理論研究了亞穩的β″-Mg5Si6相和α-Al間的界面性質,給出了界面能,晶格錯配度及應變能的定量值。當由上述不同方法得到的界面能接近時,可以認為所計算的界面能是準確的。

. 其它熱物理性能

其它熱物理性能,如熱導率和液相密度也是影響凝固過程的重要因素。例如,熱導率和熱擴散決定相變界面間的熱傳遞速率,熱傳導速率越快,凝固過程越容易實現。

很多測定熱導率的實驗方法主要集中在研究單相試樣的性能。最近,Lamvik[62]描述了一種適合研究金屬固/液界面處熱導率的方法。這種方法耦合定向凝固/溶化過程,通過獲得界面傳播速率,界面處的溫度梯度及相變潛熱等數據來估算固/液相在轉變溫度處的熱導率。

液相密度的計算主要根據Sung等提供的多元合金密度計算公式[63]。該方法不僅考慮液相中各個元素的體積隨溫度的變化,而且也考慮不同原子尺寸的各種元素混合在一起后產生的混合體積,該理論計算的結果能夠準確地反映液相密度隨溫度的變化,并且該公式計算結果與多種工業用高溫合金的實際結果較吻合,計算結果的誤差一般在±2.5%。

多元合金凝固過程的組織演變模擬

系統的穩定狀態由吉布斯自由能決定,而系統的動態微觀組織演變,除滿足熱力學條件(即能量最小原理外),還需要考慮微觀粒子的動力學行為。因而獲得精確的多元Al合金的相圖熱力學數據,各相之間的界面能,以及各組元在合金熔體中的擴散系數等熱物理數據后,可以預測多元Al合金凝固過程中的顯微組織以及這些顯微組織在時間和空間上的動態演化過程。在顯微組織演變的動態模擬中,除了需要熱力學和熱物理性能數據外,還需要選擇合適的動力學模型來描述微觀組織演化的動力學規律。

在現有的各種微觀組織演變的動力學模型中,相場動力學模型是非常成功的一種,并在多元合金凝固過程的顯微組織演變模擬方面獲得了廣泛的應用。本節將首先簡介相場方法的基本原理,然后介紹使用該方法模擬多元Al合金凝固過程中顯微組織演變的幾個實例。

. 相場方法簡介

相場動力學模型是利用與位置和時間相關的場變量來描述任意的微結構,并認為微結構演化的動力學過程是基于能量最小原理,通過相場動力學方程中總自由能的降低來實現[64]。不同的場變量,滿足不同的相場動力學方程。比如對于微結構演變過程中的守恒量,滿足Cahn-Hilliand方程[65]:

式中,c為演化過程中的守恒量,比如濃度等。t為時間。M是與原子移動性有關的參量。Ftot是指的系統總的熱力學自由能函數。而微結構演變過程中的非守恒量,則滿足Allen-Cahn方程[66]:

式中,η為演化過程中的不守恒量,比如有序度,晶粒取向等。L是與界面能等有關的參量。在上面兩個方程中,需要知道材料的熱力學與熱物理性能參數,包括體系的自由能變化,界面能,擴散系數等。這些熱力學與熱物理參數,可以通過動態的與現有的熱力學軟件,比如ThermoCalc等耦合,來獲得不同溫度與成分下所研究體系合理的熱力學與熱物理性能數據[67]。當上述參數已知的情況下,多元合金微結構的演化可看成是由相場動力學方程來控制。

根據對相場方程解法的不同,人們通常把相場模型分為連續模型和微觀模型兩大類。這兩種模型都是基于Ginzburg-Landao理論[68]派生的方法。微觀模型與連續模型的主要區別在于場變量的不同。微觀相場模型是用原子占據晶格位置的幾率來描述原子組態和相形貌,該相場模型由Khachaturyan創建[69],并由Chen等人做了發展[70]。目前應用較廣泛的多元合金連續相場模型主要有3種:WBM模型[71-72],KKS模型[73]和Steinbach模型[74]。WBM模型保持與熱力學一致,并假定固/液界面是由濃度相同的固、液相的混合,從而使得模型中會增加一個額外的雙勢阱;而KKS模型,同樣保持與熱力學一致,但假定平衡時固/液界面是由化學勢相同的固、液相的混合,該模型還假定合金熔體為理想溶液,從而忽略了溶質濃度的影響,因而主要適合于稀二元合金凝固過程的模擬。Steinbach模型則通過純物質與合金之間變量的匹配,直接將純物質的相場模型擴展為合金的相場模型,從而實現了多元多相合金微結構的動態模擬。

. 三元合金的等溫枝晶粗化

王錦成等人[75]采用相場方法研究了Ni-Al-Nb三元合金的枝晶粗化問題。通過研究不同比例的固相與液相初始條件下的微結構演化,他們發現當固/液比從64%變化到88%時,模擬得到的微結構將從連接的盤狀過渡到桿狀形貌,這與實驗觀測到的共晶反應中由于界面能的變化所引起的層狀形貌過渡到桿狀形貌的變化一致。對相場模擬結果的進一步分析表明,在相對較低的固/液比下,枝晶等溫粗化主要是由第3臂的再熔化,聚集和平滑化等效應主導。而在相對較高的固/液比下,聚集,平滑,瑞利失穩,小液滴的圓滑化與收縮等效應將主導枝晶的粗化過程。在模擬過程中還發現,增加液相的擴散系數將有利于瑞利失穩的發生,從而加速液相從扁平的盤狀形貌到圓柱狀形貌的轉變。圖8給出了Ni-Al-Nb三元合金固/液比和界面比例隨時間的變化關系。圖9給出了相應于圖8曲線上A,B,C時刻(tA=0.056 s,tB=0.124 s,tC=0.600 s)的微結構圖。

圖8 Ni-Al-Nb三元合金凝固過程的固/液比和界面比例隨時間的變化Fig.8 Varialion of solid-liquid ratio and interfacial proportion with time during solidification for ternary alloysNi-Al-Nb

圖9 圖8曲線上A,B,C時刻(tA=0.056 s,tB=0.124 s,tC=0.600 s)的微結構圖Fig.9 Selected patterns of microstructural evolution during solidification for ternary alloysNi-Al-Nb,a,b,c corresponding to timesA,B,and C in figure 8 respectively(tA=0.056 s,tB=0.124 s,andtC=0.600 s)

. 熱力學與動力學數據庫耦合相場方法模擬多元合金凝固

介萬奇等人[76]基于相場方法,通過集成熱力學與擴散動力學數據庫,模擬了Al-Cu-Mg多元合金的枝晶生長過程。在多元合金的相場模擬中,通過引入界面中固相與液相具有相同的化學勢,但成分不同的假設,特別考慮了不同類的溶質與溶質之間的相互作用。在模擬過程中,通過與ThermoCalc軟件的動態連接,在相場模擬中直接使用優化的熱力學和動力學數據庫,能夠及時準確的預測出相平衡與溶質擴散。對于二維的Al-Cu-Mg三元合金枝晶生長,模擬得到了定量的固相和液相的溶質分布與擴散矩陣。計算表明擴散系數與成分密切相關。圖10給出了二維的Al-Cu-Mg三元合金枝晶生長過程中溫度為900 K,時間為6×10-5秒時的溶質分布形貌圖。從圖中可以看出,在枝晶的主干里溶質濃度較低,而在第二枝晶臂之間則有最高的溶質分布。這一現象與實驗觀測是一致的。

圖10 Al-2%Cu-3.5%Mg合金(原子分數)枝晶生長過程中溫度為900 K,時間為6×10-5s時的溶質分布形貌圖Fig.10 Solute concentration profiles by simulation for the dendrite growth process of Al-2%Cu-3.5%Mg alloy(atomic percentage)at 900 K(t=6×10-5s):(a)Cu and(b)Mg

. 再加熱過程中枝晶的碎塊化

枝晶在凝固過程可能會發生碎塊化,從而影響微結構。這樣的微結構演化在再加熱熱處理定向生長時往往會發生。因而對于枝晶碎塊化的研究有助于加深人們對單晶的定向生長過程中微結構從柱狀向等軸狀轉變,以及缺陷生成過程的理解。Boettinger等人[64]對C2H4(CN)2-C10H16O采用相場方法,在模擬過程中將體系溫度瞬時從1 574 K加熱到1 589 K,同時將無量綱的過飽和度從0.86降低到0.25。通過施加這樣劇烈的熱起伏,他們的模擬結果表明枝晶將快速停止生長,同時枝晶臂前段有部分重新熔化而出現碎塊化。在模擬過程中如果降低熱起伏的程度,這樣的碎塊化將不會發生。圖11給出了C2H4(CN)2-C10H16O枝晶的熔化與碎塊的形成過程。其中圖a給出的是溫度為1 574 K時的微結構;圖b,c,d給出的是當體系溫度瞬時提高到1 589 K后,微結構的變化情況。從圖中明顯可以看出在施加了劇烈的熱起伏后,枝晶的再熔化與碎塊化過程。用實驗方法觀察Al合金在凝固過程中枝晶可能發生的碎塊化,現階段還很困難。

. .多元合金定向凝固后顯微組織的測定和模擬

杜勇等[77]使用所建立的熱力學及熱物理數據庫計算模擬了多元Al合金(Al356.1)在冷卻速度為2 K/s定向凝固條件下的顯微組織及成分偏析,并把模擬結果同實驗結果進行了比較。在模擬計算時,考慮了固相、液相中的擴散、過冷及二次樹枝晶的長大對凝固過程的影響。圖12所示為Al356.1合金定向凝固后樣品橫截面的金相組織。經EPMA成分分析,樣品中存在的相有(Al),(Si),α-A lMnSi和β-AlFeSi。采用BSE自動分析得到(Al),(Si),α-A lMnSi和β-AlFeSi相的體積分數分別為0.93,0.051,0.009和0.01。所測的各相體積分數分別同計算預測結果((Al):0.93,(Si):0.048,α-A lMnSi:0.007,β-AlFeSi:0.014)吻合得相當好。

圖11 C2H4(CN)2-C10H16O枝晶的熔化與碎片的形成模擬圖:(a)為1 574 K時的微結構,(b),(c),(d)為溫度瞬時提高到1 589 K后的微結構Fig.11 Simulated diagram ofmelting of dendritic structure and formation of fragments:dendritic micro structure at growth temperature of 1 574 K(a)and change of dendritic micro structure at temperature instantly elevated to 1 589 K(b),(c),(d)for C2H4(CN)2-C10H16O

圖12 合金樣品Al356.1橫截面的金相照片(冷卻速度2 K/s)Fig.12 Metallograph of transverse section of the Al 356.1 alloys directionally solidified with cooling rate of 2 K/s

. 液相擴散系數對合金(柱狀晶向等軸晶轉變)的影響

由于文獻中合金液相擴散系數的實驗報道較少,相場模擬中通常假設液相的擴散系數與溫度成分無關,即恒為10-9m2/s。然而,這種簡化處理會在一定程度上影響定量描述凝固微觀組織結構的準確性。最近,張利軍,杜勇等[78]在綜合考慮各種實驗測定和理論預測擴散系數的基礎上,建立了一套精確的Al-Ni合金液相的原子移動性數據庫。通過所建立的熱力學、固相擴散動力學及熱物理性能數據庫,張利軍,杜勇等[78]采用MICRESS相場軟件[79],系統研究了液相擴散系數對Al-Ni合金CET(即:柱狀晶向等軸晶轉變)的影響。圖13所示為所模擬的Al-1.9%Ni(質量分數)合金在溫度梯度為100 K/cm及冷卻速率為1 K/s的情況下的CET的組織形貌演變過程。圖13 a,13 b中除液相擴散系數外,其它的熱物理參數均保持一致。圖13 a中液相的擴散系數由所獲得的原子移動性數據庫提供,而圖13 b中液相的擴散系數恒為10-9m2/s。如圖13所示,液相擴散系數的不同對Al-1.9%Ni合金CET的轉變位置及組織形貌的影響較大。

圖13 相場模擬的Al-1.9%Ni合金CET的組織形貌演變(溫度梯度100 K/cm,冷卻速率1 K/s,(a)和(b)中除液相擴散系數外其它的熱物理參數均保持一致):(a)液相擴散系數由所獲得的原子移動性數據庫提供,(b)液相擴散系數恒為10-9m2/sFig.13 The phase-filed simulated micro structure evolution ofCET in the Al-1.9%Ni alloy.All the ther mophysicalparameters are kept as the same except for the liquid diffusivity.The liquid diffusivity is calculated from the atomic mobility database in(a),while that is assumed to be 10-9m2/s in(b)

總結與展望

現實世界廣泛使用的金屬材料絕大多數都是多組元多相體系。建立精確的相圖熱力學及熱物理性能數據庫是提升傳統材料性能及設計出新型材料的關鍵所在。由于多組元多相體系的復雜性,只有采用多種理論和實驗方法才有可能構筑其準確的相圖熱力學及熱物理性能數據庫。用實驗方法測定關鍵的相圖和熱力學數據,并通過相圖計算方法與第一性原理計算方法的結合是構筑多組元多相體系在寬廣成分和溫度范圍內熱力學數據庫的有效途徑。類似地可通過實測關鍵擴散系數,并通過宏觀模擬方法(如擴散控制相變的原子移動性D ICTRA模擬方法)和微觀計算方法(如分子動力學、第一性原理計算等)的結合構筑多相體系的擴散系數矩陣。同相圖熱力學及擴散動力學數據的構筑相比,界面能等其它熱物理性能數據庫的獲取任重而道遠。

目前國內外學者對重要工業材料平衡相的熱力學性質及熱物理性能所做工作相對較多。與之相比,對非平衡相或亞穩相所做工作很少。今后研究工作的重點應是獲得工業應用上重要相的熱力學及熱物理性能數據,逐步建立含亞穩相的多元多相體系豐富的相圖熱力學及熱物理性能數據庫,實現多元多相材料制備過程微結構演變的定量模擬,滿足新材料設計和制備的需求。

[1]Emmerich H.Advances of and by Phase-Field Modelling in Condensed-Matter Physics[J].Advances in Physics,2008,57(1):1-87.

[2]Moelans N,Blanpain B,Wollants P.An Introduction to Phase-Field Modeling of Microstructure Evolution[J].CALPHAD,2008,32:268-294.

[3]Steinbach I.Phase-Field Models in Materials Science[J].Modelling and Simulation inMaterials Science and Engineering,2009,17:1-31.

[4]Chen L Q.Phase-Field Models for Microstructure Evolution[J].Annu RevMater Res,2002,32:113-140.

[5]Xueshan Lu(陸學善).PhaseD iagram and Phase Transfor mation(相圖與相變)[M].Unviversity of Science and Technology of China Press,1990.

[6]Zhang Shengbi(張圣弼),Li Daozi(李道子).Phase D iagram-Principle,Calculation and Their the Applications in Metallurgy(相圖-原理、計算及在冶金中的應用)[M].Beijing:Metallurgical Industry Press,1986.

[7]Zhao J C.Methods for Phase Diagram Deter m ination[M].Oxford:Elsevier,2007.

[8]Wang Peiming(王培銘),Xu Qianwei(許乾慰).Methods for Material Research(材料研究方法)[M].Beijing:Science Press,2005.

[9]Kleppa O J.Systematic Aspects of the High-Temperature Ther mochemistry of Binary Alloys and Related Compounds[J].Journal of Phase Equilibria,1994,15:240-63.

[10]Feutelais Y,Legendre B,GuymontM,et al.Standard Enthalpy of Formation ofAl0.28Fe 0.72 at298 K[J].Journal of Alloys and Compounds,2001,322(1/2):184-189.

[11]Kaufman L,Bernstein H.Computer Calculation of Phase D iagrams[M].New York:Academic Press,1970.

[12]SaundersN,Miodownik A P.CALPHAD(Calculation of Phase Diagrams)-A Comprehensive Guide[M].Great Britain:Elsevier Science Ltd,1998.

[13]Anderson J O,Helander T,Hoglund I,et al.Ther mo-Calc and Dictra,Computational Tools forMaterials[J].CALPHAD,2002,26(2):273-312.

[14]Bale C W,Chartrand P,Degterov S A,et al.FactSage Ther moche mical Software and Databases[J].CALPHAD,2002,26(2):189-228.

[15]Chen S L,Daniel S,Zhang F,et al.The PANDAT Software Package and Its Application[J].CALPHAD,2002,26(2):175-188.

[16]Jiang C,Gleeson B.A Combined First-Principles/CALPHAD Modeling of the Al-Ir System[J].Acta Materialia,2006,54(15):4 101-4 110.

[17]Golumbfskie W J,Prins S N,Eden T J,et al.Predictions of the Al-Rich Region of the Al-Co-Ni-Y System Based upon First-Principles and Experimental Data[J].Calphad,2009,33(1):124-135.

[18]ZhangL J,WangJ,Du Y,et al.Thermodynamic Propertiesof the Al-Fe-Ni System Acquired via a Hybrid Approach Combining Calorimetry,First-Principles and CALPHAD[J],ActaMaterialia,2009,57(18):5 324-5 341.

[19]Chumak I,Richter KW,IpserH.The Fe-Ni-Al Phase Diagram in the Al-Rich(>50%Al)Corner[J].Inter metallics,2007,15:1 416-1 424.

[20]Zhang L J,Du Y.Thermodynamic Description of the Al-Fe-Ni System over theWhole Composition and Temperature Ranges:Modeling Coupled with Key Experiment[J].CALPHAD,2007,31:529-540.

[21]Bitterlich H,L?serW,SchultzL.Reassessment ofAl-Ni and Ni-Fe-Al Solidus Temperatures[J].J Pha Equilib,2002,23(4):301-304.

[22]Bradley A J.Microscopical Studies on the Iron-Nickel-Aluminium System.Part I.α+βAlloys and Isother mal Section of the Phase Equilibrium Diagram[J].J Iron Steel Inst,1949,163(1):19-30.

[23]Masahashi N,Kawazoe H H,Takasugi T,et al.Phase Relations in the Section Ni3Al-Ni3Fe of the Al-Fe-Ni System[J].Z Metallkd,1987,78(11):788-794.

[24]Jia C C,Ishida K,Nishizawa T.Phase Equilibria in Ni Rich Region of Ni-Al-Fe System[J].Metall Mater Trans A,1994,25(4):473-281.

[25]Goman'KovV I,Tret'yakova S M,Monastyrskaya E V,et al.Structural Diagrams of Quasi Binary Alloys Ni3Fe-Ni3Al,Ni3Mn-Ni3Al,and Ni3Mn-Ni3Ga[J].Russ Metall.1998,6:125-131.

[26]Himuro Y,Tanaka Y,Kamiya N,et al.Stability of Ordered L12 Phase in Ni3Fe-Ni3X(X:Si and Al)PseudobinaryAlloys[J].Inter metallics,2004,12(6):635-643.

[27]Andersson J-O,?gren J.Models for Numerical Treatment of Multicomponent Diffusion in Simple Phases[J].J Appl Phys,1992,72(4):1 350-1 354.

[28]ZhangL J,Du Y,Chen Q,et al.Atomic Mobilities and Diffusivities in the fcc,L12 andB2 Phasesof the Ni-Al System[J].Int J Mater Res,in press(2010).

[29]Hancock G F.Diffusion of Nickel in Alloys Based on the Intermetallic Compound Ni3Al[J].Phys Stat SolA,1971,7:535-540.

[30]Shi Y,Frohberg G,Wever H.Diffusion of63Ni and114mIn in theγ′-Phase Ni3Al[J].Phys Stat Sol A,1995,152(2):361-375.

[31]Hoshino K,Rothman S J,Averback R S.Tracer Diffusion in Pure and Boron-Doped Nickel-Aluminide(Ni3Al)[J].Acta Metall,1988,36:1 271-1 279.

[32]Bronfin MB,Bulatov G S,Drugova I A.Self-Diffusion of Nickel in the Nickel-Aluminum(Ni3Al)Inter metallide and in Pure Nickel[J].FizMetalMetalloved.1975,40:363.

[33]Larikov L N,Geichenko V V,Fal'schenko V M.Diffusion Processin Ordered Alloys[M].Kiev:Naukova Dumka Press,1975.

[34]Frank S,S?dervall U,Herzig Chr.Self-Diffusion of Ni in Single and Polycrystals of Ni3Al.A Study of S IMS and Radiotracer Analysis[J].Phys Stat Sol B,1995,191(1):45-55.

[35]Cserháti C,SzabóIA,Márton Z S,et al.Tracer Diffusion of 63Ni in Ni3(Al,Ge)Ternary Inter metallic Compound[J].Inter metallics,2002,10:887-909.

[36]Campbell C E.Private Communication[M].USA:NIST,2009.

[37]Yamamoto T,Takashima T,Nishida K.Inter diffusion in the Zeta-Solid Solution of a Ni-Al System[J].Trans Jpn InstMetals,1980,21:601-608.

[38]MantinaM,Wang Y,Arroyave R,et al.First Principle Calculation of Self-Diffusion Coefficients[J].PRL,2008,100:215-901.

[39]Fradin F Y,Rowland T J.Nuclear Magnetic Resonance Measurement of the Diffusion Coefficient of Pure Aluminum[J].Appl Phys Lett,1967,11(6):207-209.

[40]Allen P S,Andrew E R,Bates C A.Proceeding of the18th Ampere Congress Nottingham,England,1974[C].Amsterdam,North-Halland:This Scientific Conference Press,1975:327.

[41]Lundy T S,Murdock J F.Diffusion of Al26 and Mn54 in Aluminum[J].J Appl Phys,1962,33:1 671-1 673.

[42]Stoebe T G,Gulliver R D,Ogurtani T O,et al.Nuclear Magnetic Resonance Studies of Diffusion of 27Al in Aluminum and Aluminum Alloys[J].Acta Metall,1965,13(7),701-708.

[43]Volin T E,Balluffi R W.Annealing Kinetics of Voids and the Self-Diffusion Coefficient in Aluminum[J].Phys Status Solidi,1968,25(1):163-173.

[44]MantinaM.,Wang Y,Chen L Q,et al.First Principle I mpurity Diffusion Coefficients[J].ActaMater,2009,57:4 102-4 108.

[45]Fujikawa S I,Takada Y.Interdiffusion between Aluminum and Al-MgAlloys.Diffusion and Defect Data-Solid State Data Part A[J].Defect D iffus Forum,1997,143-147:409-414.

[46]Minamino Y,Yamane T,Miyake T,et al.Effect of High Pressure on Diffusion Reactions and Phase diagram in Aluminum-Magnesium System[J].Mater Sci Technol,1986,2(8):777-783.

[47]Rothman S J,Peterson N L,NowickiL J,et al.Tracer Diffusion of Magnesium in Aluminum Single Crystals[J].Phys Status Solidi B,1974,63(1):29-33.

[48]Fujikawa S,Hirano K.Diffusion of Magnesium-28 in Aluminum[J].Mater Sci Eng,1977,27(1):25-33.

[49]Hisayuki K,Yamane T.Diffusion of Magnesium in Commercial Al-Mg Alloy under High Pressure[J].Z MetaIlkd,1999,90:423.

[50]Moreau G,CornetJ A,CalaisD.Acceleration ofChemical Diffusion during Irradiation in the Aluminum-Magnesium System[J].J NuclMater,1971;38:197-202.

[51]Verlinden J,Gijbbels R.I mpurity Diffusion in Aluminum as Determined from Ion-Probe Mass Analysis[J].Adv Mass Spectrom A,1980,8:485-495.

[52]Fujikawa S,Hirano K-I,Fukushima Y.Diffusion of Silicon in Aluminum[J].Metall Trans A,1978,9(12):1 811-1 815.

[53]Beerwald A H.Diffusion ofVarious Metals in Aluminum[J].Z Elektrochem Angew Phys Chem,1939,45:789-795.

[54]BergnerD,Cyrener E.Diffusion of Foreign Elements in Aluminum Solid Solutions,Part 3:Investigations Into the Diffusion of Silicon in Aluminum Using the Microprobe[J].Neue Hutte,1973,18:356-361.

[55]Mehl R H,Rhines F N,Vonden Steinen KA.Diffusion in Alpha Solid Solutions of Aluminum[J].Metals and Alloys,1941,13:41-44.

[56]Van derVen V,Ceder G.First Principles Calculation of the Inter diffusion Coefficient in Binary Alloys[J].PRL,2005,94:045 901.

[57]Daw MS,Baskes MI.Embedded-Atom Method:Derivation and Application to Impurities,Surfaces,and Other Defects in Metals[J].Phys Rev B,1984,29:6 443-6 453.

[58]ZhangW,Smith J R,Evans A G.The connection between ab Initio Calculations and Interface Adhesion Measurements on Metal/Oxide Systems:Ni/Al2O3and Cu/Al2O3[J].Acta Mter,2002,50:3 803-3 816.

[59]Hu S Y,BaskesMI,Stan M,et al.Atomistic Calculations of Interfacial Energies,Nucleus Shape and Size ofθ′Precipitates in Al-Cu Alloys[J].Acta Mater,2006,54:4 699-4 707.

[60]Lu X G.An Effective Approach to Estimate the Interfacial Energy[C]//Chang YW.The Sixth Pacific Rim International Conference on AdvancedMaterials and Processing.Jeju Island,Korea:Scientific Conference Press,2007:5-9.

[61]Wang Y,Liu Z K,ChenL Q,et al.First-Principle Calculations ofβ″-Mg5Si6/α-Al Interfaces[J].Acta Mater,2007,55:5 934-5 947.

[62]Lamvik M,Zhou J M.A Novel Technique for Measuring the Ther mal Conductivity of Metallic Materials during Melting and Solidification[J].Meas Sci Technol,1995,6:660-887.

[63]WangLing(王 玲),Reseasch of Effect of Solidification Parameters on Segregation Behaviour and Stability of Pasty Region for Superalloy(凝固參數對高溫合金偏析行為和糊狀區穩定性的影響研究)[D].Beijing:University of Science and Technology Beijing,2007.

[64]BoettingerW J,Warren J A,Beckermann C,et al.Phase-Field Simulation of Solidification[J].Annu RevMater Res,2002,32:163-194.

[65]Cahn J W,Hilliard J E.Free Energy of a Nonunifor m System:Interfacial Free Energy[J].J Chem Phys,1958,28:258-267.

[66]Allen S M,Cahn J W.A Microscopic Theory for Antiphase Boundary Motion and Its Application to Antiphase Domain Coarsening[J].Acta Metall,1979,27:1 085-1 095-1 979.

[67]Suzana G Fries,Bernd Boettge,Janin Eiken,et al.Upgrading CALPHAD to Microstructure Simulation:the Phase-Field Method[J].Int J Mat Res,2009,100:128-134.

[68]GinzburgV L,Landau L D.On the Theory of Superconductivity[J].J Exp Theoret Phys,1950,20,1 064.

[69]Khachatuyran A G.Theory of Structural Transfor mations in solids[M].New York:W iley,1983:129.

[70]Chen L Q,Khachatuyran A G.Computer Simulation of Structural Transfor mation during Precipitation of an OrderedIntermetallic Phase[J].Acta MetallMater,1991,39:2 533.

[71]Wheeler A A,Boettinger W J,McFadden G B.Phase-Field Model for Isothermal Phase Transitions in Binary Alloys[J].Phys Rev A,1992,45:7 424.

[72]Wheeler A A,McFadden G B,Boettinger J B.Phase-Field Model for Solidification of a Eutectic Alloy[J].Proc R Soc London SerA,1996,452:495.

[73]Kim S G,Kim W T,Suzuki T.Phase-Field Model for Binary alloys[J].Phys Rev E,1999,60(6):7 186.

[74]Tiaden,Nestler B,Diepers H J,et al.The Multiphase-Field Model with an Integrated Concept for Modeling Solute Diffusion[J].Physica D,1998,115:73-86.

[75]Wang J C,Yang G.Phase-Field Modeling of Isothermal Dendritic Coarsening in Ternary Alloys[J].Acta Mater,2008,56:4 585-4 592.

[76]Zhang R J,Jing T,Jie W Q,et al.Phase-Field Simulation of Solidification in Multicomponent Alloys Coupled with Ther modynamic and Diffusion Mobility Databases[J].Acta Mater,2006,54:2 235-2 239.

[77]Du Y,Chang YA,Liu S H,et al.Thermodynamic Description of the Al-Fe-Mg-Mn-Si System and Investigation of Microstructure and Microsegregation during Directional Solidification of an Al-Fe-Mg-Mn-SiAlloy[J].Z Metallkd.2005,96(12):1 351-1 362.

[78]Zhang L J,Du Y,Steinbach I,et al.Diffusivities of the Al-Fe-Ni Melt and Their Effectson the Microstructure during Solidification[J].Acta Mater,2010,58(10):3 664-3 675.

[79]http://www.access.r wth-aachen.de/MICRESS/.

Thermodynam ic and Thermo physical Databases of Multi-Component Al Alloys and Simulation on Microstructural Evolution during Solidification

Du Yong,Xu Honghui,Kong Yi,Liu Shuhong,Zhang Lijun,Zhao Dongdong,Li Yiwei,W ang Peisheng,Jin Zhanpeng,Huang Baiyun

(State Key Laboratory of powder metallurgy,Central South University,Changsha 410083,China)

A hybrid approach to establish thermodynamic and thermo physical databases is described,followed by a brief summary on the establishment of thermodynamic and thermo physical databases of multi-component Al alloys.A few of case studies to demonstrate the microstructural evolution of Al alloys during solidification by means of phase field method using the knowledge of both thermodynamic and thermo physical properties are described.Future directions in the establishment of next generation of ther modynamic and thermo physical databases aswell as the quantitative microstructural evolution ofmulti-component and multi-phases alloys are addressed.

ther modynamics;thermo physical properties;phase field method;microstructural evolution;solidification

TG111.3,TG111.4

A

1674-3962(2010)06-0028-12

2009-12-09

國家自然科學基金重點項目(50831007)

杜 勇,男,1964年生,教授,博士生導師

猜你喜歡
擴散系數界面數據庫
國企黨委前置研究的“四個界面”
當代陜西(2020年13期)2020-08-24 08:22:02
基于FANUC PICTURE的虛擬軸坐標顯示界面開發方法研究
數據庫
財經(2017年2期)2017-03-10 14:35:35
人機交互界面發展趨勢研究
數據庫
財經(2016年15期)2016-06-03 07:38:02
數據庫
財經(2016年3期)2016-03-07 07:44:46
數據庫
財經(2016年6期)2016-02-24 07:41:51
基于Sauer-Freise 方法的Co- Mn 體系fcc 相互擴散系數的研究
上海金屬(2015年5期)2015-11-29 01:13:59
FCC Ni-Cu 及Ni-Mn 合金互擴散系數測定
上海金屬(2015年6期)2015-11-29 01:09:09
手機界面中圖形符號的發展趨向
新聞傳播(2015年11期)2015-07-18 11:15:04
主站蜘蛛池模板: 亚洲天堂网视频| 久久综合AV免费观看| 欧美日韩国产精品综合| 亚洲VA中文字幕| 精品人妻系列无码专区久久| 午夜免费视频网站| 国产一区二区三区在线无码| 免费观看三级毛片| 亚洲免费播放| 国产日本欧美亚洲精品视| 狠狠v日韩v欧美v| 伊人无码视屏| 久久久久免费看成人影片 | 亚洲—日韩aV在线| 国产成人久久综合777777麻豆| 人妻丝袜无码视频| 制服丝袜一区二区三区在线| 亚洲欧美在线综合一区二区三区 | 欧美日韩动态图| 亚洲国产一区在线观看| 欧美成人国产| 亚洲va视频| 激情六月丁香婷婷| 欧美亚洲一区二区三区在线| 亚洲免费黄色网| 成人免费黄色小视频| 国产流白浆视频| 岛国精品一区免费视频在线观看| 综合天天色| 精品无码国产自产野外拍在线| 天堂成人在线| 一区二区三区在线不卡免费| 亚洲首页在线观看| 日韩无码黄色| 爱做久久久久久| 无码精品国产VA在线观看DVD| 国产极品美女在线观看| 久久精品娱乐亚洲领先| 日韩精品成人在线| 国产在线拍偷自揄拍精品| 久久99国产综合精品女同| 国产精品欧美激情| 国产91麻豆免费观看| 三级国产在线观看| 日韩东京热无码人妻| 综合成人国产| 免费无码又爽又黄又刺激网站| 一本久道久综合久久鬼色| 国产精女同一区二区三区久| 国产成人综合日韩精品无码首页| 毛片免费在线| 成人福利在线视频| 国产网站免费| 欧美国产日产一区二区| 日韩高清成人| 色婷婷综合在线| 激情六月丁香婷婷四房播| 亚洲精品无码专区在线观看| 国产成人亚洲欧美激情| 久久精品丝袜高跟鞋| 国产一区自拍视频| 亚洲综合精品第一页| 久久久久亚洲AV成人网站软件| 人妻免费无码不卡视频| 欧美日韩v| 亚洲成人动漫在线观看| 九色91在线视频| 91麻豆精品国产高清在线| 国产男女免费视频| 天天综合网在线| 伊人色在线视频| 午夜激情婷婷| 久久情精品国产品免费| 国产高清国内精品福利| 九色最新网址| 日韩美一区二区| 亚洲高清无在码在线无弹窗| 国产成人麻豆精品| 国产成人高清精品免费5388| 99久视频| 99re精彩视频| 免费av一区二区三区在线|