型砂壓實(shí)過程的非線性有限元數(shù)值模擬

李文珍，吳浚郊，周克林

（1.清華大學(xué)材料學(xué)院，北京 100084；2.上海汽車集團(tuán)，上海 200041）

在濕型砂鑄造生產(chǎn)中，為了保證鑄件的幾何尺寸和力學(xué)性能符合設(shè)計(jì)和使用要求，高質(zhì)量砂型的生產(chǎn)起決定作用。型砂緊實(shí)過程是生產(chǎn)砂型的關(guān)鍵工藝環(huán)節(jié)，利用現(xiàn)代模擬仿真技術(shù)對型砂緊實(shí)過程進(jìn)行模擬，預(yù)測砂型緊實(shí)度分布并進(jìn)而預(yù)測砂型質(zhì)量，是提高造型生產(chǎn)效率獲得高質(zhì)量鑄型的有效手段。但是黏土濕型砂屬于離散材料體系，與一般的連續(xù)材料體系相比具有如下特點(diǎn)：

1）材料非線性即砂型由大量非連續(xù)的外層包覆黏土-水膜的細(xì)小砂粒組成，砂粒間相互接觸；

2）幾何非線性即型砂在緊實(shí)過程中砂型的變形量很大，可達(dá)30%以上；

3）接觸非線性即在緊實(shí)過程中型砂與砂箱和模樣之間存在非線性接觸。因此在型砂緊實(shí)過程中同時(shí)存在著材料非線性、幾何非線性和接觸非線性，這類問題一般用非線性有限元法進(jìn)行求解，如圖1所示。目前所能查到的濕砂造型數(shù)值模擬領(lǐng)域最早的研究是Leone和R.L.Lewis[1]使用有限元法對壓實(shí)造型中型砂密實(shí)過程進(jìn)行的數(shù)值模擬。

1 用ANSYS軟件求解非線性問題方法簡介

ANSYS軟件是一個(gè)商品化的大型通用有限元分析軟件，目前廣泛應(yīng)用于土木工程、水利工程、機(jī)械工程、航空航天和汽車交通等領(lǐng)域。ANSYS具有一些突出的技術(shù)特點(diǎn)：將前處理模型的信息（單元、節(jié)點(diǎn)、材料），邊界條件（載荷、約束）以及后處理信息（求解結(jié)果）集成在統(tǒng)一的數(shù)據(jù)庫中，極大地方便了用戶的使用。ANSYS具有強(qiáng)大的非線性分析功能，可以處理幾何非線性、材料非線性、狀態(tài)非線性及單元非線性問題。ANSYS提供了100余種非線性材料模型如橡膠、泡沫、土壤和巖石等。另外ANSYS軟件能夠靈活地在前處理中進(jìn)行網(wǎng)格劃分，確保了單元的形態(tài)和精度。用ANSYS軟件求解非線性問題需要一系列帶校正的線性近似來實(shí)現(xiàn)，如圖2.

圖1 型砂壓實(shí)過程中非線性問題的求解流程

圖2 純粹增量近似與牛頓－拉普森近似的關(guān)系[2]

ANSYS軟件使用牛頓－拉普森平衡迭代進(jìn)行計(jì)算，它迫使在每一個(gè)載荷增量的末端解達(dá)到平衡收斂（在某個(gè)容限范圍內(nèi)）。圖2b）描述了在單自由度非線性分析中牛頓－拉普森平衡迭代的使用。在每次求解前，NR方法估算出殘差矢量，這個(gè)矢量是回復(fù)力（對應(yīng)于單元應(yīng)力的載荷）和所加載荷的差值。程序然后使用非平衡載荷進(jìn)行線性求解，且核查收斂性。如果不滿足收斂準(zhǔn)則，重新估算非平衡載荷，修改剛度矩陣，獲得新解。持續(xù)這種迭代過程直到問題收斂[2]。

