Al-Cu合金凝固微觀組織的三維模擬及優(yōu)化

卜曉兵, 李落星, 張立強, 朱必武, 王水平

(1. 湖南大學(xué) 汽車車身先進設(shè)計與制造國家重點實驗室，長沙 410082；2. 湖南大學(xué) 材料科學(xué)與工程學(xué)院，長沙 410082)

Al-Cu合金凝固微觀組織的三維模擬及優(yōu)化

卜曉兵1,2, 李落星1,2, 張立強2, 朱必武2, 王水平2

(1. 湖南大學(xué) 汽車車身先進設(shè)計與制造國家重點實驗室，長沙 410082；2. 湖南大學(xué) 材料科學(xué)與工程學(xué)院，長沙 410082)

采用CA-FE模型，對同一鑄件不同尺寸的Al-2%Cu合金凝固微觀組織進行三維模擬及預(yù)測，并引入空位形成能對固相擴散系數(shù)進行了優(yōu)化。結(jié)果表明，當(dāng)鑄件直徑分別為20、40和60 mm時，柱狀晶占總晶體比率分別為65.2%、52.1%和21.4%，逐漸減少，而等軸晶所占比率逐漸增加；當(dāng)換熱系數(shù)由500 W/(m2·K)增大到5 000 W/(m2·K)時，同一鑄件中不同直徑的鑄件凝固組織中柱狀晶組織比例顯著增大。模擬結(jié)果與實驗結(jié)果吻合較好，能夠較為準(zhǔn)確地反映等軸晶和柱狀晶的分布位置、比例和大小。

Al-Cu合金；3D元胞自動機；固相，擴散系數(shù)；凝固組織

鑄件中的細(xì)小晶粒能夠提高材料的強度與韌性，因而預(yù)測微觀組織形成是控制晶粒尺寸、形貌和鑄件產(chǎn)品質(zhì)量的重要前提。自O(shè)LDFIELD等對微觀組織模擬開展研究以來，微觀組織數(shù)值模擬的研究經(jīng)歷了由定性模擬到定量模擬，再由單純金屬凝固組織模擬到多元合金凝固組織模擬的歷程，并形成了不同模擬方法：確定性模擬、隨機性模擬、相場法等[1-4]。

元胞自動機法是一種隨機性模擬，該方法以形核過程的物理機制和晶體生長動力學(xué)理論為基礎(chǔ)，可以得到晶粒的尺寸、大小及分布，也可以描述柱狀晶、等軸晶的形成及柱狀晶向等軸晶的轉(zhuǎn)變過程。該方法模擬的微觀組織不依賴計算過程中的單元網(wǎng)格劃分結(jié)構(gòu)。RAPPAZ等[2]把晶粒的形核與長大的物理機制引入到模型中，提出采用基于高斯分布的連續(xù)性形核方法分別處理型壁和液體中形核，并考慮了晶核的晶體學(xué)位向關(guān)系和枝晶尖端的生長動力學(xué)，成功預(yù)測了從柱狀晶到等軸晶的轉(zhuǎn)變并得到了實驗驗證；GANDIN等[5-8]將宏觀溫度場與微觀組織形成相結(jié)合來模擬在宏觀溫度場下微觀組織的形成與演化。康秀紅等[9]、許慶彥等[10]也開展了相應(yīng)的工作，考慮了凝固過程中溫度場的分布、凝固過程溶質(zhì)再分配以及對流等因素的影響，實現(xiàn)了對金屬凝固組織的計算機模擬。ZHU等[11-21]改進了CA方法，將枝晶尖端生長所需的過冷度用 KGT模型來描述，修正了枝晶尖端生產(chǎn)僅依賴于熱過冷的模型。

本文作者采用元胞自動機耦合有限差分的方法(CA-FE)，模擬了同一鑄件不同半徑Al-2%Cu合金凝固組織的形貌，探究并模擬同一鑄件中不同壁厚對晶粒大小、密度形核參數(shù)等對柱狀晶向等軸晶轉(zhuǎn)變的影響；預(yù)測金屬凝固過程中柱狀晶等軸晶的轉(zhuǎn)變；計算合金的固相擴散系數(shù)時，通過引入了空位形成能提高了模擬的準(zhǔn)確性。

