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

Al-Ti-B 細化工業純鋁凝固組織演變過程數值模擬*

2021-05-06 01:03:46宋巖江鴻翔趙九洲何杰張麗麗李世欣
物理學報 2021年8期
關鍵詞:生長

宋巖 江鴻翔? 趙九洲? 何杰 張麗麗 李世欣

1) (中國科學院金屬研究所, 師昌緒先進材料創新中心, 沈陽 110016)

2) (中國科學技術大學材料科學與工程學院, 沈陽 110016)

1 引 言

晶粒細化對提高鑄件強度、韌性及后續加工性能, 減輕合金成分偏析及熱裂、疏松傾向至關重要[1-3].目前, 常通過向合金熔體中添加Al-Ti-B[4,5],Al-Ti-C[6]等含有異質形核質點的中間合金的方法來細化鋁合金鑄件的晶粒組織.經過細化處理的鋁合金凝固組織演變包括α-Al 形核、晶粒枝晶化以及α-Al 晶粒長大等過程[7], 是熔體內異質核心粒子數量密度與尺寸分布[8]、合金成分[9]、凝固過程中溶質再分配[10]、熔體冷卻條件及熔體內溫度場、濃度場等因素共同作用的結果.僅通過實驗手段難以明確解析凝固組織演變機理及各影響因素的作用.近年來, 人們利用數值模擬的方式對合金凝固過程開展了大量的研究.目前, 鋁合金凝固過程的數值模擬基本上可以分為兩類, 即宏觀尺度上的模擬和微觀尺度上的模擬.宏觀尺度上的模擬主要是指對溫度場、流場及應力場等的模擬; 微觀尺度上的模擬主要是指在晶粒尺度上對凝固過程進行仿真, 通常涉及晶粒的長大、柱狀晶與等軸晶的轉變等.

材料顯微組織對材料性能有重要影響.開展微觀組織演變模擬對于深入理解晶粒形核和長大/粗化過程, 有效控制合金凝固組織具有重要意義.在合金微觀組織模擬方面, 目前采用的方法主要有:相場法(PF)[11-13]和元胞自動機法(CA)[14,15].其中,CA 方法能夠充分地考慮影響枝晶生長的各種物理機理, 與具有明確物理意義的參數相結合, 算法簡單, 計算效率較高.1994 年, Gandin 和Rappaz[16]首先建立了研究凝固組織形成過程的二維CA 模型,計算了晶粒在模具壁或熔體內部的異質形核, 再現了柱狀晶向等軸晶的轉變(CET)以及柱狀晶的競爭生長, 但該模型只適用于模擬純金屬等溫凝固過程.此后, CA 模擬方法得到了迅速的發展.2001 年, Zhu 和Hong[17]在考慮固、液相之間溶質再分配的基礎上提出了改進的CA 模型, 模擬了Al-Cu 合金凝固過程中柱狀晶和等軸晶的形貌演變.石玉峰等[18]建立了適用于三元鋁合金枝晶生長的CA 模型, 模擬了Al-Si-Mg 合金的枝晶組織演變.隨后, 潘詩琰和朱鳴芳[19]又將CA 模型從二維擴展到了三維.張顯飛等[20]建立了模擬合金樹枝晶生長的三維元胞自動機模型, 考慮了合金元素之間相互作用對溶質擴散的影響, 給出了耦合枝晶生長動力學和合金凝固熱力學計算模擬凝固組織形成過程的方法, 并用所建模型模擬了二元Al-Cu合金和三元Al-Cu-Mg 合金的枝晶生長過程.迄今為止, CA 模型的發展經歷了從二維到三維、從純金屬到合金的過程, 并逐漸引入溫度場、濃度場等對晶粒微觀組織形貌演變的影響, 能夠很好地呈現界面處溶質再分配、成分過冷等對生長的影響.然而, 以往的CA 模型通常依據實驗結果來推測合金凝固時的形核過程, 不能真正地實現晶粒形核與枝晶生長過程的耦合模擬.