ANSYS程序提供了一系列命令來增強(qiáng)問題的收斂性，如自適應(yīng)下降，線性搜索，自動(dòng)載荷步，及二分等，可被激活來加強(qiáng)問題的收斂性，如果不能得到收斂，那么程序依據(jù)輸入的命令或者繼續(xù)計(jì)算下一個(gè)載荷前或者終止[2]。

2 型砂材料本構(gòu)關(guān)系建立

型砂本構(gòu)關(guān)系即應(yīng)力應(yīng)變關(guān)系的建立是對砂型緊實(shí)過程進(jìn)行有限元分析的基礎(chǔ)。在土力學(xué)中已提出了很多測定和建立砂土本構(gòu)關(guān)系的方法，主要實(shí)驗(yàn)方法為有側(cè)限壓縮，它是目前最常見的用來測量土的壓縮性能的試驗(yàn)方法[3]。

2.1 型砂本構(gòu)關(guān)系測試原理[4]

圖3所示為型砂靜態(tài)本構(gòu)關(guān)系測試原理。具體方法是：將混好的型砂加入試樣筒中，隨后利用液壓式萬能試驗(yàn)機(jī)進(jìn)行慢速加載，通過試樣上的壓板作用在松散的型砂表面。在加載過程中，壓桿上的BLR-1載荷傳感器測量壓板載荷P的變化，位移傳感器測量壓板相應(yīng)的位移值U。這兩路傳感器信號經(jīng)Y6DL-1型應(yīng)變儀放大后，由LZ3-304函數(shù)記錄儀以P為縱坐標(biāo)，U為橫坐標(biāo)輸出P-U曲線。經(jīng)下列轉(zhuǎn)換得到在靜態(tài)載荷作用下，有側(cè)限情況下的型砂應(yīng)力－應(yīng)變（本構(gòu)關(guān)系）曲線。

式中：S－壓板面積；

h－試樣筒的高度。

圖3 測量型砂靜態(tài)本構(gòu)關(guān)系實(shí)驗(yàn)原理圖[4]

2.2 試驗(yàn)結(jié)果及分析

取型砂緊實(shí)率40%，膨潤土含量為10%，初始密度分別為 0.9g/cm3、1.0g/cm3、1.1 g/cm3、和 1.2g/cm3，測得的型砂應(yīng)力應(yīng)變關(guān)系如圖4所示。

從圖4可以看出型砂應(yīng)力應(yīng)變關(guān)系為非線性，呈下凹形狀，這點(diǎn)與大多數(shù)連續(xù)材料的特性不同。未緊實(shí)的松散型砂在整個(gè)壓縮過程中具有非常大的變形量。其原因如下：型砂在開始受壓時(shí)的變形主要表現(xiàn)為砂粒重新排列填補(bǔ)空隙，因此應(yīng)力較??；當(dāng)型砂緊實(shí)到一定階段，砂粒之間的接觸變得更加緊密，砂粒運(yùn)動(dòng)不僅需要克服黏土膜之間的黏滯力，還要克服砂粒之間的摩擦力；在緊實(shí)的最后階段，砂粒之間石英內(nèi)核發(fā)生相互作用的正應(yīng)力，因?yàn)槭㈩w粒非常堅(jiān)硬，所以緊實(shí)后期應(yīng)力急劇增大。

圖4 型砂本構(gòu)關(guān)系曲線

3 型砂壓實(shí)過程的數(shù)值模擬

在獲得型砂側(cè)限條件的本構(gòu)關(guān)系后，本文在有限元軟件ANSYS的平臺上，采用APDL進(jìn)行二次開發(fā)，編寫了相應(yīng)的計(jì)算程序并耦合到現(xiàn)有ANSYS軟件中，對型砂壓實(shí)過程進(jìn)行了數(shù)值模擬。

3.1 材料模型的建立