1 固相溶質(zhì)擴散系數(shù)的優(yōu)化

凝固過程中，溶質(zhì)分布是影響枝晶生長和形貌的重要因素，溶質(zhì)擴散的模擬是進行微觀組織模擬的一個重要方面。

模擬溶質(zhì)擴散時，可以采用下面兩種方法。

1) 溶質(zhì)擴散采用Scheil模型：假設(shè)固相中不存在擴散，液相完全混合。

Scheil模型表達(dá)式為

固相擴散系數(shù)是 CA模擬中重要的晶體生長參數(shù)，影響濃度場的計算。物質(zhì)的遷移可通過對流和擴散兩種方式進行。但實際凝固過程中固相中是存在擴散的，固體中不發(fā)生對流，擴散是唯一的物質(zhì)遷移方式。為了提高模擬的準(zhǔn)確性，采用下面的方法。

2) 用液相擴散系數(shù)和固相擴散系數(shù)描述液相和固相擴散過程。在固相晶體中存在空位，這些空位的存在使原子遷移變得更容易，從而產(chǎn)生原子空位跳躍遷移能，該遷移能等于原子的擴散激活能Q，通過菲克第一定律可以計算出固相擴散系數(shù)。

而實驗測得的Al-Cu合金固相擴散系數(shù)和菲克定律計算的值差別很大，這個差別可以用銅在鋁中以空位機制擴散解釋，除了考慮空位遷移能外，還必須考慮擴散原子近旁空位的形成能。銅原子周圍空位形成的幾率可用空位平衡濃度(ne/N)表示，所以真正的擴散系數(shù)應(yīng)為

2 基于CA方法的宏-微觀耦合模擬

2.1 宏觀溫度場的計算

溫度場的模擬計算是微觀組織模擬的前提。液體金屬充滿鑄型以后，金屬和鑄型之間的導(dǎo)熱主要以不穩(wěn)定的導(dǎo)熱方程進行。宏觀溫度場的計算采用有限元方法(FE)求解，三維宏觀溫度場不穩(wěn)定導(dǎo)熱方程為

式中：ρ為密度，g/cm3；cp為定壓比熱容，J/(kg·℃)；θ為溫度，℃；t為時間，s；λ為熱傳導(dǎo)率，W/(m·℃)；Q為熱源項，；L為熔化潛熱，J/kg；φs為固相率；x，y和z為坐標(biāo)值，m。

輻射換熱按照Stefen-Boltzman定律： q =εσT4，0T為表面的絕對溫度，ε為輻射黑度；σ0為Stefen-Boltzman常數(shù)。

2.2 形核模型

形核采用 RAPPAZ等[22-23]提出的基于高斯分布的連續(xù)形核模型，利用統(tǒng)計方法，假設(shè)形核現(xiàn)象發(fā)生在一系列不同的形核位置上(忽略形核本身所需要的時間，即形核是瞬時出現(xiàn))，這些形核位置由連續(xù)非離散的分布函數(shù)dn/d(ΔT)來描述。在某一過冷度ΔT所形成的晶粒密度n(ΔT)即可由對該分布曲線的積分求得：

式中：ΔTN為平均形核過冷度；ΔTσ為標(biāo)準(zhǔn)曲率過冷度；nmax為異質(zhì)形核襯底的數(shù)目。

2.3 生長模型

生長模型采用由 KURZ等[24]提出的枝晶尖端生長動力學(xué)模型，即KGT模型

式中：ΔTc為溶質(zhì)擴散引起的成分過冷度；ΔTt為熱擴散引起的過冷度；ΔTk為動力學(xué)過冷度；ΔTr為固/液界面引起的曲率過冷度。

在正常凝固條件下，后三項數(shù)值較小，近似有ΔT=ΔTc，?為枝晶尖端溶質(zhì)過飽和度，R為枝晶尖端半徑，Gc為枝晶前沿液相中的溶質(zhì)濃度梯度，G為溫度梯度，Γ為Gibbs-Thompson系數(shù)，m為液相線斜率，Pe為溶質(zhì)的 Peclet數(shù)的函數(shù)，在低速生長時取 1，D為液相內(nèi)的溶質(zhì)擴散系數(shù)，v為枝晶尖端生長速度。

3 實驗

