周學(xué)君，陳　丁

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

?

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

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

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

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

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

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

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

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

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

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

1.1SPH方法的近似

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

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

(1)

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

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

(2)

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

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

(3)

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

(4)

1.2控制方程的SPH離散

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

(5)

(6)

(7)

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

(8)

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

(9)

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

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

(10)

(11)

(12)

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

(13)

其中

(14)

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

1.3邊界條件

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

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

(15)

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

2　鑄造充型模擬

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

2.1弓型件充型

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

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

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

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

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

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

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

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

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

2.2環(huán)型件充型

本算例模型是一個圓盤型中孔零件，實驗中著色水以18m·s-1的速度從底部入口進入型腔，鑄模幾何尺寸如圖3所示。

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

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

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

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

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

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

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

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

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

[1]帕坦卡SV. 傳熱與流體流動的數(shù)值計算[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]曹文炅. 鑄造充型過程SPH方法建模及數(shù)值模擬[D]. 廣州：華南理工大學(xué)，2011年.

[10]強洪夫，韓亞偉，王坤鵬，等.基于罰函數(shù)SPH新方法的水模擬充型過程的數(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é)君，男，湖北蘄春人，講師，博士，主要研究方向為力學(xué)中的數(shù)值方法等。

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