基于SPH方法的鑄造充型過(guò)程數(shù)值研究

周學(xué)君，陳　丁

(1. 黃岡師范學(xué)院 數(shù)理學(xué)院，湖北 黃州 438000；2. 河海大學(xué) 力學(xué)與材料學(xué)院，江蘇 南京 210098)

?

基于SPH方法的鑄造充型過(guò)程數(shù)值研究

周學(xué)君1,2，陳丁2

(1. 黃岡師范學(xué)院 數(shù)理學(xué)院，湖北 黃州 438000；2. 河海大學(xué) 力學(xué)與材料學(xué)院，江蘇 南京 210098)

數(shù)值方法模擬鑄造充型過(guò)程一直是鑄造學(xué)科的研究熱點(diǎn)，考慮利用光滑粒子流體動(dòng)力學(xué)(Smoothed Particle Hydrodynamics，SPH)方法研究鑄造充型過(guò)程。SPH方法是一種成熟的無(wú)網(wǎng)格粒子數(shù)值計(jì)算方法，特別適合模擬大變形問(wèn)題。在建立充型過(guò)程的計(jì)算模型的基礎(chǔ)之上，通過(guò)弓形件和環(huán)形件兩個(gè)充型的水模擬算例，并與實(shí)驗(yàn)和文獻(xiàn)中結(jié)果比較，展示了本文SPH算法在處理鑄造充型問(wèn)題中的優(yōu)勢(shì)。

光滑粒子流體動(dòng)力學(xué)；鑄造充型；無(wú)網(wǎng)格；數(shù)值模擬

鑄造是基礎(chǔ)工業(yè)中重要的工藝方法，它利用液態(tài)金屬的流動(dòng)性，在特定形狀的鑄模的約束下，把金屬材料制造成各種用途的鑄件。充型過(guò)程是鑄件成型的第一個(gè)階段，過(guò)程中伴隨著一系列的復(fù)雜的物理和化學(xué)變化，易發(fā)生鑄件缺陷對(duì)鑄件的產(chǎn)品質(zhì)量產(chǎn)生重要影響；特別是在充型過(guò)程中，型腔始終是密閉狀態(tài)，很難直接觀測(cè)金屬流在型腔內(nèi)的運(yùn)動(dòng)規(guī)律。數(shù)值模擬為鑄造充型過(guò)程的研究提供了一種方便廉價(jià)的思路，不僅可以幫助工程師直接觀察到金屬流在型腔內(nèi)的流動(dòng)狀態(tài)，而且再現(xiàn)性好，縮短產(chǎn)品開(kāi)發(fā)周期和降低成本。

鑄造充型過(guò)程的模擬是通過(guò)求解一系列的守恒方程，得出腔內(nèi)流體的運(yùn)動(dòng)規(guī)律；主要困難在于守恒方程中壓力的計(jì)算和自由表面的位置、形狀的判斷。自上世紀(jì)60年代開(kāi)始，人們開(kāi)始利用數(shù)值方法來(lái)模擬鑄造充型過(guò)程，大多是基于網(wǎng)格類方法。針對(duì)守恒方程中壓力場(chǎng)的計(jì)算，出現(xiàn)了SIMPLE方法(1972)[1]和PISO 方法(1986)[2]等；針對(duì)自由表面的追蹤方面，提出了MAC方法(1965)[3]，VOF方法(1981)[4]，PIC方法(1987)[5]和Level Set方法(1988)[6]等。以這些方法為基礎(chǔ)，人們也開(kāi)發(fā)一些鑄造充型過(guò)程的仿真和分析的商業(yè)軟件，如Magmasoft，F(xiàn)low3d，A DSTEFAN，A nycasting，華鑄 CAE，F(xiàn)T-STAR等。眾所周知，基于網(wǎng)格的數(shù)值方法在處理極度大變形問(wèn)題時(shí)，會(huì)出現(xiàn)因網(wǎng)格扭曲或者纏繞而引發(fā)的計(jì)算精度降低和失敗等情況。無(wú)網(wǎng)格方法擺脫了網(wǎng)格的限制，在處理大變形、自由表面、運(yùn)動(dòng)交界面等網(wǎng)格類算法難以解決的問(wèn)題時(shí)，展現(xiàn)了獨(dú)特的優(yōu)勢(shì)。