鑄錠材料選擇Al-2%Cu合金。鑄型形狀及尺寸如圖1所示。為了減少縱向熱對流對微觀組織的影響，鑄型材料選用鋼板，鑄型底部和四周用保溫石棉包好，產(chǎn)生絕熱效果。冷卻過程采取空冷方式。圖1中A、B、C位置沿橫切面剖開，分別對同一鑄件中直徑為20、40和60 mm的橫切面進行粗磨、細(xì)磨和拋光后采用氫氟酸進行腐蝕 5 s，便可以觀察到不同直徑的晶粒組織。

圖1 鑄型形狀示意圖Fig.1 Shape and dimensions of experimental casting mold(unit: mm)

4 模擬結(jié)果及分析

采用CAFé模型，計算尺寸與實驗鑄件一致。鑄件尺寸為600 mm×100 mm×150 mm，宏觀網(wǎng)格劃分為1 mm，微觀計算時網(wǎng)格再細(xì)分為100×100×100。在模擬宏觀溫度場時，鑄件和鑄型之間的換熱系數(shù)設(shè)置為2 kW/(m2·K)，鑄型頂部與空氣的換熱系數(shù)設(shè)置為100 W/(m2·K)。形核參數(shù)如表1所列。

表1 所使用材料的性能和模擬參數(shù)Table1 Material properties and model parameters used in simulations

4.1 不同鑄造半徑下凝固組織的模擬分析

鑄件尺寸是影響凝固組織的重要因素，由圖2可以看出，鑄件微觀組織由表層細(xì)晶區(qū)，柱狀晶區(qū)和中心等軸晶區(qū)3部分組成。

圖2 鑄件橫切面三維凝固組織和縱切面三維凝固組織模擬圖Fig.2 Transverse (a) and upright section (b) of threedimensional results

液態(tài)Al-Cu合金注入鑄型后，由于型壁溫度低，與型壁接觸的一層很薄熔液產(chǎn)生強烈過冷，同時型壁還可作為非均勻形核的基底，導(dǎo)致立刻形成大量的晶核，并且這些晶核迅速長大，互相接觸，形成細(xì)小、方向雜亂的等軸晶粒。這些晶核很快就因彼此碰撞而無法繼續(xù)生長，由圖 2(a)可以看出，在不同鑄件尺寸的表面均形成一層很薄很細(xì)的等軸晶區(qū)。

隨著細(xì)晶區(qū)外殼的生成，鑄件的型壁被熔液加熱而不斷升溫，使剩余液體的冷卻速度變慢，并且由于結(jié)晶釋放潛熱，細(xì)晶區(qū)前沿液體的過冷度減小，形核變得困難，只有細(xì)晶區(qū)中現(xiàn)有的晶體向液體中生長，且只有一次軸垂直于型壁(散熱最快的方向)的晶體才能得到優(yōu)先生長，而其他取向的晶粒，由于受鄰近的限制而不能發(fā)展。因此，這些與散熱相反方向的晶體擇優(yōu)生長而形成柱狀晶區(qū)。隨著鑄造尺寸的增加，柱狀晶占總晶體比率減小，分別為 65.2%、52.1%和21.4%，而等軸晶比率增加。這是由于柱狀晶生長到一定程度，越小尺寸的前沿液體越遠(yuǎn)離型壁，散熱困難，冷卻速度變慢，而且熔液中的溫差隨之減小，這將阻止柱狀晶的快速生長，當(dāng)整個熔液溫度降到熔點以下時，熔液中出現(xiàn)許多晶核并沿著各個方向生長，形成等軸晶區(qū)。

如果鑄型不加熱，采用室溫，則熔體冷卻速度快，凝固過程中，已形核的晶粒根據(jù)其不同的相對位置和擇優(yōu)取向，鄰近枝晶較早開始相互接觸形成晶界并制約了對方的生長，使晶粒沒有足夠的時間長大，故形成的晶粒較為細(xì)小，在鑄造尺寸較小時(20 mm)，以等軸晶為主。而實驗中鑄型加熱至250 ℃，鑄造冷卻速度慢，形核率低，晶核具有充足的時間生長，因此形成的晶粒比較粗大。在鑄件尺寸較小時，以柱狀晶為主。

