粉末共注射成形充模流動過程前沿位置及場分布的數值模擬

封 娟，何 浩，李益民，王光耀

(1. 中南大學 數學科學與計算技術學院，長沙 410075；2. 中南大學 粉末冶金研究院，長沙 410083)

粉末共注射成形充模流動過程前沿位置及場分布的數值模擬

封 娟1，何 浩2，李益民2，王光耀2

(1. 中南大學 數學科學與計算技術學院，長沙 410075；2. 中南大學 粉末冶金研究院，長沙 410083)

對粉末共注射充模流動過程進行數值分析和實驗驗證。采用實驗和數值擬合的方法確定芯/殼層界面的厚度，并運用改進的控制體積法對芯、殼層喂料前沿進行追蹤。采用有限元和有限差分法對控制方程組進行數值求解，用Matlab進行程序開發，獲得芯、殼層充模過程中的熔體前沿分布以及溫度場和壓力場的分布情況。將模擬結果與實驗結果進行對比分析，發現在充填初期，模擬的喂料前沿位置與實驗較為吻合，但隨著充填的進行，兩者偏差增大，其原因可能是在模擬過程中沒有考慮注射坯的收縮。

粉末共注射成形；數值模擬；控制體積法；充模流動

粉末共注射成形是由傳統粉末冶金技術與現代塑料共注射技術相結合而產生的金屬零部件近凈成形新技術[1]。粉末共注射成形通過不同方式將多種組分在模具內復合注射成形，該技術整合了不同組分的優越性能，具有高性能、低成本、短流程、凈成形的工藝優勢。在零部件的開發設計及生產方面具有重要的應用價值。因此成為國內外學術界研究的熱點。目前，國內外研究機構如英國的Cranfield大學[2]、德國第三研究院、日本Taisei-Kogyo公司、德國的博世集團以及中國的中南大學粉末冶金國家工程研究中心[3]等都對其產品開發和過程工藝理論進行了研究。充模過程是粉末共注射成形的關鍵環節，其模擬研究為優化工藝參數、改進注射模具、預測產品質量以及指導生產具有重要意義。但由于加入了大量粉末，粉末共注射流變行為比塑料共注射要復雜得多，同時喂料的微小變化也會改變其流動行為和產品最終性能，均加大了粉末共注射成形研究和控制的難度。因此，對該過程的研究顯得非常必要。

目前，國內在這方面的研究起步較晚，關于共注射的研究大都是塑料共注射成形的充模模擬研究，對粉末共注射成形充模流動的研究基本是空缺的。對粉末共注射充模流動的模擬可以在注塑共注射充模流動模擬的研究上開展。關于注塑充模流動，LI等[4?5]建立了塑料共注射成形充填過程的物理模型，并給出了數值解法。周國發等[6]和匡唐清等[7]將塑料共注射成形充填過程視為多組分分層流動，采用控制體積法追蹤熔體芯、殼層前沿。塑料共注射成形模型推導出的控制方程涉及壓力和溫度在三維方向上的梯度或散度，場的求解需要解決芯層熔體在充填方向以及厚度方向上界面變化的問題，但目前在充填方向上研究較多，對厚度方向研究卻較少。由于共注射成形涉及兩種喂料且存在芯/殼層界面，這個界面位置對熔體充填的影響在壓力控制方程中有明顯體現，仍采用上述忽略厚度方向的方法顯然是不合理的。

因此，直接采用傳統的塑料共注射成形模擬方法對粉末共注射成形數值模擬問題研究存在欠缺，關鍵在于芯/殼層界面厚度位置無法確定。本文作者通過實驗與分析提出一種確定芯/殼層界面厚度位置的方法，基于改進的控制體積法，對共注射芯、殼層界面進行同時追蹤，從而實現粉末共注射成形過程的動態模擬，獲得芯、殼層熔體在型腔中壓力場、溫度場的分布。最后將模擬結果和實驗結果進行對比分析。

1 數學模型

圖1 共注射型腔示意圖Fig. 1 Geometry of plat for co-injection molding (mm)

