1 3 D Numerical Modelling of Mould Filling of a Coat Hanger Distributer and Rectangular Cavity

Filling processes occur in a wide range of industries, ranging from packaging of consumer products to manufacturing processes for making polymeric, metal and ceramic components. These processes involve the complex interplay of extrusion of a viscous liquid into a mould or container where it displaces a gas phase. Numerical modelling based on computational fluid dynamics can be useful for understanding the filling process. However, complexity arises in that the fluid dynamics in both the viscous liquid and gas phase must be resolved while concurrently determining the location of the fluid-gas interface and the interaction of this interface with the solid surface, i.e., the wetting behaviour. Determining the free surface location and wetting behaviour is an integral part of the numerical method.


Introduction
Filling processes occur in a wide range of industries, ranging from packaging of consumer products to manufacturing processes for making polymeric, metal and ceramic components.These processes involve the complex interplay of extrusion of a viscous liquid into a mould or container where it displaces a gas phase.Numerical modelling based on computational fluid dynamics can be useful for understanding the filling process.However, complexity arises in that the fluid dynamics in both the viscous liquid and gas phase must be resolved while concurrently determining the location of the fluid-gas interface and the interaction of this interface with the solid surface, i.e., the wetting behaviour.Determining the free surface location and wetting behaviour is an integral part of the numerical method.
Numerical methods have been applied to bottle and container filling for consumer products where the rheology can include shear thinning and viscoelastic effects and instabilities such as buckling and coiling may be prevalent [Tome, et al., 2001;Oishi et al., 2008;Roberts & Rao, 2011;Ville et al., 2011].In metal casting simulations, the fluids are generally Newtonian, but complexity arises from the high injection rates leading to turbulent flow and temperature-dependent behaviour such as solidification [Ilinca & Hetu, 2000;Cross et al, 2006].For injection moulding of polymers, time-and temperature-dependent effects are seen in conjunction with non-Newtonian rheology [Ilinca & Hetu, 2001;Kumar & Ghoshdastidar, 2002].In powder injection moulding for ceramic and metal forming, a suspension of particles is injected into a mould to create a green part, which later sees further processing steps to produce the final part [Hwang & Kwon, 2002;Ilinca & Hetu, 2008].Numerical methods for these problems range from finite difference, to finite volume, and finite element.
General classes of algorithms for determining the location of the free surface include Eulerian, Lagrangian and arbitrary Lagrangian-Eulerian (ALE) descriptions.Eulerian methods use a fixed-grid with an interface capturing technique such as the volume of fluid [Hirt & Nichols, 1981] or level-set method [Sethian, 1999] to determine the location of the free surface.Traditional Lagrangian methods use a moving mesh as a material interface that advects with the fluid.These methods often require multiple remeshing steps to avoid mesh distortion and tangling [see for instance, Bach & Hassager, 1985;R. Radovitzky & M. Ortiz, www.intechopen.comNumerical Modelling 4 1998; Zhang & Khayat, 2001].To avoid these problems, Lagrangian mesh-free methods have been developed using smooth particle hydrodynamics or the material point method to avoid meshing issues [Kulasegaram et al, 2006;Kauzlaric et a;., 2011;Love & Sulsky, 2006].These methods work well for being able to capture a moving interface with topological changes, but have difficulty in accurately solving the base physics, e.g.viscosity, and applying boundary conditions such as surface tension.ALE methods are hybrid techniques that seek to exploit the benefits of both the Eulerian and Lagrangian description in a hybrid manner, to determine the location of the interface [i.e.Sackinger et al, 1996;Lewis et al., 1998;Nithiarasu, 2005].
Injection loading of a ceramic paste is a high-rate process used to create green ceramic parts, which subsequently experience binder burnout and sintering to produce the final ceramic part.In the process of interest, the mould is a rectangular cavity, with an inflow from a coat hanger die that should distribute the flow evenly across the mould inflow.The cavity is small with dimensions of 1.3 cm by 3.6 cm in plane and a height of 0.4 cm.The inflow tube to the distributer has a diameter of 0.5cm.Figure 1 shows short shots (or incomplete filling of a mould) for the injection loading process, illustrating some of the defects that can occur when a fluid with complex shear and temperature-dependent rheology interacts with a high rate process.Fig. 1.Short shots of injection loading of a ceramic paste in a part are shown [Rao et al, 2006].The photo on the left is injection loading using a slow injection speed, while the photo on the right uses a higher injection speed.At low filling speed, the paste acts like a solid material.Even at high filling rates, when the paste begins to act as a fluid, pooling at the centre of the mould is seen and the desired mould shape is not achieved.
Temperature control and high-rate processing can limit the folding instability seen on the left of Figure 1.However, the pooling phenomenon shown on the right of Figure 1 occurs when optimal processing conditions are used, indicating that the design of the distributer may be responsible for material building up in the centre of the die.The complex shear-rate and temperature-dependent rheology of the ceramic paste was determined to follow a power-law dependence on shear rate and a Williams-Landau-Ferry temperature dependence [Rao et al., 2006].The material shear thinned quickly to a constant viscosity value at moderate shear rates.Thus, because of the high shear rates in the injection loader, it was determined that the rheology was essentially Newtonian as long as the temperature remained constant at the processing temperature.Therefore, a study was undertaken to better understand the filling dynamics and reduce pooling by changing the distributer design using a Newtonian fluid.