圖3 直徑20 mm鑄件的凝固組織模擬圖Fig.3 Predicted grain structures of d20 mm casting: (a) h=500 W/(m2·K); (b) h=1 000 W/(m2·K); (c) h=2 000 W/(m2·K); (d) h=5 000 W/(m2·K)

4.2 不同冷卻強度下凝固組織的預(yù)測

冷卻強度也是影響鑄件凝固組織的主要因素之一，常被用來控制凝固組織的形成。實驗中，鑄件壁面與鑄型的換熱系數(shù)設(shè)置為2 kW/(m2·K)，并且分別在不同換熱系數(shù)下對直徑20、40和60 mm的鑄件進行了模擬。

對比圖 3～5可以看出，隨著換熱系數(shù)(h)從 500 W/(m2·K)增大到 5 000 W/(m2·K)，不同半徑的鑄件凝固組織中柱狀晶比例顯著增大，而同樣半徑的鑄件中晶粒度變化并不大。這是由于當(dāng)換熱系數(shù)較小時，液體金屬在鑄型壁上形成穩(wěn)定凝固殼層的時間較長，晶體從鑄型壁處游離出來的過程也較長，游離晶粒增多，而游離到液體中的晶粒又可以阻止柱狀晶的形成，進而使柱狀晶區(qū)域變?。环粗?，當(dāng)換熱系數(shù)較大時，液態(tài)金屬在鑄型壁上形成穩(wěn)定的凝固殼層時間較短，一旦穩(wěn)定的凝固殼層形成，柱狀晶就直接由表面細(xì)等軸晶凝固層中的晶粒為基底向內(nèi)生長，發(fā)展成由外向內(nèi)生長的柱狀晶區(qū)，有利于柱狀晶的形成。

圖4 直徑40 mm鑄件的凝固組織模擬圖Fig.4 Predicted grain structures of d40 mm casting: (a) h=500 W/(m2·K); (b) h=1 000 W/(m2·K); (c) h=2 000 W/(m2·K);(d) h=5 000 W/(m2·K)

圖5 直徑60 mm鑄件的凝固組織模擬圖Fig.5 Predicted grain structures of d60 mm casting: (a) h=5 00 W/(m2·K); (b) h=1 000 W/(m2·K); (c) h=2 000 W/(m2·K);(d) h=5 000 W/(m2·K)

4.3 凝固組織模擬實驗驗證

鑄件尺寸如圖1所示，由3個不同壁厚和半徑的圓柱組成：1) 壁厚150 mm，直徑60 mm；2) 壁厚450 mm，直徑40 mm；3) 壁厚750 mm，直徑20 mm。鑄件的澆注溫度720 ℃，鑄型溫度250 ℃, 鑄件三維凝固組織模擬結(jié)果如圖2所示，截面厚度均為4 mm。比較圖6～8可以發(fā)現(xiàn)，采用CAFé模型進行模擬的結(jié)果與實驗結(jié)果能較好吻合，能夠較為準(zhǔn)確的反應(yīng)出等軸晶和柱狀晶的分布位置、比例和大小等。如圖6所示，對直徑60 mm的橫截面進行微觀組織模擬時，橫截面為圖1中C位置(高125 mm處，直徑60 mm和直徑40 mm的交界面)斷面。此橫截面外側(cè)20 mm厚的部分直接鑄型接觸，形成較粗的等軸晶；而與直徑40 mm鑄件直接接觸的部分，先形成柱狀晶，再形成較細(xì)的等軸晶。所以，三維CAFé模型是模擬鑄件凝固組織的一種有效模型，應(yīng)用CAFé模型研究了澆注溫度和冷卻強度對晶粒組織的影響。

圖6 直徑60 mm鑄件凝固組織的模擬結(jié)果和實驗結(jié)果Fig.6 Predicted (a) and experimental(b) solidification grain structures of d60 mm casting

圖7 直徑40 mm鑄件凝固組織的模擬結(jié)果和實驗結(jié)果Fig.7 Predicted (a) and experimental (b) solidification grain structures of d40 mm casting

圖8 直徑20 mm鑄件凝固組織的模擬結(jié)果和實驗結(jié)果Fig.8 Predicted (a) and experimental (b) solidification grain structures of d 20 mm casting

5 結(jié)論