由于共注射成形過程以及共注射熔體流變性質的復雜性，本研究采用塑料共注射成型過程的研究方法，將粉末共注射喂料充填過程視為Hele-Shaw流動[8?9]。在傳統單組分注射成形充填過程分析所作假設的基礎上，對三維薄壁型腔充填過程進行考慮時，針對共注射具有芯、殼層界面的特點作如下假設：1) 將實際復雜的型腔簡化為薄板狀型腔，型腔在x和y方向上的厚度尺寸遠大于方向上的厚度尺寸，如圖1所示。由于z方向上的尺寸較小，可以忽略喂料熔體在充填過在厚度方向上的速率；2) 喂料熔體中粉末與黏結劑混合均勻，喂料在充填過程中不發生兩相分離，并且將喂料視為不可壓縮的連續介質，即忽略密度的變化，且物性參數不變；3) 喂料在充填過程中，充填方向以對流傳熱為主，忽略充填方向的熱傳導，在厚度方向上以熱傳導為主，忽略厚度方向上的對流傳熱；4) 忽略慣性力、彈性力和重力喂料在充填過程的影響， 假定壓力沿厚度方向不變，殼層前沿熔體的壓力為零；5) 假設喂料的流動是層流，且流動是充分發展的，無滑移邊界，忽略流動前沿處的“噴泉”效應；6) 忽略喂料充模過程中的凝固現象；7) 忽略物性參數隨時間的變化；8) 熔體在芯殼層處有明顯的界面存在，在界面處無明顯的混合，即混合可以忽略。結合流體力學的基本方程[10?11]，得到連續性方程：

壓力控制方程：

式中：p、T與S代表熔體的壓強、溫度及流導；δ為芯層厚度的一半；ui和vi分別表示在界面處、芯層或者殼層熔體在x和y方向上的速度分量；b為型腔厚度的一半；ηc和ηs分別為芯、殼層的黏度；T為喂料熔體的溫度；iρ、ci、λi、ηi和i分別為熔體i的密度、比熱容、熱傳導系數、黏度和剪切速率。

根據流動模型的簡化及注射成形充填過程中的實際情況，設邊界條件如下：1)殼層流動前沿處，p=0；2)芯/殼層界面處，3)沿型腔側壁，滿足不滲透條件，在型腔入口處，Q=Qi，T=T0。其中，Q、Qi和T0分別表示流率、澆口處的流率和注射溫度。

壓力場、溫度場的求解交替進行，兩者通過喂料熔體黏度模型耦合。考慮熔體為非牛頓流體，其黏度模型采用Cross-WLF黏度模型[12]：

2 共注射成形前沿界面和芯/殼層熔體界面位置的確定方法

2.1 芯/殼層界面厚度方向的分析

前面推出的控制方程組，表面上是一個三維方程組，但由于所作出的簡化在推導過程中忽略了厚度方向上的速度，從而忽略了在厚度方向上的傳導，因此該方程組實際是一個介于二維和三維的方程組。雖然這種2.5維的公式推導在一定程度上簡化了控制方程，但實際上并沒有忽略三維空間中厚度方向上的位置函數δ的確定。對于傳統單組份注射成形不存在這個厚度位置函數δ，因而可以對三維型腔采用三角單元平剖，在厚度方向分成等高的w層，各層之間的流動通過能量交換聯系起來，然后采用傳統的控制體積法就可以實現前沿界面的追蹤。與其相比，共注射成形由于芯層和殼層組分的不同，在芯/殼層處有一個明顯的界面也需要進行追蹤，如圖2所示。由于芯/殼層界面厚度方向的位置對于充填情況的影響不可忽略，如果采用上述傳統研究方法處理共注射成形，就無法確定位置函數δ在任意時刻的值，因此是不合理的。

圖2 共注射橫剖、縱剖示意圖Fig. 2 Sketch map of cross and vertical section of co-injection molding

為了分析芯/殼層界面厚度隨時間的變化情況，采用316L(40%)不銹鋼喂料為殼層、316L(60%)不 銹鋼喂料為芯層，選擇不同的芯層注射時間進行實驗，觀測芯層熔體厚度的變化情況。圖3所示為實驗觀測到的穿伸位置分別為26.1、37.8、47.0、60.0和76.1 mm處芯層熔體厚度隨著注射時間的變化情況。由圖3可知，殼層的預填充量是影響芯層熔體厚度方向分布情況的首要因素，并且在某一固定的穿伸位置 處，芯層熔體在厚度方向上的充填時間僅需0.2～0.4 s，與芯層熔體在注射水平方向上的充填時間相比，充填時間較短。通過圖3還可以發現，在固定的穿伸位置處，隨著充填過程的進行，在開始注射的0.2～0.4 s內芯層熔體在厚度方向上的變化與注射時間大致呈線性關系，在厚度達到最大值后，不再隨著注射時間延長發生變化。