Dynamic wetting models and mould filling
Understanding dynamic wetting, or the interaction between the free surface and the mould walls, has been the subject of numerous experimental and theoretical studies and is still an outstanding research topic [see for instance Blake, 2006;Ren et al, 2011].The difficulty arises from the contradiction of a moving contact line at the fluid-gas interface and the no slip boundary conditions traditionally applied at solid surfaces.How can the contact line advance when the velocity vanishes at the solid surface?Highly viscous materials such as polymers and particle suspensions have large capillary numbers (a measure of the ratio of viscous forces to capillary forces) and are often hypothesized to obey a rolling motion condition with an 180 o dynamic contact angle.Numerically, this approach is difficult to apply and other methods have been proposed.The simplest models use a Navier slip condition to allow either slip on the entire solid surface, slip for the gas phase only, or slip only at the dynamic contact line.These models are ad hoc and ignore any thermodynamics considerations such as the static contact angle and surface energy.More advanced wetting models generally give dynamic contact angle as a function of the local Ca, the static contact angle and other material properties of the fluid and the solid surface [see Schunk et al., 2006 for a brief review of this work].Hoffman [1975] used experimental measurements to develop a universal correlation for the dynamic contact angle as a function of static contact angle and a local capillary number, while Cox [1986] developed a competing model using asymptotic analysis.Kistler [1983] used a linear model that was easy to implement in numerical computer codes.Shikhmurzaev [1994] used hydrodynamic theory and included a surface phase as part of the wetting model.Blake [1969] developed a molecular kinetic theory that reduces to a linear model for small contact angles.Blake noted that the advancing angle is a monotonically increasing function of Ca.The degree of velocity dependence will however increase steeply as viscosity increases or surface tension decreases.
To model dynamic contact for the filling process, the dynamic angle is tied to the balance of forces at the advancing wetting line, namely the tangential wetting line force, liquid-gas surface tension force, and fluid viscous force.Here, a version of the Blake model is used that is straightforward to populate with experiments.The dynamic contact angle is measured in the laboratory as a function of the velocity of the wetting line for the fluids and surface used in our experiments as input to the wetting model.The Blake model is also easy to implement numerically, since the wetting speed can be written as a function of the dynamic and static contact angles.The performance of the Blake model at the high capillary number limit may be suspect, since it exhibits an unbounded dynamic contact angle while more physical models such as Cox and Hoffman reach a limit of 180 o for large capillary numbers [Schunk et al, 2006].

Mould design
Because the initial design of the mould exhibited pooling of the fluid in the centre of the cavity and this would lead to poor filling (see Figure 1), we tried two minor redesigns to the distributor to see if we could improve the flow into the mould and reduce pooling.Ideally, we would like to see a flat profile coming out of the distributor and a more one-dimensional front shape.Figure 2 shows the original mesh, Mesh 1, and a variation, Mesh 2, with a longer distributor, and Mesh 3 with a longer-taller distributor.The idea behind Mesh 2 was to give a longer length for the flow to develop a flat profile and fill up the distributor before entering the main cavity.Mesh 3 kept this longer distributor and made it wider on inflow to ease the fluid entering the cavity.These ideas were inspired by discussions in Sartor [1990] about die design.The meshes themselves all have the same cavity size and a similar amount of refinement, though Mesh 2 and Mesh 3 have more elements and unknowns due to the longer distributor.
Mesh 1: Original Geometry Mesh 2: Longer Distributor Mesh 3: Longer-Taller Distributor Fig. 2. Geometries to be investigated: Mesh 1 is the original design, Mesh 2 incorporates a longer distributer, and Mesh 3 incorporates a longer distributer with a wider inflow to the cavity.