本文在耦合合金熱力學計算的基礎上, 基于群體動力學方法描述α-Al 形核過程及球形生長過程, 基于CA 方法描述α-Al 晶粒枝晶生長.實現了晶粒形核與枝晶生長過程的耦合計算.通過與實驗結果比較, 驗證了該模型的準確性.利用所建模型研究了不同凝固條件下Al-5Ti-1B 中間合金細化處理工業純鋁過程中的凝固組織演變行為.

2 模型描述

經細化處理的鋁合金熔體凝固過程可以分為兩個階段: 1) 當熔體的溫度降低到臨界形核溫度時,α-Al 開始在有效異質核心粒子(粒子半徑大于α-Al 臨界形核半徑r?)上形核, 并進行球形生長; 2) 當α-Al 球形晶粒的半徑超過 2 1r?時[21], 球形晶粒轉變為樹枝狀, 直到凝固結束.

2.1 α-Al 形核及晶粒球形長大的群體動力學模型

定義函數f(r,t)來描述α-Al晶粒的尺寸分布.則t時刻, 熔體內半徑處于r—r+ dr范圍內α-Al晶粒的數量密度為f(r,t)dr.α-Al 晶粒的數量密度n, 平均半徑〈r〉, 體積分數φ可以通過下式計算:

在考慮α-Al 形核, 擴散長大/溶解的共同作用下,f(r,t)滿足如下控制方程[22]:

式中,I是α-Al 的形核率,v是α-Al 球狀晶粒擴散長大/溶解的速度, 通過下式計算[22]:

式 中,D是溶質在Al 熔體中的擴散系數.Cl是遠離固液界面處熔體中溶質的摩爾濃度,CI=Ce(α/r)為固液界面處溶質在熔體中的摩爾濃度(Ce為熔體中溶質的平衡濃度,α為毛細長度),Cs是α-Al 晶粒內部溶質的摩爾濃度.

當Al 熔體中有充足數量的有效異質形核粒子時,α-Al 的形核率通過經典形核理論計算[23]:

式中,N0是單位體積熔體內Al 原子的數量,nc為具有臨界晶核半徑r*的α-Al 晶粒內的原子數量,為形核能壘,σS/L為固液界面能, ΔGV是形核體積驅動力,δ為Al 原子的平均跳躍距離,kB是玻爾茲曼常數,T為溫度.f(θ) = (2 — 3 cosθ+ cos3θ)/4,θ是α-Al 與異質形核粒子間的接觸角.

當熔體中有效異質形核粒子數量不足時, 在t—t+ Δt時間段內,α-Al 晶粒形核過程受該時間間隔內熔體中有效異質形核粒子數量密度ΔN(t)的限制, 該時間段內,α-Al 的形核率表示為

2.2 α-Al 枝晶生長的CA 模型

當晶粒不能維持球形生長形態時, 會轉變為樹枝狀的形貌繼續長大.在枝晶生長過程中, 液相中溶質擴散的濃度場控制方程為[17]

在Δt時間內, 界面元胞固相分數的增量Δfs通過下式計算[20]:

式中,Lφ為沿界面法向方向經過元胞中心的線段長度.

基于固液界面處溶質守恒原則, 界面的法向生長速度vn為[24]

式中,為界面法向量.

枝晶生長過程中, 枝晶表面的過冷度ΔT可以表示為[18]

式中,T*為枝晶表面實際溫度,Tmelt為熔體實際溫度; ΔTc=Teq—TL是成分過冷,Teq是界面處局部平衡液相線溫度,TL是合金初始的液相線溫度;ΔTK是動力學過冷, 通常較小, 可以忽略; 曲率過冷ΔTr可表示為[25]

式中,Γˉ 為Gibbs-Thomson 系數,K為界面平均曲率, 由下式給出[24]:

式中, Δs和Z分別是元胞尺寸和鄰胞的數量.對于三維計算區域, 取26 鄰居構型(Z= 26), Δs可以根據計算要求取不同的尺寸.是界面能各向異性函數, 由下式計算[19]:

式中,ε為界面能各向異性強度值.

2.3 合金熱力學計算

基于置換溶液模型, 合金體系的吉布斯自由能由下式計算[26]:

式中,Gφ為φ相的吉布斯自由能(l 為液相, s 為固相),xi為組元i的摩爾分數,為純組元i的標準摩爾吉布斯自由能,R= 8.314 J·mol—1·K—1, 為氣體常數.為Redlich-Kister 參數.熱力學計算所用參數來自于COST 507 數據庫[27].

Al 熔體與α-Al 之間的界面能σS/L可以通過下式計算[28]:

式中, ΔHf,m為合金的摩爾熔化焓變, ΔSf,m為合金的摩爾熔化熵變.NA為Avogadro 常數,Vm為Al熔體的摩爾體積.對于面心立方鋁合金體系(fcc-Al), 無量綱界面能β= 0.561[28].

ΔSf,m可以通過熱力學第二定律導出式計算[29]:

式中,Sl(T,Cl)和Ss(T,Cs)分別為α-Al 液固兩相的熵值.Gl和Gs分別為液固兩相的摩爾吉布斯自由能.

ΔHf,m通過下式計算[29]:

式中,Hl(T,Cl)和Hs(T,Cs)分別為α-Al 形核時液固兩相的焓值.

2.4 溫度場控制方程

凝固過程中, 合金熔體的溫度變化通過下式計算:

式中,T(t)為熔體t時刻溫度,T0為澆鑄前熔體初始溫度,vcool為熔體的冷卻速度, ΔHv和Cp分別是Al 熔體的凝固潛熱和比熱.

3 結果與討論

為了驗證所建模型準確性, 模擬了Al-5Ti-1B(質量分數%, 下同)中間合金細化處理工業純鋁凝固微觀組織演變過程, 并與Greer 等[30]的實驗結果進行了對比(實驗過程中所用工業純鋁成分如表1 所列).溶質對凝固過程中晶粒形核與長大的影響通過生長限制因子Q=mi(1 —ki)C0,i表征[31].其中,mi為溶質i的液相線斜率,ki為平衡分配系數,C0,i為溶質初始成分.計算發現: 工業純鋁中溶質Ti 的Q值遠大于其他元素, 因此模擬過程中將工業純鋁簡化為二元Al-Ti 合金.考慮到TiB2粒子在鋁溶體中的穩定性較好, 計算過程中忽略其溶解/粗化過程, 假設TiB2粒子在熔體中穩定存在[30].模擬使用的熱物性參數如表2 所列.

表1 工業純鋁成分[30]Table 1.The composition of solute in the CP-Al alloys used by Greer et al.[30].

表2 模擬使用的工業純鋁參數[22,30,32]Table 2.The parameters of CP-Al alloys used in the simulations[22,30,32].

加入熔體中的TiB2粒子尺寸分布通過下式描述:

式中,NTiB2是直 徑為d的TiB2粒子數量,d0=0.72 μm 為粒子的平均直徑,N0為加入到熔體中的TiB2粒子總的數量.

3.1 細化處理條件下鋁合金凝固微觀組織演變

通過擬合計算確定α-Al 與TiB2粒子之間的接觸角θ.即取不同的接觸角θ, 模擬0.3%的Al-5Ti-1B 中間合金添加到工業純鋁熔體中, 并以3.5 K/s 的冷卻速度凝固后工業純鋁的最終凝固組織及晶粒尺寸.對比模擬結果與實驗結果, 確定θ為4.8°.