圖3 不同穿伸位置處共注射充填過程中芯層厚度的變化情況Fig. 3 Change of thickness of core melt during filling process

因此，基于以上實驗研究和分析，對型腔進行平剖在x和y方向上追蹤熔體的充填情況時，通過實驗先確定芯層熔體在厚度方向上的最大值，再在數值模擬分析時，假定芯、殼層充填量為定值的情況下，厚度方向上芯殼層界面位置是勻速變化的，從而芯殼層界面厚度方向的位置(δ)可以通過式(9)確定

將δ代入到控制方程中，按照數值求解的方法即可獲得任意時刻的芯、殼層熔體前沿的分布情況以及場的分布情況。由此便可確定任意時刻(t)厚度方向上芯/殼層界面的位置。

在殼層充填量一定的情況下，通過實驗數據擬合可得芯/殼層界面在厚度方向上的最大值bmax，理論注射時間ttheory=0.3 s，則在數值模擬中假定芯層界面位置δ隨時間的變化率為k=bmax/ttheory。

2.2 改進的控制體積法追蹤前沿界面

在共注射成形充填過程中，除了判斷控制體積是否被充填滿之外，還需要判斷控制體積被何種熔體充填，為此對傳統的控制體積法[13?14]做出一些改進。設f為填充因子，其定義為

式中：?Vi為節點i的控制體內已填充喂料熔體的體積；Vi為控制體的體積。根據填充因子f值的大小判斷節點的屬性，即f=0，空節點(控制體內無流體)；0＜f＜1，前沿節點(控制體未被殼層熔體充滿)；f=1，內部節點或入口節點(控制體完全被殼層熔體充滿)；0＜f＜2，內部節點(控制體被芯、殼層體充滿)；f=2，內部節點或入口節點(控制體完全被芯層熔體充滿)。

首先計算所有前沿節點控制體內的體積流率，由于所有前沿節點和空節點壓力為零，因此其體積流率(qi)的計算[15]如下：

需要指出的是，在注射時先注入一定量的殼層熔體，再將芯殼層兩種熔體一起注入。當注射殼層熔體時，最小步長的計算采用式(12a)，當注射芯層熔體時，由于型腔中存在兩個前沿，即殼層前沿、芯/殼層界面。因此，最小步長的取值需要取殼層前沿、芯/殼層界面充填的最小值，計算需采用式(12a)和式(12b)。

當一個前沿節點充滿而變為內部節點時，與之相鄰的所有空節點成為新的前沿節點，同時更新所有其他前沿節點的填充因子：

即可進入下一步壓力場的計算，重復這個步驟便可動態地跟蹤流動前沿的發展。

這樣，在數學模型、邊界條件以及熔體前沿界面追蹤的問題解決后，獲得任意時刻熔體在型腔中的流動情況并進行分析。

3 數值求解

在充模流動過程中，對場的求解以及熔體前沿的追蹤作如下處理：對時間步長的選取采用的是使每個時間剛好有一個前沿節點被充滿。此時，與其相鄰的所有空節點都變成前沿節點。假定流動開始時第一個控制體積已被完全充滿，此時可認為熔體處于等溫狀態，這樣便獲得初始時刻的前沿位置和溫度場。由此可計算出此時的剪切黏度η、流通率s和平均體積流量Qa，然后迭代求解壓力場。在獲得壓力場數據后計算下一時間間隔和熔體的前沿位置，以便計算下一時刻的溫度場。獲得溫度場數據后，程序返回計算新時刻的壓力場，如此循環直至型腔被熔體充滿。

首先對流動區域進行剖分，沿型腔中面將其劃分成一系列的三角形單元的結構性網格，在厚度方向上由中性面向模壁分成高度相同的w層。采用有限元法對壓力場進行求解，設三角形單元內的壓力按線性分布，單元向其節點控制體積所提供的流率qj(j=1, 2, 3)，則

采用Galerkin法[16]得到單元的局部剛度矩陣

式中：kij為線性插值函數矩陣。對每個單元的局部剛度矩陣累加，得到總體剛度矩陣，通過迭代求解確定該時刻的壓力場。