Chapter organization
In this chapter, we investigate the design of an injection loading mould using flow visualization experiments and numerical models.The finite element method is used to understand the interaction of the inflowing viscous liquid with the geometry and the displaced gas phase.Filling dynamics are determined with a diffuse-interface implementation of the level set method [Sethian 1999].The Blake wetting model is used to represent the interaction of the free surface with the mould surface at the dynamic contact line [Blake & Haynes, 1969;Blake, 2006].The flow visualization experiments are carried out under isothermal conditions using an acrylic mould and a viscous Newtonian fluid.Three different moulds are examined in two different orientations with gravity.Simulation results are given for these six cases.
The chapter is organized in the following manner.First, the equations and numerical method are presented.Next, the experimental methods used to provide input parameters for the models and flow visualization studies to better understand the filling dynamics and provide confidence in the numerical method are discussed.In the subsequent section, the results for injection moulding process are given for a Newtonian fluid into a coat hanger die distributer and a rectangular cavity, where the 3D level set simulations are compared to experiments.We conclude by summarizing the results and discussing future efforts.

Equations of motion
We can write the equations of motion for a single-phase fluid and then generalize them for our multiphase flow problem, where a viscous fluid displaces a gas.The fluids of interest are assumed to be incompressible and have a constant density, meaning that the velocity field, u, will be solenoidal and the continuity equation contains no density or pressure variables.(Note, this is a good assumption for the viscous liquid but a simplification for the gas phase, which is actually compressible.) Conservation of momentum for a Newtonian fluid takes into account gradients in the fluid stress tensor, T, defined as the product of the viscosity and the shear rate tensor, () t uu  , gradients in the pressure, p, as well as gravitational effects and inertial terms that can be dependent on time, t, and the fluid density, .Note that gravity, g, can be an important body force in filling processes and most filling processes fill counter to gravity.

()
Because we have a viscous fluid displacing a gas phase, the location of the interface between fluids is unknown a priori.To determine the location of the free surface as it evolves in time, we use an Eulerian interface-capturing scheme based on the level set method, the details of which are included in the following section.

Interface capturing
We use the level set method of Sethian [1999] to determine the evolution of the interface with time.The level set is a scalar distance function, the zero of which coincides with the free surface or fluid-gas interface, e.g.
We initialize this function to have a zero value at the fluid-gas interface, with negative distances residing in the fluid phase and positive distances in the gas phase.An advection equation is then used to determine the location of the interface over time.
Derivatives of the level set function can give us surface normal, n, and curvature, , at the interface, which is useful for applying boundary conditions.
We use the equations of motion described in the previous section, but vary the material properties across the phase interface.This variation is handled using a smooth Heaviside function that modulates material properties to account for the change in phase.

()
( 1 s i n ( ) ) ; 2 () 1 () This is a diffuse interface implementation of the level set method, which allows for an interfacial zone of length 2.This zone is usually chosen to be four to six elements wide.
Equation averaging is done using a Heaviside for the gas and viscous fluid equations.
Because the properties are linear, this process results in Heaviside-averaged properties in a single momentum equation and an unchanged continuity equation: The properties have fluid properties in some regions and gas properties in other regions as modulated by the numerical Heaviside.In the diffuse interface region, the properties are averaged between fluid and gas values.Figure 3 shows a schematic representation of the Heaviside function.The regularized Dirac delta function, which is defined as is used to apply surface tension and capillary boundary conditions via a continuous surface force approach [Brackbill et al, 1992].This applies surface tension as a volumetric body force on the momentum equation, which is distributed throughout the interfacial zone region through the regularized Dirac delta function.