1) 采用有限元和元胞自動機結(jié)合的CA-FE模型，通過加入空位形成能的固相擴散系數(shù)，對Al-Cu合金的凝固組織進行了三維模擬和預(yù)測，其結(jié)果能夠準(zhǔn)確地反映等軸晶和柱狀晶的分布位置、比例和大小等。

2) 隨著鑄件尺寸的增加，鑄件凝固組織中柱狀晶占總晶體比率減小，分別為65.2%、52.1%和21.4%，而等軸晶比率增加。

3) 當(dāng)換熱系數(shù)由從 500W /(m2·K)增大到 5 000 W/(m2·K)時，鑄件凝固組織中柱狀晶比例顯著增大，而同樣鑄造半徑的晶粒度變化并不大。

REFERENCES

[1] OLDFIELD W. A quantitative approach to casting solidification:Freezing of cast iron[J]. ASM Trans, 1966, 59(2): 945-960.

[2] RAPPAZ M, GANDIN C A. DESBIOLLES J L, THEVOZ P H.Prediction of grain structures in various solidification processes[J]. Metallurgical and Materials Transactions A, 1996,27(3): 695-705.

[3] ZHU P P, SMITH R W. Dynamic simulation of crystal growth by Monte-Carlo method[J]. Acta Metall, 1992, 40(12):3369-3379.

[4] PENROSE O, FIFE P C. Thermodynamically consistent models of phase field type for the kinetics of phase transitions[J].Physical D, 1990, 43(10): 44-62.

[5] GANDIN C A, DESBIOLLES J L, RAPPAZ M, THEVOZ P H.A three-dimensional cellular-finite element model for the prediction of solidification grain structures[J]. Metall Master Trans A, 1999, 30(12): 3153-3165.

[6] GUILLEMOT G, GANDIN, C A, BELLET M. Modeling of casting, welding and advanced solidification processes Ⅺ[M].Warrendale, PA: TMS, 2006: 307-314.

[7] GUILLEMOT G, GANDIN C A, COMBEAU H. Modeling of macrosegregation and solidification grain structures with a coupled cellular automation finite element model[J]. ISIJ Int,2006, 46: 880-895.

[8] GUILLEMOT G, GANDIN C A, BELLET M. Interaction between single grain solidification and macroseg-regation:Application of a cellular automaton-finite element model[J].Crystal Growth, 2007, 303(1): 58-68.

[9] 康秀紅, 杜 強, 李殿中, 李依依. 用元胞自動機與宏觀傳輸模型耦合方法模擬凝固組織[J]. 金屬學(xué)報, 2004, 40(5):452-456.

KANG Xiu-hong, DU Qiang, LI Dian-zhong, LI Yi-yi.Modeling of the solidification microstructure evolution by coupling cellular automaton with macrotransport model [J]. Acta Metallrugica Sinica, 2004, 40(5): 452-456.

[10] 許慶彥, 馮偉明, 柳百成, 熊守美. 鋁合金枝晶生長的數(shù)值模擬[J]. 金屬學(xué)報, 2002, 38(8): 799-803.

XU Qing-yan, FENG Wei-ming, LIU Bai-cheng, XIONG Shou-mei. Numerical simulation of dendrite growth of aluminum alloy[J]. Acta Metallrugica Sinica, 2002, 38(8): 799-803.

[11] ZHU M F, HONG C P. A modified cellular automation model for the simulation of dendritic growth in solidification of alloys[J]. ISIJ Int, 2001, 41(5): 436-445.

[12] 陳 晉, 朱鳴芳, 孫國雄. 用 CA 方法模擬過冷熔體中自由樹枝晶的生長[J]. 金屬學(xué)報, 2005, 41(8): 799-803.

CHEN Jin, ZHU Ming-fang, SUN Guo-xiong. Numerical simulation on free dendrite growth in undercooled melt using cellular automaton method[J]. Acta Metallurgica Sinica, 2005,41(8): 799-803.

[13] 朱鳴芳, 于 金, 戴 挺. 金屬凝固過程數(shù)值模擬的最新進展[J]. 鑄造, 2005, 54(2): 115-120.

ZHU Ming-fang, YU Jin, DAI Ting. Recent progress in numerical modeling of solidification processes[J]. Foundry, 2005,54(2): 115-120.

