夏伯雄，杜 軍，汪明玲，葛曉健，張 玲

（武漢軟件工程職業(yè)學院機械制造工程系，湖北 武漢 430205）

熔模鑄造對于生產(chǎn)凈形零件是一個重要的生產(chǎn)過程。這個過程基本上是在一個預熱陶瓷沙箱內(nèi)凝固的。熱分析是此系統(tǒng)設(shè)計的關(guān)鍵領(lǐng)域[1]。雖然凝固時間的估計是車間的第一手需求，然而當微觀結(jié)構(gòu)與散熱路徑密切相關(guān)時，熱隨時間的變化決定著產(chǎn)品的質(zhì)量。熱從熔模鑄造鑄型的轉(zhuǎn)移受輻射和對流的雙重控制。研究人員試圖提出一次性模型來預測凝固時間，這是以簡單的假設(shè)為基礎(chǔ)的，如恒壁溫。隨著數(shù)字化計算機的問世和商業(yè)工程軟件包的激增，在設(shè)計優(yōu)化和原型開發(fā)的過程中，該行業(yè)已經(jīng)成熟到能足以使用這樣的代碼。值得注意的是今天的模擬已達到了成熟的水平，造型既進行到了宏觀尺寸，又擴展到了更精細的枝狀水平[2]。這種戰(zhàn)略尤其重要的是消除了鑄件的缺陷。

從文獻調(diào)查來看，全面計算流體動力學（CFD）和集中參數(shù)都曾經(jīng)嘗試過用于分析熔模鑄造。然而，這兩種方法從未比較過。設(shè)計師第一手的預測總是運用集中計算，這樣可能導致誤差。因此，在本文中，將全相能量數(shù)值計算法與集中參數(shù)公式法的求值進行了比較。為此，以一個圓柱形幾何圖形為對象考慮分析。鑄型假設(shè)是充滿了處于908 K時的A356合金，即20℃的過熱。壁上的熱傳導既有對流又有輻射。對流換熱系數(shù)取固定值5 W/m2K，這是對靜態(tài)空氣的合理估計。輻射部分的計算不妨以15元左右的人民幣價值來描述。該材料性能列于表1。模具填充物的熱效應(yīng)被忽略，換句話說，已經(jīng)在最初的過熱質(zhì)量要求中予以了考慮。

表1 材料特性

1 方法介紹

1.1 數(shù)值模擬方法

凝固過程由能量方程控制：

在這里，ρ為密度，H為焓，T是溫度，t表示時間。等式右邊的s項是一個原始項，當沒有內(nèi)熱產(chǎn)生時其值為零。速度部分的影響可被忽略，當尺寸是幾cm的中型鑄件時，速度沒有影響到鑄件凝固前沿的熱放出。

液體部分糊狀區(qū)與其溫度的關(guān)系如下：

對于 T＞Tliq，β=1；而當 T＜Tsol時，β=0

Tliq為液相線溫度；Tsol為固相線溫度。

液體部分非線性方程也有報道，公式固然接近于物理學。然而，Bakhtiyarov等[3]表明，對于不同類型的方程，如Sceil模型，鋁合金等不會有太多的液體方程的變化，建議在工程實踐中可直接采用線性規(guī)律。

能量方程是一個線性方程，當適當?shù)倪吔鐥l件已知時可進行數(shù)值求解。邊界條件可能是Dirrichlet型或Neumann型或它們的組合。在Dirrichlet型中，可提供溫度值。在Neumann型中，可使用熱通量值。絕熱，即絕緣邊界是其中特殊的一類，熱通量為零。

焓的定格法是建模凝固非常流行的一個方法，這要歸因于Voller和Prakash[4]。這種方法不是將液體和固體區(qū)域進行分開處理，而是將液體和固體區(qū)域當成一個共同的網(wǎng)格用于整個區(qū)域。作為焓方程的部分解法是跟蹤凝固前沿。糊狀區(qū)的液體部分可從方程（2）描述的線性關(guān)系中得到。使用上述方程是要建立一個模型來運用FLUENT商業(yè)CFD代碼[5]。流體體積的方法已被用于捕捉邊界[6]。

邊界條件是：

在坐標軸上：偏導數(shù)趨進于零，即 U/ri和T/ri=0，其中ri指坐標，U和T分別為速度和焓場。

在鑄型壁上：傳導熱系數(shù)值h是已知的。鑄型壁的傳熱率由q=-h（Tw-Tamb）給出，其中T為鑄型壁溫度；Tamb為環(huán)境溫度。

1.2 集中參數(shù)分析法

Upadhya等人[7]提出的集中分析法，可以用以下方式進行。鑄造的熱傳導率由下式給出：

其中h是傳熱系數(shù)，A為表面積，Ts和T∞分別是外殼和周圍環(huán)境的溫度。

如果L*是金屬溶融的有效熱，總凝固過程的傳熱時間為tsol，則總傳導熱量為：

其中，ρ為密度，V為體積。

其中，L是潛熱，ΔTsuper是過熱度。

凝固所需時間：

我們將物體的模數(shù)定義為體積與面積的比率，即M=V/A。模數(shù)值低就表示物體相對體積具有較高的表面積。因此，模數(shù)接近零，標志著物體具有很高的表面面積。

現(xiàn)在可以看出：

因此，我們觀察到凝固時間與模數(shù)成正比。

熱傳導包括輻射和對流。熱傳導系數(shù)的輻射部分可從Stephan-Boltzmann輻射方程獲得，給出如下：