Finite element implementation
The equations of motion (7) and the level set advection (4) were solved using a finite element method as implemented in ARIA [Notz et al., 2006].Bilinear shape functions were used for the three velocity components, pressure, and level set.The LBB requirement on the velocity and pressure space was circumvented using Dohrmann-Bochev pressure stabilized pressure-projection (PSPP) [Dohrmann and Bochev, 2004] to allow for this equal order, bilinear, interpolation of all variables.(LBB compliant elements have the velocity space higher than the pressure space [Hughes, 2000].)The velocity vector and pressure unknowns were solved in the same matrix, while the level set equation was solved in a separate matrix, but at the same time step intervals.The level set equation was stabilized using a Taylor-Galerkin method [Donea, 1985].The PSPP stabilization method greatly improved the condition number of the discretized matrix equations when compared to LBB elements or other stabilization methods, allowing for the use of an ILU preconditioner with a BiCGStab Krylov iterative solver.Further details of the modelling approach and equations, the numerical methods used and the finite element implementation can be found in Rao et al. [2011].

Contact-line wetting model
Boundary conditions for the dynamic contact line where the free surface and wall intersects are handled with a Blake wetting condition [Blake and Haynes, 1969;Blake, 2006].
Parameters for the model are informed by goniometer experiments that determine the wetting speed, v wet as a function of the dynamic contact angle, θ, and the static contact angle, θ s for the various fluids and surfaces of interest [Mondy et al., 2007].
sinh( (cos cos )) The dynamic contact angle can be calculated from the level set function.
When we integrate the stress, T, by parts in the finite element implementation of the momentum equation, a surface term is created as a natural boundary condition.We exploit the surface stress term and add on a Navier slip condition that includes the wetting speed from the Blake model.
The value of the tangential wall velocity ramps from zero at a level set length scale away from the contact line to v wet from equation ( 11) at the contact line.Away from the contact line, we revert to no slip for the tangential wall velocity boundary condition.The normal velocity is enforced as no penetration everywhere.The shape of this ramp must be smooth in order to get a realistic wetting line that shows a smooth transition from the contact line to the bulk flow.For a sharp transition from slip to Blake, we ended up with unphysical looking cusps in the interface shape near the wall.The transient terms introduce dynamics into the wetting line motion and allow smooth movement of the contact line.The  parameter is taken to be approximately the time step.The  parameter is generally small so that equation ( 13) functions almost like a Dirichlet requirement on the fluid velocity.

Experiments
It was decided to visually record the flow of a simple, single phase, Newtonian liquid through transparent moulds to build confidence in the front tracking and wetting models used in the computations.Details of the geometries, experimental conditions, and properties of the materials used in each test will be given in the following sections.

Geometries
For these tests, we built transparent acrylic moulds identical to the three mesh geometries (Figure 2).The strength of the acrylic material making up the transparent moulds dictated that we inject at a lower injection pressure and lower operating temperature than the actual injection loading process.To mimic the actual injection process, we used a pressure driven syringe held at a constant pressure of 29.95 ± 0.10 psig during injection.The syringe of liquid was degassed in a vacuum chamber prior to the experiment.The syringe was modified to prevent leakage around the plunger and subsequent bubble formation.However, this resulted in more friction and a flow rate t h a t v a r i e d s o m e w h a t f r o m e x p e r i m e n t t o experiment.Hence, the time to fill the moulds varied from test to test, but was determined and recorded for each test.Reported below is the median and spread of the measured filling time for four repeated experiments in each geometry.
The tests were conducted at a room temperature that ranged from 23 to 24 o C. To minimize the effect of temperature on the viscosity, the liquid was held in a water bath set to 23.5 o C after degassing and before being used in an experiment.These conditions resulted in an injection rate that was approximately ten times slower than the actual process.However, the Reynolds number for the actual process and the validation experiments were similar and in the Stokes regime of much smaller than one.For an average fill time of 20 s, a cavity volume of 1.87 cm 3 , and characteristic length scale from the inflow port of 0.5cm, the Reynolds number was 0.0006 and the capillary number was 4.5.

Materials and properties
We chose a liquid, UCON 75-H-90,000 (oxyethylene/oxypropylenes from Dow Chemical) as our model fluid.This UCON was chosen since its viscosity is an order of magnitude lower than the ceramic paste at processing conditions.Combined with an order of magnitude decrease in the injection rate compared to the real system, this gives a similar Reynolds number.The liquid surface tension and wetting properties of the UCON lubricant were determined at 23.5°C.The viscosity of the lubricant was measured with a Rheologica™ constant stress rheometer and was equal to 390 Poise over a shear rate ranging from 0.1 to 10 sec -1 , indicating Newtonian rheological behaviour.The density of the liquid was measured with a densitometer to be 1.09 g/cm 3 .The surface tension measured with a Du Noüy ring (mean circumference of 5.935 cm) was 42.40.1 dyne/cm.
The dynamic contact angle on acrylic was measured with a feed-through goniometer [Mondy et al., 2007], in which liquid can be continuously injected to achieve "high" velocities, or the sessile drop can be allowed to relax to obtain "low" velocities.Figure 4 shows a schematic of the experimental apparatus and the results of this wetting test.During each experiment, the angle of contact and the location of the triple point of contact was recorded.The wetting line speed was then determined from the location of the triple point with time.The data fit to the Blake model (equation 9) gives a wetting constant v o of 0.00130cm/s and scale factor  of 2.29 with a static contact angle θ of 37.3°.