求解溫度場時需要考慮型腔被充填的具體情況，對于分別被芯、殼層熔體充填的控制體，能量守恒方程(4)中的物性參數iρ、ci、λi和ηi，分別取為芯、殼層熔體對應的物性系數。而對于控制體被芯、殼層熔體充填的區域，即1＜f≤2處，該處控制體遵守的能量守恒方程的物性參數值取芯、殼層熔體物性參數的加權平均值，即

能量方程的物性參數確定后，采用有限元與有限差分混合法對溫度場進行求解。在能量方程中熱對流項和黏性項按體積作加權平均，其中熱對流項采用上風法進行差分，以保證數值計算的穩定性。能量方程的具體差分形式為

式中：?z為z坐標變化量；λ為熱傳導系數；F1為熱對流項差分矩陣；F2為黏性項差分矩陣。

4 數值模擬的實現與驗證

4.1 芯/殼層界面在厚度方向上最大值的擬合

通過2.1的分析可知，為了獲得芯/殼層界面δ的值，需要先確定芯/殼層界面在厚度方向上最大值，本文作者選用316L(40%)為殼層、316L(60%)為芯層進行實驗，其中殼層充填量為60%。實驗和模擬選用的型腔為薄板型腔，具體尺寸如圖1所示，實驗的工藝參數見表1所示。實驗用注射成型機是海天塑機集團有限公司生產的HTSJ160B共注射成型機。研究假設喂料在充填過程中，充填方向以對流傳熱為主，對流傳熱系數通過黏度測定獲得。由實驗獲得不同芯層穿伸距離處芯層的厚度，采用Origin軟件對實驗數據進行6次多項式函數擬合，其形式如式(21)所示。多項式參數值見表2，實驗擬合曲線圖如圖4所示。擬合結果與實驗結果的相關系數R=0.992 95。

式中：x表示芯層穿伸距離，z表示芯層厚度的最小值，即bmax。通過式(21)可獲得不同芯層穿伸距離處芯層的厚度最大值bmax，然后通過本文作者所提出的芯層熔體在厚度方向上的變化與注射時間呈線性關系的假設，可得任意時刻在任意芯層穿伸距離處芯殼層界面的δ值。

4.2 界面的追蹤實現以及溫度場、壓力場的實現

將4.1節求得的芯殼層界面δ值代入控制方程中，按照數值求解的方法進行編程獲得任意時刻芯殼層熔體前沿的分布情況以及場的分布情況。所開發的Matlab程序在惠普公司生產的HPJ400工作站上運行[17]。圖5所示為共注射成形充填過程，芯、殼層喂料在不同注射時間下的前沿運動分布圖，其中，A、B、C、D、E、F為實驗曲線，a、b、c、d、e、f為模擬曲線。由圖5可以看出，模擬的結果大致反映出芯、殼層喂料充填過程中前沿的運動情況，模擬結果與實驗結果吻合較好。對共注射充填過程的模擬，隨著充填的進行，模擬結果與實驗結果偏差增大。原因是模擬過程中沒有考慮注射坯的收縮，造成模擬獲得的注射坯與實驗獲得的注射坯在尺寸上存在差異(見曲線D、E、F和d、e、f)。

表1 模擬和實驗工藝參數Table 1 Injection parameters of simulation and experiment

表2 多項式擬合參數值Table 2 Fitted values of polynomial

圖4 實驗擬合曲線圖Fig. 4 Fitting curve of experiment data

圖5 芯、殼層熔體在不同芯層注射時間下的前沿分布圖Fig. 5 Front distribution of skin and core melt at various time (A, B, C, separately, is front position of core melt at 0.8, 1.5, 2.0 s through experiment; a, b, c, separately, is front position of core melt at 0.902, 1.653, 2.194 s through simulation; D, E, F, separately, is front position of skin melt at 0.8, 1.5, 2.0 s through experiment; d, e, f, separately, is front position of skin melt at 0.902, 1.653, 2.194 s through simulation)