[14] 陳 晉, 朱鳴芳, 孫國雄. 過冷熔體中的樹枝晶分枝機制模擬[J]. 鑄造, 2005, 54(9): 895-898.

CHEN Jin, ZHU Ming-fang, SUN Guo-xiong. Numerical simulation on dendrite branching in undercooled melt[J].Foundry, 2005, 54(9): 895-898.

[15] 蘇 敏, 孫東科, 朱鳴芳. Lattice Boltzmann方法在枝晶生長模擬中的應(yīng)用[J]. 機械工程材料, 2007, 31(2): 75-78.

SU Min, SUN Dong-ke, ZHU Ming-fang. Application of Lattice Boltzmann method in modeling of dendritic growth[J]. Materials For Mechanical Engineering, 2007, 31(2): 75-78.

[16] ZHU M F, CAO W, CHEN S L, HONG C P, CHANG Y A.Modeling of microstructure and microsegregation in solidification of multi-component alloys[J]. Phase Equilibria and Diffusion, 2007, 28(1): 130-138.

[17] ZHU M F, STEFANESCU D M. Virtual front tracking model for the quantitative modeling of dendritic growth in solidification of alloys[J]. Acta Mater, 2007, 55(5): 1741-1755.

[18] ZHU M F, HONG C P, STEFANESCU D M, CHANG Y A.Computational modeling of microstructure evolution in solidification of aluminum alloys[J]. Metall Mater Trans B, 2007,38(4): 517-524.

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

[20] 楊朝蓉, 孫東科, 潘詩琰, 戴 挺, 朱鳴芳. CA-LBM模型模擬自然對流作用下的枝晶生長[J]. 金屬學(xué)報, 2009, 45(1): 43-50.

YANG Chao-rong, SUN Dong-ke, PAN Shi-yan, DAI Ting,ZHU Ming-fang. CA-LBM model for the simulation of dendritic growth under natural convection[J]. Acta Metallurgica Sinica,2009, 45(1): 43-50.

[21] PAN S Y, ZHU M F. A three-dimensional sharp interface model for the quantitative simulation of solutal dendritic growth[J].Acta Materialia, 2010, 58(1): 340-352.

[22] RAPPAZ M, GANDIN C A. Probabilistic modeling of microstructure formation in solidification processes[J]. Acta Materialia, 1993, 41(12): 345-360.

[23] GANDIN C A, RAPPAZ M. A coupled finite element-celluar automaton model for the prediction of dendritic grain structure in solidification processes[J]. Acta Metall Mater, 1994, 42(7):2233-2246.

[24] KURZ W, GIOVANOLA B, TRIVEDI R. Theory of microstructural development during rapid solidification[J]. Acta Metall Mater, 1986, 34(5): 823-830.

Three-dimensional numerical simulation and optimization of solidification structure of Al-Cu alloy

BU Xiao-bing1,2, LI Luo-xing1,2, ZHANG Li-qiang2, ZHU Bi-wu2, WANG Shui-ping2
(1. State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,Hunan University, Changsha 410082, China;2. College of Materials Science and Engineering, Hunan University, Changsha 410082, China)

A 3D cellular automaton-finite element (CA-FE) method was used to predict the solidification structure of Al-Cu alloy at different casting radii, and optimized the solution diffusion coefficient in solid phase. The results show that when the casting radius changes from 10 mm to 30 mm, the proportion of columnar grains changes from 65.2% to 21.4%.With increasing heat transfer coefficient from 500 W/(m2·K) to 5 000 W/(m2·K), the proportion of columnar grains also gradually increases. The simulated results are in accord with the experimental ones well, and can accurately reflect the distribution, proportion, size of equiaxed and columnar grains.

Al-Cu alloy; 3D cellular automation method; solid phase; solution diffusion coefficient; solidification structure

TG292

A

1004-0609(2011)09-2195-07

國家自然科學(xué)基金資助項目(51075132)；湖南省杰出青年基金資助項目(09JJ1007)；高等學(xué)院博士點專相科研基金資助項目(20090161110027)；湖南大學(xué)汽車車身先進設(shè)計制造國家重點實驗室自主研究課題項目(60870005)

2010-07-19；

2010-12-28

李落星；教師，博士；電話：0731-88821445; E-mail：llxly2000@163.com

(編輯 何學(xué)鋒)