Flow tests
Figures 5 through 7 are representative frames from the video recordings of the fill process using a vertical alignment of the mould as in the proposed ceramic injection process.The time it took to completely fill the original geometry, Mesh 1, (defined by the time when the liquid began to exit along the entire length of the vent) was 24.6 ± 1.2 s, to fill the geometry of Mesh 2 was 26.9 ± 2.0 s, and to fill the geometry of Mesh 3 was 24.6 ± 2.3 s.Because the fill rates varied, the time is shown in these figures in a nondimensional form.The initial time is taken to be when the front passes a line in the square entry channel 0.16 cm from the entrance of the distributor.Here, one can see the effects of changing the distributor geometry.In Figure 5 the liquid has just filled the distributor.All of the modified geometries help flatten the leading front, especially Mesh 2. Mesh 2 fills the distributor with the fastest relative time.

Mesh 3
Time/total time=0.24 Figure 6 shows the times at which the leading front hits the wall opposite the injection port.
In this case, the front in Mesh 2 takes the longest relative time to reach this stage.The other geometries once again follow the same pattern as before, with Mesh 2 giving the flattest flow profile and the original mesh the displaying the most curvature of the interface.At this point, the fluid has wetted more of the side walls with the Mesh 2 geometry and the front is flatter.Because Mesh 2 has a flatter profile, it takes the longest time to reach the top wall compared to Mesh 1 and Mesh 3.

Mesh 3
Time/total time=0.82Fig. 6.Comparison of the effect of distributor geometry on the time it takes to reach the wall farthest from the injection port for vertical mould orientation.
Figure 7 shows the locations of the voids remaining in the mould once filled but without any over pressure (a relative time of 1).Small bubbles away from the corners are artefacts of the syringe loading process.The voids in geometries with redesigned distributors remain in the same locations as those seen in the original geometry.The relative areas of the bubbles on the images, which reflect the volume of air trapped, were determined and compared quantitatively in Table 1.One pixel resolution represents about 1×10 -5 cm 2 .All redesigned distributors result in smaller bubbles in the upper corners than those in the original mesh.
The lower bubbles of Mesh 2 are also smaller, whereas those in Mesh 3 are approximately the same as those in the original geometry.

Mesh 1 Mesh 2 Mesh 3
Fig. 7. Voids remain in the front upper and lower corners of each geometry for the vertical mould orientation.Each frame consists of two views of the mould (top view and side view).These voids can be seen more easily from the side view.Next, we experimented with the orientation of the mould.Moulds were turned so that the distributor was perpendicular to gravity on the lower surface and the vent was moved to be on the upper surface.In other words, the flow direction and gravity are now perpendicular, with gravity acting in the thinnest cavity direction.Results for the original mesh showed that orientation with respect to gravity had a large impact on the likelihood of voids remaining in the corners of the mould.Figure 8 shows representative video frames at three stages of the fill: 1) when the distributor was completely filled, 2) when the front hit the far wall, and 3) when the liquid began to exit through the vent and the sides had completely wetted.Because the fill rates varied, the time is shown in a nondimensional form.The time it took to completely fill the original geometry in the horizontal position was 23.3 ± 0.9 s, to fill the geometry of Mesh 2 was 28.9 ± 1.5 s, and to fill the geometry of Mesh 3 was 27.3 ± 1.8 s.
The top views of the horizontal orientation show that the liquid wets the top later than the bottom.The oval area of the middle image in Figure 8 indicates where the liquid has wetted the upper surface.The side views are more difficult to interpret because of the lighting challenges accompanying trying to see through the thickest section of the mould.However, the dark areas of the side view are where the liquid has wetted the sides.
It is interesting to note that the horizontal alignment causes the front to hit the far wall before even wetting the sides at all.In other words, the front was less "flat" entering the mould from the distributor.Nevertheless, the results also show that the horizontal orientation, with the vent on top and the distributor entrance on bottom, resulted in bubbles only in the upper corners nearest the distributor.When oriented vertically, bubbles are trapped in corners both opposite the distributor and opposite the vent.The bubbles nearest the distributor with the horizontal orientation are approximately the same size as the bubbles trapped in the upper corners opposite the vent in the vertical orientation.
The effect of distributor geometry on the filling in the horizontal orientation is shown in Figure 9 and Figure 10, which correspond to Figure 5 and Figure 6 in the vertical orientation.Again, the geometry of Mesh 2 results in the flattest front entering the mould from the distributor (Figure 9).The front reaches the back wall at approximately the same time using the distributors of Mesh 2 and 3.By the time the front reaches the back wall both modified geometries result in more liquid filling the mould than with the original geometry; however, Mesh 2 results in somewhat more than Mesh 3 (Figure 10).Sizes of the bubbles left in the various geometries in the horizontal orientation are listed in Table 2. Comparing these values with the measurements in Table 1, one can see that the amount of gas left in the vertical orientation is roughly the same as that in the horizontal orientation, although in the horizontal orientation there are fewer bubbles.Bubbles observed in other locations than the corners are almost always artefacts of the syringe loading process.