光滑粒子流體動(dòng)力學(xué)是目前發(fā)展最成熟的無(wú)網(wǎng)格方法之一，是由Lucy(1977)[7]、Gingold和Monaghan(1977)[8]并行提出的一種Lagrangian形式的無(wú)網(wǎng)格粒子計(jì)算方法。SPH方法最初用于解決天體物理學(xué)相關(guān)問(wèn)題，而后被廣泛應(yīng)用于流體力學(xué)和固體力學(xué)的諸多領(lǐng)域。在SPH方法中，將求解域離散成一系列粒子，粒子間通過(guò)核函數(shù)來(lái)建立相互作用，對(duì)于自由表面的追蹤能夠通過(guò)粒子的位置自動(dòng)確定，從而方便地獲取充型過(guò)程中流體的流動(dòng)軌跡，在鑄造充型過(guò)程模擬中具有很大優(yōu)勢(shì)。曹文炅在其博士學(xué)位論文[9]較系統(tǒng)地闡述了SPH方法在鑄造充型中的應(yīng)用；強(qiáng)洪夫等[10]采用罰函數(shù)來(lái)處理流體與壁面之間的作用，對(duì)水模擬充型過(guò)程進(jìn)行SPH數(shù)值分析。

本文以鑄造充型水模擬模型為算例，并與相關(guān)文獻(xiàn)和實(shí)驗(yàn)結(jié)果比較，研究該SPH算法在充型仿真中的可靠性。

1　SPH方法的近似和控制方程的離散

1.1SPH方法的近似

在SPH方法中，將連續(xù)體離散成有限個(gè)粒子，通過(guò)對(duì)某個(gè)粒子的支持域內(nèi)其他粒子的加權(quán)求和獲得該點(diǎn)的數(shù)值解。SPH方法近似分兩步進(jìn)行，第一步是積分近似；第二步是粒子近似。

場(chǎng)函數(shù)f(x)在問(wèn)題域Ω內(nèi)的積分近似表示式為

(1)

式(1)中x，x′為包含在問(wèn)題域Ω內(nèi)點(diǎn)的坐標(biāo)向量；dx′表示x處無(wú)窮小體元；h是光滑長(zhǎng)度；W為光滑函數(shù)(核函數(shù))，這里采用三次樣條函數(shù)[11]為光滑函數(shù)。

積分表示式(1)可寫(xiě)成離散化的粒子近似式：

(2)

這里mj，ρj分別表示粒子j的質(zhì)量和密度，j=1,2,…，N,N為在x處粒子的支持域內(nèi)的粒子總數(shù)。

于是，在粒子i處的場(chǎng)函數(shù)的粒子近似式可寫(xiě)成

(3)

其中Wij=W(xi-xj,h)，且f(x)的導(dǎo)數(shù)的粒子近似式為

(4)

1.2控制方程的SPH離散

Lagrangian描述下的流體動(dòng)力學(xué)控制方程包括連續(xù)性方程、動(dòng)量方程和運(yùn)動(dòng)方程：

(5)

(6)

(7)

其中α表示Cartesian坐標(biāo)分量x,y和z，重復(fù)下標(biāo)代表Einstein求和法則；ρ，v，p分別為密度、速度和壓力；g是體力；D/Dt代表物質(zhì)導(dǎo)數(shù)。另外，本文研究水模擬鑄造充型問(wèn)題，故壓力p是由弱可壓縮流體的狀態(tài)方程[12]得出：

(8)

式中ρ0取流體初始密度，參數(shù)γ取7，且參數(shù)B由

(9)

確定，其中vmax為流體的最大速度。

將上述流體控制方程應(yīng)用SPH核近似和粒子近似后，離散化的SPH公式為

(10)

(11)

(12)

這里式(11)和(12)中包含對(duì)稱結(jié)構(gòu)，可以減少由于粒子不一致問(wèn)題產(chǎn)生的誤差[13]，而文獻(xiàn)[9]中采用移動(dòng)最小二乘法來(lái)處理此問(wèn)題；另外，式(11)中Πij為人工粘度項(xiàng)，它的作用是消除SPH方法在模擬流體動(dòng)力學(xué)問(wèn)題時(shí)產(chǎn)生的數(shù)值不穩(wěn)定性，本文采用的人工粘度[14]為

(13)

其中

(14)

式中αΠ和βΠ是常數(shù)，取1.0左右； c為聲速，為保證流體的弱可壓縮性，本文取c=10vmax；vij=vi-vj，xij=xi-xj。

1.3邊界條件

鑄造充型模擬的邊界為固定邊界，對(duì)于固定邊界條件的施加，一般有三種方法：邊界排斥力法、鏡像虛粒子法和耦合邊界法。針對(duì)充型模型中所用到的鑄模幾何形狀復(fù)雜，采用邊界排斥力法來(lái)處理固定邊界。