從上述型砂本構(gòu)關(guān)系試驗(yàn)結(jié)果可以看出，型砂的本構(gòu)關(guān)系是一種非線性應(yīng)力應(yīng)變關(guān)系，其應(yīng)力對應(yīng)變的一階導(dǎo)數(shù)即彈性模量在變形過程中隨著應(yīng)變的增加而遞增，對于這樣的模型，通常的商業(yè)分析軟件無法直接予以處理，必須對實(shí)際本構(gòu)關(guān)系進(jìn)行一定的簡化處理，并采用APDL語言進(jìn)行編程來完成計(jì)算，以實(shí)現(xiàn)對這種復(fù)雜應(yīng)力應(yīng)變關(guān)系的求解。

本文根據(jù)型砂材料的實(shí)際特點(diǎn)，對型砂材料模型作了如下假定：

1）型砂是均勻、連續(xù)、各向同性的介質(zhì)；

2）由于型砂緊實(shí)過程相對較快，因此，計(jì)算過程中不考慮回彈；

3）應(yīng)力加載過程為從零單調(diào)增加到最大，卸載過程為無位移（很小）回彈。根據(jù)這些假定，對試驗(yàn)測得的型砂非線性本構(gòu)關(guān)系（圖4）采用多線性本構(gòu)關(guān)系進(jìn)行簡化，如圖5所示。在曲線上，每兩點(diǎn)之間為一直線段，彈性模量ΔE=Δσ/Δε為常數(shù)，型砂表現(xiàn)為線彈性材料。本文選取了初始密度為0.9 g/cm3，膨潤土含量為8%，緊實(shí)率為40%的型砂作為計(jì)算對象，線性段階躍點(diǎn)的應(yīng)力應(yīng)變值如表1所示。

圖5 型砂的多線性本構(gòu)關(guān)系

表1 初始密度為0.9 g/cm3的型砂各線性段應(yīng)力應(yīng)變

3.2 型砂幾何模型

計(jì)算中，型砂幾何模型如圖6所示。設(shè)砂箱和余砂框的總高度為A，壓實(shí)后型砂頂面距離分型面的高度為B，未壓實(shí)前成型壓頭距離模樣頂面的距離為a，壓實(shí)后型砂定面距離模樣頂面的距離為b.為了使型砂壓實(shí)后型砂更加均勻，則成型壓頭尺寸滿足以下等式條件：

圖6 幾何模型尺寸

3.3 加載過程和邊界條件的處理

在型砂壓實(shí)過程中，型砂與砂箱之間存在很大的摩擦力，這對型砂壓實(shí)過程有著不可忽略的影響，因此在計(jì)算過程中，在砂箱和型砂之間采用接觸條件，并參考文獻(xiàn)[3]中的數(shù)據(jù)，取平均摩擦系數(shù)為0.5.

由于型砂和砂箱之間存在摩擦力，型砂采用單面壓實(shí)后，緊實(shí)度分布情況一般為上緊下松，這對于處在分型面處的型腔在充填金屬液體時(shí)是極為不利的。因此，本文分別計(jì)算了兩種加載情況：

1）采用上加壓單面壓實(shí)；

2）采用上下加壓雙面壓實(shí)，且按照上面加載量為下面加載量的3倍進(jìn)行加載，并比較了兩種加載情況的計(jì)算結(jié)果。本文計(jì)算了緊實(shí)率為40%，膨潤土含量為8%，初始密度為0.9 g/cm3的型砂。采用位移加載，且按照壓縮40%的位移量加載；由于型砂的本構(gòu)關(guān)系特殊，將型砂分為幾個(gè)線性段進(jìn)行近似模擬，采用分步加載。在加載過程中，根據(jù)每個(gè)線性段的彈性應(yīng)變量進(jìn)行分步位移加載。當(dāng)前一個(gè)分步加載求解完畢后，提取計(jì)算結(jié)果并改變模型中各個(gè)節(jié)點(diǎn)的位置，重新劃分網(wǎng)格，然后改變彈性模量，在進(jìn)行下一個(gè)分步加載求解時(shí)，將前一個(gè)分步計(jì)算的結(jié)果以預(yù)應(yīng)力載荷加在重新劃分網(wǎng)格的模型上，重復(fù)該計(jì)算過程，直到計(jì)算求解到最后一個(gè)線性段結(jié)束為止。