Simulations
T h e s i m u l a t i o n s r u n s w e r e m a d e o n 6 4 p r ocessors of Thunderbird (Sandia National Laboratories capacity computing platform) and ran in less than 3.5 hours.This allowed us to do real time design and sensitivity calculations for parameters such as wetting speed, second phase viscosity, level set length scale and inflow pressure.Properties for the liquid used for the validation simulations were the measured values discussed in the experimental section above.The density and viscosity of the displaced gas phase were taken as fictitious values of one thousand times smaller than the liquid phase density and viscosity.These values are summarized in Table 3.At 25 o C, the viscosity of air is 2.0x10 -4 Poise and the density of air is 0.0012 g/cm 3 , thus our second phase properties are very close for density, but three orders of magnitude too high for viscosity.The numerical method fails to converge for values of the liquid/gas viscosity ratio of more than 1000 for a diffuse interface implementation of the level set equations, so this is a necessary expedient requiring us to adhere to this fictitious viscosity value.

Material Property Value
For the inflow condition, we used a constant inflow pressure of 1.0x10 6 dyne/cm 2 .This value was chosen to match the horizontal fill time of 23 seconds in the original mesh and then used for all other meshes and geometries.A shooting method was used, where different values of the inflow pressure were used and the solution was examined to see if it filled in the correct time.This required many simulations to be run.It is believed that the actual boundary condition for the experiments is somewhere between a constant velocity and constant pressure condition, but this is hard to replicate numerically.From the simulations, we found that the velocity changed quickly in the beginning for the pressure inflow boundary condition and subsequently reached a steady value.
The initial 3D mesh and boundary conditions are given Figure 11 for the mould filling simulations.We assume symmetry about the centreline and only solve half the problem to improve the computational efficiency.The mesh contains 6744 8-Node hexahedral elements giving 41300 total degrees of freedom for bilinear velocity/bilinear pressure interpolation.This mesh was shown to be adequate, as a more refined version of this mesh gave the same fill times and meniscus shapes [Rao et al., 2006].The time to fill the mould for the vertical orientation for Mesh 1, Mesh 2, and Mesh 3 are 15.2s, 17.5s, and 13.2s, which was much faster than the experimental values of 24.6s, 26.9s, and 24.6s.We defined our fill time as when the vent had filled completely to a distance of the half the level-set length-scale or two elements.Unfortunately, the numerical fill times are difficult to obtain accurately as the gas phase viscosity makes a large difference in how fast a simulation will fill for the same value of pressure.For instance, when we reduced the second phase viscosity by a factor of 10 we got an increase in the inflow velocity when keeping all other parameters constant.Thus, the fact that the gas phase is harder to push out than it is for the experiments, adds a great deal of uncertainty to the fill times.However, we hope to be able to predict trends.Figure 12 shows the effect of the distributor geometry on the meniscus shape for Mesh 1, Mesh 2, and Mesh 3 in a vertical orientation.Comparing Figure 12 and Figure 5, the numerical and experimental version of this profile, we can see that the simulations are exhibiting the physically correct trends.The original mesh takes the longest fractional time to fill the distributor, 42%, and gives the most bulging front shape.Mesh 2 is an improvement, taking 18% of the time to fill the distributor, while Mesh 3 is somewhere in between at 24%.The values for the experimental distributor dimensionless fill times are 32%, 13%, and 24%, so the simulations are also capturing the correct trends for fill time though they are not quantitative.The shape of the meniscus for Mesh 1 has more of a bulge at the edge of the distributor than the experimental meniscus, which looks as if the front is pinned at the distributor.
We can also look at the profiles and dimensionless time to hit the back wall.These results are given in Figure 13.Comparing Figure 13 to Figure 6, for the simulation versus experiment, we can see differences in the meniscus shape.The numerical interface reaches the back wall for Mesh 1 and Mesh 3, before it wets the sidewall and Mesh 2 has a flatter profile in the experiments than the simulations.The percentage time to reach the back wall for the simulations on Mesh 1, Mesh 2, and Mesh 3 are 60%, 70%, and 55% compared to 80%, 86% and 82% for the experiments.Again, we capture the correct trends, but are still unable to match the data quantitatively.Figure 14 shows the full meniscus shape and void locations for the simulations on Mesh 1, Mesh 2, and Mesh 3. The void in the corner near the distributor does eventually fill in, since its size is less than the level-set length scale and we are using a diffuse interface method.The larger void at the vent never fills in as the viscous gas phase is trapped away from the vent by the fluid.For the numerical solutions it is hard to make any predictions about void size, though we can say that for similar values of the dimensionless time the voids for Mesh 2 will be smaller than Mesh 1, with Mesh 3 being somewhere in between, which does follow the experimental trend.
We also examined the horizontal orientation numerically, which gave fill times for Mesh 1, Mesh 2, and Mesh 3 of 21.4s, 23.1s, and 22.2s compared to 23.3s, 28.9s, and 27.3s for the experiments.Again, we follow the trends of the experiment, but do not match quantitatively.Mesh 2 seems to take a longer time to fill for the experiment than one would predict numerically.Figure 15 shows a comparison of the profiles for Mesh 1 in a vertical and horizontal orientation.Comparing Figure 15 to the experimental equivalent, Figure 8, we can see that we have quantitative differences but do match some trends.The numerical meniscus shape leaving the distributor looks similar to the experimental profile as it bulges more in the centre for the horizontal orientation, though the simulation is less dramatic.The numerical profiles when the fluid first hits the back wall are flatter for the vertical orientation than the horizontal, though the vertical should be even flatter to match the data.The numerical solutions predict two voids for each orientation, though the experiments do not show a second void for the horizontal orientation near the vent.However, this void may just be difficult to see experimentally.Conversely, the horizontal void at the outflow may be an artefact of the numerical method as we have a difficult balance at the outflow between wetting forces, gravity, the gas phase viscosity, and the material flowing out the vent.Also, our numerical vent is not identical to the experimental one and exhibits a slightly different area and shape.