圖1 為添加0.01% Al-5Ti-1B 中間合金后, 熔體以3.5 K/s 的速度冷卻時, 工業純鋁凝固過程中熔體的過飽和度、液-固相變的形核驅動力以及α-Al形核率與凝固時間的關系.隨著熔體溫度的降低,熔體逐漸變為過飽和狀態.熔體的過飽和度和形核驅動力隨著熔體冷卻時間的延長逐漸增大.當熔體的溫度降低到某一臨界溫度時,α-Al 在Al 熔體中開始形核.形核過程非常短暫, 大約為0.07 s.由于受到有效TiB2粒子數量的影響,α-Al 的形核過程可分為兩個階段: 階段I 為形核剛開始發生時, 熔體中有效TiB2粒子數量充足,α-Al 形核率由(6)式確定; 階段II 為隨著熔體過飽和度、α-Al 形核驅動力的增大, 有效TiB2粒子數量增加速度不能滿足熔體形核時, 有效TiB2粒子數目將限制熔體中α-Al 的形核(圖1 中階段II).熔體冷卻過程中, 晶粒形核和長大釋放凝固潛熱趨向于降低熔體過冷度, 熔體與外界熱交換引起的熔體溫度降低將增加熔體過冷度.在凝固的初始階段, 熔體中α-Al晶粒的形核率及數量比較低, 因此熔體過冷度隨著凝固時間的增加逐漸增加.高的過冷度使得熔體內部晶核形核率增加, 形成較多新的α-Al 晶粒.當凝固潛熱釋放的速度超過熔體與外界換熱速度時,再輝發生, 熔體過冷度減小, 不能進一步提供新的有效異質核心, 形核過程停止.

圖1 0.01% Al-5Ti-1B 中間合金細化處理工業純鋁凝固過程中, 熔體過飽和度S = Cl — Ce, α-Al 形核驅動力ΔGV,異質形核率I 隨時間的變化.熔體冷速為3.5 K/sFig.1.The supersaturation of the melt S = Cl — Ce (dash line), the driving force of nucleation ΔGV (dot line) and the heterogeneous nucleation rate I (solid line) of α-Al during cooling the CP-Al melt inoculated by 0.01% Al-5Ti-1B master alloy at the cooling rate of 3.5 K/s.

α-Al 異質形核期間, 晶粒數量隨著凝固時間增加逐漸增加, 尺寸分布變寬(如圖2 所示).形成的α-Al 晶粒在形核早期以球狀方式快速生長, 當晶粒半徑超過球狀生長的臨界值( 2 1r?)時, 界面失穩, 晶粒開始以樹枝狀的方式進行生長.圖3 給出了兩種不同生長方式的晶粒數量密度在α-Al異質形核期間隨凝固時間的變化.可見, 球狀α-Al晶核數量在形核初期快速增加, 隨后逐漸地演變成為樹枝晶.最終, 晶核都以樹枝晶的方式生長.從圖4(a)—圖4 (c)可以看出, 生長初期, 樹枝晶距離較遠, 以多重對稱的方式沿著各自優先生長取向生長.隨著樹枝晶周圍溶質富集程度的不斷增加,固-液界面的不穩定性加劇, 一次枝晶臂粗化, 二次枝晶臂開始顯現.當不同的樹枝晶逐漸長大并互相接觸時, 枝晶周圍的溶質場發生碰撞, 枝晶臂生長互相阻礙, 樹枝晶逐漸失去對稱性.在凝固末期,不同的樹枝晶相互接觸, 枝晶臂不斷生長和粗化,最終形成具有不規則形貌的凝固組織(如圖4(d)所示).

圖2 0.01% Al-5Ti-1B 中間合金細化處理工業純鋁凝固過程中, α-Al 晶核半徑r 分布隨時間的變化.熔體冷速為3.5 K/sFig.2.The radius distribution of α-Al nucleus at different time during cooling the CP-Al melt inoculated by 0.01% Al-5Ti-1B master alloy at the cooling rate of 3.5 K/s.

