型砂壓實過程的非線性有限元數值模擬

李文珍，吳浚郊，周克林

（1.清華大學材料學院，北京 100084；2.上海汽車集團，上海 200041）

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

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

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

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

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

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

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

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

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

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

2 型砂材料本構關系建立

型砂本構關系即應力應變關系的建立是對砂型緊實過程進行有限元分析的基礎。在土力學中已提出了很多測定和建立砂土本構關系的方法，主要實驗方法為有側限壓縮，它是目前最常見的用來測量土的壓縮性能的試驗方法[3]。

2.1 型砂本構關系測試原理[4]

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

式中：S－壓板面積；

h－試樣筒的高度。

圖3 測量型砂靜態本構關系實驗原理圖[4]

2.2 試驗結果及分析

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

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

圖4 型砂本構關系曲線

3 型砂壓實過程的數值模擬

在獲得型砂側限條件的本構關系后，本文在有限元軟件ANSYS的平臺上，采用APDL進行二次開發，編寫了相應的計算程序并耦合到現有ANSYS軟件中，對型砂壓實過程進行了數值模擬。

3.1 材料模型的建立

從上述型砂本構關系試驗結果可以看出，型砂的本構關系是一種非線性應力應變關系，其應力對應變的一階導數即彈性模量在變形過程中隨著應變的增加而遞增，對于這樣的模型，通常的商業分析軟件無法直接予以處理，必須對實際本構關系進行一定的簡化處理，并采用APDL語言進行編程來完成計算，以實現對這種復雜應力應變關系的求解。

本文根據型砂材料的實際特點，對型砂材料模型作了如下假定：

1）型砂是均勻、連續、各向同性的介質；

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

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

圖5 型砂的多線性本構關系

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

3.2 型砂幾何模型

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

圖6 幾何模型尺寸

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

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

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

1）采用上加壓單面壓實；

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

單元采用8節點平面等參元，初始網格劃分結果如圖7所示。由于型砂在壓實過程中具有流動性特征，若以彈性材料進行求解，用ANSYS求解后會在成型壓頭附近出現一個缺口，為了解決此問題，將成型壓頭下面的節點在Y方向的位移進行耦合求解。此外，成型壓頭兩邊的型砂變形量大，單元若太小，在計算過程中容易發生不收斂而導致計算中止，因此特意將成型壓頭兩邊的型砂單元劃分的較大，由于我們感興趣的區域為分型面附近的應力分布，因此在型砂頂面將單元劃大對結果不會有很大的影響[5]。

圖7 模型網格劃分結果

4 數值模擬結果及分析

4.1 單面壓實和雙面壓實數值模擬結果的比較

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

表2 沿砂箱邊緣高度方向的應力分布各點的值

表3 沿砂箱中心線高度方向的應力分布各點的值

圖8 數值模擬應力場的分布

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

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

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

4.2 數值模擬結果和實驗測量值的比較

由于實驗條件所限，本文參考了文獻[4]中的試驗數據（如圖10所示，圖中應力單位為MPa），并將模擬結果與之作了數量級上的比較以驗證模擬結果。

圖10 參考文獻[4]的實測應力分布圖

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

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

5 結論

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

1）“拱橋效應”現象；

2）“雙面加壓”要優于“單面加壓”，從而可以提高型砂分型面的強度；

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

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

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

［2］王國強.實用工程數值模擬技術及其在ANSYS上的實踐［M］.西安：西北工業大學出版社，1999.

［3］張翼飛.射壓造型的數值模擬和試驗研究［D］.北京：清華大學機械工程系.2001.

［4］謝濱.壓實造型和空氣沖擊造型砂型應力場研究［D］.北京：清華大學機械工程系，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.