Mesh1
Figure 16 shows the meniscus profiles for Mesh 1, Mesh 2, and Mesh 3 after filling the distributor for the horizontal orientation.Comparing Figure 16 and Figure 9, the numerical and experimental version of this profile, we can see that the simulations are again exhibiting the correct trends of the physical situation.Mesh 1 shows the most pooling at the centre of the mould, Mesh 2 has a flatter profile as does Mesh 3. The experiments predict filling times for Mesh 3 to be in between Mesh 1 and Mesh 2, and the simulations follow this trend.The dimensionless times to fill the distributor numerically for Mesh 1, Mesh 2, and Mesh 3 are 29%, 15%, and 21% compared to 26%, 13% and 22% for the experiments.In general, the simulations predict the vertical filling to be faster overall than the horizontal by several seconds for each of the geometries, whereas the experiments are faster in the vertical for Mesh 2 and 3, but slower for Mesh 1.This could have resulted from some experimental errors or from the uncertainty in injection rates and flow profiles from the experiments to the simulations.
Figure 17 shows the free surface profile as the fluid hits the back wall for Mesh 1, Mesh 2, and Mesh 3 in the horizontal orientation.The dimensionless times it takes to hit the back wall for Mesh 1, Mesh 2, and Mesh 3 are 41%, 55%, and 52% compared to values of 82%, 71%, and 70% for the experiments seen in Figure 10.For the simulations, Mesh 1 hits the back wall for the smallest dimensionless time while Mesh 2 and Mesh 3 take about the same time.For the experiments, Mesh 1 takes the longest time, while Mesh 2 and Mesh 3 do take about the same time.Thus for this case, we are capturing one trend, but not the differences between Mesh 1 and Mesh 2. From Table 4 we can see that we capture some of the correct trends, especially for the vertical orientation, though some features elude us like the time to fill to the back wall for Mesh 1 in the horizontal orientation.Sources of uncertainty in the simulations include: 1) Lack of clarity of inflow conditions from the experiment to the simulation, since it is somewhere in between constant pressure and constant velocity, 2) Possible poor performance of the Blake wetting model at moderate capillary number and large dynamic contact angles, 3) Possible poor performance of the wetting model at wetting speed higher than experiments used to populate the model.Sources of error in the simulations include: 1) Numerically expedient of high gas phase viscosity, 2) Lack of compressibility for the gas dynamics, 3) The use of a diffuse interface model that smears out material property jumps, allowing viscous bleed through of the liquid phase into the gas phase.For future work, we will try to reduce these uncertainties and errors by using some of the advanced features in ARIA, which should be available soon, to allow for smaller values of the gas phase viscosity, such as sharp integration and a compressible gas phase.The optimal choice of wetting model for moderate to high capillary numbers continues to be an ongoing focus of our research.