Monaghan(1994)[12]在固定邊界上布置一組虛粒子用于對(duì)鄰近邊界的實(shí)粒子產(chǎn)生排斥力，從而阻止鄰近邊界的這些實(shí)粒子非物理穿透邊界。排斥力的表達(dá)式與計(jì)算分子力所使用的Lennard-Jones方程相似，即當(dāng)邊界虛粒子位于鄰近邊界的實(shí)粒子的影響域內(nèi)，則會(huì)在沿著兩粒子的中心線處對(duì)實(shí)粒子產(chǎn)生一個(gè)作用力：

(15)

其中，參數(shù)n1和n2滿足n1>n2，且取值一般分別為12和4；系數(shù)D的取值一般由具體問(wèn)題而定，一般取與速度最大值的平方相等的量級(jí)；截止半徑r0的取值一般等于粒子的初始間距。

2　鑄造充型模擬

本節(jié)所選取的兩個(gè)算例來(lái)自Schmid和Klein(1995)[15]的鑄造充型水模擬實(shí)驗(yàn)，鑄模都由透明材料制成，并用高速攝像機(jī)記錄充型過(guò)程。利用SPH方法來(lái)模擬這兩個(gè)算例，均視為二維問(wèn)題，時(shí)間積分采用跳蛙法，算例中流體都為水，密度為998kg·m-3，動(dòng)力粘度為1.003×10-3Pa·s。

2.1弓型件充型

鑄模幾何尺寸如圖1所示，整體上來(lái)看該型腔為“弓”字型，含有5個(gè)直角轉(zhuǎn)角，2個(gè)弧形轉(zhuǎn)角，底部入口處有一垂直延伸段以便水能夠被腔壁約束順利入腔。實(shí)驗(yàn)中水從型腔底部入口以8.7m·s-1的速度充入腔內(nèi)，不計(jì)重力的影響。

SPH模擬中，在弓型件邊界上布置4 290個(gè)虛粒子，水粒子有27 000個(gè)，粒子的初始間距為1mm，時(shí)間步長(zhǎng)為2×10-7s；充型過(guò)程中粒子的最大速度大致為42m·s-1，則聲速取420m·s-1，邊界排斥力表達(dá)式(15)中D取值為1 764，截止半徑r0=1mm。

本文SPH模擬文獻(xiàn)[9]中模擬結(jié)果和實(shí)驗(yàn)結(jié)果如圖2所示，共選取4個(gè)時(shí)間節(jié)點(diǎn)的結(jié)果進(jìn)行比較，分別是7.15ms，25.03ms，39.34ms和53.64ms。

圖1　弓型腔幾何尺寸示意圖

圖2　圖(a)～(d)為弓形件充型模擬的實(shí)驗(yàn)結(jié)果[15]、圖(e)～(h)為文獻(xiàn)[9]中模擬結(jié)果、圖(i)～(l)為本文SPH仿真結(jié)果

在t=7.15ms時(shí)刻，水流通過(guò)第1個(gè)弧形轉(zhuǎn)角，流向由垂直變?yōu)樗健煞NSPH模擬結(jié)果基本都和實(shí)驗(yàn)保持一致，水頭的位置和形態(tài)較相符。

在t=25.03ms時(shí)刻，水流到達(dá)第4個(gè)直角轉(zhuǎn)角，并在第2個(gè)轉(zhuǎn)角處形成空腔。從水的流態(tài)來(lái)看，兩種模擬結(jié)果都與實(shí)驗(yàn)吻合得較好，比較精確地預(yù)測(cè)到水流在經(jīng)過(guò)第4個(gè)轉(zhuǎn)角后的水頭所呈現(xiàn)出來(lái)的形態(tài)，但對(duì)于第2個(gè)轉(zhuǎn)角處的空腔的大小的預(yù)測(cè)，都比實(shí)驗(yàn)值略小；另外本文結(jié)果中水頭前端上升的高度比文獻(xiàn)中要大，更符合實(shí)驗(yàn)結(jié)果。

在t=39.34ms時(shí)刻，水流經(jīng)過(guò)第6個(gè)轉(zhuǎn)角， 兩種SPH模擬準(zhǔn)確預(yù)測(cè)了水頭撞擊腔壁后所形成的水花的形狀。對(duì)第3個(gè)轉(zhuǎn)角到第4個(gè)轉(zhuǎn)角之間的空腔的大小和形狀的預(yù)測(cè)，本文結(jié)果與實(shí)驗(yàn)結(jié)果更相符；在第4個(gè)轉(zhuǎn)角到第6個(gè)轉(zhuǎn)角之間的空腔的預(yù)測(cè)上，文獻(xiàn)中結(jié)果更客觀；但兩種SPH模擬結(jié)果中第2個(gè)轉(zhuǎn)角處空腔都已經(jīng)消失，可能的原因是SPH模擬中沒(méi)有考慮氣流的影響。

