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

        ?

        Al-4%Cu凝固過程枝晶生長的數(shù)值模擬

        2016-08-16 00:57:20徐藹彥李露露
        材料工程 2016年6期
        關(guān)鍵詞:方向生長

        張 敏,徐藹彥,汪 強(qiáng),李露露

        (西安理工大學(xué) 材料科學(xué)與工程學(xué)院,西安 710048)

        ?

        Al-4%Cu凝固過程枝晶生長的數(shù)值模擬

        張敏,徐藹彥,汪強(qiáng),李露露

        (西安理工大學(xué) 材料科學(xué)與工程學(xué)院,西安 710048)

        改進(jìn)了模擬枝晶生長常用的二維元胞自動機(jī)和有限差分(CA-FD)模型,新模型引入擾動函數(shù)來控制二、三次枝晶的生長;在枝晶生長過程中,將溶質(zhì)濃度明確地分為液相溶質(zhì)濃度和固相溶質(zhì)濃度兩部分;并在溶質(zhì)再分配與擴(kuò)散過程中采用八鄰位差分以減少網(wǎng)格形狀導(dǎo)致溶質(zhì)擴(kuò)散的各向異性。模擬了Al-4%Cu二元合金過冷熔體中,單個和多個等軸晶沿不同擇優(yōu)方向生長及單方向和多方向柱狀樹枝晶競爭生長過程中的枝晶形貌、液相溶質(zhì)濃度和固相溶質(zhì)濃度分布情況。模擬結(jié)果表明:擾動的引入能夠促使枝晶產(chǎn)生分支,并控制二、三次枝晶的生長速率;液/固相溶質(zhì)計(jì)算模型能夠準(zhǔn)確地模擬出枝晶生長過程中液/固相溶質(zhì)分布;此外改進(jìn)后的模型實(shí)現(xiàn)了枝晶沿任意方向的競爭生長。

        數(shù)值模擬;元胞自動機(jī)法;枝晶形貌;溶質(zhì)濃度

        金屬材料的性能除與合金元素含量緊密相關(guān)外,還受到其微觀組織的制約,有效掌握過冷熔體凝固過程的特征是分析和控制凝固后微觀組織的最佳途徑[1]。過冷金屬在凝固過程中常常以樹枝晶的形式形核并生長,因此,長期以來對枝晶生長規(guī)律的研究受到國內(nèi)外專家學(xué)者的廣泛關(guān)注[2-4]。

        然而影響合金凝固微觀組織的因素很多,僅僅通過實(shí)驗(yàn)的方法研究結(jié)晶過程不僅耗時、耗力和耗資,而且很難直接觀察到枝晶生長的整個過程。因此,通過數(shù)值模擬技術(shù)研究凝固規(guī)律,預(yù)測顯微組織具有巨大的研究潛力,近幾十年來,通過材料領(lǐng)域科研工作者的共同努力,建立了各種各樣的微觀組織模擬方法[5-8],其中最著名的有相場法(Phase Field Method,PFM)[5]和元胞自動機(jī)法(Cellular Automaton Method,CAM)[6]。雖然元胞自動機(jī)法運(yùn)用于材料領(lǐng)域相對于相場法較晚,但其以高效和較強(qiáng)的工程應(yīng)用能力而發(fā)展迅猛。在國外,Minoru等[9]利用元胞自動機(jī)法模擬出了Fe-0.7%C二元合金鑄造過程中枝晶的生長。在國內(nèi),付振南等[10]建立了1種基于概率捕獲的元胞自動機(jī)模型,并順利模擬出鎂合金的枝晶生長過程;占小紅等[11]采用對角線模擬角度法實(shí)現(xiàn)了枝晶沿特定方向生長,但能夠模擬的枝晶角度有限,且未能實(shí)現(xiàn)枝晶全方位地?fù)駜?yōu)生長。

        本工作對現(xiàn)有的元胞自動機(jī)和有限差分(CA-FD)模型做進(jìn)一步改進(jìn),使之能夠模擬出Al-4%Cu(質(zhì)量分?jǐn)?shù))合金凝固過程中等軸晶和柱狀樹枝晶沿任意方向生長的枝晶形貌,同時能夠再現(xiàn)枝晶生長過程中液相溶質(zhì)濃度和固相中溶質(zhì)濃度的分布狀態(tài)。

        1 數(shù)學(xué)物理模型

        1.1枝晶生長模型

        過冷熔體的二維凝固區(qū)域劃分為均勻的正方形網(wǎng)格,并將溫度、溶質(zhì)濃度等多物理量賦予網(wǎng)格;另外,網(wǎng)格分為3種狀態(tài),即液相、界面和固相。界面向前推進(jìn)過程采用Moore鄰域(八鄰胞)捕捉液相網(wǎng)格[12]。

        枝晶生長過程伴隨著溫度和溶質(zhì)濃度的變化,包括枝晶尖端的穩(wěn)態(tài)生長和形成分支時的非穩(wěn)態(tài)生長兩部分。穩(wěn)態(tài)生長時,已形核的枝晶在熱過冷ΔTt、成分過冷ΔTc和曲率過冷ΔTr的作用下不斷捕捉液相而長大,則生長界面的總過冷度為:

        (1)

        根據(jù)Gibbs-Thompson方程[13],可建立tn時刻固/液界面的平衡關(guān)系:

        (2)

        (3)

        采用傳統(tǒng)的尖銳界面模型[14],界面生長速度與總過冷度之間的關(guān)系為:

        (4)

        式中:μk(θ)為界面動力學(xué)系數(shù)。

        另外,枝晶生長時不僅會受到生長動力學(xué)各向異性的影響,還會受到界面能各向異性的制約。界面生長動力學(xué)系數(shù)μ(θ)和Gibbs-Thompson系數(shù)Γ(θ)可表示為:

        (5)

        (6)

        (7)

        固相率的增長是影響過冷熔體中液/固轉(zhuǎn)變最主要的變量之一,它與界面速率成正比。用式(8)進(jìn)行計(jì)算:

        (8)

        式中:Δx為網(wǎng)格尺寸;Δt為步長時間;G為鄰位網(wǎng)格狀態(tài)參數(shù);A為擾動因子;rand()能夠在[0,1]產(chǎn)生一個隨機(jī)數(shù)。

        網(wǎng)格的形狀對枝晶生長的方向有很大的影響,為了更好地再現(xiàn)枝晶生長擇優(yōu)方向的任意性,引進(jìn)鄰位網(wǎng)格狀態(tài)參數(shù)[15],其值由式(9)確定:

        (9)

        式中:b為恒定的經(jīng)驗(yàn)參數(shù);s′和s″分別為最鄰近的4個網(wǎng)格和次鄰近的4個網(wǎng)格的狀態(tài)參數(shù)。如果相鄰的元胞為固相時,s′和s″值為1,為液相或界面元胞時,其值為0。

        界面元胞從被捕捉開始,固相率不斷增加。時間tn內(nèi),某一標(biāo)記的界面網(wǎng)格j內(nèi)固相率為:

        (10)

        1.2溶質(zhì)再分配與擴(kuò)散控制模型

        溶質(zhì)再分配過程是枝晶生長狀態(tài)與溶質(zhì)場耦合的樞紐,合理有效的溶質(zhì)分配機(jī)制是枝晶生長研究者追求的目標(biāo)。本文在保證整個凝固場溶質(zhì)守恒的前提下,針對界面元胞固液共存現(xiàn)象,提出1種修正的溶質(zhì)再分配方法。

        假定凝固過程中,液/固界面的液相和固相溶質(zhì)濃度滿足:

        (11)

        (12)

        (13)

        從界面元胞排出的溶質(zhì)導(dǎo)致枝晶周圍液相溶質(zhì)濃度升高,使液相元胞間出現(xiàn)較大的濃度梯度,這必然加劇溶質(zhì)的擴(kuò)散。對于二維非穩(wěn)態(tài)溶質(zhì)擴(kuò)散,采用如下控制方程:

        式中:DL,DS分別為液相、固相中溶質(zhì)擴(kuò)散系數(shù)。在數(shù)值計(jì)算過程中,為了消除網(wǎng)格形狀因素對溶質(zhì)擴(kuò)散過程的影響,文中采用八鄰胞的溶質(zhì)擴(kuò)散模型,即:

        (16)

        (17)

        (18)

        式中:υmax為枝晶在當(dāng)前時間步長內(nèi)的最大生長速率。

        本文涉及的Al-Cu合金熱物性參數(shù)見表1,其中,Cu為溶質(zhì),且溶質(zhì)濃度為質(zhì)量分?jǐn)?shù)。

        表1 Al-Cu合金熱物性參數(shù)

        1.3模型的基本假設(shè)

        (2)材料的熱物性參數(shù)不變,為一常數(shù);不考慮熔體流動引起的換熱及傳質(zhì),且假定溶質(zhì)的擴(kuò)散只在同種狀態(tài)下進(jìn)行,即固相與液相之間不存在溶質(zhì)擴(kuò)散;過冷熔體凝固過程中,假定液相溶質(zhì)擴(kuò)散系數(shù)和固相溶質(zhì)擴(kuò)散系數(shù)恒定不變。

        2 模擬結(jié)果與分析

        2.1單個等軸晶生長

        為了分析過冷熔體中枝晶生長時溶質(zhì)濃度的變化,將計(jì)算區(qū)域劃分為400×400個正方形網(wǎng)格,網(wǎng)格尺寸為0.5μm。假設(shè)整個計(jì)算區(qū)域溫度均勻且恒定,Al-4%Cu合金在過冷度ΔT=20K的作用下開始凝固,在計(jì)算區(qū)域中心放置一個具有擇優(yōu)生長方向的晶粒。

        圖1為0.002s時枝晶生長形貌及溶質(zhì)濃度分布,擇優(yōu)取向?yàn)?°。從圖1可以看出,一次枝晶臂沿4個方向的生長形貌均勻,一次枝晶臂上有微小的突起,這是二次枝晶臂的雛形。本模型還實(shí)現(xiàn)了液相溶質(zhì)濃度和固相溶質(zhì)濃度的分離,如圖1(b),(c);圖1(d)為0.002s時枝晶生長沿一次枝晶臂的溶質(zhì)濃度分布曲線,可以看出,在過冷度一定的情況下,枝晶在生長初期處于非穩(wěn)態(tài)生長階段,固相中溶質(zhì)濃度逐漸增高,且增加的趨勢減弱,當(dāng)生長一段時間后,固相溶質(zhì)濃度趨于一穩(wěn)定值;而液相溶質(zhì)濃度沿枝晶生長方向形成一個圓滑的過渡區(qū)域,這與金屬凝固基本原理一致。

        圖1 單個等軸晶生長0.002s模擬結(jié)果(a)枝晶生長形貌;(b)液相溶質(zhì)濃度分布;(c)固相溶質(zhì)濃度分布;(d)枝晶沿主干臂方向的溶質(zhì)濃度分布曲線Fig.1 Simulation results of single equiaxed growth after 0.002s(a)dendrite morphology;(b)distribution of liquid solute concentration;(c)distribution of solid solute concentration;(d)the curve of solute distribution along the primary arm

        枝晶在生長過程中,除形成粗大的一次枝晶臂外還會生長出參差不齊的二、三次枝晶臂。為了更好地展現(xiàn)枝晶生長過程中二、三次枝晶生長位置的隨機(jī)性,施加一個隨機(jī)的擾動函數(shù),其擾動因子在(0,1)區(qū)間內(nèi)取值。圖2為擇優(yōu)取向?yàn)?°時不同擾動因子作用下0.005s時枝晶的生長形貌,可以看出,隨著擾動因子的增加,二次枝晶臂數(shù)量不斷增多,且隨著擾動因子的增大而變粗;還可以看出,枝晶形貌的對稱性幾乎不隨著擾動因子的改變而失穩(wěn),即模型不會因?yàn)閿_動的施加而失衡。

        圖2 不同擾動因子作用下枝晶生長形貌(a)A=0.05;(b)A=0.3Fig.2 Dendrite morphology at different noise(a)A=0.05;(b)A=0.3

        圖3為擾動因子為0.3時,不同擇優(yōu)方向枝晶的生長形貌和液相溶質(zhì)濃度分布,擇優(yōu)方向依次為0°,30°,45°和60°。從圖3(a-1)~(d-1)可以看出,本文的模型能夠很好地實(shí)現(xiàn)枝晶沿不同方向的生長,在枝晶生長過程中,一次枝晶生長優(yōu)勢突出,二次枝晶大小沿一次枝晶臂方向程階梯狀變化。從圖3(a-2)~(d-2)可得,在溶質(zhì)濃度再分配與擴(kuò)散過程中,由于二次枝晶的阻礙,溶質(zhì)富集現(xiàn)象十分明顯,生長方向不同的二次枝晶臂生長到一定程度時,臂與臂間會形成一個封閉的區(qū)域,該區(qū)域中的溶質(zhì)因無法排出而濃度急劇增加,形成多個高溶質(zhì)濃度的空腔,在空腔中極易孕育形成二次相。

        圖3 不同擇優(yōu)取向時枝晶生長形貌(1)及液相中溶質(zhì)濃度分布(2)(a)0°;(b)30°;(c)45°;(d)60°Fig.3 Dendrite morphology (1) and distribution of liquid solute concentration (2) at different preferred orientation(a)0°;(b)30°;(c)45°;(d)60°

        2.2多等軸晶生長

        為了研究過冷熔體中Al-4%Cu合金多枝晶的生長行為。采用500×500個均勻的正方形網(wǎng)格劃分整個過冷熔體,網(wǎng)格的尺寸為0.5μm,初始過冷度ΔT=15K,溶質(zhì)濃度為4%。凝固開始時,在過冷熔體中設(shè)定14個大小相等,擇優(yōu)方向不同的晶粒。在凝固過程中,假定整個熔體溫度均勻,冷卻速率為20K/s。

        圖4 多個等軸晶生長形貌及溶質(zhì)濃度分布(a),(b),(c)t=0.005s;(a′),(b′),(c′)t=0.01sFig.4 Equiaxed dendrite morphology and distribution of solute concentration at different time(a),(b),(c)t=0.005s;(a′),(b′),(c′)t=0.01s

        圖4為Al-4%Cu合金凝固過程中不同生長階段內(nèi)枝晶生長形貌和溶質(zhì)濃度分布。由圖4可看出,枝晶有0°,30°,45°和60°4種不同的擇優(yōu)方向。晶粒1,5,8和14的擇優(yōu)方向?yàn)?°,晶粒2,7,10和13的擇優(yōu)方向30°,晶粒3和9擇優(yōu)方向?yàn)?5°,晶粒4,6,11和12的擇優(yōu)方向?yàn)?0°。在枝晶生長初期,晶粒沿著各自不同的擇優(yōu)方向生長,相互影響較小,如圖4(a)所示;而在生長后期后,枝晶相互阻礙形成不對稱的等軸晶,如圖4(a′)中的晶粒1,這是因?yàn)樯L初期,枝晶生長排出的溶質(zhì)有限,擴(kuò)散所影響的區(qū)域較??;隨著凝固的進(jìn)行,排除溶質(zhì)不斷累積,擴(kuò)散加強(qiáng),多個溶質(zhì)場相互疊加,致使枝晶間相互影響程度加劇。從圖4還可以看出,凝固前期,二次枝晶臂相對于一次枝晶臂很小,而在凝固后期,由于枝晶間的相互阻礙,一次枝晶臂生長受限,幾乎停止生長,此時,處于有利生長方向的二、三次枝迅速生長,部分二次枝晶臂大小甚至超過了一次枝晶臂,如圖4(a′)中的晶粒2。從圖4(c)可以看出,從枝晶中心沿一次枝晶臂向外,固相溶質(zhì)濃度依次升高,二次枝晶臂內(nèi)溶質(zhì)濃度高于該臂生長位置處一次枝晶臂內(nèi)的溶質(zhì)濃度,枝晶邊緣處固相溶質(zhì)濃度最高,且各個晶粒外圍固相溶質(zhì)濃度差異較小。

        2.3定向凝固單方向及多方向柱狀樹枝晶競爭生長

        在凝固過程中,當(dāng)溫度梯度G滿足特定條件時,枝晶會沿一定的方向生長形成柱狀樹枝晶。假設(shè)定向凝固溫度梯度G=1000K/m,且溫度梯度垂直于熔體底部。將枝晶生長區(qū)域劃分為400×400個正方形網(wǎng)格,網(wǎng)格尺寸為1μm。模擬單方向柱狀樹枝晶生長時,在過冷熔體底部設(shè)定4個擇優(yōu)方向?yàn)?°晶核;模擬多方向柱狀樹枝晶生長時,假設(shè)過冷熔體底部為弧狀,在熔體底部預(yù)設(shè)11個不同擇優(yōu)生長方向的晶核,相鄰兩晶核間擇優(yōu)生長方向相差7.5°。

        模擬結(jié)果如圖5所示。從圖5(a),(c),(d)可以看出,單方向柱狀樹枝晶生長時,枝晶間激烈地競爭生長,二次枝晶臂很難擴(kuò)展;液相中溶質(zhì)濃度呈山峰狀分布,兩枝晶距離越近處溶質(zhì)濃度富集越嚴(yán)重;固相中溶質(zhì)濃度分布規(guī)律與等軸晶生長過程中的相似。從圖5(b)可以看出,模擬得到多方向柱狀樹枝晶均垂直于熔體底部,實(shí)現(xiàn)了柱狀樹枝晶沿任意方向競爭生長。

        圖5 多個柱狀晶生長形貌及溶質(zhì)濃度分布(a),(b)枝晶生長形貌;(c)液相溶質(zhì)濃度;(d)固相溶質(zhì)濃度Fig.5 Columnar dendrite morphology and distribution of solute concentration(a),(b) columnar dendrite morphology;(c)distribution of liquid solute concentration;(d)distribution of solid solute concentration

        3 結(jié)論

        (1)基于二維枝晶生長的新模型,有效地再現(xiàn)了單個和多個等軸晶沿不同擇優(yōu)方向生長及單方向和多方向柱狀樹枝晶競爭生長過程中枝晶形貌變化、液相溶質(zhì)濃度和固相溶質(zhì)濃度分布情況。

        (2)等軸晶生長過程中,擾動的引入能夠促使枝晶產(chǎn)生分支,并控制二、三次枝晶生長;不同擇優(yōu)方向的多個等軸晶在生長前期互不干涉,隨著凝固排出溶質(zhì)的累積,枝晶間相互影響加劇,形成不規(guī)則的枝晶形貌。

        (3)新模型實(shí)現(xiàn)了柱狀樹枝晶垂直于弧狀邊界向熔體中心的競爭生長。

        [1]HE Y Z, DING H L, LIU L F, et al. Computer simulation of 2D grain growth using a cellular automata model based on the lowest energy principle [J]. Materials Science and Engineering:A, 2006, 429(1-2):236-246.

        [2]陳晉.基于元胞自動機(jī)方法的凝固過程微觀組織數(shù)值模擬[D].南京: 東南大學(xué),2005.

        CHEN J. Numerical simulation on solidification microstructures using cellular automaton method[D]. Nanjing: Southeast University,2005.

        [3]HONG C P, ZHU M F, LEE S Y. Modeling of dendritic growth in alloy solidification with melt convection [J]. Tech Science Press, 2006, 2(4):247-260.

        [4]LI Q, GUO Q Y, LI R D. Numerical simulation of dendrite growth and microsegregation formation of binary alloys during solidification process [J]. Chinese Physics, 2006, 15(4):778-791.

        [5]LI J J, WANF J C, YANG G C. Phase-field simulation of microstructure development involving nucleation and crystallographic orientations in alloy solidification [J]. Journal of Crystal Growth, 2007, 309(1):65-69.

        [6]LI D, LI R, ZHANG P W. A cellular automaton technique for modelling of a binary dendritic growth with convection [J]. Applied Mathematical Modelling, 2007,31(6):971-982.

        [7]YIN H, FELICELLI S D. Dendrite growth simulation during solidification in the lens process [J]. Acta Materialia, 2010, 58(4):1455-1465.

        [8]陳守東,陳敬超,呂連灝. 基于CA-FE的雙輥連鑄純鋁凝固組織模擬[J]. 材料工程, 2012, (10):48-53.

        CHEN S D, CHEN J C, LU L H. Numerical simulation of solidified microstructure of Twin-roll continuous casting pure aluminum based on CA-FE method[J]. Journal of Materials Engineering, 2012, (10):48-53.

        [9]MINORU Y, YUKINOBU N, HIROSHI H, et al. Numerical simulation of solidification structure formation during continuous casting in Fe-0.7mass%C alloy using cellular automaton method[J]. ISIJ International, 2006, 46(6):903-908.

        [10]付振南,許慶彥,熊守美. 基于概率捕獲模型的元寶自動機(jī)方法模擬鎂合金枝晶生長過程[J]. 中國有色金屬學(xué)報(bào), 2007,17(10): 1567-1573.

        FU Z N, XU Q Y, XIONG S M. Numerical simulation on dendrite growth process of Mg alloy using cellular automaton method based on probability capturing model[J]. The Chinese Journal of Nonferrous Metals, 2005, 41(6):583-587.

        [11]ZHAN X H, WEI Y H, DONG Z B. Cellular automaton simulation of grain growth with different orientation angles during solidification process[J]. Journal of Materials Processing Technology, 2008, 208(1-3):1-3.

        [12]HAKAN H, MATTI R. Microstructure evolution influenced by dislocation density gradients modeled in a reaction diffusion system [J]. Computational Materials Science, 2013, 67:373-383.

        [13]MOHSEN A Z, HEBI Y. Comparison of cellular automaton and phase field models to simulate dendrite growth in hexagonal crystals [J]. Materials Science Technology, 2012, 28(2):137-146.

        [14]ZHU M F, DAI T, LEE S Y, et al. Modeling of solutal dendritic growth with melt convection [J]. Computers and Mathematics with Applications, 2008, 55(7):1620-1628.

        [15]LUO S, ZHU M Y. A two-dimensional model for the quantitative simulation of the dendritic growth with cellular automaton method [J]. Computational Materials Science, 2013, 71:10-18.

        Numerical Simulation on Dendrite Growth During Solidification of Al-4%Cu Alloy

        ZHANG Min,XU Ai-yan,WANG Qiang,LI Lu-lu

        (School of Material Science and Engineering,Xi’an University of Technology,Xi’an 710048,China)

        A new two-dimensional cellular automata and finite difference (CA-FD) model of dendritic growth was improved, which a perturbation function was introduced to control the growth of secondary and tertiary dendrite, the concentration of the solute was clearly defined as the liquid solute concentration and the solid-phase solute concentration in dendrite growth processes, and the eight moore calculations method was used to reduce the anisotropy caused by the shape of the grid in the process of redistribution and diffusion of solute. Single and multi equiaxed dendrites along different preferential direction, single and multi directions of columnar dendrites of Al-4% Cu alloy were simulated, as well as the distribution of liquid solute concentration and solid solute concentration. The simulation results show that the introduced perturbation function can promote the dendrite branching, liquid/solid phase solute calculation model is able to simulate the solute distribution of liquid/solid phase accurately in the process of dendritic growth, and the improved model can realize competitive growth of dendrite in any direction.

        numerical simulation;cellular automaton method;dendritic morphology;solute concentration

        張敏(1967-),男,博士,教授,博士生導(dǎo)師,主要從事焊接成型過程的力學(xué)行為及其結(jié)構(gòu)質(zhì)量控制焊接和凝固過程的組織演變行為及其先進(jìn)焊接材料的研究工作,聯(lián)系地址:陜西省西安市金花南路5號西安理工大學(xué)(710048),E-mail:zhmmn@xaut.edu.cn.

        10.11868/j.issn.1001-4381.2016.06.002

        TG292

        A

        1001-4381(2016)06-0009-08

        國家自然科學(xué)基金資助項(xiàng)目(51274162);國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)(2013AA031303);陜西省教育廳專項(xiàng)科研計(jì)劃項(xiàng)目(14JK1539)

        2014-05-23;

        2015-07-23

        猜你喜歡
        方向生長
        2022年組稿方向
        2022年組稿方向
        2021年組稿方向
        碗蓮生長記
        小讀者(2021年2期)2021-03-29 05:03:48
        2021年組稿方向
        2021年組稿方向
        共享出行不再“野蠻生長”
        生長在哪里的啟示
        華人時刊(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
        亚洲国产成人av二区| 麻豆国产乱人伦精品一区二区| 欧美精品一本久久男人的天堂| 国产午夜精品综合久久久| 免费观看91色国产熟女| 亚洲av无码专区首页| 久久久久亚洲女同一区二区| 国产一区二区三区视频了| 国产精品国三级国产a| 内射欧美老妇wbb| 国产普通话对白视频二区| 女同另类激情在线三区| 少妇高潮久久蜜柚av| 国产色xx群视频射精| 色偷偷88888欧美精品久久久 | 国产成人综合久久大片| 高清毛茸茸的中国少妇| 少妇寂寞难耐被黑人中出| 亚洲人成影院在线高清| 中文日本强暴人妻另类视频 | 暖暖视频在线观看免费| 国产系列丝袜熟女精品视频| 中文字幕日韩精品亚洲精品| 丝袜美腿亚洲一区二区| 亚洲成av人片在线观看无码| 美女高潮流白浆视频在线观看| 精品国产免费一区二区久久| 夹得好湿真拔不出来了动态图| 色偷偷一区二区无码视频| 97碰碰碰人妻视频无码| 91精品久久久中文字幕| 中文字幕aⅴ人妻一区二区| 97久久精品人人妻人人| 一本久道在线视频播放| 久久亚洲中文字幕精品一区| 中文字幕日韩一区二区三区不卡 | 久久日日躁夜夜躁狠狠躁| 东京热加勒比无码少妇| 中出高潮了中文字幕| 亚洲女同av在线观看| 亚洲精品乱码8久久久久久日本 |