Conclusion
A diffuse interface finite element/level-set algorithm has been used to investigate filling behaviour for injection loading using a Blake wetting model.The modelling has been successful in matching experimental data qualitatively, but quantitative agreement is still lacking especially for the wetting dynamics and meniscus shape.For future work, we will investigate an advanced version of the level set method termed the conformal decomposition finite element method (CDFEM).CDFEM is a hybrid moving boundary algorithm, which uses a level set field to determine the location of the fluid-fluid interface and then dynamically adds mesh on the interface to facilitate the resolution of discontinuous material properties and fields, as well as the application of boundary conditions such as capillarity.This is a sharp interface method, where it is possible to apply jumps in material properties, material models, and field variables [Noble et al, 2010].We believe this algorithm, which should be available soon, will lead to better agreement with experiments and should allow for straightforward inclusion of a compressible gas phase.

Acknowledgments
This work was supported by the Laboratory Directed Research and Development program at Sandia National Laboratories.We would like to thank Dr. Pin Yang, the project principal investigator, for his support of this study.We would also like to thank our Sandia reviewer Daniel Guildenbecher and P. Randall Schunk for their insightful editorial comments.
Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.

Fig. 5 .
Fig. 5. Comparison of the effect of distributor geometry on the shape of the fluid front entering the mould for vertical mould orientation.

Fig. 8 .Fig. 9 .
Fig. 8. Original geometry (Mesh 1) oriented horizontally (A) and vertically (B).Images on the left show when the distributor fills completely, middle images show when the front first hits the back wall, and images on the right show when the part is filled to the point that the fluid leaves the vent area through the entire length of the vent.Bubbles left in the corners are circled in red.

Fig. 11 .
Fig. 11.Initial mesh and boundary conditions for 3D level set simulations.The outflow vent is on the same side as the distributor for vertical simulations and opposite the distributor for horizontal simulations.

Fig. 12 .
Fig. 12. Free surface profile after filling the distributor for Mesh 1, Mesh 2, and Mesh 3 for vertical mould orientation.
Fig. 13.Free surface profile after hitting the back wall for Mesh 1, Mesh 2, and Mesh 3 for vertical mould orientation.
. Mesh 1 oriented horizontally (A) and vertically (B).Leftmost pictures show profiles when the distributor is filled, middle shows profile when the fluid hits the back wall, and rightmost pictures show final profile.Both front and side views are given to highlight void location.A B www.intechopen.com Fig. 16.Free surface profile after filling the distributor for Mesh 1, Mesh 2, and Mesh 3 for horizontal mould orientation.
Fig. 17.Free surface profile after hitting the back wall for Mesh 1, Mesh 2, and Mesh 3 for horizontal mould orientation.

Table 1 .
Bubble sizes remaining in the corners for vertical mould orientation.

Table 3 .
Material properties used for validation simulations.

Table 4
summarizes the fill times to reach the distributor, back wall, and complete mould for the simulations and experiments for vertical and horizontal orientations on all three meshes.

Table 4 .
Summary of fill times for experiment and simulations to reach the distributor, the wall and completely full.