對流部分可以用不同的方式加以分析。一個簡單的方法是建議取 hconv=c（Ts-T∞）1/3，其中 c 是一個依表面幾何形狀而定的常數(shù)，hrad為輻身傳熱系數(shù)；hconv為對流傳熱系數(shù)。另一種方法是，以表面的方向性為基礎(chǔ)來選擇自然對流的相關(guān)性，見參考文獻[8]。然而，對于目前的計算，我們使用了以前討論過的hconv固定值。

2 分析與討論

模擬是在模數(shù)M的幾個不同層次上進行的。包括使用了以前所提到的約10000個節(jié)點的軸對稱方程。通過比較三種格層次的結(jié)果來確保格的獨立性。格獨立凝固時間的推導值與目前網(wǎng)格配置的計算值相差僅約0.5%。

控制傳導方程可利用已經(jīng)提及的商業(yè)CFD求解器FLUENT進行數(shù)值求解。方程式基本上整理成以下形式，以利于用一種標準形式進行方程的離散化：

按照SIMPLE算法運用以壓力為基礎(chǔ)的有限體積法使方程離散化[9]。其中φ表示任何守恒的變量，S是原始邊界。在計算流體力學（CFD）過程中，上述類型的微分方程可轉(zhuǎn)換成一個由數(shù)字計算機解決的代數(shù)方程。對于控制方程的解法，規(guī)定了一個帶有1×10-7偏差的收斂判斷準則。作為離散處理，使用了具有較高準確性級別輸出的二階展開方法。詳細數(shù)值實施過程可參見[5]，對模具進行注射成形工藝的填充工作過程可參見[10]。

通過監(jiān)測最內(nèi)側(cè)的液體部分的凝固端點，可估算出凝固時間。

如果已知物體的模數(shù)，運用方程（2）也可估計出凝固時間。在這里應(yīng)當指出的是，Upadhya等沒有以物體對象作為整體來確定模數(shù)[7]，而是將物體對象離散成若干個小體積，顯然這樣做會更好地預測凝固時間。在[7]中提到的這個方案稱為節(jié)模數(shù)方法，并發(fā)現(xiàn)這種方法在精度和計算資源方面非常具有競爭性[11]。

圖1顯示了模數(shù)為0.0125的物體對象的溫度隨時間變化的曲線。這個M值對應(yīng)一個長度為5.0cm，直徑10.0cm的圓柱體。從圖1可觀察到，最初有一個平緩區(qū)，這說明過熱材料的冷卻約21 s。一旦凝固開始，潛熱釋放到環(huán)境。液體部分隨著時間逐漸減少，直到凝固停止在固相線的溫度。

現(xiàn)在定義一個術(shù)語Tr，Tr是通過全相的分析獲得的凝固時間與集中參數(shù)公式法預測的凝固時間的比率。Tr將隨著材料內(nèi)較大梯度的大體積的模數(shù)增加而變化。

這可從圖2中觀察到，時間比Tr隨著模數(shù)增加的二次方而增加。Tr的表達式為：

從以上函數(shù)中可以看出，當模數(shù)為零時，Tr=2.36。極限值零模數(shù)對應(yīng)于非常高的表面積情形。零模數(shù)的極限情況是研究人員特別感興趣的，因為它涉及到的能源的最高可能傳導區(qū)域。

圖1 液體部分隨時間的變化（M=0.0125）

圖2 模數(shù)對凝固時間的影響

3 總 結(jié)

本文表明，可以從能量方程的完整解法得到的凝固時間和集中分析法得到的凝固時間之間建立的一個函數(shù)關(guān)系。這樣的函數(shù)關(guān)系一旦建立，可以通過計算對一定模數(shù)物體的凝固時間做出快速的預測。分析還表明，對于模數(shù)接近零值的高表面體積比，實際凝固時間大約為集中分析法預測時間的2.36倍。一般認為是低模數(shù)高價值的小件物品如黃金飾品往往走熔模鑄造的生產(chǎn)路線。

［1］G.H.Geiger.and D.R.Poirier.Transport Phenomena in Metallurgy［J］.Addesson Wesley，1980.

［2］Hamilton RW，See D，Butler S et al.Multiscalemodeling for the prediction of casting defects in investment castaluminum alloys［J］.Mater Sci Eng，2003，A 343：290-300.

［3］Bakhtiyarov SI，Overfet RA，Teodorescu SG.Fractionalsolid measurements of solidifying melt［J］.J Fluids Eng 2004，126:193-197.

［4］Voller VR，Prakash C.A fixed-grid numerical modelingmethodology for convection-diffusion mushy region phase-change problems［J］.Int J Heat Mass Transfer 1987，30：1709-1720.

［5］Fluent Inc.FLUENT Users Manual[M].Lebanon，2004.

［6］Hirt W，Nichols BD.method for thedynamics of free boundaries［J］.J Comput Phys，1981，39：202-225.

［7］Upadhya GK，Das S，Chandra U et al.Modeling theinvestment casting process：a novel approach for view factorcalculations and defect prediction［J］.Appl Math Model，1995，19：354-362.

［8］Bejan A.Convective heat transfer[M].Wiley,New Jersey，2004.

［9］Patankar SV.Numerical heat transfer and fluid flow［J］.Hemisphere Publication，Washington，1980.

［10］Samanta SK，Chattopadhyay H，Godkhindi MM et al.A numerical study of solidification inpowder injection molding process［J］.Int J Heat Mass Transfer，2008，5：672-682.

［11］Kotschi R，Plutshak L.An easy and inexpensive techniqueto study solidification of castings in three dimensions ［J］.AFS Trans，1981，89：601-608.