圖3 0.01% Al-5Ti-1B 中間合金細化處理工業純鋁凝固過程中, α-Al 晶核總數量nall, 球狀晶nspherical 以及樹枝晶數目ndendrities 隨時間的變化.熔體冷速為3.5 K/sFig.3.The number density of all nucleus nall in the Al melt, the number density of spherical nucleus nspherical and the number density of dendrities ndendrities during cooling the CP-Al melt inoculated by 0.01% Al-5Ti-1B master alloy at the cooling rate of 3.5 K/s.

圖4 0.01% Al-5Ti-1B 中間合金細化處理工業純鋁凝固微觀組織演變過程 (a) 固相分數 = 0.1%; (b) 固相分數 = 1.0%; (c) 固相分數 = 25.0%; (d) 固相分數 = 70.0%.熔體冷速為3.5 K/s.計算區域尺寸為600 μm × 600 μm ×600 μmFig.4.Solidification microstructure evolution during cooling the CP-Al melt inoculated by 0.01% Al-5Ti-1B master alloy at a cooling rate of 3.5 K/s: (a) Solid fraction =0.1%; (b) solid fraction = 1.0%; (c) solid fraction = 25.0%;(d) solid fraction = 70.0%.The size of computational domain is 600 μm × 600 μm × 600 μm.

3.2 中間合金添加量對鋁合金凝固過程的影響

圖5 為工業純鋁熔體中添加不同數量Al-5Ti-1B 中間合金凝固時,α-Al 異質形核率與熔體過冷度隨時間的變化關系.可以看出, 當Al-5Ti-1B 中間合金加入量為0.01%時,α-Al 在熔體過冷度約為0.4 K 時開始形核.形核過程中, 階段I 的形核主要受熔體過冷度的控制, 形核率隨熔體過冷度的增大不斷增加.異質形核過程進入階段II 時, 熔體中有效TiB2粒子數量不足,α-Al 形核受到熔體中有效TiB2粒子數量限制.隨著凝固潛熱釋放速度的增加, 熔體過冷度開始降低, 再輝發生, 形核過程停止.形核過程中階段I 和II 持續時間分別占總形核時間的50%左右.當Al-5Ti-1B 加入量為0.4%和1.0%時, 由于Al-Ti 合金體系富鋁角液相線斜率為正值, 額外的溶質Ti 含量的增加使液相線溫度提高, 因此熔體形核時達到的最大形核過冷度較大, 熔體中有效TiB2粒子數量隨之增加.當熔體過冷度分別約為0.6 和1.0 K 時,α-Al 開始形核.形核過程中, 階段I 發生的時間占整個形核過程的90%左右.階段II 發生較短時間后, 熔體中形核及晶粒生長帶來的固相分數增加釋放凝固潛熱速率足夠快, 熔體過冷度開始降低, 再輝發生,形核過程停止.

圖5 不同數量Al-5Ti-1B 中間合金細化工業純鋁凝固時異質形核率I, 熔體過冷度ΔT 隨時間的變化.熔體冷速為3.5 K/sFig.5.Calculated heterogeneous nucleation rate I of α-Al and the undercooling of the melt ΔT for the CP-Al inoculated by different amount of Al-5Ti-1B master alloy at a cooling rate of 3.5 K/s.

圖6 為模擬計算的工業純鋁最終凝固組織.可以看到, 隨著中間合金添加量的增加, 凝固組織中含有的晶粒數目逐漸增多, 晶粒尺寸逐漸減少.當凝固組織晶粒尺寸較大時, 枝晶的枝晶臂均較發達, 晶粒大小分布較分散.這類凝固組織通常具有較差力學能與后續加工性能.隨著晶粒尺寸的降低, 枝晶臂的生長受到限制, 凝固組織主要以均勻細小的等軸晶粒組成.圖7 給出了不同數量細化劑細化工業純鋁晶粒尺寸的預測值與實驗值.可以發現, 模擬結果與實驗結果符合得較好.