圖5中芯、殼層熔體在不同芯層注射時間下的前沿分布實驗曲線A、B、C、D、E、F是通過如下方法獲得的：設置芯層注射時間分別為0.8、1.5和2.0 s，在HTSJ160B夾層塑料注射機上進行實驗，得到3個不同芯層注射時間的注射坯。由于研究假定熔體在型腔中的流動為層流，對獲得的注射坯進行平剖，并沿注射方向將平剖面垂直平分N等份，如圖6所示，任取其中相鄰的兩層，記為第i層和第i+1層，其中第i層更靠近對稱軸，然后在坐標紙上相應標記出該層的芯、殼層熔體前沿位置，當所有層都標記完畢時，用曲線連接即可得到此芯層注射時間下的芯、殼層熔體前沿位置，其中Y軸方向只考慮了一半。

通過數值模擬，還可獲得芯、殼層熔體在充填過程中壓力場、溫度場的分布情況。圖7和8所示分別為芯、殼層喂料在注射時間為0.902 s時的壓力場和溫度場分布圖。由圖7和8可以看出，在充填過程中，殼層熔體在澆口處的壓力和溫度比較大，隨著穿伸位置的增大，殼層熔體的壓力和溫度隨之減小，這與共注射矩形模腔的生產實踐相符合。

圖6 共注射樣品平剖面垂直平分圖Fig. 6 Midnormal drawing of athwartship plane of co-injected part

圖7 芯、殼層熔體在時間為0.902 s時的壓力場分布圖Fig. 7 Pressure field of core and skin feed stocks at 0.902 s

圖8 芯、殼層熔體在時間為0.902 s時的溫度場分布圖Fig. 8 Temperature field of core and skin feedstocks at 0.902 s

通過對芯、殼層注射過程芯/殼層界面、殼層前沿追蹤的模擬結果與實驗結果的對比，可以看出，模擬獲得的芯/殼層界面、殼層前沿的結果與實驗結果符合程度相對較好。需要指出的是，利用數值模擬技術除了能追蹤前沿界面以外，還獲得熔體在充填過程中的壓力場、溫度場等物理場的分布情況，但是，這些場的分布在實驗中較難驗證，需要通過間接的方法進行驗證。例如，可以嘗試利用注射坯的裂紋來驗證壓力場的分布情況，利用熔接線來驗證溫度場的分布等。這需要在以后的工作中進行進一步的驗證。

5 結論

1) 在建立聚合物熔體充填三維薄板型腔的數學模型基礎上，對傳統的控制體積法進行改進并追蹤芯、殼層熔體前沿的運動情況，闡述了粉末共注射成形與普通注射成形充填過程數值模擬的不同，基于芯層在厚度方向上的充填時間較短，和殼層預填充量是影響芯層厚度變化的首要因素的實驗結論，提出了在固定的穿伸距離處芯層厚度在未達到最大厚度前，芯層厚度呈線性變化的假設。

2) 在充填過程中利用Malta軟件進行程序開發，實現了共注射過程中芯、殼層前沿界面的追蹤，以及壓力場、溫度場的分布的可視化，模擬結果與實驗結果吻合很好。

REFERENCES

[1]ALCOCK J R, LOGAN P M, STEPHENSON D J. Surface engineering by co-injection moulding[J]. Surface and Coatings Technology, 1998, 105: 65?71.

[2]朱 計, 類彥威, 張 杰. 共注射成型技術及其發展[J]. 工程塑料應用, 2007, 35(4): 31?35. ZHU Ji, LEI Yan-wei, ZHANG Jie. Co-injection moulding and its development[J]. Engineering Plastics Application, 2007, 35(4): 31?35.

[3]ROSOF B H. The metal injection molding process comes of age[J]. Journal of Mater, 1989, 41(8): 13?16.

[4]LI C T, ISAYEV A I. Interface development and encapsulation in simultaneous co-injection molding of disk (Ⅰ): Twodimensional modeling and formulation[J]. Journal of Applied Polymer Science, 2003, 88: 2300?2309.

[5]LI C T, LEE D J, ISAYEV A I. Interface development and encapsulation in simultaneous co injection molding of disk (Ⅱ): Two-dimensional simulation and experiment[J]. Journal of Applied Polymer Science, 2003, 88: 2310?2318.

[6]周國發, 孫 懋, 柳和生, 李紳元, 羅忠民. 順序共注成型的實驗研究[J]. 中國塑料, 2003, 17(4): 53?57. ZHOU Guo-fa, SUN Mao, LIU He-sheng, LI Shen-yuan, LUO Zhong-min. Experimental study of sequential co-injection molding[J]. China Plastics, 2003, 17(4): 53?57.