單元采用8節(jié)點(diǎn)平面等參元，初始網(wǎng)格劃分結(jié)果如圖7所示。由于型砂在壓實(shí)過程中具有流動(dòng)性特征，若以彈性材料進(jìn)行求解，用ANSYS求解后會(huì)在成型壓頭附近出現(xiàn)一個(gè)缺口，為了解決此問題，將成型壓頭下面的節(jié)點(diǎn)在Y方向的位移進(jìn)行耦合求解。此外，成型壓頭兩邊的型砂變形量大，單元若太小，在計(jì)算過程中容易發(fā)生不收斂而導(dǎo)致計(jì)算中止，因此特意將成型壓頭兩邊的型砂單元?jiǎng)澐值妮^大，由于我們感興趣的區(qū)域?yàn)榉中兔娓浇膽?yīng)力分布，因此在型砂頂面將單元?jiǎng)澊髮Y(jié)果不會(huì)有很大的影響[5]。

圖7 模型網(wǎng)格劃分結(jié)果

4 數(shù)值模擬結(jié)果及分析

4.1 單面壓實(shí)和雙面壓實(shí)數(shù)值模擬結(jié)果的比較

數(shù)值模擬結(jié)果如圖8a）和圖b）分別是單面壓實(shí)和雙面壓實(shí)后的計(jì)算結(jié)果。為了比較兩種加載情況的計(jì)算結(jié)果，每一種加載情況都提取了兩條線上的應(yīng)力分布（如圖6）：一條靠近砂箱邊緣，一條則為模型的對稱中心線，兩條線方向均由下向上。兩條線上的應(yīng)力值如表2和表3.將表中的應(yīng)力制成圖，如圖9所示。

表2 沿砂箱邊緣高度方向的應(yīng)力分布各點(diǎn)的值

表3 沿砂箱中心線高度方向的應(yīng)力分布各點(diǎn)的值

圖8 數(shù)值模擬應(yīng)力場的分布

圖9 沿砂箱邊緣和模型中心兩條線上的應(yīng)力分布

從圖8和圖9的模擬計(jì)算結(jié)果可以看出，由于型砂砂粒之間摩擦力及型砂與砂箱壁和模型表面之間摩擦力的影響，在靠近砂箱的附近應(yīng)力沿砂箱高度方向增加（如圖9），而在砂箱中心線附近應(yīng)力則是沿砂箱高度降低（如圖9），模擬結(jié)果符合實(shí)際緊實(shí)情況。同時(shí)值得注意的是，在模樣拐角處和型砂狹窄區(qū)域入口處，應(yīng)力場為顏色較深，表明在該處出現(xiàn)一個(gè)低應(yīng)力區(qū)域，應(yīng)力低于其上方區(qū)域的應(yīng)力，模擬出了文獻(xiàn)[4]試驗(yàn)中的“搭橋”現(xiàn)象。

此外從圖9還可以看出，在高度較低的地方，雙面壓實(shí)砂箱邊緣應(yīng)力曲線的值要比單面壓實(shí)砂箱邊緣應(yīng)力曲線的值要大；而在高度較高的地方，雙面壓實(shí)砂箱邊緣應(yīng)力曲線的值反而要比單面壓實(shí)砂箱邊緣應(yīng)力曲線的值要小。由此說明，雙面壓實(shí)后，可以提高分型面處型砂的緊實(shí)度，且整個(gè)型砂應(yīng)力分布要更加均勻，符合雙面壓實(shí)要優(yōu)于單面壓實(shí)的規(guī)律。

4.2 數(shù)值模擬結(jié)果和實(shí)驗(yàn)測量值的比較

由于實(shí)驗(yàn)條件所限，本文參考了文獻(xiàn)[4]中的試驗(yàn)數(shù)據(jù)（如圖10所示，圖中應(yīng)力單位為MPa），并將模擬結(jié)果與之作了數(shù)量級上的比較以驗(yàn)證模擬結(jié)果。