圖6 不同數量Al-5Ti-1B 中間合金細化工業純鋁凝固組織模擬結果 (a) 0.005%; (b) 0.01%; (c) 0.4%; (d) 1.0%.熔體冷速為3.5 K/s.計算區域尺寸為900 μm × 900 μm ×900 μmFig.6.Simulated solidification microstructure for the CP-Al melt inoculated by different amount of Al-5Ti-1B master alloy at a cooling rate of 3.5 K/s: (a) 0.005%; (b) 0.01%;(c) 0.4%; (d) 1.0%.The size of computational domain is 900 μm × 900 μm × 900 μm.

圖7 不同數量Al-5Ti-1B 中間合金細化工業純鋁晶粒尺寸預測結果與實驗結果.熔體冷速為3.5 K/s.誤差棒為標準偏差Fig.7.Predicted and measured grains size vs.the additive amount of Al-5Ti-1B master alloy for the inoculated CP-Al at a cooling rate of 3.5 K/s.The error bars represent the standard deviations.

圖8 不同冷卻速度下經0.5% Al-5Ti-1B 中間合金細化處理, 工業純鋁凝固時α-Al 形核率隨時間的變化 (a) 1.0 K/s;(b) 2.0 K/s; (c) 5.5 K/s; (d) 10.0 K/sFig.8.Calculated heterogeneous nucleation rate of α-Al during cooling the CP-Al inoculated by 0.5% Al-5Ti-1B master alloy at the different rate: (a) 1.0 K/s; (b) 2.0 K/s; (c) 5.5 K/s; (d) 10.0 K/s.

3.3 冷卻速度對細化處理條件下鋁合金凝固過程的影響

圖8 給出了不同冷速條件下α-Al 異質形核率隨時間的變化.可以看出, 當熔體的冷卻速度為1.0, 2.0 和5.5 K/s 時, 第I 階段形核占據了α-Al異質形核的大部分, 第II 階段形核只持續很短的時間后, 再輝發生, 形核停止(如圖8(a)—圖8(c)).當熔體的冷卻速度為10 K/s 時, 形核階段I 持續時間所占比例降低, 之后階段II 的形核過程受熔體中異質核心粒子數量的限制, 直到再輝發生, 形核停止, 如圖8(d).

圖9 為不同冷速條件下, 0.5% Al-5Ti-1B 中間合金細化處理工業純鋁凝固組織的模擬結果.可以看到, 隨著熔體冷卻速度的增加, 凝固組織中含有的晶粒數目逐漸增多, 晶粒尺寸逐漸減少.當熔體冷卻速度較低時, 凝固過程中形核數目較少, 凝固組織晶粒較大.當熔體冷卻速度較高時, 凝固過程中短時間內產生的晶核數目較多, 晶粒沒有充足的時間長大, 枝晶的枝晶臂生長受到周圍晶粒的阻礙, 合金最終凝固組織主要以均勻細小的等軸晶粒組成.圖10 給出了不同冷速條件下經0.5% Al-5Ti-1B 中間合金細化, 工業純鋁晶粒尺寸的預測值與實驗值.可以發現, 模擬結果與實驗結果符合得較好.

圖9 0.5%Al-5Ti-1B 中間合金細化工業純鋁, 在不同冷卻速度凝固后凝固組織的模擬結果 (a) 1.0 K/s; (b) 2.0 K/s;(c) 5.5 K/s; (d) 10 K/s.計算區域尺寸為900 μm × 900 μm ×900 μmFig.9.Simulated solidification microstructure for the CP-Al melt inoculated by 0.5% Al-5Ti-1B master alloy at the different cooling rate of the melt: (a) 1.0 K/s; (b) 2.0 K/s;(c) 5.5 K/s; (d) 10 K/s.The size of computational domain is 900 μm × 900 μm × 900 μm.

圖10 0.5%Al-5Ti-1B 中間合金細化工業純鋁, 在不同冷卻速度下凝固后晶粒尺寸預測結果與實驗結果.誤差棒為標準偏差Fig.10.Predicted and measured grains size vs.the cooling rate of the melt for the CP-Al inoculated by 0.5%Al-5Ti-1B master alloy.The error bars represent the standard deviations.