在t=53.64ms時(shí)刻，水流到達(dá)型腔頂部左邊的腔壁，兩種SPH模擬客觀地捕捉第4個(gè)轉(zhuǎn)角與第6個(gè)轉(zhuǎn)角之間的空腔的位置；但第7個(gè)轉(zhuǎn)角到頂部左邊的腔壁之間的空腔，文獻(xiàn)中SPH模擬與實(shí)驗(yàn)結(jié)果更相符。

2.2環(huán)型件充型

本算例模型是一個(gè)圓盤(pán)型中孔零件，實(shí)驗(yàn)中著色水以18m·s-1的速度從底部入口進(jìn)入型腔，鑄模幾何尺寸如圖3所示。

圖3　環(huán)形腔幾何尺寸示意圖

SPH模擬中，在環(huán)型件邊界上布置1 550個(gè)虛粒子，水粒子有30 910個(gè)，邊界虛粒子的初始間距為0.4mm，水粒子的初始間距為0.8mm，時(shí)間步長(zhǎng)為2×10-7s；充型過(guò)程中粒子的最大速度大致為40m·s-1，則聲速取400m·s-1，邊界排斥力表達(dá)式中D取值為1 600，截止半徑r0=0.4mm。

圖4給出了不同時(shí)刻實(shí)驗(yàn)和兩種SPH仿真結(jié)果，其中SPH結(jié)果中還包含粒子軸向速度云圖信息。

圖4　圖(a)～(c)為環(huán)形件充型模擬的實(shí)驗(yàn)結(jié)果[15]、圖(d)～(f)為文獻(xiàn)[9]中模擬結(jié)果、圖(g)～(i)為本文SPH仿真結(jié)果

在t=8.82ms，水流通過(guò)下方垂直流道進(jìn)入型腔沖擊中心型芯，被分為兩股支流，分別向型芯兩側(cè)流動(dòng)繼續(xù)充填型腔，在沖擊型腔外輪廓后又再次分叉，分別沿壁面流向型腔頂部和底部，形成四股支流。兩種SPH模擬與實(shí)驗(yàn)結(jié)果吻合程度都較高，水流形態(tài)和流向預(yù)測(cè)準(zhǔn)確；但本文模擬結(jié)果中的流體粒子更有序。

在t=11.76ms，型腔頂部的兩股支流匯合，沖向中心型芯上部，而底部的支流則繼續(xù)沿型腔外輪廓流動(dòng)，并與垂直流道進(jìn)入型腔的水流匯合，形成左右兩個(gè)空腔。兩種SPH模擬對(duì)于水流走向和空腔形成預(yù)測(cè)基本準(zhǔn)確，但對(duì)于型腔底部匯流后的水流沖向型芯上部后的形態(tài)的預(yù)測(cè)，本文模擬結(jié)果與實(shí)驗(yàn)更符合。

在t=16.17ms，型腔內(nèi)出現(xiàn) 4個(gè)空腔，隨著后續(xù)水流不斷沖入型腔，空腔逐漸縮小并完成充型過(guò)程。本文SPH模擬捕捉到空腔的形成過(guò)程，且空腔形狀飽滿，粒子分布有序；但在本文模擬結(jié)果中上部?jī)蓚€(gè)空腔的大小比實(shí)驗(yàn)值略小，可能與氣壓有關(guān)。

從弓形件和環(huán)形件的充型模擬結(jié)果來(lái)看，本文SPH模擬結(jié)果比較客觀地反映了水模擬充型實(shí)驗(yàn)進(jìn)程，與文獻(xiàn)[9]中的結(jié)果對(duì)比也說(shuō)明了本文所采用對(duì)稱結(jié)構(gòu)的SPH離散公式，能夠較成功地處理粒子不一致問(wèn)題。