[7]匡唐清, 柳和生, 周國發, 賴佳美. 共注射成型充填過程的運動界面追蹤[J]. 塑料工業, 2005, 33(3): 32?35. KUANG Tang-qing, LIU He-sheng, ZHOU Guo-fa, LAI Jia-mei. Tracing of moving interface during polymer mold filling in co-injection molding process[J]. China Plastics Industry, 2005, 33(3): 32?35.

[8]QU Xuan-hui. Numerical simulation of feedstock melt filling in a cylindrical cavity with solidification in powder injection molding[J]. Transactions of Nonferrous Metals Society of China, 1998, 8(4): 544?545.

[9]曲選輝, 李益民, 黃伯云. 金屬粉末注射成形技術[J]. 粉末冶金材料科學與工程, 1996, 1(2): 33?37. QU Xuan-hui, LI Yi-min, HUANG Bo-yun. The technology of metal powder injection molding[J]. Materials Science and Engineering of Powder Metallurgy, 1996, 1(2): 33?37.

[10]LI C T, LEE D J, ISAYEV A I. Interface and encapsulation in simultaneous co-injection molding of disk: Two-dimensional simulation and experiment[C]//Annual Conference 2002, Society of Plastics Engineers, 2002, 48: 465?469.

[11]HIEBER C A, SHEN S F. A finite-element/finite-difference simulation of the injection-molding filling process[J]. Journal of Non-Newtonian Fluid Mechanica, 1987, 7: 1?32.

[12]KOSZKUL J, NABIALEK J. Viscosity models in simulation of the filling stage of the injection molding process[J]. Journal of Materials Processing Technology, 2004, 157/158: 183?187.

[13]HIRT C W, NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39: 201?225.

[14]RIDER W J, KOTHE D B. Reconstructing volume tracking[J]. Journal of Computational Physics, 1998, 141: 112?152.

[15]劉 全, 柳和生, 水鴻壽, 張曉軼. 數值模擬界面流方法進展[J]. 力學進展, 2002, 32(2): 259?274. LIU Quan, LIU He-sheng, SHUI Hong-shou, ZHANG Xiao-yi. Advances in numerical simulation of interfacial/free-surface flows[J]. Advances in Mechanics, 2002, 32(2): 259?274.

[16]李慶揚. 數值分析[M]. 武漢: 華中理工大學出版社, 1982. LI Qing-yang. Numerical analysis[M]. Wuhan: Huazhong University of Science and Technology Press, 1982.

[17]孫 祥, 徐留美, 吳 青. MATLAB7.0基礎教程[M]. 北京:清華大學出版社, 2005. SUN Xiang, XU Liu-mei, WU Qing. MATLAB7.0 based tutorial[M]. Beijing: Tsinghua University Press, 2005.

(編輯 龍懷中)

Numerical simulation of melt front and filed profile in powder co-injection molding filling process

FENG Juan1, HE Hao2, LI Yi-min2, WANG Guang-yao2
(1. School of Mathematical Science and Computing Technology, Central South University, Changsha 410075, China; 2. State Key Laboratory of Powder Metallurgy, Central South University, Changsha 410083, China)

The powder co-injection molding filling process was analyzed and verified by the theory of numerical analysis and experiments. The thickness of core/shell interface was determined by experiments and numerical fitting. The melt front was traced by a modified control volume. The pressure equation and energy equation were solved by hybrid finite element/finite difference approaches. The simulation program was developed by means of Matlab, as a result the tracking of movement of skin front and core/skin interface and profile of temperature and pressure fields were obtained. Finally, simulation results were compared with the experimental ones. The simulated results agree well with the movement of feed’s forefront. With the filling of the feed, the deviation between simulation results and the test results increases. The reason may be that the contraction of injected parts is ignored in the process of numerical simulation.

powder co-injection molding; numerical simulation; control volume method; filling flow

TF124；O242.1

A

國家高技術研究發展計劃資助項目(2007AA03Z114)；國家自然科學基金資助項目(51075405)；中南大學博士后基金資助項目(7604130008)

2011-08-01；

2011-10-26

何 浩，助理研究員，博士；電話：13973114617；E-mail: he_hao555@yahoo.com.cn

1004-0609(2012)08-2333-07