- Email: [email protected]
The analysis of mould filling in castings using the finite element method
Journal of Materials Processing Technology, 38 (1993) 291-302
291
Elsevier
The analysis of mould filling in castings using the finite element method A.S.Usmani J.T.Cross R.W.Lewis* Institute for Numerical Methods in Engineering University College of Swansea Singleton Park, Swansea, SA2 8PP, UK
Abstract
For an accurate modelling of the process of solidification in castings, it is necessary in most cases to include the modelling of the process of filling of the mould. The flow associated with the pouring of molten metal into moulds has been modelled using the Finite Element Method. This paper gives a brief outline of how an FE based Navier-Stokes equation solver has been used to analyse the flow and the velocities obtained used to advect a pseudo-concentration function for nmdelling the fluid front. A heat transfer algorithm is then used Io predict the temperatures throughout the mould and the moving fluid. Some examples have been included, which show the application of these techniques to real casting problems. Examples include the solution of the plain mould filling of a complex practical problem, mould filling with heat transfer and a problem in which mould filling and solidification are coinbined.
1
Introduction
The primary benefit of a mould filling analysis preceding a thermal analysis of a casting is that instead of the usual method of starting a thermal analysis by assuming a constant temperature in the mould and melt, a realistic (non-uniform) initial temperature will be available. Similarly, an initial (solenoidal) velocity field will be available if the effects of inertial/natural convection are to be taken into account. Examination of the flow patterns during filling can reveal a wide range of potential difficulties with the height, dimensions and positioning of risers or the gating system. By including phase change effects during filling it is possible to predict freeze-off, where the flow through a narrow section freezes before reaching its intended destination, and cold shuts, where molten metal meets already *Correspondance to: Professor it.W.Lewis
0924-0136/93/$06.00 © 1993 Elsevier SciencePublishers B.V. All rights reserved.
292
solidified m e t a l causing a fault liue. W i t h evcu siulph' dies costing tens of tt!.,m~ands o! p o u n d s the ability to expose flaws iu a die design b~:tore ',he die is madu c a i , - , d i i~e considerable cost reduclion. The Inethod used %r the modelling of mould tilling is based upon the vo/v;~u , , j j i u , [ m e t h o d (VOF)[]] (also called the pseudo concentration m e t h o d [2] o," sat uraih~, iuu~ ion, m e t h o d [3, 4]) which is. in turn, based upoll the w(!I[ knowll Hmv~:(~" o.nd ,'f// } ~ t . \ ( ' technique. The M A C iechnique involves using a large uumber of ula.rkcr par![~ h> -:!n~ ad all over the fluid occupied region and ea,,:h parlicle is tracked as ]1 movc~s wiI!~ } i , flow velocity at. its p a r t i c u l a r location, li extracts a high price iu terms of ('ou~i)ul.u~ -:?(,ra,~,u and C P U time. especially for 3 D at)plicaiious. F h e \ O F method r e t a h > ih~ , < ~ ~.pi ~i~i s i m p l i c l y of l.he M A C technique but d e m a n d s c(msidcrab[ 3 less (omt)uicr c[[,, '. ['}..c m a r k e r particles a.r(~ replaced by a nodal variable which is advecicd u.-ing t}lc ,. '!,~<~ i,.~; from the solutiol~ of lhe NaMer Stokes equat.ions. The main a,lvautag, d ~}..-: ~o!},.i is t h a t the m o v i n g t ' r o n i does not. have to be d('tcrnfiuedexplicitly, uei!h,'r &~X}, ;<~i b o u n d a r y conditions havo to t)e satisfied oxplici!lx. T h e flow o[ molten metal into a ~asliug is iuvar}a.btv iurbulenl., ~-t tii,,~ . ~u~iii} cation for modelling it, as lalninar. Mii,aic ct. al.[5] negle(i lurt)ulenc(' o]I Ihc g,~,iiui> :ha! the t i m e taken t\)r a mould to till is generally very short and ihu uc('d [or {imu a,-iagi~e: in t u r b u l e n c e modelling tileal|s t h a t the tlow must haw., a fin't(' tim~ lo devciop. \i~!i [\ui m a n ' s laws of the wall have beer, used to obtain bouudary {ractions oi) the \~a[{,~ (~l i h ( cavity and this lorces the fluid velocity profiles to mimic those of turbulenl flow. l ii~v,c~(.~. the chief m o t i v a t i o n for modelling mould filling is io obtaii~ a realistic *cu,per, t u~,,, iieh! at the end of filling. This will be possible if the time ,)f filling from a laminar iio'v }ms.'d model is close to reality.
2
E q u a t i o n s and formulation
The Navier Stokes equations for incon,t)ressible tlow iu two-dimeusions arc used i , ruo,.hfi the ttow. The equations include the mass col,scrva.tiou (co~ltitmiiy equal ioil) w}iici~ }. zl~ ,i to enforce i n c o m p r e s s i b i l i b ~. (7onsc.,'valiol~
oi"ma.~s V-v
I)
'}:
in which v is velocity. Conservatiov of mornentum Dv
,oh7
= V-
r ~ f-,g
(Q
Here, r is the stress tensor, p and g are density and gravitational acceleration. ~md D represents the total or substantial derivative. 'Fhese equatiolls have been used cxit!usive]/ for m o d e l l i n g incompressibh, flow and further details (;an be [ound iu a b u n d a m ' c in the literature[6, 7].
293 Different material properties are used for the elements in the fluid filled region and the empty region. Fictitious values of density, and viscosity are used for the fluid to achieve a low Re and keep the flow in the laminar regime. Density and viscosity close to that of air are used for the elements in the empty region. For the elements containing the fluid front the properties are averaged according to the full and empty fractions of lhe element by interpolation at the integration points. 1)irichlet or essential boundary conditions for the Navier -Stokes equations are specified velocities at the boundaries. Pressure may not be specified at the bo,mda.ry as it is an implicit variable in an incompressible flow[8] which 'adjusts' itself to deliver a solenoidal velocity field. Nenmann or natural boundary, conditions (:an be the normal and t.angenlial traction forces,
f..
=
~?Ln
- P + 2~ t),-~
"f~ = # \ O r
(3)
+ O,z J
(1)
where n and v are the unit normal and tangent vectors respectively, P is 1he pressure, 11 is the dynamic viscosity. The equations are discretized by using a conventional Galerkin weighted resid~la.I technique and a mixed velocity pressure formulation is used to solve the equations. \"elocities (v), pressure (P), and pseudo-concentration (F) are approximated over each elemenl by the equations r~
i=1
l)~' = E
'i Pi
i=1
i=1
For an element e, rz represents the number of velocity and pseudo-concentration nodes and n' the number of pressure nodes. Similarly N represenls the quadratic basis functions used for velocity and pseudo-concentration and A,." the linear basis functions for pressure. This difference in interpolation functions prevents mesh locking caused by over-constrainment[9, page 208]. Substituting the shape functions (5) into equations (1) (4) yiehts a set of fully coupled first order ordinary differential equations with respect to time which can be written M0+K0=F (6) where 0 represents all the variables. The details of the matrices in Eqn. (6) can be found in references, such as [8]. Equations (6) are discretized in time using the implMt backward difference method which results in,
[Mn+l -~ AIK,~+I] (0n+l) = [Mn+l] (On) 4- (,~tF..+l)
(7)
294
3
Fluid front tracking
As mentioned in Section 1, a pseudo-concentration flm<:tion is used to track the ire( fluid surface. If this l.unction is represented by F ( x , y . t ) for 2-I) flow. the firm ordel pure advection equation which conserves the function F ( x , Y, t) is written as. DF
l)t This function moves with the flow, and so, if a particular value of this functio> t:' is associated with the free fluid surf,ace [2, 10], it can be tracked over the timeslei)!~ by simply plotting the contour of Fc at each timestep. In reality l" is a slep function: hut, here a continuous function has been used to avoid tim numerical diffictlllies inh,'ra.q~, i~i advecting a step function. The value of F ( x , y , i ) - 0.0 I;; is used 4<) m a r k ihe fr,,, ~ surface, while F(x, y,l) > 0.0 indicates the fluid regi
(+ii
This is a first order ordinary differential equation which is discretized in tiln~' u.4ng !,he backward Euter method. After solving the Navier Stokes equations as described ,~a[liet:. Eqn. (9) is solved at. each timestep using the velocities obtained, which restl!{.~ in ,n u p d a t e d pseudo-concentration field. The position of the free fluid front is coimid,:ut vvit}a the contour of/+'+ or ]"(x, y, 1) = 0.0. The material properties of the elements a++ :~,,)ditied according to the new position of the fluid l,ronl. Von K a r m a n ' s laws of the wall are used to oblaiII the ttow houn'.tar 5 (t,!idiiioi{s at the wMls of the cavity. These laws are normally applied to turbulenl flow i>~oblem>, however here ~s in [10], they are used to obtain l,angenlial tractions at the wall. ~+l i;h(~ mould cavity. The velocity normal to the walls is specified as zero. l ' h c i~.~++genliat tractions are calculated from
where u~, is the tangential velocity to the wall. and ~+ = 5.5 + 2.5 In .V+ where an appropriate value of y + > 30 must be assuined. In fact, the tra ctio~s obtained fl-om Eqn. (10) are very low and the value of y+ does not seem to affect the final analysis significantly. The outflow boundary conditions are specified as traction fr('{~ The basic algorithm as described above was successfully tested using a. be~tchmark problem of flow between two parallel plates [11]. "lb keep the boundaries not parallel to either of the cartesian axes, impervious, it is necessary to use rolation transtbrmation techniques. Pinder and Gray [12, page 281] and Enge/man eL. a/.[13] have described the m e t h o d in detail. The details of the procedure with respect to the presenl problem cem also be found in [11].
295
4
Filling of a real casting
A 2-D representation of a real casting as shown in Fig. 1 including the mesh and a typi<'al velocity field, was solved for simulating the fill using the code developed. The results appear in Fig. 2 as fluid front profiles at different times.
5
Improved algorithms for filling and coupled heat transfer calculations
Owing to the central difference nature of conventional Galerkin finite element approximations, oscillations occur in the solution when the finite element mesh is of insufficient refinement to model the variations of the field variable. In most cases such spurious oscillations may be eliminated by judicious mesh refinement in the zones where regions of high gradient are expected. In transient problems a timestep size larger than the limits dictated by the physical properties and the mesh size also results in an oscillatory solution, which may well become meaningless if the magnitude of the oscillations masks the real values. Spatial and temporal refinement will remedy this for most diffusion dominated problems, but will often cease to be practical for problems of pure advection or advection dominated transport. In the present advection dominated problem of filling and coupled heat transfer, the mathematically robust Taylor-Galerkin method forwarded by Donea[14] has been used. lie developed these methods for different time marching schemes and was guided by the work of Morton and Parrott[15] who showed that for each particular timestepping method there is a corresponding optimal form of Petrov-Galerkin weighting functions. This led Donea to discretize the pure advection equation firstly in time, improving the difference approximation of the time derivative by including higher order Taylor series terms, and secondly in space, by using the conventional Bubnov-Galerkin finite element approach, thus leading to the 'Taylor-Oalerkin' method. This method can be seen as the finite element counterpart of the Lax-Wendroff schemes used in finite differences [16]. '['he details of the Taylor Galerkin procedure for this problem can again be found in reference [11].
5.0.1
A test example of coupled filling and thermal analysis
Fig. 3 shows a simple problem of gravity filling between two vertical walls of a mould, which absorb heat from the molten metal as it flows into the void. For this problem realistic thermal properties of mould, metal and air haw~ been used in the appropriate regions. The properties of a typical aluminium alloy have been used for the full region, air properties have been used for the empty regions and steel properties for the mould regions. These are shown in Fig. 3 The properties used for the Navier Stokes solver correspond to an Re of approximately 20.0. When this pr~,blem was solved for these properties, using GFEM and the
r~
J
i r
~
II................... r............•
........
©
c ~
,H
,H ~i~i
~'
~~ ~ii~i!!~ii!i!ii!i!ii!iii~ ~....
i
~iiiii~!~!i!i~i~i~i!~
iiii~i~i~i~i~
~
~ii~i~
iiii~i~i
H
298
,1 ....
-
I!.tJ
T= 70~C I"
o.o2m ~'-r
•
o.o3m
.
!iiii iTiiiiTiiii iiiiiiii771 | ................
li.:: T = 25~C
~,~.~
!
:.i:i:iiiit
~
li::::!i:::ii:iT:il
:: " , ":::: .............. :::Mould::
:
~
i
,
k=200mJ~-,i~
_.
~L
i:!
i •
. ::
L-:
=900 kglK
,{
"
li
P = 2700 ~ ' l
~
o.o2m
:
'lP :: 28' C
.
i i MOuldi.i
'::::~i2~5~Ciii~
+
:: /"
:
' ) ) T.=i25~iCi i i "-
.lm
Air To= 25~' C
.. . . .. . .. .. . . .. . .. . li:i: -
:i
.
i
.:
.
i
.
.
.
.
.
.
.
p=LO
iiliL
- ............
.
5-.
.
.
. :: ~:-7soo!~:! 5
_ :: e a-:450 i:~ c":
.
lt2g 71'~,--
kg m3
r: _.
k = .025mJY-,, K .
.
.
.
.
.
.
.
.
.
.
_
........
i
:
. . . . .
/.,
c = lO00kg~
-.-i =•.
.
.
.
. . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. .
:
•
-
0.0
F i g u r e 3: C o n f i g u r a t i o n for vertical filling with coupled heat transfer.
299
i
/ •
t ~k,S
' • ~
Min. value =
25.0000 Interv~ = 3 2 , 1 4 2 9
Max. value= 700.0000
Time =
0.02
MirLv , a u c = 2 5 . 0 0 0 0 Max. value = 7 0 0 . 0 0 0 0
ln~r,,al=
32.1429
Time = 0 . 0 6
r
~i ~'~ --1 Min. value= 25.12000 I n t e ~ = 32.1429
Min. value=
Ms,x. value= 700.0000 Time= 0.08
M ~ . value = 7 0 0 . 0 0 0 0
2 5 . ( X ) 0 0 interval = 3 2 . 1 4 2 9 Time = 0 . 1 0
ti
;I ~I
i
k~C
Min. v N u e = 2 5 . 0 0 0 0 M~.value=
700.0000
Intern = 32.1429
Min. vMue =
Time= 0.12
Max. v a l u e =
2 5 . 0 0 0 0 Interva.l = 3 2 . 1 4 2 9 7130.0000 T i m e = 0 . 1 4
F i g u r e 4: Explicit solution for vertical filling with coupled heat transfer: isotherms.
©
©
11
o
~z
~o
Q
J
~i~i~i~~ ~ ~
o
~:~ ~
~'~~ "
....i ......
~'
301
implicit backward Euler (BE) timestepping scheme, the resulls were riddled with oscillations, in addition to being overdiffuse. Results obtained fl'om the explicit Taylor Galerkin analysis, are shown in Fig. 4. It is clear from the isotherm plots that a sharp thermal front is advected faithfully without oscillation or excessive dissipation. Being explicit this solution was, however, expensively obtained from 140 timesteps. The timestep size was 0.001 corresponding to a Courant number (C) of roughly 0.5.
6
S i m u l a t i o n of solidification during filling
To simulate solidification during filling several additional problems must be addressed. The first and most important is the introduction of an algorithm to model the latent heat evolution. The second problem is of simulating the freezing of the fluid, which has been achieved by arresting the flow at the boundaries and increasing the viscosity of the fluid away from the boundary at the freezing point. A preliminary example of this procedure appears in Fig. 5.
7
Conclusions
From the number of problems solved here, it can be concluded that a finite element flow and heat transfer code such as that described in this work, when used with a front tracking algorithm such as the VOF method or the pseudo-concentration method based on the pure adw.'ction equation, provides an effective means ot' soving problems such as that of mould filling. This method can be used in other applications as well, such as forming processes involving creeping viscous flows, multiphase flow etc. The main advantage of this method is that the moving front does not have to be explicitly determined, neither do the boundary conditions on the front have to be explicitly satisfied. Furthermore, it is necessary to use special techniques, such as the Taylor Galerkin nmthod to deal with the problems arising from transporting a scalar when advection is the dominant mode. This is particularly important for obtaining meaningful temperature fields associated with filling.
A c k n o w l e d g e m e n t : The authors wish to thank the SERC who supported this work via ACME grant no. G R / F 32790.
References [1] C.W.Hirt and B.D.Nichols. Volume of fluid (vof) inethod for the dynamics of free boundaries. Journal of Computational Phgsics , 39:201 225, 1981.
302 [2] E.Thompson. Use of pseudo-concentrations to follow creeping viscous flows during transient analysis. International Journal for Numerical Methods i7~.Fluids, 6:749 7(i l, 1986. [3] T.J.Smith and D.B.Welbourn. The integration of geometric modelling with iir~ii~'(ei~ ment analysis for tile computer-aided design of castings. Applied Scientific t~i,,~archo 44:139- 160, 1987. [4] R.W.Lewis, K.Morgan, E.D.L.Pugh, and T.,I.Smith, A now on discontinuous ~m~m::r ical solutions of the kinematic wave equation. Inlernatio~al Journal for Vum,:'~-,',l Methods in Engineering, 20:555 563, 198'1. [5] B.Minaie, K.A.Stelson, and V.R.Volter. Fluid flow and solidification modellim< o~ die castings. ASME, Modelling of Materials Processing, MD, 3:35 50, 1987. [6] H.Sehlichting. Boundary-Layer 7'heorg. McGraw-lfill Book Company, i968. [7] D.K.Edwards, V.E.Denny, and A.F.Mills. 1)'ansfer Proc~.~,~es - A~z h~trodn,~o, to Diffusion, ConvectioT~ and Radiation. McGraw-Hill Book ('.ompany, 197:1. [8] P.M.Gresho, R.b.Lee, and R..L.Sani. On the time-dependent solution ~,1 the ~,~ou~ pressible Navier-Stokes equations in two and three dimensions. In Rece~d Ad~a~c¢,~ in Numerical Methods in Fluids, volume I. Pineridge Press l&nited. Swansea I;180. [9] T..l.R.Hughes. The Finite Element Method - Linear Static and Dynamic /:i,lt~ t;l~ment Analysis, Prentice-Hall international, Inc., Englewood Cliffs, New ,Jersey 07632, 1987. [10] D.M.Gao, G.Dhatt, 3.Belanger, and A,B.Cheikh. A finite element simulation of met.ai flow in moulds. [n Sixth International Co~@rence for Numerical Methods in 7 ht:~"mal Problems, Swansea, U.K., July 1989. Pineridge Press, Swansea. [11] A.S.Usmani, 3.T.Cross, and R.W.Lewis. A finite element model h)r the simuia!i(m of mould filling in metal casting and the associated heat transl~r. International ,t,-~t~l~taJ for Numerical Methods in Engineering, to be published in 1992. [12] G.F.Pinder and W.G.Gray. Finite h;lement Simulation i7~ Surface aT~d ,5'ub.~tlJ),:~ Hydrology. Academic Press, 1977. [13] M.S.Engelman, R.L,Sani, and P.M.Gresho. The implementation of normal ami/t~I tangential boundary conditions in finite element codes ior incompressible ttuid ttow International dournal for Numerical Methods in Fluids: 2:225 238, 1982. [14] J.Donea. A Taylor-Galerkin method for convective transport problems, lnter~ational Journal for Numerical Methods in Engineering, 20:1(11 119, t984. [15] K.W.Morton and A.K.Parrott. Generalized Galerkin methods for first order ~!yl)er bolic equations. Journal of Computational Physics, 36:249 270, 1980. [16] P.J.Roache. Computational Fluid Mechanics. U.S.A., 1976.
Hermosa Publishers, Albuquerque,
291
Elsevier
The analysis of mould filling in castings using the finite element method A.S.Usmani J.T.Cross R.W.Lewis* Institute for Numerical Methods in Engineering University College of Swansea Singleton Park, Swansea, SA2 8PP, UK
Abstract
For an accurate modelling of the process of solidification in castings, it is necessary in most cases to include the modelling of the process of filling of the mould. The flow associated with the pouring of molten metal into moulds has been modelled using the Finite Element Method. This paper gives a brief outline of how an FE based Navier-Stokes equation solver has been used to analyse the flow and the velocities obtained used to advect a pseudo-concentration function for nmdelling the fluid front. A heat transfer algorithm is then used Io predict the temperatures throughout the mould and the moving fluid. Some examples have been included, which show the application of these techniques to real casting problems. Examples include the solution of the plain mould filling of a complex practical problem, mould filling with heat transfer and a problem in which mould filling and solidification are coinbined.
1
Introduction
The primary benefit of a mould filling analysis preceding a thermal analysis of a casting is that instead of the usual method of starting a thermal analysis by assuming a constant temperature in the mould and melt, a realistic (non-uniform) initial temperature will be available. Similarly, an initial (solenoidal) velocity field will be available if the effects of inertial/natural convection are to be taken into account. Examination of the flow patterns during filling can reveal a wide range of potential difficulties with the height, dimensions and positioning of risers or the gating system. By including phase change effects during filling it is possible to predict freeze-off, where the flow through a narrow section freezes before reaching its intended destination, and cold shuts, where molten metal meets already *Correspondance to: Professor it.W.Lewis
0924-0136/93/$06.00 © 1993 Elsevier SciencePublishers B.V. All rights reserved.
292
solidified m e t a l causing a fault liue. W i t h evcu siulph' dies costing tens of tt!.,m~ands o! p o u n d s the ability to expose flaws iu a die design b~:tore ',he die is madu c a i , - , d i i~e considerable cost reduclion. The Inethod used %r the modelling of mould tilling is based upon the vo/v;~u , , j j i u , [ m e t h o d (VOF)[]] (also called the pseudo concentration m e t h o d [2] o," sat uraih~, iuu~ ion, m e t h o d [3, 4]) which is. in turn, based upoll the w(!I[ knowll Hmv~:(~" o.nd ,'f// } ~ t . \ ( ' technique. The M A C iechnique involves using a large uumber of ula.rkcr par![~ h> -:!n~ ad all over the fluid occupied region and ea,,:h parlicle is tracked as ]1 movc~s wiI!~ } i , flow velocity at. its p a r t i c u l a r location, li extracts a high price iu terms of ('ou~i)ul.u~ -:?(,ra,~,u and C P U time. especially for 3 D at)plicaiious. F h e \ O F method r e t a h > ih~ , < ~ ~.pi ~i~i s i m p l i c l y of l.he M A C technique but d e m a n d s c(msidcrab[ 3 less (omt)uicr c[[,, '. ['}..c m a r k e r particles a.r(~ replaced by a nodal variable which is advecicd u.-ing t}lc ,. '!,~<~ i,.~; from the solutiol~ of lhe NaMer Stokes equat.ions. The main a,lvautag, d ~}..-: ~o!},.i is t h a t the m o v i n g t ' r o n i does not. have to be d('tcrnfiuedexplicitly, uei!h,'r &~X}, ;<~i b o u n d a r y conditions havo to t)e satisfied oxplici!lx. T h e flow o[ molten metal into a ~asliug is iuvar}a.btv iurbulenl., ~-t tii,,~ . ~u~iii} cation for modelling it, as lalninar. Mii,aic ct. al.[5] negle(i lurt)ulenc(' o]I Ihc g,~,iiui> :ha! the t i m e taken t\)r a mould to till is generally very short and ihu uc('d [or {imu a,-iagi~e: in t u r b u l e n c e modelling tileal|s t h a t the tlow must haw., a fin't(' tim~ lo devciop. \i~!i [\ui m a n ' s laws of the wall have beer, used to obtain bouudary {ractions oi) the \~a[{,~ (~l i h ( cavity and this lorces the fluid velocity profiles to mimic those of turbulenl flow. l ii~v,c~(.~. the chief m o t i v a t i o n for modelling mould filling is io obtaii~ a realistic *cu,per, t u~,,, iieh! at the end of filling. This will be possible if the time ,)f filling from a laminar iio'v }ms.'d model is close to reality.
2
E q u a t i o n s and formulation
The Navier Stokes equations for incon,t)ressible tlow iu two-dimeusions arc used i , ruo,.hfi the ttow. The equations include the mass col,scrva.tiou (co~ltitmiiy equal ioil) w}iici~ }. zl~ ,i to enforce i n c o m p r e s s i b i l i b ~. (7onsc.,'valiol~
oi"ma.~s V-v
I)
'}:
in which v is velocity. Conservatiov of mornentum Dv
,oh7
= V-
r ~ f-,g
(Q
Here, r is the stress tensor, p and g are density and gravitational acceleration. ~md D represents the total or substantial derivative. 'Fhese equatiolls have been used cxit!usive]/ for m o d e l l i n g incompressibh, flow and further details (;an be [ound iu a b u n d a m ' c in the literature[6, 7].
293 Different material properties are used for the elements in the fluid filled region and the empty region. Fictitious values of density, and viscosity are used for the fluid to achieve a low Re and keep the flow in the laminar regime. Density and viscosity close to that of air are used for the elements in the empty region. For the elements containing the fluid front the properties are averaged according to the full and empty fractions of lhe element by interpolation at the integration points. 1)irichlet or essential boundary conditions for the Navier -Stokes equations are specified velocities at the boundaries. Pressure may not be specified at the bo,mda.ry as it is an implicit variable in an incompressible flow[8] which 'adjusts' itself to deliver a solenoidal velocity field. Nenmann or natural boundary, conditions (:an be the normal and t.angenlial traction forces,
f..
=
~?Ln
- P + 2~ t),-~
"f~ = # \ O r
(3)
+ O,z J
(1)
where n and v are the unit normal and tangent vectors respectively, P is 1he pressure, 11 is the dynamic viscosity. The equations are discretized by using a conventional Galerkin weighted resid~la.I technique and a mixed velocity pressure formulation is used to solve the equations. \"elocities (v), pressure (P), and pseudo-concentration (F) are approximated over each elemenl by the equations r~
i=1
l)~' = E
'i Pi
i=1
i=1
For an element e, rz represents the number of velocity and pseudo-concentration nodes and n' the number of pressure nodes. Similarly N represenls the quadratic basis functions used for velocity and pseudo-concentration and A,." the linear basis functions for pressure. This difference in interpolation functions prevents mesh locking caused by over-constrainment[9, page 208]. Substituting the shape functions (5) into equations (1) (4) yiehts a set of fully coupled first order ordinary differential equations with respect to time which can be written M0+K0=F (6) where 0 represents all the variables. The details of the matrices in Eqn. (6) can be found in references, such as [8]. Equations (6) are discretized in time using the implMt backward difference method which results in,
[Mn+l -~ AIK,~+I] (0n+l) = [Mn+l] (On) 4- (,~tF..+l)
(7)
294
3
Fluid front tracking
As mentioned in Section 1, a pseudo-concentration flm<:tion is used to track the ire( fluid surface. If this l.unction is represented by F ( x , y . t ) for 2-I) flow. the firm ordel pure advection equation which conserves the function F ( x , Y, t) is written as. DF
l)t This function moves with the flow, and so, if a particular value of this functio> t:' is associated with the free fluid surf,ace [2, 10], it can be tracked over the timeslei)!~ by simply plotting the contour of Fc at each timestep. In reality l" is a slep function: hut, here a continuous function has been used to avoid tim numerical diffictlllies inh,'ra.q~, i~i advecting a step function. The value of F ( x , y , i ) - 0.0 I;; is used 4<) m a r k ihe fr,,, ~ surface, while F(x, y,l) > 0.0 indicates the fluid regi
(+ii
This is a first order ordinary differential equation which is discretized in tiln~' u.4ng !,he backward Euter method. After solving the Navier Stokes equations as described ,~a[liet:. Eqn. (9) is solved at. each timestep using the velocities obtained, which restl!{.~ in ,n u p d a t e d pseudo-concentration field. The position of the free fluid front is coimid,:ut vvit}a the contour of/+'+ or ]"(x, y, 1) = 0.0. The material properties of the elements a++ :~,,)ditied according to the new position of the fluid l,ronl. Von K a r m a n ' s laws of the wall are used to oblaiII the ttow houn'.tar 5 (t,!idiiioi{s at the wMls of the cavity. These laws are normally applied to turbulenl flow i>~oblem>, however here ~s in [10], they are used to obtain l,angenlial tractions at the wall. ~+l i;h(~ mould cavity. The velocity normal to the walls is specified as zero. l ' h c i~.~++genliat tractions are calculated from
where u~, is the tangential velocity to the wall. and ~+ = 5.5 + 2.5 In .V+ where an appropriate value of y + > 30 must be assuined. In fact, the tra ctio~s obtained fl-om Eqn. (10) are very low and the value of y+ does not seem to affect the final analysis significantly. The outflow boundary conditions are specified as traction fr('{~ The basic algorithm as described above was successfully tested using a. be~tchmark problem of flow between two parallel plates [11]. "lb keep the boundaries not parallel to either of the cartesian axes, impervious, it is necessary to use rolation transtbrmation techniques. Pinder and Gray [12, page 281] and Enge/man eL. a/.[13] have described the m e t h o d in detail. The details of the procedure with respect to the presenl problem cem also be found in [11].
295
4
Filling of a real casting
A 2-D representation of a real casting as shown in Fig. 1 including the mesh and a typi<'al velocity field, was solved for simulating the fill using the code developed. The results appear in Fig. 2 as fluid front profiles at different times.
5
Improved algorithms for filling and coupled heat transfer calculations
Owing to the central difference nature of conventional Galerkin finite element approximations, oscillations occur in the solution when the finite element mesh is of insufficient refinement to model the variations of the field variable. In most cases such spurious oscillations may be eliminated by judicious mesh refinement in the zones where regions of high gradient are expected. In transient problems a timestep size larger than the limits dictated by the physical properties and the mesh size also results in an oscillatory solution, which may well become meaningless if the magnitude of the oscillations masks the real values. Spatial and temporal refinement will remedy this for most diffusion dominated problems, but will often cease to be practical for problems of pure advection or advection dominated transport. In the present advection dominated problem of filling and coupled heat transfer, the mathematically robust Taylor-Galerkin method forwarded by Donea[14] has been used. lie developed these methods for different time marching schemes and was guided by the work of Morton and Parrott[15] who showed that for each particular timestepping method there is a corresponding optimal form of Petrov-Galerkin weighting functions. This led Donea to discretize the pure advection equation firstly in time, improving the difference approximation of the time derivative by including higher order Taylor series terms, and secondly in space, by using the conventional Bubnov-Galerkin finite element approach, thus leading to the 'Taylor-Oalerkin' method. This method can be seen as the finite element counterpart of the Lax-Wendroff schemes used in finite differences [16]. '['he details of the Taylor Galerkin procedure for this problem can again be found in reference [11].
5.0.1
A test example of coupled filling and thermal analysis
Fig. 3 shows a simple problem of gravity filling between two vertical walls of a mould, which absorb heat from the molten metal as it flows into the void. For this problem realistic thermal properties of mould, metal and air haw~ been used in the appropriate regions. The properties of a typical aluminium alloy have been used for the full region, air properties have been used for the empty regions and steel properties for the mould regions. These are shown in Fig. 3 The properties used for the Navier Stokes solver correspond to an Re of approximately 20.0. When this pr~,blem was solved for these properties, using GFEM and the
r~
J
i r
~
II................... r............•
........
©
c ~
,H
,H ~i~i
~'
~~ ~ii~i!!~ii!i!ii!i!ii!iii~ ~....
i
~iiiii~!~!i!i~i~i~i!~
iiii~i~i~i~i~
~
~ii~i~
iiii~i~i
H
298
,1 ....
-
I!.tJ
T= 70~C I"
o.o2m ~'-r
•
o.o3m
.
!iiii iTiiiiTiiii iiiiiiii771 | ................
li.:: T = 25~C
~,~.~
!
:.i:i:iiiit
~
li::::!i:::ii:iT:il
:: " , ":::: .............. :::Mould::
:
~
i
,
k=200mJ~-,i~
_.
~L
i:!
i •
. ::
L-:
=900 kglK
,{
"
li
P = 2700 ~ ' l
~
o.o2m
:
'lP :: 28' C
.
i i MOuldi.i
'::::~i2~5~Ciii~
+
:: /"
:
' ) ) T.=i25~iCi i i "-
.lm
Air To= 25~' C
.. . . .. . .. .. . . .. . .. . li:i: -
:i
.
i
.:
.
i
.
.
.
.
.
.
.
p=LO
iiliL
- ............
.
5-.
.
.
. :: ~:-7soo!~:! 5
_ :: e a-:450 i:~ c":
.
lt2g 71'~,--
kg m3
r: _.
k = .025mJY-,, K .
.
.
.
.
.
.
.
.
.
.
_
........
i
:
. . . . .
/.,
c = lO00kg~
-.-i =•.
.
.
.
. . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. .
:
•
-
0.0
F i g u r e 3: C o n f i g u r a t i o n for vertical filling with coupled heat transfer.
299
i
/ •
t ~k,S
' • ~
Min. value =
25.0000 Interv~ = 3 2 , 1 4 2 9
Max. value= 700.0000
Time =
0.02
MirLv , a u c = 2 5 . 0 0 0 0 Max. value = 7 0 0 . 0 0 0 0
ln~r,,al=
32.1429
Time = 0 . 0 6
r
~i ~'~ --1 Min. value= 25.12000 I n t e ~ = 32.1429
Min. value=
Ms,x. value= 700.0000 Time= 0.08
M ~ . value = 7 0 0 . 0 0 0 0
2 5 . ( X ) 0 0 interval = 3 2 . 1 4 2 9 Time = 0 . 1 0
ti
;I ~I
i
k~C
Min. v N u e = 2 5 . 0 0 0 0 M~.value=
700.0000
Intern = 32.1429
Min. vMue =
Time= 0.12
Max. v a l u e =
2 5 . 0 0 0 0 Interva.l = 3 2 . 1 4 2 9 7130.0000 T i m e = 0 . 1 4
F i g u r e 4: Explicit solution for vertical filling with coupled heat transfer: isotherms.
©
©
11
o
~z
~o
Q
J
~i~i~i~~ ~ ~
o
~:~ ~
~'~~ "
....i ......
~'
301
implicit backward Euler (BE) timestepping scheme, the resulls were riddled with oscillations, in addition to being overdiffuse. Results obtained fl'om the explicit Taylor Galerkin analysis, are shown in Fig. 4. It is clear from the isotherm plots that a sharp thermal front is advected faithfully without oscillation or excessive dissipation. Being explicit this solution was, however, expensively obtained from 140 timesteps. The timestep size was 0.001 corresponding to a Courant number (C) of roughly 0.5.
6
S i m u l a t i o n of solidification during filling
To simulate solidification during filling several additional problems must be addressed. The first and most important is the introduction of an algorithm to model the latent heat evolution. The second problem is of simulating the freezing of the fluid, which has been achieved by arresting the flow at the boundaries and increasing the viscosity of the fluid away from the boundary at the freezing point. A preliminary example of this procedure appears in Fig. 5.
7
Conclusions
From the number of problems solved here, it can be concluded that a finite element flow and heat transfer code such as that described in this work, when used with a front tracking algorithm such as the VOF method or the pseudo-concentration method based on the pure adw.'ction equation, provides an effective means ot' soving problems such as that of mould filling. This method can be used in other applications as well, such as forming processes involving creeping viscous flows, multiphase flow etc. The main advantage of this method is that the moving front does not have to be explicitly determined, neither do the boundary conditions on the front have to be explicitly satisfied. Furthermore, it is necessary to use special techniques, such as the Taylor Galerkin nmthod to deal with the problems arising from transporting a scalar when advection is the dominant mode. This is particularly important for obtaining meaningful temperature fields associated with filling.
A c k n o w l e d g e m e n t : The authors wish to thank the SERC who supported this work via ACME grant no. G R / F 32790.
References [1] C.W.Hirt and B.D.Nichols. Volume of fluid (vof) inethod for the dynamics of free boundaries. Journal of Computational Phgsics , 39:201 225, 1981.
302 [2] E.Thompson. Use of pseudo-concentrations to follow creeping viscous flows during transient analysis. International Journal for Numerical Methods i7~.Fluids, 6:749 7(i l, 1986. [3] T.J.Smith and D.B.Welbourn. The integration of geometric modelling with iir~ii~'(ei~ ment analysis for tile computer-aided design of castings. Applied Scientific t~i,,~archo 44:139- 160, 1987. [4] R.W.Lewis, K.Morgan, E.D.L.Pugh, and T.,I.Smith, A now on discontinuous ~m~m::r ical solutions of the kinematic wave equation. Inlernatio~al Journal for Vum,:'~-,',l Methods in Engineering, 20:555 563, 198'1. [5] B.Minaie, K.A.Stelson, and V.R.Volter. Fluid flow and solidification modellim< o~ die castings. ASME, Modelling of Materials Processing, MD, 3:35 50, 1987. [6] H.Sehlichting. Boundary-Layer 7'heorg. McGraw-lfill Book Company, i968. [7] D.K.Edwards, V.E.Denny, and A.F.Mills. 1)'ansfer Proc~.~,~es - A~z h~trodn,~o, to Diffusion, ConvectioT~ and Radiation. McGraw-Hill Book ('.ompany, 197:1. [8] P.M.Gresho, R.b.Lee, and R..L.Sani. On the time-dependent solution ~,1 the ~,~ou~ pressible Navier-Stokes equations in two and three dimensions. In Rece~d Ad~a~c¢,~ in Numerical Methods in Fluids, volume I. Pineridge Press l&nited. Swansea I;180. [9] T..l.R.Hughes. The Finite Element Method - Linear Static and Dynamic /:i,lt~ t;l~ment Analysis, Prentice-Hall international, Inc., Englewood Cliffs, New ,Jersey 07632, 1987. [10] D.M.Gao, G.Dhatt, 3.Belanger, and A,B.Cheikh. A finite element simulation of met.ai flow in moulds. [n Sixth International Co~@rence for Numerical Methods in 7 ht:~"mal Problems, Swansea, U.K., July 1989. Pineridge Press, Swansea. [11] A.S.Usmani, 3.T.Cross, and R.W.Lewis. A finite element model h)r the simuia!i(m of mould filling in metal casting and the associated heat transl~r. International ,t,-~t~l~taJ for Numerical Methods in Engineering, to be published in 1992. [12] G.F.Pinder and W.G.Gray. Finite h;lement Simulation i7~ Surface aT~d ,5'ub.~tlJ),:~ Hydrology. Academic Press, 1977. [13] M.S.Engelman, R.L,Sani, and P.M.Gresho. The implementation of normal ami/t~I tangential boundary conditions in finite element codes ior incompressible ttuid ttow International dournal for Numerical Methods in Fluids: 2:225 238, 1982. [14] J.Donea. A Taylor-Galerkin method for convective transport problems, lnter~ational Journal for Numerical Methods in Engineering, 20:1(11 119, t984. [15] K.W.Morton and A.K.Parrott. Generalized Galerkin methods for first order ~!yl)er bolic equations. Journal of Computational Physics, 36:249 270, 1980. [16] P.J.Roache. Computational Fluid Mechanics. U.S.A., 1976.
Hermosa Publishers, Albuquerque,