4 結 論

采用耦合群體動力學方法與元胞自動機方法建立了描述細化處理條件下工業純鋁凝固微觀組織演變的數值模型.該模型使用群體動力學方法計算α-Al 的非均勻形核過程及α-Al 晶粒的初始球形生長過程, 使用元胞自動機方法來計算晶粒的枝晶演變過程.利用建立的模型, 通過輸入中間合金內TiB2粒子尺寸分布、添加量以及合金熔體冷卻速度等工藝參數, 可描述不同凝固條件下工業純鋁凝固微觀組織演變過程, 預測細化處理后工業純鋁凝固微觀組織形貌.

猜你喜歡
生長
野蠻生長
碗蓮生長記
小讀者(2021年2期)2021-03-29 05:03:48
生長的樹
少兒美術(2020年3期)2020-12-06 07:32:54
自由生長的家
現代裝飾(2020年11期)2020-11-27 01:47:48
美是不斷生長的
快速生長劑
共享出行不再“野蠻生長”
生長在哪里的啟示
華人時刊(2019年13期)2019-11-17 14:59:54
野蠻生長
NBA特刊(2018年21期)2018-11-24 02:48:04
生長
文苑(2018年22期)2018-11-19 02:54:14
主站蜘蛛池模板: 亚洲成人高清无码| 精品视频免费在线| 中文字幕在线观看日本| 日韩经典精品无码一区二区| 黄色网站不卡无码| 999精品色在线观看| 91九色视频网| 国内精品免费| 久久人搡人人玩人妻精品| 成人午夜亚洲影视在线观看| 欧美日韩91| 久久男人资源站| 青青草一区| 一级不卡毛片| 中文字幕天无码久久精品视频免费| 亚洲有无码中文网| 亚洲成aⅴ人在线观看| 国产精品网拍在线| 欧美日韩专区| 999在线免费视频| 日韩欧美在线观看| 午夜视频日本| 欧美第一页在线| 亚洲欧洲自拍拍偷午夜色| 日韩午夜福利在线观看| 亚洲侵犯无码网址在线观看| 精品国产91爱| 欧洲精品视频在线观看| 日韩精品一区二区三区中文无码| 亚洲天堂高清| 国产精品55夜色66夜色| 国产在线无码一区二区三区| 美女潮喷出白浆在线观看视频| 日韩一区二区三免费高清| 在线观看的黄网| 亚洲一区毛片| 欧美va亚洲va香蕉在线| 91网站国产| 无码乱人伦一区二区亚洲一| 一区二区在线视频免费观看| 国内丰满少妇猛烈精品播| 日本免费新一区视频| 亚洲欧美日韩天堂| 色AV色 综合网站| 国产网站免费观看| 日韩无码黄色| 91视频日本| 亚洲欧美日韩成人高清在线一区| 亚洲午夜片| 日韩激情成人| 亚洲无码高清视频在线观看| 亚洲 日韩 激情 无码 中出| 国产日韩欧美视频| a亚洲天堂| 欧美精品影院| 国产激情无码一区二区APP| 91综合色区亚洲熟妇p| 色综合狠狠操| 中文字幕2区| jizz在线免费播放| 永久免费av网站可以直接看的 | 免费人成黄页在线观看国产| 亚洲日韩在线满18点击进入| 又黄又爽视频好爽视频| 天天摸夜夜操| 午夜激情福利视频| 1769国产精品视频免费观看| 国产中文在线亚洲精品官网| 制服丝袜无码每日更新| 国产精品密蕾丝视频| 免费观看国产小粉嫩喷水| 亚洲精品另类| 日本成人精品视频| 亚洲无码视频图片| 欧美另类一区| 天天综合亚洲| 日本不卡在线播放| 成年人国产视频| 国产精品美女网站| 国产精品综合色区在线观看| 亚洲资源站av无码网址| 幺女国产一级毛片|