本文將光滑粒子流體動(dòng)力學(xué)(SPH)方法應(yīng)用到鑄造充型問(wèn)題上，通過(guò)構(gòu)建弱可壓縮流體的計(jì)算模型，利用SPH算法仿真水模擬鑄造充型過(guò)程，并與文獻(xiàn)中的模擬和實(shí)驗(yàn)結(jié)果比較，發(fā)現(xiàn)本文SPH算法結(jié)果不僅能夠較準(zhǔn)確地模擬流體自由表面的形態(tài)，而且對(duì)于流體走向和空腔的形成過(guò)程預(yù)測(cè)也比較吻合。在SPH算法設(shè)計(jì)中，控制方程采用具有對(duì)稱結(jié)構(gòu)的離散形式，動(dòng)量方程壓力項(xiàng)中加入人工粘度，利用邊界排斥力法施加固定邊界條件，這些技術(shù)的應(yīng)用對(duì)于逼真的模擬效果起到顯著作用。結(jié)果表明，本文所采用SPH方法及其技術(shù)處理對(duì)于模擬鑄造充型過(guò)程是有效的。

[1]帕坦卡SV. 傳熱與流體流動(dòng)的數(shù)值計(jì)算[M]. 郭寬良，譯. 合肥：安徽科學(xué)技術(shù)出版社，1984：132-142.

[2]IssaRI.Solutionoftheimplicitlydiscretizedfluidflowequationbyoperatorsplitting[J].JournalofComputationalPhysics,1986,62：40-65．

[3]HarlowFH,WelchJE.Numericalstudyoflarge-amplitudefree-surfacemotions[J].PhysicsofFluid, 1966,9：842-851．

[4]HirtCW,NicholsBD.Volumeoffluid(VOF)methodforthedynamicsoffreeboundaries[J].JournalofComputationalPhysics,1981,39(1)：201-225．

[5]HarlowFH.PICanditsprogeny[C].WorkshoponParticleMethodsinFluidDynamicsandPlasmaPhysics.LosAlamos,NM,USA,13-15April1987．

[6]OsherS,FedkiwR.Levelsetmethods:Anoverviewandsomerecentresults[J].JournalofComputationalPhysics, 2001,169：463-502．

[7]LucyLB.Numericalapproachtotestingthefissionhypothesis[J].AstronomicalJournal，1977， 82：1013-1024.

[8]GingoldRA,MonaghanJJ.Smoothedparticlehydrodynamics:Theoryandapplicationtonon-sphericalstars[J].RoyalAstronomicalSociety，MonthlyNotices1977，181：375-389.

[9]曹文炅. 鑄造充型過(guò)程SPH方法建模及數(shù)值模擬[D]. 廣州：華南理工大學(xué)，2011年.

[10]強(qiáng)洪夫，韓亞偉，王坤鵬，等.基于罰函數(shù)SPH新方法的水模擬充型過(guò)程的數(shù)值分析[J].工程力學(xué)，2011，28(1)：245-250.

[11]LiuMB,LiuGR.SmoothedParticleHydrodynamics(SPH):anOverviewandRecentDevelopments[J].ArchComputMethodsEng， 2010，(17): 25-76.

[12]MonaghanJJ.SimulationfreesurfaceflowswithSPH[J].JournalofComputationalPhysics, 1994, 110: 399-406.

[13]MonaghanJJ.AnintroductiontoSPH[J].ComputerPhysicsCommunications, 1988, 48: 89-96.

[14]LattanzioJC,MonaghanJJ,PongracicH, et al.Controllingpenetration[J].SIAMJournalonScientificandStatisticalComputing, 1986, 7(2): 591-598.

[15]SchmidM,KleinF.Fluidflowindiecavities-experimentalandnumericalsimulation[C].Indianapolis:NADCA18.InternationalDieCastingCongressandExposition, 1995: 93-99.

責(zé)任編輯王菊平

A numerical study on foundry filling process based on SPH method

ZHOU Xue-jun1,2, CHEN Ding2

(1.College of Mathematics and Physics, Huanggang Normal University, Huangzhou 438000, Hubei, China;2. College of Mechanics and Materials, Hohai University, Nanjing 210098, Jiangsu, China)

Numerical simulation for foundry filling process is always a hot research topic of casting. Smoothed Particle Hydrodynamics (SPH) is presented to simulate foundry filling process in this paper. A mature meshless particle numerical method, the SPH method is good for mimicking large deformation. By establishing the model of filling process, the SPH method is applied to simulate water filling process of bow-shaped and annular models. Compared with the results of the literature and experiments, the SPH algorithm presented in this paper shows its advantages in solving foundry filling problems.

smoothed particle hydrodynamics; foundry filling; meshless; numerical simulation

TG21

A

1003-8078(2016)03-0012-06

2016-04-25

10.3969/j.issn.1003-8078.2016.03.04

周學(xué)君，男，湖北蘄春人，講師，博士，主要研究方向?yàn)榱W(xué)中的數(shù)值方法等。

湖北省教育廳科學(xué)技術(shù)研究項(xiàng)目(B2015218)。