圖10 參考文獻(xiàn)[4]的實(shí)測應(yīng)力分布圖

文獻(xiàn)[4]中的試驗(yàn)條件如下：砂箱平面尺寸為390 mm×350 mm，型砂初始密度為0.9 g/cm3，緊實(shí)率40%，膨潤土含量為8%。這些試驗(yàn)條件除了砂箱尺寸和本文計(jì)算模型有一定的差別（本文模型尺寸為400 mm×300 mm），其他條件均相等，因而有一定的參考意義。

圖9中計(jì)算出的沿砂箱高度方向的應(yīng)力曲線和圖10b）實(shí)測應(yīng)力曲線非常相似，趨勢相同，應(yīng)力均隨著高度的增加而增加，且都在砂箱和模樣之間存在一塊應(yīng)力比較小的區(qū)域，也即“拱橋效應(yīng)”。在應(yīng)力值的數(shù)量級上，在“拱橋效應(yīng)”區(qū)域，計(jì)算出的應(yīng)力值大約在0.34 MPa～0.9 MPa之間，實(shí)測的應(yīng)力值在0.36 MPa～0.7MPa之間；而在高出模樣的區(qū)域，計(jì)算出的應(yīng)力值大約在0.9 MPa～1.9 MPa之間，而試驗(yàn)測得的值也在0.8 MPa以上。此外，從計(jì)算結(jié)果和試驗(yàn)結(jié)果可以看出，在模樣上方都有一塊比較緊實(shí)的區(qū)域。上述結(jié)果表明數(shù)值模擬的結(jié)果和試驗(yàn)結(jié)果比較吻合。

5 結(jié)論

利用商品化有限元分析軟件ANSYS，通過APDL二次開發(fā)編程對黏土砂造型的壓實(shí)過程進(jìn)行了數(shù)值模擬。通過對兩種加載方法的計(jì)算結(jié)果比較以及計(jì)算結(jié)果與相應(yīng)的試驗(yàn)數(shù)據(jù)比較，結(jié)果表明ANSYS可以很好地解決型砂緊實(shí)這樣的非線性問題。濕型砂壓實(shí)過程的數(shù)值模擬結(jié)果可以較好地解釋如下現(xiàn)象：

1）“拱橋效應(yīng)”現(xiàn)象；

2）“雙面加壓”要優(yōu)于“單面加壓”，從而可以提高型砂分型面的強(qiáng)度；

3）在砂箱邊緣應(yīng)力分布從下至上逐漸變大。

通過ANSYS軟件對型砂緊實(shí)過程進(jìn)行數(shù)值模擬取得了較好的效果。但是，由于實(shí)驗(yàn)數(shù)據(jù)來自參考文獻(xiàn)，因而以后需要通過進(jìn)一步的實(shí)驗(yàn)來驗(yàn)證模擬結(jié)果的精度。此外，本文模擬仍局限在二維模型上，對于三維模型模擬有待進(jìn)一步的開發(fā)。

［1］Leone M，Lewis R L.Numerical Modeling of Green Sand During Compaction［J］.AFSTransactions，1988，136：763-774.

［2］王國強(qiáng).實(shí)用工程數(shù)值模擬技術(shù)及其在ANSYS上的實(shí)踐［M］.西安：西北工業(yè)大學(xué)出版社，1999.

［3］張翼飛.射壓造型的數(shù)值模擬和試驗(yàn)研究［D］.北京：清華大學(xué)機(jī)械工程系.2001.

［4］謝濱.壓實(shí)造型和空氣沖擊造型砂型應(yīng)力場研究［D］.北京：清華大學(xué)機(jī)械工程系，1993.

［5］Li Wenzhen，Wu Junjiao.Numerical Simulation of Compacting Process of Green Sand Molding Based on Sand Filling［J］，Materials Science Forum，2007,561-565:1879-1882.