亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        高Cr 鑄鐵中M7C3 碳化物與奧氏體共生長的元胞自動機模擬*

        2021-11-19 05:16:22張山張紅偉苗淼馮苗苗雷洪王強
        物理學報 2021年21期
        關鍵詞:元胞溶質碳化物

        張山 張紅偉? 苗淼 馮苗苗 雷洪 王強

        1) (東北大學,材料電磁過程研究教育部重點實驗室,沈陽 110819)

        2) (東北大學冶金學院,沈陽 110819)

        高鉻鑄鐵中M7C3 碳化物大小適中、彌散均勻分布,有利于提高合金的耐磨性.為分析凝固過程中M7C3 碳化物晶粒在基體中的形貌及分布、M7C3 碳化物與奧氏體晶粒生長的相互作用、引起的溶質偏聚對最終M7C3 碳化物粒徑分布的影響,本文開發(fā)了Fe-C-Cr 三元合金小面晶M7C3 碳化物與奧氏體晶粒共生長的二維微觀元胞自動機模型,模型中加入潛熱釋放對凝固過程溫度場的影響,由C,Cr 兩溶質界面擴散共同確定晶體生長速度,由凝固路徑數(shù)據(jù)表插值獲取液相元胞的溶質平衡濃度,設定M7C3 碳化物鄰胞結構并優(yōu)化形狀因子來保持M7C3 碳化物小面晶形貌,模擬了Fe-4%C-17%Cr 三元合金(C 和Cr 的質量分數(shù)分別為4%和17%)初生M7C3 碳化物和共晶奧氏體晶粒的生長演變過程.研究表明,M7C3 碳化物和奧氏體晶粒各自的生長速度隨著界面液相中C,Cr 溶質的超飽和度和貝克列數(shù)的增大而增大;隨著奧氏體的析出和晶粒生長,M7C3 碳化物晶粒的生長速度明顯增快;當奧氏體晶粒逐漸接觸并包圍M7C3 碳化物晶粒時,兩相晶粒生長速度逐漸降低.凝固過程中,奧氏體晶粒生長向外排出C,Cr 溶質,與吸收C,Cr 溶質生長的M7C3 碳化物晶粒互補,致使二者生長互相促進,最終奧氏體晶粒包圍M7C3 碳化物晶粒生長.預測的冷卻曲線與實驗冷卻曲線變化趨勢相符;最終凝固組織形貌和M7C3 碳化物體積分數(shù)與實驗相符;剩余液相、奧氏體中C,Cr 溶質濃度演變也與Gulliver-Scheil,Partial Equilibrium,Lever Rule 模型預測結果相符.

        1 引言

        高鉻鑄鐵被廣泛應用于發(fā)電廠、礦山、煉鋼廠的高耐磨零部件中[1],其優(yōu)異的耐磨性主要來自增強相M7C3碳化物的高硬度[2-5].M7C3碳化物大小適中,彌散均勻分布于基體中,有利于提高基體的韌性與合金的耐磨性.學者們[6-8]針對高鉻鑄鐵凝固過程中Cr/C 比、冷卻速率等對M7C3碳化物、奧氏體的數(shù)量、尺寸、分布等展開了大量實驗研究,本文作者[9]也預測了高鉻鑄鐵的凝固路徑,以期改進材料的耐磨性.對于凝固過程中M7C3小面晶碳化物在基體中的形貌分布受何影響,M7C3碳化物與奧氏體枝晶生長的相互作用,引起的溶質偏聚對最終M7C3碳化物粒徑分布作用規(guī)律的相關數(shù)值模擬未見他人報道.

        晶粒生長及形貌演變的數(shù)值模擬方法中,微觀元胞自動機 (microscopic-scale cellular automata,MCA)方法能夠耦合溫度場、溶質場、流場等各物理場演變過程,有效描述凝固過程中枝晶內部形貌的演化特征,相比于相場法,MCA 法能夠模擬較大尺寸材料的組織演變,計算耗時較短,因而被廣泛應用于模鑄、連鑄、定向凝固、增材制造及焊接等過程的組織模擬中[10,11].

        目前,采用微觀元胞自動機方法進行單相枝晶生長的模擬已發(fā)展至三元以上合金[12-18],為確定多元合金中枝晶的生長速度,一種做法是采用界面動力學系數(shù)與過冷度的乘積來計算,其界面平衡固相成分、平衡液相線溫度通過熱力學計算引擎(如PanEngine、Thermo-Calc 等)來獲得.基于此,Zhu等[12]預測了Al-Cu-Mg 合金中單一枝晶相的微觀結構;戴挺[13]等進行了Al-Cu-Mg-Si 四元合金枝晶生長模擬.此法未體現(xiàn)多元合金溶質間擴散對生長速度的相互作用關系.再一種做法是,生長速度采用局部溶質平衡濃度與實際濃度的差值求得.如Zhu 等[19]針對Al-Si 二元合金共晶組織演變的模型,以及張蕾等[20]針對球墨鑄鐵(Fe-C 二元合金)中石墨與奧氏體離異共晶生長的模型.此方法多針對二元合金.第三種做法是基于界面生長速度與多溶質生長貝克列(Pe)數(shù)對應關系,通過耦合各溶質的超飽和度方程和界面溶質守恒方程求得界面生長速度,可明確給出多元合金各溶質間擴散對生長速度的相互作用.其中界面溶質濃度對應的液相線溫度及液相線斜率等熱力學參數(shù)由熱力學平衡計算獲得.基于此,本文作者模擬了Fe-C-Cr三元合金凝固過程中奧氏體[14]及鐵素體+奧氏體包晶相[15]的生長速度,石玉峰等[16]模擬了Al-Si-Mg 三元合金枝晶生長形貌.另外,Michelic 等[17]等基于元胞自動機法改進的VFT (virtual fronttracking)模型,引入外插函數(shù)求解界面兩側溶質平衡濃度,再利用界面溶質守恒方程求得生長速度,模擬了Fe-C-Si-Mn-P 五元合金枝晶的生長過程.在枝晶和第二相晶粒共生長的模擬中,Zhu等[19]模擬了Al-Si 二元合金共晶組織演變;張蕾等[20]模擬了球墨鑄鐵(Fe-C 二元合金)中石墨與奧氏體依循離異共晶模式的協(xié)同競爭生長,并將石墨晶粒處理成球狀.本課題組張晛韋等[18]模擬了亞共晶球墨鑄鐵(Fe-C-Si 三元合金)中石墨與奧氏體的離異共晶生長,其石墨和奧氏體的生長速度采用上述第三種做法.

        本研究中,M7C3碳化物為密排六方(HCP)小面晶結構[21,22],在二維截面中呈現(xiàn)六邊形形狀.對于小面晶生長過程中形狀的保持,Stefan-Kharicha等[23]采用的前沿跟蹤模型,分別計算生長面法向和切向速度,并通過保持晶角位置隨鄰面生長移動來保持晶面形狀,模擬了熱浸鍍鋅過程中六邊形面晶(Fe2Al5)的生長.付振南等[24]通過特別規(guī)定元胞鄰居關系及界面單元捕獲規(guī)則,模擬了AZ91D鎂合金中密排六方結構枝晶的生長.董祥雷等[25]在相場模型中引入小面晶各向異性修正方程,模擬了各向異性穩(wěn)態(tài)六方GaN 螺旋結構形貌.綜上,目前僅有學者針對其他體系中相近的六方小面晶形貌進行了模擬,未見針對Fe 基三元合金中M7C3碳化物生長,尤其是M7C3碳化物與奧氏體共生長模擬的報道.

        本課題組前期工作[26]中,采用微觀元胞自動機模型結合傳熱傳質過程,引入界面溶質貝克列數(shù)和超飽和度方程來求解界面多溶質平衡濃度,模擬了高Cr 鑄鐵(Fe-C-Cr 合金)凝固過程中M7C3碳化物和奧氏體晶粒的協(xié)同生長過程.本文進一步加入了潛熱釋放影響.將液相元胞溶質平衡濃度更新,由原采用的Fe-C 二元合金液相溶質平衡濃度公式,替換為通過Gulliver-Scheil (GS)模型結合熱力學平衡計算預測的Fe-C-Cr 三元合金凝固路徑獲得;將相變界面溶質擴散系數(shù)由線性插值改為調和插值,更合理反映相變界面擴散特性;并優(yōu)化了M7C3碳化物形狀修正因子.模擬了Fe-4%C-17%Cr合金(C 和Cr 的質量分數(shù)分別為4%和17%,下同)凝固過程中M7C3碳化物和奧氏體晶粒形貌演變及濃度遷移,深入分析了M7C3碳化物和奧氏體共生長過程中溶質貝克列數(shù)、超飽和度與界面溶質平衡濃度、生長速度之間關系,預測的冷卻曲線和最終凝固組織形貌分別與實驗進行了對照,并將剩余液相和奧氏體中C,Cr 溶質濃度演變規(guī)律與凝固路徑預測結果進行了對比.

        2 數(shù)學模型

        建立Fe-4%C-17%Cr 三元合金M7C3碳化物和奧氏體共生長演變過程二維微觀元胞自動機模型.

        2.1 基本假設

        1)不考慮熔體的流動;

        2)忽略動力學過冷;

        3)忽略溶質間互擴散;

        4)固/液界面始終處于平衡狀態(tài);

        5)忽略M7C3碳化物在奧氏體中的生長;

        6)熱擴散的尺度比溶質擴散高3—4 個數(shù)量級,全場溫度均勻變化.

        2.2 溫度場

        由于熱量傳遞速率遠大于溶質擴散速率,認為溫度場均勻變化.考慮凝固過程中潛熱釋放的影響,溫度場控制方程如下:

        2.3 溶質擴散方程

        在計算區(qū)域內引入勢函數(shù)P(等效液相成分),把整個計算區(qū)域處理為單相.統(tǒng)一的溶質擴散方程為

        方程(2)離散時,當兩個相鄰元胞有1 個為M7C3碳化物或奧氏體時界面擴散系數(shù)Di/e相應取DM,i或Dγ,i,否則通過調和插值計算界面擴散系數(shù):

        式中,下標e 為兩元胞界面位置,下標P 和E 為中心元胞和其鄰居元胞;每個元胞中的溶質擴散系數(shù)為其中各相的加權平均值,即Di=fγDγ,i+fMDM,i+fLDL,i;fe為界面e 距E 元胞中心距離與P 和E元胞中心間距之比,這里網格均勻,fe=1/2.

        為平均溶質濃度:

        式中,fL為液相質量分數(shù),CM,i為對應的M7C3碳化物溶質濃度.濃度場采用無擴散的邊界條件:.

        2.4 形 核

        模型中,為與實驗[22]結果對照,M7C3碳化物、奧氏體晶粒的形核位置和形核數(shù)量參照實驗形貌設定.實際形核數(shù)量為形核密度、形核過冷度等的函數(shù)關系[27],還受鑄錠尺寸影響,同時由形核Gaussian 分布可知,增大形核過冷度會減小形核數(shù)量.實驗分析表明[28],增大冷卻速率,一次、二次枝晶臂尺寸減小,說明形核數(shù)量隨冷卻速率增大而增大.由GS 模型凝固路徑預測,Fe-4%C-17%Cr合金中M7C3碳化物和奧氏體晶粒的析出溫度分別是1304 ℃和1266 ℃,形核過冷度隨機設定.當元胞局部過冷度大于形核過冷時,該元胞形核,單元狀態(tài)即刻由液態(tài)轉變?yōu)榻缑嬖?并被賦予1 個隨機的晶體學取向θ.基于奧氏體為四重對稱的枝晶狀,奧氏體晶粒的取向范圍設為[—π/4,π/4];而M7C3碳化物晶粒在二維截面具有六重對稱性,其晶粒取向范圍設為[—π/6+2,π/6+2],此處+2 是為了區(qū)分M7C3碳化物與奧氏體晶粒而設置.

        2.5 生長動力學

        2.5.1 界面熱力學平衡

        界面元胞滿足局部熱力學平衡條件,局部溫度與液相線溫度之間滿足(5)式:

        式中,T是界面胞溫度;TL,k是液相線溫度(下標k指析出M 或γ相); ΔTR和 ΔTC分別是曲率過冷和成分過冷.

        2.5.2 曲率過冷

        曲率過冷表達式為

        式中,Γ為固液界面Gibbs-Thomson 系數(shù),M7C3碳化物和奧氏體分別對應ΓM,Γγ;K為界面曲率,采用計數(shù)法計算[29];f(φ,θ) 為界面各向異性函數(shù).根據(jù)Gibbs-Thomson 方程,界面各向異性函數(shù)如下:

        式中,k為晶體對稱相關數(shù),對于面心立方晶體,如奧氏體取k=4,對于密排六方晶體,如M7C3碳化物,取k=6;ε為界面能各向異性系數(shù),界面各向異性大于臨界值后會在固液界面上出現(xiàn)晶向缺失現(xiàn)象,因此界面能各向異性系數(shù)ε須滿足,否則會出現(xiàn)晶向缺失現(xiàn)象[30],導致凝固組織處于熱力學不穩(wěn)定狀態(tài),對于奧氏體,ε=0.04,對于M7C3碳化物,ε=0.025.θ為晶體取向;φ為晶體界面法線方向與參考坐標軸(x軸)的夾角,由法向定義計算[26].

        2.5.3 成分過冷

        界面前沿的成分過冷與界面實際溶質濃度及平衡溶質濃度有關,Fe-C-Cr 三元合金需考慮C和Cr 兩種溶質濃度對于成分過冷的影響:

        式中,mL/k,i為液相線斜率(下標中k指M,γ相),為溶質C 的界面平衡液相濃度,為溶質Cr 的界面平衡液相濃度.

        2.5.4 界面前沿溶質超飽和度

        對于溶質C 和Cr 而言,不考慮溶質之間的互擴散[31],界面前沿溶質超飽和度滿足

        設奧氏體/液相界面溶質分配系數(shù)為常數(shù),則

        M7C3碳化物晶粒中溶質濃度CS,i設為常數(shù),則M7C3碳化物/液相界面的液相平衡濃度:

        Pei為溶質i貝克列(Peclet)數(shù),表達為

        式中,Vn為晶粒生長速度,R為晶粒生長半徑.

        在相界面處,生長速度唯一,PeC和PeCr存在如下關系:

        通過(8)—(13)式,可確定成分過冷 ΔTC與Cr 溶質貝克列數(shù)PeCr間單一函數(shù)關系.

        2.6 界面法向速度

        M7C3碳化物和奧氏體晶粒生長過程中,由固液界面處溶質質量守恒,得到界面元胞的法向生長速度:

        界面速度分量為

        式中μk為動力學系數(shù)(下標k表示M/L,γ/L 界面).Burbelko 等[32]采用動力學系數(shù)進行生長速度的計算,μγ/L=10—3,μgr/L=10—8.Zhu 等[12]通過比較其模擬與GS 模型模擬的固相分數(shù)隨溫度變化差異確定界面動力學系數(shù)μγ/L=3×10—4.孫玉成[33]引入 [1-(fS-π/6)]n(n為經驗常數(shù))對共晶團晶粒生長速度進行修正,以解決在模擬計算過程中晶粒長大后其晶粒度預測失準的問題.本研究通過與實驗觀測到M7C3碳化物晶粒的大小進行對比,對M7C3碳化物和奧氏體晶粒的生長速度進行動力學修正,動力學系數(shù)分別為μM/L=4.8×10—4,μγ/L=2.6×10—4.

        2.7 更新固相分數(shù)

        晶體的固相分數(shù)增量與固-液界面推進速度有關,可用(17)式表示:

        式中 Δs為元胞單元尺寸(單位m),Δx=Δy=Δs.

        對于奧氏體晶粒,為體現(xiàn)界面擾動的因素,將固相分數(shù)增量乘上1 個擾動函數(shù),即1+A×[1-2Rand(·)],其中A代表隨機擾動振幅,取值0.1,Rand(·)為[0,1]之間的隨機數(shù).

        對于M7C3碳化物晶粒,為了減小網格各向異性對M7C3碳化物晶粒形貌的影響,以及減少M7C3碳化物晶粒內不完全凝固點的問題,將M7C3碳化物晶粒的固相分數(shù)增量乘G,G是與鄰位網格狀態(tài)有關的形狀因子,其表達式如下:

        圖1 是不同系數(shù)a所對應的同一時刻同一位置下碳化物晶粒質量分數(shù)圖.可以看出隨著系數(shù)a的不斷增加,M7C3碳化物晶粒內不完全凝固點的數(shù)量逐漸減小,當a=0.50 時,M7C3碳化物晶粒內完全凝固;然而當a≥ 0.50,隨著系數(shù)a的不斷增加,M7C3碳化物晶粒的六邊形外輪廓偏向弧度化.因此本模擬最終選取a=0.40.

        圖1 t=99.99 s 時,不同系數(shù)a 條件下模擬的M7C3 質量分數(shù)圖 (a) a=0.25;(b) a=0.40;(c) a=0.50;(d) a=0.60Fig.1.M7C3 morphology in form of mass fraction with several values of coefficient a at t =99.99 s:(a) a=0.25;(b) a=0.40;(c) a=0.50;(d) a=0.60.

        在用(17)式求出元胞內的固相分數(shù)增量后,更新元胞新時刻的固相分數(shù):

        2.8 更新固相濃度

        對于奧氏體晶粒內溶質固相濃度,用(20)式更新:

        由GS 凝固路徑預測模型計算,凝固過程中C,Cr溶質在M7C3碳化物中濃度變化較小,因此,本研究將M7C3碳化物晶粒內溶質固相濃度設為常數(shù).

        2.9 M7C3 碳化物與奧氏體元胞鄰居關系

        通過實驗觀測分析,Fe-4%C-17%Cr 合金中的M7C3碳化物晶粒在二維平面上理想狀態(tài)具有六重對稱的特征,但實際合金凝固過程中受到一些凝固缺陷等因素的影響,二維平面觀測的結果有的為規(guī)則六邊形,有的為近似六邊形等狀態(tài)[34-36].據(jù)其形貌特征[24],本研究定義了CA 方法中M7C3碳化物晶粒的二維鄰居單元如圖2 所示,其中黑色方框為某一CA 單元,周圍緊鄰的16 個單元為其鄰居.奧氏體晶粒采用Moore 型鄰居單元,包括最近鄰的8 個元胞,如圖3 所示.

        圖2 M7C3 元胞鄰居單元Fig.2.Neighborhood relations for M7C3 grain.

        圖3 γ 元胞Moore 型鄰居單元Fig.3.Moore neighborhood relations for austenite grain.

        2.10 晶粒捕獲規(guī)則

        模擬區(qū)域中每個元胞共有3 種狀態(tài):液相、固相及界面(固液相區(qū))狀態(tài).捕獲規(guī)則定義如下:液相元胞形核后,此元胞為界面狀態(tài);當該元胞生長至固相分數(shù)為1 時,此元胞為固相狀態(tài),周圍鄰居元胞即時轉變?yōu)榻缑嬖?當某一元胞同時被奧氏體和M7C3碳化物元胞捕獲時,該元胞變?yōu)閮上喙泊娴慕缑姘?且變?yōu)楣滔嘣笾荒芨髯圆东@周圍8 個鄰居元胞;其余界面元胞成長為固相元胞后,按奧氏體或M7C3碳化物各自的鄰居關系捕獲周圍鄰居元胞,逐漸生長.

        2.11 溶質再分配

        在生長界面,M7C3碳化物晶粒從界面液相中吸收C,Cr 溶質(M/L 界面分配系數(shù)均大于1);奧氏體晶粒生長向界面液相中排出C,Cr 溶質(γ/L界面分配系數(shù)小于1),進而引起界面前沿液相溶質重新進行分配.

        奧氏體界面胞液相的溶質濃度為

        2.12 時間步長

        式中 Δt0是初始時間步長.

        2.13 計算步驟

        1)所有元胞均賦予相同的初始溫度和初始濃度,計算時間步長;

        2)溫度場考慮潛熱釋放,溫度緩慢下降((1)式);

        3)當達到M7C3碳化物或奧氏體的形核過冷度,該元胞轉變?yōu)榻缑嬖?

        4)計算M7C3碳化物和奧氏體界面各向異性函數(shù)f(φ,θ)((7)式),求得曲率過冷 ΔTR((6)式);

        5)通過(5)式,計算界面元胞成分過冷 ΔTC;

        6)通過(8)—(13) 式,求解界面元胞PeCr,再求得PeC和界面平衡濃度;

        7)液相元胞的平衡濃度按對應溫度,通過GS凝固路徑預測模型計算的數(shù)據(jù)表插值獲得;

        8)求解界面法向生長速度((1)—(16) 式)并更新固相分數(shù)((17)—(19) 式);

        9) M7C3碳化物和奧氏體界面胞進行局部溶質再分配((21)—(22)式);

        10)計算擴散引起的濃度場變化((2)—(4)式);

        11)再重復步驟2)—10)進行時間層循環(huán),直到全場凝固結束,計算中止.

        表1 列出了 Fe-4%C-17%Cr 三元合金模擬所用的物性參數(shù)[22,37,38].

        表1 Fe-4%C-17%Cr 三元合金模擬所用的物性參數(shù)(單位%是指質量分數(shù))Table 1.Physical properties used for Fe-4%C-17%Cr ternary alloy.Unit of % represents mass fraction (wt%).

        3 模擬結果與分析

        本研究利用微觀元胞自動機方法研究二維空間Fe-4%C-17%Cr 三元合金中初生相M7C3碳化物和奧氏體晶粒共生長過程.選取的計算域參照實驗[22]顯微組織圖片尺寸,確定為700 μm × 520 μm,劃分為350 × 260 個均勻正方形網格,網格尺寸為2 μm × 2 μm.形核位置和形核數(shù)量參照實驗[22]結果進行設置.熔體初始溫度參照實驗冷卻曲線的起始溫度設定為1334.32 ℃,冷卻速率參照石墨鑄型[22]實驗冷卻曲線選取,分為4 個階段:M7C3碳化物析出且獨立生長階段(T> 1266 ℃),冷卻速率為4.6 K/s;奧氏體析出初始階段(1220 ℃ <T≤1266 ℃),冷卻速率為1.4 K/s;奧氏體和M7C3碳化物潛熱釋放引起模擬區(qū)域溫度上升階段(1194 ℃ <T≤ 1220 ℃),冷卻速率為2.4 K/s;接近凝固結束潛熱作用較小階段(T≤ 1194℃),冷卻速率為2.0 K/s.隨著溫度的不斷下降,M7C3碳化物和奧氏體晶粒依次形核生長,直至模擬區(qū)域達到最終凝固狀態(tài).

        本文模擬的Fe-4%C-17%Cr 合金凝固冷卻曲線,將與文獻[22]中同成分添加1.5%Ti (質量分數(shù))合金及文獻[39]中Fe-3.23%C-23.8%Cr 合金(C 和Cr 的質量分數(shù)分別為3.23%和23.8%)添加4%Ti與否的凝固冷卻曲線進行對比.

        3.1 合金凝固路徑

        表2 列出了上述二種成分合金由GS 模型預測和在Fe-C 偽二元相圖中的析出相(M7C3碳化物、奧氏體、滲碳體)及對應析出溫度.可以看出,同一成分GS 模型預測和Fe-C 偽二元相圖給出的各相析出趨勢一致:Fe-4%C-17%Cr 合金中,添加Ti 元素后,M7C3碳化物析出溫度均變低、奧氏體析出溫度變高;二者給出的各相析出溫度差較小.由表2 可知,添加Ti 元素后,至固相線或Cementite(CEM)相析出,Fe-4%C-17%Cr-1.5%Ti 體系,TiC體積分數(shù)為3.00%,奧氏體體積分數(shù)為63.04%,M7C3體積分數(shù)為26.38%;Fe-3.23%C-23.8%Cr-4%Ti 體系,TiC 體積分數(shù)為7.92%,奧氏體體積分數(shù)為74.28%,M7C3體積分數(shù)為17.80%.因而冷卻曲線的拐點和溫度回復主要反映奧氏體和M7C3相的析出及共生長演變.

        表2 GS 模型和Fe-C 偽二元相圖中二種成分Fe-C-Cr 合金析出相及析出溫度Table 2.Phase type and precipitation temperature in two Fe-C-Cr alloys by GS prediction and in Fe-C pseudo binary phase diagram.

        由圖4 中Fe-C 偽二元合金相圖可知,添加Ti元素后,共晶點右移,導致Fe-4%C-17%Cr 合金(圖4(c)—(d))中M7C3碳化物析出溫度降低、奧氏體析出溫度升高.但從凝固路徑來看,未改變M7C3→γ析出順序,即由無Ti 添加的L→M7C3→γ→CEM,至添加Ti 的L→TiC→M7C3→γ→CEM.但是,對于文獻[39]實驗合金成分Fe-3.23%C-23.8%Cr,由于離共晶點較近,添加Ti 后,相變路徑已移至共晶點左側,析出順序由無Ti 添加的L→M7C3→γ→CEM,變?yōu)樘砑覶i 后,奧氏體先于M7C3相析出的L→TiC→γ→M7C3→Solidus.但合金內主要析出相仍為M7C3和γ.對比分析圖4(a)與圖4(b)可知,Fe-4%C-17%Cr 合金與Fe-3.23%C-23.8%Cr 均為過共晶三元合金,析出相M7C3與奧氏體及順序均相同,冷卻曲線拐點和溫度回復應相近.

        圖4 Fe-C 偽二元合金相圖(A 區(qū),液相;B 區(qū),液相+TiC;C 區(qū),液相+M7C3;D 區(qū),液相+M7C3+γ;E 區(qū),M7C3+γ;F 區(qū),液相 +TiC+M7C3;G 區(qū),液相+TiC+γ;H 區(qū),液相+TiC+M7C3+γ;I 區(qū),TiC+γ+M7C3) (a) Fe-C-17%Cr;(b) Fe-C-23.8%Cr;(c) Fe-C-17%Cr-1.5%Ti;(d) Fe-C-23.8%Cr-4%TiFig.4.Fe-C pseudo binary phase diagram:(a) Fe-C-17%Cr;(b) Fe-C-23.8%Cr;(c) Fe-C-17%Cr-1.5%Ti;(d) Fe-C-23.8%Cr-4%Ti.A zone,Liquid;B zone,Liquid+TiC;C zone,Liquid+M7C3;D zone,Liquid+M7C3+γ;E zone,M7C3+γ;F zone,Liquid+TiC +M7C3;G zone,Liquid+TiC+γ;H zone,Liquid+TiC+M7C3+γ;I zone,TiC+γ+M7C3.

        3.2 冷卻曲線和相分數(shù)演變

        圖5 給出了本模擬得到的Fe-4%C-17%Cr 合金冷卻曲線(黑色).可以看出,在溫度高于1266 ℃階段,僅有M7C3碳化物析出,相分數(shù)較少,潛熱釋放未對冷卻曲線產生明顯影響,冷卻曲線基本以恒定冷卻速率下降;至1266 ℃奧氏體析出,溫度場受到M7C3碳化物與奧氏體凝固釋放潛熱的共同作用,冷卻曲線斜率產生改變.直到82.51 s 時(見圖7),大部分M7C3碳化物與奧氏體晶粒相互接觸,生長相互促進,使得凝固過程潛熱釋放量高于散熱量,模擬域冷卻曲線溫度回復明顯.82.51—120.19 s,冷卻曲線產生明顯溫度回復,溫度從1200.39 ℃升高到1210.03 ℃,此時間段M7C3與奧氏體晶粒兩相體積分數(shù)之和從24.72%升高到68.62%.隨后,M7C3碳化物與奧氏體晶粒生長速度逐漸減小,隨之潛熱作用減小,冷卻曲線恢復下降趨勢,直至最終凝固.

        圖5 預測的Fe-4%C-17%Cr 合金冷卻曲線及添加1.5%Ti合金石墨型實驗冷卻曲線[22]對比Fig.5.Comparison of predicted cooling curve for Fe-4%C-17%Cr alloy and experimental cooling curve with 1.5%Ti addition in graphite mold [22].

        圖6 Fe-3.23%C-23.8%Cr 合金添加4%Ti 與否砂型實驗冷卻曲線[39]Fig.6.Experimental cooling curves in sand mold for Fe-3.23%C-23.8%Cr alloy with or without 4%Ti addition[39].

        圖7 本模擬冷卻曲線及各相質量分數(shù)演變Fig.7.Predicted cooling curves and evolution of phase mass fraction.

        本文模擬的Fe-4%C-17%Cr 合金的實驗冷卻曲線[40]僅監(jiān)測到1200 ℃溫度以下的相變行為,未能反映奧氏體和M7C3相的析出.為此,采用文獻[39]中同為過共晶成分的相近體系(即Fe-3.23%C-23.8%Cr 合金添加4%Ti 與否)的凝固路徑及文獻砂型實驗冷卻曲線(如圖6 所示)進行對照分析,間接驗證本模擬冷卻曲線的正確性.

        將圖5 中模擬得到的Fe-4%C-17%Cr 合金冷卻曲線(黑色)與添加1.5%Ti 后實驗冷卻曲線(紅色+實心方框)和圖6 添加Ti 與否的實驗[39]冷卻曲線進行對比.本模擬合金添加Ti 前后冷卻曲線變化趨勢與文獻合金添加Ti 前后實驗冷卻曲線變化趨勢相同,從而驗證了本模擬得到的冷卻曲線的正確性.另外,圖5、圖6 還分別給出GS 凝固路徑計算的兩成分合金添加Ti 與否的析出相順序,可以看出,Fe-4%C-17%Cr 合金添加1.5%Ti后,合金的析出相順序不變;Fe-3.23%C-23.8%Cr合金添加4%Ti 后M7C3與奧氏體析出順序發(fā)生了交換,但仍以M7C3和奧氏體為主要相.圖5 和圖6中黑色曲線的趨勢相近,拐點和溫度回復都發(fā)生在GS 凝固路徑預測的相析出溫度以下較低溫度處,溫度滯后明顯,砂型冷卻凝固時間比本模擬石墨型冷卻時間更長.

        圖7 給出了本模擬過程中奧氏體、M7C3碳化物和液相體積分數(shù)的演變過程.開始時,熔體處于完全液相,初始溫度高于M7C3碳化物形核溫度.在t=7.03 s 時,熔體冷卻至M7C3碳化物形核過冷度,M7C3碳化物形核,生長初始階段M7C3碳化物相分數(shù)增長緩慢,冷卻曲線斜率未發(fā)生改變.在t=15.22 s 時,熔體冷卻至奧氏體形核過冷度,溫降曲線出現(xiàn)第一個拐點,表明奧氏體晶粒形核并開始生長.在15.22—82.51 s 時大部分M7C3碳化物與奧氏體晶粒獨自生長、未接觸.當時間t> 82.51 s時,M7C3碳化物與奧氏體晶粒相互接觸共同生長,各相分數(shù)增長較快.冷卻曲線發(fā)生溫度回復階段,與模擬過程中M7C3碳化物、奧氏體相體積分數(shù)快速增長階段相對應;冷卻曲線溫度回復結束之后,M7C3碳化物、奧氏體相體積分數(shù)增長緩慢,接近最終凝固狀態(tài),冷卻曲線又恢復下降趨勢.通過比較M7C3碳化物晶粒獨自生長階段(t< 15.22 s)與M7C3碳化物和奧氏體晶粒共同生長階段(t≥15.22 s),可以看出合金凝固過程中M7C3碳化物與奧氏體晶粒的生長主要集中在二者相互接觸、共生長階段.當t=200.00 s 時,模擬區(qū)域中合金幾乎完全凝固,M7C3碳化物、奧氏體體積分數(shù)之和為97.744%,剩余液相分數(shù)為2.256%.此后,模擬區(qū)域內M7C3碳化物、奧氏體相分數(shù)及液相分數(shù)隨著時間的推移幾乎沒有改變,本文以此刻狀態(tài)為最終凝固狀態(tài).此時已達到Cementite 相形核溫度,本模擬僅考慮M7C3碳化物和奧氏體晶粒的生長,因此,模擬結束狀態(tài)時奧氏體枝晶臂間仍存在少量液相.

        圖8 給出了本模擬與GS 模型預測的凝固路徑中各相分數(shù)隨溫度的變化曲線,其中1304 ℃對應圖7 中7.03 s;1266 ℃對應15.22 s;1194 ℃對應147.68 s;1200.39 ℃為模擬冷卻曲線溫度回復起點;1151 ℃為模擬結束溫度.為與實驗晶粒尺寸進行比對分析,本模擬中晶粒生長速度乘了動力學系數(shù)(10—4數(shù)量級).因此,圖8 中1266 ℃時模擬的M7C3碳化物晶粒所占的模擬區(qū)域體積分數(shù)為0.0025%,比GS 模型預測值小;1200.39 ℃時,M7C3碳化物晶粒所占的模擬區(qū)域體積分數(shù)為3.15%,奧氏體晶粒所占的模擬區(qū)域體積分數(shù)為21.57%.此后,M7C3碳化物和奧氏體晶粒相互接觸共生長,快速生長階段與GS 凝固路徑快速生長階段相同,各相分數(shù)的整體變化趨勢相同,說明模擬過程各相生長的合理性.在1200.39—1210.03 ℃區(qū)間,各相分數(shù)曲線非單調增加,這是由于相分數(shù)大幅增加致使?jié)摕後尫盼醇皶r散出導致溫度回復所引起的,該階段M7C3碳化物、奧氏體相分數(shù)增長較快,潛熱釋放作用持續(xù)增強,冷卻曲線呈現(xiàn)上升趨勢.相比較而言,圖7 中各相分數(shù)隨時間的變化曲線呈單調變化,冷卻曲線溫度回復與圖8 中潛熱釋放各相體積分數(shù)在溫度回復階段繼續(xù)增長相對應.

        圖8 本模擬與GS 模型中各相質量分數(shù)隨溫度變化Fig.8.Evolution of phase mass fraction with temperature by present model and GS model.

        3.3 界面平衡濃度、晶粒生長速度隨溶質貝克列數(shù)演變

        隨著凝固進行,界面前沿的成分過冷不斷增大,溶質超飽和度將不斷增大.圖9 給出了M7C3碳化物/液相和奧氏體/液相界面液相中溶質C,Cr 的貝克列數(shù)以及溶質界面平衡濃度隨溶質Cr超飽和度的變化,均取模擬區(qū)域中所有界面胞數(shù)據(jù)的平均值作圖.可以看出,隨Cr 超飽和度增大,M7C3碳化物/液相和奧氏體/液相界面液相中溶質貝克列數(shù)均增大;在相同超飽和度下,Cr 溶質貝克列數(shù)大于C 溶質貝克列數(shù),這與界面液相中溶質貝克列數(shù)與溶質的擴散系數(shù)成反比((12)式)、而液相中溶質Cr 擴散系數(shù)小于C 擴散系數(shù)(見表1)的規(guī)律相一致.

        M7C3碳化物/液相界面C,Cr 溶質液相平衡濃度隨Cr 超飽和度呈下降趨勢(圖9(a)),與之相反,奧氏體/液相界面C,Cr 溶質液相平衡濃度隨Cr 超飽和度呈上升趨勢(圖9(b)).這是由于M7C3碳化物/液相界面溶質C,Cr 的分配系數(shù)均大于1,而奧氏體/液相界面C,Cr 溶質的分配系數(shù)均小于1,由(10)式,碳化物/液相界面溶質液相平衡濃度與超飽和度成反比,奧氏體界/液相界面溶質液相平衡濃度與超飽和度成正比,導致隨著Cr 超飽和度增大,奧氏體/液相界面溶質液相平衡濃度不斷增大,M7C3碳化物/液相界面溶質液相平衡濃度不斷減小.這與奧氏體生長向外排出C,Cr 溶質,M7C3碳化物生長吸收C,Cr 溶質的特性相符.

        圖9 M7C3 碳化物和奧氏體界面液相溶質貝克列數(shù)、平衡濃度演變 (a) M7C3 界面;(b)奧氏體界面Fig.9.Evolution of solute Peclet number and equilibrium concentration in liquid at M7C3 and austenite interface cell:(a) M7C3/liquid interface;(b) austenite/liquid interface.

        圖10 進一步給出M7C3碳化物/液相和奧氏體/液相界面生長速度隨Cr 溶質貝克列數(shù)的變化曲線,界面生長速度分別為M7C3/液相或奧氏體/液相界面各元胞生長速度的平均值.在M7C3/液相界面(圖10(a)),PeCr=0.0116 時,奧氏體析出,此后M7C3碳化物晶粒的生長速度明顯增快,說明二相晶粒相互促進生長.在PeCr< 0.0413 (交點1)區(qū)間,晶粒生長速度整體隨Cr 溶質貝克列數(shù)的增加逐漸增大,此處(交點1)對應冷卻曲線上溫度回復起始位置;PeCr=0.0366 (交點2)為冷卻曲線溫度回復達到最大值位置.同樣,在奧氏體/液相界面處(圖10(b)),在PeCr< 0.2688 (交點3)區(qū)間內,晶粒生長速度隨Cr 溶質貝克列數(shù)的增加逐漸增大,并在PeCr=0.2688 (交點3)時,生長速度達最大,此處對應冷卻曲線上溫度回復起始位置;PeCr下降至0.2114 (交點4)為冷卻曲線回復溫度達到最大值位置.在冷卻曲線溫度開始回復后,兩相晶粒生長速度逐漸降低,這是由于該時刻已有部分奧氏體晶粒與M7C3碳化物晶粒相接觸,其后奧氏體晶粒不斷包圍M7C3碳化物晶粒,接觸部位越來越多,晶粒沒有繼續(xù)生長空間,導致M7C3碳化物、奧氏體晶粒的平均生長速度降低.若冷卻速率增大,會加快模擬區(qū)域溫度的下降,由(5)式可知,會增加成分過冷和曲率過冷;通過圖9 和圖10 可知,隨著成分過冷度的不斷增大,貝克列數(shù)逐漸增大,進而會增大M7C3碳化物和奧氏體獨立生長階段界面的生長速度,縮短彼此獨立生長時間,相對來說奧氏體晶粒生長速度較快,會快速包圍M7C3碳化物晶粒,從而減小M7C3碳化物尺寸.

        圖10 M7C3 碳化物/液相和奧氏體/液相界面生長速度演變:(a) M7C3 界面;(b)奧氏體界面Fig.10.Evolution of growth velocity with PeCr at M7C3/liquid and austenite/liquid interface:(a) M7C3/liquid interface;(b) austenite/liquid interface.

        3.4 奧氏體、M7C3 碳化物晶粒共生長

        圖11 為凝固過程中不同時刻奧氏體晶粒的形貌演變,以奧氏體相質量分數(shù)顯示;圖12 為相應時刻晶粒取向圖,顯示出奧氏體晶粒(綠-藍顏色晶粒)與M7C3碳化物晶粒(粉-紅顏色晶粒)共同生長過程.t=70.83 s 時刻(圖11(a)和圖12(a)),奧氏體晶粒A,B,C 與M7C3碳化物晶粒彼此獨立生長;t=95.83 s 時刻(圖11(b)和圖12(b)),奧氏體晶粒A 的右側與M7C3碳化物晶粒左側接觸,奧氏體晶粒A 外輪廓右側形貌發(fā)生改變;t=120.83 s時刻(圖11(c)和圖12(c)),奧氏體晶粒B,C 接觸到M7C3碳化物晶粒,此時奧氏體晶粒A,B,C 圍繞M7C3碳化物晶粒生長,呈現(xiàn)出包圍趨勢;t=145.83 s 時刻(圖11(d)和圖12(d)),M7C3碳化物晶粒被奧氏體晶粒完全包圍.

        圖11 不同時刻奧氏體晶粒形貌 (a) t=70.83 s;(b) t=95.83 s;(c) t=120.83 s;(d) t=145.83 sFig.11.Morphologies of austenite grains at different moment:(a) t=70.83 s;(b) t=95.83 s;(c) t=120.83 s;(d) t=145.83 s.

        圖12 不同時刻晶粒取向 (a) t=70.83 s;(b) t=95.83 s;(c) t=120.83 s;(d) t=145.83 sFig.12.Crystallographic orientations of M7C3 and austenite grains at different moment:(a) t=70.83 s;(b) t=95.83 s;(c) t=120.83 s;(d) t=145.83 s.

        圖13 為圖11、圖12 中y=180 μm 處紅色線位置奧氏體晶粒A 和M7C3碳化物晶粒中C,Cr溶質濃度的變化曲線.當t=70.83 s 時,奧氏體晶粒A 與M7C3碳化物晶粒未接觸,此時奧氏體晶粒和M7C3碳化物晶粒中間尚有液相區(qū)域(圖13(a)對應1—2 段;圖13(b)對應 1′— 2′段)與圖12(a)黑色框內奧氏體晶粒和M7C3碳化物晶粒之間的液相區(qū)域相對應;當t=95.83 s 時,奧氏體晶粒和M7C3碳化物晶粒接觸;當t=120.83,145.83 s時M7C3碳化物晶粒左側邊界幾乎未發(fā)生改變,僅M7C3碳化物晶粒的右側向外生長.隨著時間的推移,M7C3碳化物晶粒逐漸長大,并最終被奧氏體晶粒完全包圍.由圖13 還可知,M7C3碳化物晶粒內C,Cr 溶質的濃度高于界面液相濃度,奧氏體晶粒內的C,Cr 溶質溶度低于界面液相溶質濃度,這與M7C3/液相及奧氏體/液相界面溶質分配規(guī)律相符.

        圖13 奧氏體晶粒A 和鄰近M7C3 碳化物晶粒周圍溶質濃度的演變:(a) C 濃度;(b) Cr 濃度Fig.13.Evolution of solute concentration around austenite grain A and adjacent M7C3 carbide grain:(a) C concentration;(b) Cr concentration.

        3.5 最終凝固組織形貌及體積分數(shù)與實驗對比

        圖14 給出Fe-4%C-17%Cr 合金在石墨型冷卻條件下顯微組織[22](圖14(a))與相同條件下本模擬的最終凝固結果(圖14(b)—(d)).其中,圖14(b)為以奧氏體質量分數(shù)顯示的晶粒形貌,圖中藍色區(qū)域為M7C3碳化物晶粒,白色區(qū)域為奧氏體晶粒,紅色區(qū)域為奧氏體晶臂處剩余液相區(qū);圖14(c)和圖14(d)為模擬的最終凝固時的C,Cr 溶質濃度分布.通過Image-Pro Plus 6.0 軟件測得實驗顯微組織圖中M7C3碳化物體積分數(shù)為12.00%,本模擬的M7C3碳化物體積分數(shù)為11.94%,模擬結果與實驗結果中M7C3碳化物晶粒大小形貌均相近.

        圖14 Fe-4%C-17%Cr 合金實驗[22]和預測的凝固組織 (a)實驗凝固形貌;(b)本模擬奧氏體質量分數(shù);(c)本模擬C 濃度場;(d)本模擬Cr 濃度場Fig.14.Experimental[22] and predicted solidification microstructure for Fe-4%C-17%Cr alloy:(a) experimental microstructure;(b) predicted austenite mass fraction;(c) predicted C concentration field;(d) predicted Cr concentration field.

        3.6 剩余液相與奧氏體中溶質濃度

        3.6.1 剩余液相中C,Cr 溶質濃度演變

        圖15 為本模擬(黑色曲線)與GS (各溶質在固相無擴散,紅色曲線)、PE (C 溶質在固相充分擴散、Cr 無擴散,綠色曲線)、LR (各溶質在固相充分擴散,藍色曲線)模型預測的剩余液相中溶質C,Cr 濃度的變化趨勢對比.本模擬的剩余液相溶質濃度由所有元胞中剩余液相溶質總和除以剩余總液相分數(shù)得到.凝固初期,M7C3碳化物先行析出,其生長過程中會吸收液相中的C,Cr 溶質,所以C,Cr 溶質都為下降趨勢;γ相析出后,γ生長過程中向液相中排出C,Cr 溶質(kp,γ/L,C

        圖15 剩余液相中C,Cr 溶質濃度演變Fig.15.Evolution of C and Cr concentration in residual liquid.

        3.6.2 奧氏體中C,Cr 溶質濃度演變

        圖16 給出了本模擬的奧氏體中C,Cr 溶質濃度隨奧氏體體積分數(shù)的變化曲線及與GS,PE,LR 模型曲線的對比.本模擬中C,Cr 溶質在奧氏體中有限擴散,為對比GS 結果,亦給出C,Cr 溶質在奧氏體中Dγ=0 結果.圖16(a)中初期C 溶質有擴散曲線在無擴散曲線上方,在奧氏體相分數(shù)約為0.42 時交叉,之后變?yōu)闊o擴散曲線在有擴散曲線上方.這種趨勢是由于,交叉點處為奧氏體和M7C3碳化物快速生長階段,奧氏體中溶質C 的擴散系數(shù)比Cr 大4 個數(shù)量級,奧氏體晶粒和M7C3碳化物晶粒接觸、包圍M7C3晶粒后,緩慢的凝固過程為奧氏體晶粒中C 溶質擴散提供了更長的時間,進而導致C 濃度低于無擴散結果.圖16(b)中Cr 溶質無擴散曲線在有擴散曲線上方,二者濃度值幾乎重合,這是因為Cr 溶質擴散系數(shù)較小(3.67×10—14m2/s),接近0 的緣故.總體來看,奧氏體中C溶質濃度隨凝固的進行呈逐漸增大趨勢,Cr 溶質濃度呈持續(xù)下降趨勢,C,Cr 溶質濃度值雖與GS,PE,LR 凝固路徑曲線存在一定差異,但C,Cr 整體變化趨勢與GS,PE,LR 模型曲線演變趨勢相同,C,Cr 溶質有無擴散曲線的上下位置與GS 模型、LR 模型曲線位置對應吻合.

        圖16 奧氏體中C,Cr 溶質濃度演變與凝固路徑模型對比 (a) C 濃度;(b) Cr 濃度Fig.16.Comparison of C and Cr solute concentration evolution in austenite with GS,PE and LR solidification path prediction:(a) C concentration;(b) Cr concentration.

        4 結論

        1)開發(fā)了Fe-C-Cr 三元合金小面晶M7C3碳化物與枝晶共生長的二維微觀元胞自動機模型.模型中加入潛熱釋放對凝固過程溫度場的影響,由C,Cr 兩溶質界面擴散共同確定晶體生長速度,由凝固路徑數(shù)據(jù)表插值獲取液相元胞的溶質平衡濃度,通過設定M7C3碳化物鄰胞結構并優(yōu)化形狀因子來保持M7C3小面晶形貌,模擬了Fe-4%C-17%Cr三元合金初生M7C3碳化物和共晶奧氏體晶粒的生長演變過程.

        2) M7C3碳化物和奧氏體晶粒在液相中各自的生長速度隨界面液相中C,Cr 溶質的超飽和度和貝克列數(shù)的增大而不斷增大.隨著奧氏體的析出和晶粒生長,M7C3碳化物晶粒的生長速度明顯增快,二者生長相互促進.當奧氏體晶粒逐漸接觸并包圍M7C3碳化物晶粒,對應于冷卻曲線溫度回復開始后,晶粒沒有繼續(xù)生長空間,兩相晶粒生長速度遂逐漸降低.

        3) Fe-4%C-17%Cr 合金凝固過程中,奧氏體晶粒生長向外排出C,Cr 溶質,與吸收C,Cr 溶質生長的M7C3碳化物晶粒互補,致使二者生長互相促進,最終奧氏體晶粒包圍M7C3碳化物晶粒生長,奧氏體枝晶形貌趨近M7C3碳化物晶粒外輪廓.

        4) 預測的Fe-4%C-17%Cr 合金冷卻曲線及添加1.5%Ti 后石墨鑄型實驗冷卻曲線[22]和Fe-3.23%C-23.8%Cr 合金添加4%Ti 與否的砂型實驗冷卻曲線[39]的變化趨勢相符,最終凝固組織形貌和M7C3碳化物體積分數(shù)與實驗[22]相符,剩余液相、奧氏體中C,Cr 溶質濃度演變與GS,PE,LR模型均相符.

        5)實際工程應用中,通過增大M7C3形核過冷度(如進行熔體潔凈化處理),增加形核密度(如采取合金化手段,添加Ti,V 等合金元素增加M7C3形核質點)或增大冷卻速率,可實現(xiàn)減小M7C3碳化物尺寸和均勻分布.

        猜你喜歡
        元胞溶質碳化物
        有關溶質質量分數(shù)的計算
        改善高碳鉻軸承鋼碳化物均勻性研究
        上海金屬(2022年6期)2022-11-25 12:24:20
        滴水成“冰”
        溶質質量分數(shù)考點突破
        Cr12Mo1V1鍛制扁鋼的共晶碳化物研究
        模具制造(2019年3期)2019-06-06 02:11:04
        基于元胞自動機下的交通事故路段仿真
        智富時代(2018年5期)2018-07-18 17:52:04
        “溶質的質量分數(shù)”計算歸類解析
        Nb微合金鋼中碳化物高溫溶解行為研究
        上海金屬(2016年4期)2016-11-23 05:38:50
        基于元胞數(shù)據(jù)的多維數(shù)據(jù)傳遞機制
        北京測繪(2016年2期)2016-01-24 02:28:28
        基于AIS的航道移動瓶頸元胞自動機模型
        中國航海(2014年1期)2014-05-09 07:54:25
        亚洲中文字幕第一页免费 | 射精区-区区三区| 国产特级毛片aaaaaaa高清| 国产91网址| 人妻熟女妇av北条麻记三级| 日韩精品一区二区三区在线视频| 最爽无遮挡行房视频| 女同亚洲女同精品| 国产天堂av手机在线| 91精品国产福利在线观看麻豆| 美女av一区二区三区| 日本一本久道| 色哟哟精品中文字幕乱码| 国产亚洲精品一区二区无| 全免费a级毛片免费看网站| 日本一区二区啪啪视频| 在线观看中文字幕不卡二区| 久久婷婷五月综合97色直播| 性大片免费视频观看| 丰满熟妇人妻无码区| 国产一级黄色片在线播放| 亚洲国产精品无码中文字 | 中文字字幕人妻中文| 18精品久久久无码午夜福利 | 亚洲色成人网站www永久四虎| 无码一区二区三区AV免费换脸| 亚洲区一区二区三区四| 东北女人啪啪对白| 女厕厕露p撒尿八个少妇| 91青草久久久久久清纯 | 日韩人妻不卡一区二区三区| 国产精品午夜无码av天美传媒| 午夜精品久视频在线观看| 国产精品毛片av毛片一区二区| 国产熟妇另类久久久久| 久久AV中文一区二区三区| 精品蜜臀国产av一区二区| 中文人妻av久久人妻水蜜桃| 亚洲深深色噜噜狠狠爱网站| 日韩精品一区二区av在线| 亚洲最大水蜜桃在线观看|