1 Introduction
Among the various processes used by the industry to produce a finished part, die casting is an important manufacturing process in which liquid metal is injected into a die and solidified. Usually die casting is used for aluminum and magnesium alloys with steel molds. Die cast products are more commonly used in automotive and housing industries. Several complex processes involving a large number of process parameters affect the final product quality. Due to the recent advances in computing hardware and software, it is now possible to simulate the physics of these processes using numerical simulations. These simulations provide detailed flow and temperature histories and can be used to estimate the final product quality. Frequently, it is difficult to measure and tightly control the process parameters like initial melt temperature, mold temperature and alloy component concentration. However, they can have a significant impact on the process as well as the predicted product strength.
Temperature evolution and velocity distribution during solidification of pure metals and metal alloys has been analyzed by many researchers. Bennon and Incropera (1987a) used a continuum mixture model to solve the momentum and energy equations and applied it to the solidification problem of a rectangular cavity filled with a binary aqueous solution Bennon and Incropera (1987b). From their work, it is clear that although the continuum formulation has a limitation of smearing the interface, it is efficient than the multiple region method in which the governing equations are solved in each phase separately with appropriate interface conditions. Implementing multiple region method for a practical complex geometry with irregular interfaces is quite difficult and computationally expensive. Voller and Prakash (1987)
used the Darcy’s law of drag and the enthalpy method to model the mushy zone during solidification of a square cavity. Two vertical walls were held at below and above the melting point and the horizontal walls were thermally insulated. The damping of velocities in the mushy zone due to Darcy’s law is clearly seen in the velocity vector plot.
Vynnycky and Kimura (2007) solved a similar problem on a rectangular enclosure analytically. They used an asymptotic approach to solve the nondimensional momentum and energy equations for the case of the Rayleigh and Stefan numbers much larger and smaller than unity, respectively. Comparison with finite element transient numerical simulation showed that such an asymptotic simplification is possible. Bennett (2007) used adaptive grid with local rectangular refinement (LRR) to model solidification with fluid flow in multiple rectangular geometries. The adaptive LRR grid needed half the storage and computational effort compared to the nonadaptive grid. Wang et al. (2010) discussed a numerical model for melting in a rectangular cavity with natural convection at high Rayleigh number (). They used the consistent update technique (CUT) algorithm with a multigrid solver. Plotkowski et al. (2015) simulated analytically and numerically heat transfer and fluid flow for solidification in a rectangular cavity. They used scaling analysis to simplify the governing equations of the mixture model followed by an analytical solution and comparison with a finite volume solution for an AlCu binary alloy. Hu et al. (2017) studied the threedimensional phase change problem with natural convection using the Lattice Boltzmann method. They simulated melting in a cubical cavity and in cavities with inner rectangular cylinders and sphere. Cleary et al. (2010) used smoothed particle hydrodynamics to model filling and solidification in high pressure die casting. They simulated practical industrial case studies like differential cover, an electronic housing and a door lock plate. Gau and Viskanta (1986) performed experiments to understand the importance of natural convection during solidification and melting of pure metal. Quillet et al. (2007) and Hachani et al. (2012) reported the temperature distribution and macrosegregation of solute with an experimental solidification of a Tin based alloy in an ingot. Such experimental work is useful for calibration and validation of the numerical methods.Solidification phenomena of die casting involves interplay between heat transfer and flow due to natural convection. All the references discussed above simulated solidification on rectangular or cuboidal geometries and hence, structured Cartesian grids were used. Practical die casting geometries are complex and thus, unstructured hexahedral elements are utilized in this work. Multiple algorithms are discussed in the literature to simulate fluid flow on unstructured grids Davidson (1996); Haselbacher (1999); Kim and Choi (2000); Ferziger and Peric (2012); Mathur and Murthy (1997); Muzaferija and Gosman (1997). Strategy implemented in this research is an extension of the work by Mathur and Murthy (1997) and Muzaferija and Gosman (1997). Solution of the fluid flow and the pressure Poisson equation involves a significant computational effort. As the solidification problem proceeds, the volume of liquid zone deceases continuously and it is not necessary to solve for fluid flow in the solid zone. We therefore, have developed a consistent procedure to modify the coefficients of the discretized momentum and pressure Poisson equations in order to satisfy the mass continuity equation. Algebraic multigrid Yang and Henson (2002); Ruge and Stüben (1987) and Krylov subspace solvers from the open source library HYPRE HYP (2017) are used to solve the linear systems of equations.
For die casting, the microstructure parameters like grain size and dendritic arm spacing are important as they affect the final product quality. Phase field modeling Provatas et al. (1999); Karma and Rappel (1998); Ramirez et al. (2004) is a popular method used to study the evolution of the microstructure during solidification. The phase field method simulates the growth of each dendrite and thus, it is computationally expensive at the length scale of die cast products. In this research, an empirical relation from the work of Backer and Wang (2007) is used to estimate secondary dendritic arm spacing. There are various models Maxwell and Hellawell (1975); Greer et al. (2000); Desnain et al. (1990) suggested for grain size estimation during solidification. Here, the isothermal crystal growth model Greer et al. (2000) is used.
Use of deterministic simulations alone to analyze the engineering systems is incomplete due to the lack of precisely defined input data. Thus, there has been a growing interest Xiu and Karniadakis (2002a, 2003); Maitre et al. (2001); Carnevale et al. (2013); Marepalli et al. (2014); Ganapathysubramanian and Zabaras (2007); Fezi and Krane (2017); Hosder et al. (2006); Kumar et al. (2016); Fajraoui et al. (2017) in coupling uncertainty propagation techniques with the deterministic numerical simulations to estimate the effects of stochastic variations in the input process parameters on the outputs. The polynomial chaos expansion is a popular method used to estimate the relation between input and output parameters. Stochastic Galerkin projection Xiu and Karniadakis (2002a, 2003); Maitre et al. (2001) and collocation Carnevale et al. (2013); Marepalli et al. (2014); Ganapathysubramanian and Zabaras (2007)
are two strategies to estimate the coefficients of the polynomial chaos expansion. Stochastic Galerkin method is an intrusive method since it requires solution of a new set of equations and thus, modification of the underlying deterministic code which becomes a significant additional effort. Hence, recently nonintrusive stochastic collocation methods have been developed which need multiple evaluations of the deterministic simulation at predefined collocation points obtained by sampling from the probability distribution function of the input parameters. Values of outputs estimated at these samples are then used to estimate the coefficients of the polynomial chaos expansion.
In this paper, the traditional fractional step approach is modified to deal with the additional terms in the NavierStokes equation due to solidification. The discretized system of equations is altered in a way so as to simulate only in the liquid zone thus, increasing computational efficiency. In order to incorporate the effects of stochastic variations in the input process parameters, a parameter uncertainty propagation module has been developed in conjunction with the deterministic simulations. The method of polynomial chaos expansion is used to estimate the relation between input and output parameters and stochastic collocation is used as a wrapper on the underlying deterministic simulation. The framework has been validated against published experimental results followed by demonstration for solidification of two complex geometries.
2 Numerical Model Description
2.1 Governing Equations
Solidification, heat transfer and fluid flow due to natural convection are modeled. It is assumed that there is no macrosegregation during solidification and the metal is solidified at nominal composition. It is further assumed that the solutes do not contribute to buoyancy and the flow is incompressible. Thus, the set of governing equations consists of the standard NavierStokes equations with additional terms for solidification Plotkowski et al. (2015). They are:
(1) 
(2) 
(3) 
Here, u is the velocity vector, is density, is time, is dynamic viscosity, g is gravity vector, is coefficient of thermal expansion, is pressure, is isotropic permeability of the dendritic array, is dendrite arm spacing and is solid fraction.
To model the effects of natural convection, the Boussinesq approximation is used. This is a valid assumption for problems with moderate density variations in the domain. The fluid is modeled as a constant density fluid except for the additional buoyancy term in the momentum equation (2) Spiegel and Veronis (1960).
The Darcy drag term () represents increased resistance to the flow in the mushy zone. We have used the BlakeKozeny model (Eq. (3)) which estimates the isotropic permeability () of the dendritic array. In the liquid region, solid fraction is zero and permeability tends to infinity making the Darcy drag term to go to zero. When solid fraction is unity, permeability tends to zero and thus the coefficient of Darcy drag term goes to infinity. For stability, this coefficient is added to the diagonal term of the discretized momentum equations. As a result, the velocities in the solid region go to zero. In the mushy zone, the drag term reduces the velocities compared to the liquid zone.
The energy equation is written in terms of temperature as:
(4) 
where,
(5) 
2.2 Solution Algorithm
We have developed a new software OpenCast in an object oriented C++ environment. A finite volume method on a collocated grid is used to discretize the governing equations. The fractional step method Harlow and Welch (1965), modified to account for the solidification is used to integrate the equations. First, the momentum equations (6) without the pressure gradient term is solved to estimate an intermediate velocity (u*) field. The buoyancy and Darcy drag terms are included in this step. Second order accurate CrankNicolson scheme for the diffusion terms and the AdamsBashforth scheme for the convection terms are used for temporal discretization. The coefficient of the drag term () in the mushy zone is treated fully implicitly as:
(6) 
where, the operators , and represent the discretized convection, diffusion and buoyancy terms respectively. The full momentum equation is similarly discretized with an implicit pressure gradient term, given as:
(7) 
Subtracting equation (6) from (7) gives the velocity correction equation.
(8) 
Taking divergence of the velocity correction equation (8) and invoking the continuity constraint gives the equation for pressure:
(9) 
The overall solution algorithm to advance from timestep to is as follows:

Solve for u* using equation (6). Since the diffusion term is implicit, solution is obtained iteratively

Solve the pressure Poisson equation (9) iteratively to estimate

Correct the velocities () using equation (8)

Estimate microstructure parameters such as grain size and yield strength using the empirical relations
2.3 Discretization on an Unstructured Grid
Practical die casting geometries are quite complex. Cartesian grids introduce high staircasing errors near the boundaries. Thus, OpenCast uses unstructured grids with tetrahedral and hexahedral finite volumes. First, a tetrahedral mesh is generated using the open source software GMSH Geuzaine and Remacle (2009)
. Tetrahedral elements have higher cross diffusion terms due to mesh skewness compared to hexahedral elements. Thus, we further use an open source tool TETHEX
Artemiev (2015) to subdivide the tetrahedrons to hexahedrons.The governing equations (1), (2) and (4) can be written as a transport equation for a general scalar :
(10) 
where, is any scalar field, is the diffusion coefficient, and is the source term.
Figure 1 shows two adjacent hexahedral control volumes sharing a common face with vertices , , and . and are cell centers and is the face center. is the unit vector normal to face and in an outward direction with respect to cell . is the distance vector from to . We use a collocated finite volume formulation with all the field variables stored at cell centers.
The surface integral of the diffusion term is approximated as a summation over all the six faces of the cell. The inner product of the normal and the face centered gradient at each face is split into two terms Mathur and Murthy (1997):
(11)  
The first and second terms of equation (11) are direct and cross diffusion terms, respectively. For a structured grid, is parallel to and the cross diffusion term is identically zero as the direct diffusion term reduces to the central difference approximation of first derivative at face center.
In order to estimate the face centered gradient, the strategy used by Mathur and Murthy (1997) for two dimensional grids is extended. A local coordinate system is defined with , and (Fig. 1) as the three axes. , and
are the axes of the global frame of reference. The gradients in both these frames are related by the chain rule of differentiation.
(12) 
where, the subscripts denote derivatives. The Jacobian matrix entries come from the coordinates of the cell centers and the vertices. Value of at each vertex is estimated by averaging from the neighboring cells of the vertex. Thus, the face centered gradient is estimated by inverting the Jacobian matrix in equation (12) and multiplying by .
The surface integral of the convection term is approximated as a summation over all the six faces of the cell. The face value of the field
is estimated by interpolating from the two neighboring cells which share the face. The volume flux passing through the face (
) satisfies the discrete continuity equation. The cross diffusion term has to be accounted for in the computation of the volume flux. The details are given in Section 2.4.3.2.4 Special Modifications in Solidification Regions
For solidification problems, some additional steps are needed in order to handle the extra terms such as the Darcy drag and the latent heat terms in the momentum and energy equations respectively. The velocities in the solid region should go to zero and in the mushy zone, velocities should be significantly lower than the fully liquid region. Simultaneously, the continuity equation has to be satisfied by the face velocities for each control volume. Thus, special care has to be taken in the solution process of the pressure Poisson equation and the velocity correction step.
Figure 2 shows a typical distribution of phases during the solidification process. For ease of visualization, a two dimensional schematic is shown and the same idea has been generalized to three dimensions. The dotted line shows a solidmushy zone interface. The control volumes (cells) on the left of the line are solid and on the right are either liquid or in mushy zone. Cells are labeled with tags S:Solid, LM: Liquid or Mushy. The S or LM tag is assigned to each cell based on the temperature at the previous time step and the liquidus and solidus temperatures of the alloy.
2.4.1 Momentum equation
As described in section 2.2, the modified momentum equations (6) are solved to estimate the intermediate velocities (u*). After discretization, the Darcy drag coefficient () is added to the diagonal term of the linear equations. In the pure liquid region, this coefficient is zero and thus, it does not have any effect. In the mushy zone, it is finite and nonzero and thus, it acts like a resistance to the flow. In the fully solidified region, it is a large number and thus, the u* tends to zero. From computational efficiency point of view, it is not necessary to solve for u* in the solidified cells since it is zero. Therefore, at each time step, before solving the linearized system of equations, the matrix rows corresponding to solidified cells are removed. As these solidified cells are connected to the neighboring mushy or liquid zone cells, the rows corresponding to the neighboring cells have to be modified in a consistent manner. Consider the row corresponding to cell number 1 in Fig. 2:
(13) 
where, is any component of . Originally, is connected to all the neighboring cells from 2 to 10. But since cells 6, 7, 8 and 9 are solidified and their velocity is zero, those rows and columns are deleted from equation (13). Thus, the reduced row becomes:
(14) 
2.4.2 Pressure Poisson Equation
Similar to the momentum equations, the matrix rows corresponding to solidified cells are deleted from the discrete pressure Poisson equation. But in this case, the rows of the neighboring liquid or mushy cells cannot be updated by just deleting the connections of the solid cells as Neumann boundary conditions have to be applied. For any solidified cell, incoming or outgoing flow through all of its faces should be made zero. This is achieved as follows. Each face is shared by exactly 2 cells (face owners). All the cells which share a common vertex with a face are known as its neighbors. The face centered gradient is computed using the values at all connected neighboring cells. For example, cell numbers 1 and 2 are the owners of face F whereas, cells 1, 2, 3, 4, 5, 9 and 10 are its neighbors. The following cases arise for each face:

None of the neighbor cells is solidified: no change in the face centered gradient coefficient (eg. all the faces of cell number 3) is required

None of the owners are solidified but at least one neighbor cell is solidified: flow through the face is allowed but the face centered gradient coefficient has to be modified as the solidified cells are removed from the linear set of equation (eg. faces F and F)

At least one owner is solidified: flow through the face is blocked; thus, contribution of this face in the integrated diffusion term () is zero (eg. faces F and F)
Consider the face F for modification of face centered gradient coefficient. Original coefficients for gradient computation at face F which are valid if none of its neighbor cells are solid are given by:
(15) 
Since cell numbers 6 and 7 have solidified, their contribution has to be removed from equation (15). Thus, the last 2 columns are deleted and those coefficients are smeared equally in the remaining columns for ex., is modified to and to . Since there are 5 cells remaining, the division by 5 is required. After modification, equation (15) becomes:
(16) 
The steps for modification of the discretized pressure Poisson equation are as follows:

Identify the faces with none of the owners solidified but at least one neighbor cell is solidified and smear the coefficients as described above

If at least one owner of the face is solidified, set all of its coefficients to zero

Loop over all the cells:

If none of its face coefficients are modified, its coefficients do not change

If at least one of its face coefficients is modified, reassemble its coefficients

These steps remove the contribution of the solidified cells carefully and reduce the computational effort significantly.
2.4.3 Face Volume Flux Computation and Velocity Correction
The collocated finite volume formulation uses the face centered volume fluxes () in the continuity equation so as to avoid the checkerboarding of pressure. If there is an inconsistency in the numerical formulation of the pressure Poisson equation and the flux computations, there can be a gain or loss of mass and convergence problems. This section describes a consistent method used in the current code to handle solidification.
The volume flux is obtained by taking inner product of the velocity correction equation (8) at the face center with face normal and multiplying by face area:
(17) 
is estimated by averaging the cell values from the two owner cells of the face. is computed exactly in the same way as the regular diffusion term by splitting it into direct diffusion and cross diffusion terms (Eq. (11) with ). The face centered pressure gradient required in the cross diffusion term is estimated by the modified coefficients (Eq. (16)). This volume flux estimate satisfying the discrete continuity equation to a specified tolerance is used in the convection term (Eq. (6)).
The cell centered velocities do not satisfy the discrete continuity equation. They are computed from equation (8) and cell centered pressure gradient is estimated by averaging the face centered gradients.
2.4.4 Latent Heat Term
The GulliverScheil equation (5) which relates temperature with solid fraction is a nonlinear model. The easiest way to numerically couple this with the energy equation is to model the GulliverScheil equation fully explicitly as a source term. The problem with an explicit approach is that the source term destabilizes the discretized energy equation due to high magnitude of the latent heat coefficient. Thus, we use the source term linearization concept discussed by Patankar (1980). The nonlinear term is split into a linear term and a remainder. The linear term is modeled implicitly and the remainder term is treated explicitly. From numerical stability point of view, the coefficient of the linear term which is added to the diagonal of the matrix should be positive. If the coefficient and the remainder term are functions of the unknown (temperature in this case), the equation has to be solved iteratively till convergence.
The latent heat term of the energy equation (4) when integrated over time and control volume gives:
(18) 
where, superscripts and denote last timestep value and iteration number respectively. The value of solid fraction in the subsequent iteration () can be estimated from its latest value () by a first order Taylor expansion:
(19) 
Substituting equation (19) in (18) gives:
(20) 
and are functions of last iteration and last timestep values and thus can be computed first. Note that is always negative and when taken to the left hand side of the equation, it becomes positive and is thus added to the diagonal of the linear system matrix. Adding a positive term to the diagonal helps in stabilizing the system and speeds up convergence. Hence, this approach is found to be much better than the fully explicit method.
The GulliverScheil equation (5) plotted in Fig. 2(a) for a typical aluminum alloy shows that there is a discontinuity at the solidus temperature. Thus, the derivative cannot be computed. To deal with this difficulty, the original equation is modified by smearing the discontinuity near the solidus temperature:
(21) 
where, and is the width of linear smear which can be set to a reasonable value like 2 K. Thus, the derivative can be computed analytically. Figures 2(b) and 2(c) plot the modified solid fraction relation and its derivative respectively.
The overall iterative procedure to obtain the variables at the new time step from values at the old time step can be summarized as:

Initialize: and

Compute and using last iteration values ( and ) by equation (20) and solve the linear system of equations to estimate next iteration value

Update the solid fraction: where, is an under relaxation parameter

Estimate the relative change between the successive iteration values of temperature and solid fraction
Repeat steps until the relative change drops below a desired threshold. For the aluminum alloy used here, it is found that under relaxation is not required i.e., and the solution converges in iterations.
2.5 Solver for Linear Systems
Typical die cast geometries have high aspect ratios i.e., thin cross sections compared to the lateral dimensions. It is found that single grid iterative solvers for the elliptic pressure Poisson equation converge slowly for such geometries. Hence, in this work a multigrid solver is used. The central idea of a multigrid solver is to solve the equations on multiple coarse grids and couple the corrections from all the grids through prolongation and relaxation. The high frequency component of the residual converges fast on the fine grids while the coarse grids are used to accelerate convergence of the low frequency residual. Thus, the coarse grid solutions are used to accelerate the convergence while maintaining the discretization accuracy of the solution at the finest level.
Geometric multigrid is a technique in which multiple levels of grids are generated physically and the matrix vector system is estimated by discretizing the governing equations at each level. The main benefit of this approach is that the matrices at all the levels are obtained directly from the governing equations and thus, good convergence is observed. The main drawback is that generating coarse grids for a complex geometry with unstructured elements is nontrivial. Algebraic multigrid (AMG) tries to address this problem by coarsening the matrix using heuristics based algorithms. This is a black box approach which does not need any physical grids at coarse levels.
The BoomerAMG routine along with Krylov solvers of the open source library HYPRE HYP (2017) developed at the Lawrence Livermore National Laboratory is used in our work. The AMG solver is found to solve the modified momentum equations (6) and the energy equation (4) without difficulty. However, some consistency issues have been observed during the solution of the pressure Poisson equation (9) when extensive solidification happens. In complex geometries, as solidification proceeds, there can be disjoint pockets of metal which are yet to solidify. As these pockets are far enough from each other, they are decoupled numerically in the discrete reduced pressure Poisson equation (section 2.4.2). Hence, in some cases, there are zeros on the diagonals of the coarsest grid level generated by the heuristic coarsening algorithm of BoomerAMG, which creates problems in the convergence of the solution. Thus, the AMG solver has been coupled with a single grid BiCGSTAB solver, also from HYPRE, to be used when AMG is unable to solve the pressure Poisson equation. OpenCast has conditions programmed which automatically switch to BiCGSTAB when the AMG solver diverges. Note that this issue arises towards the end of the simulation when only a few cells are liquid (for instance 20%). Thus, the reduced system is much smaller in size compared to the original problem and a single grid solver is reasonably well convergent.
2.6 Grain Growth and Mechanical Properties Models
Grain size and Secondary Dendrite Arm Spacing (SDAS) are two important parameters used to characterize the microstructure. OpenCast uses empirical relations from the literature for estimation of microstructure parameters and mechanical properties such as yield strength.
SDAS is predicted based on the empirical relationship proposed by Backer and Wang (2007),
(22) 
The model parameters and are chosen to be 39.4 and 0.317, respectively based on the model for microstructure in aluminum alloys Backer and Wang (2007). Material behavior of the die cast alloy is predicted in terms of 2% yield strength () using empirical relationship proposed by Okayasu et al. (2015).
(23) 
Here, is in MPa, (SDAS) is in , and Okayasu et al. (2015).
Grain size estimation is based on the work of Greer et al. (2000). The grain growth rate is given by:
(24) 
where, is the grain size, is the solute diffusion coefficient in the liquid and is the time. The parameter is obtained by invariant size approximation:
(25) 
is given by
(26) 
where, is solute content in the liquid, is solute content in the solid at the solidliquid interface and is the nominal solute concentration. Thus, using the prescribed partition coefficient () and estimated solid fraction (), equations (24)–(26) are solved to get the final grain size.
2.7 Parameter Uncertainty Quantification
The final product quality in die casting is influenced by process parameters like alloy material properties, interface conditions at the mold, thermal boundary conditions etc. Due to the complexity of the process, accurate measurement and control of these parameters is difficult. Conventional deterministic simulations alone are not sufficient to estimate the effect of any stochastic variations on the product quality, thus parameter uncertainty quantification is important. From modeling point of view, parameter uncertainty quantification is a set of stochastic partial differential equations with variation in initial conditions, coefficients and boundary conditions. Stochastic variables are considered as dimensions of the problem in addition to space and time.
To estimate the relation between stochastic process parameters and output parameters, various methods have been proposed in the literature. A popular method is to use a linear combination of polynomial basis functions in the stochastic dimension to expand the output variables. Since orthogonality helps in convergence, orthogonal polynomials are used as basis functions. Wiener’s polynomial chaos Wiener (1938) with Askey family of orthogonal polynomials leads to optimal convergence of the error. Xiu and Karniadakis (2002b)
determined which polynomial leads to exponential convergence depending on the underlying probability distribution function that the stochastic variable follows. For instance, the Hermite polynomials are orthogonal with respect to the standard normal distribution as the weighting function. Hence, it is recommended to use Hermite polynomials as basis functions if the stochastic variable follows normal distribution. Using generalized polynomial chaos expansion developed by
Xiu and Karniadakis (2002b), a second order random field () can be expanded by polynomial basis (Eq. (27)). For all practical purposes, the series is truncated to order .(27) 
In this work, stochastic collocation method is used to estimate the deterministic coefficients of the polynomial chaos expansion. Collocation is a nonintrusive method and thus, modification of deterministic software is not necessary. It acts as a wrapper around existing deterministic software. The deterministic simulation is run at sample points () and a condition is imposed. The right hand side comes from each deterministic simulation and left hand side from polynomial expansion. This gives constraints written in a matrix vector form Smith (2013). ensures that the Vandermonde system (Eq. (28)) is overdetermined. Solving in the leastsquares sense gives .
(28) 
Sample points (
) have to be chosen wisely for successive implementation of stochastic collocation method. For instance, uniformly distributed samples can lead to highly oscillatory basis functions and thus poor convergence. Thus, for one dimensional stochastic problems
Smith (2013), it is popular to choose roots of the basis orthogonal polynomial as sample points. For multiple stochastic dimensions, a simple idea is to use tensor product of the single dimensional samples. The problem with tensor products is that the number of samples grows exponentially with stochastic dimensions. Each sample corresponds to a deterministic simulation and thus the computational expense grows exponentially.
Smolyak (1963) came up with an algorithm to reduce number of samples in multidimensional space maintaining the accuracy of the interpolation. In this research, sparse grid nodes are taken from the work of Heiss and Winschel (2008). A MATLAB based tool UQLab developed by Marelli and Sudret (2014) is used for post processing the simulation outputs to estimate polynomial chaos coefficients and generate the response surfaces.The strategy described above is quite general and can be applied to any numerical solution framework. In this case, the output variables ( of Eq. 27) are temperature and microstructure parameters like grain size and dendritic arm spacing. The stochastic input parameters are boundary conditions, initial conditions and alloy and material properties. The polynomial chaos expansion with stochastic collocation is combined with the deterministic computational fluid flow and heat transfer solver to estimate the sensitivity and uncertainty propagation.
3 Validation with Uncertainty
Due to the complexity of the die casting process, controlled experiments with accurate temperature measurements inside the casting during solidification are difficult. To validate our code, we have therefore used the experimental results of ingot solidification made of SnBi alloy reported by Quillet et al. (2007). In their experiments, a mm ingot is solidified by cooling one of the mm faces at a constant rate (5 K/min) and the remaining 5 faces are thermally insulated. Twenty five thermocouples on the ingot were used to measure temperatures during solidification. Due to the low cooling rate, natural convection velocities were significant and their effect can be seen from the curved contours of temperature. Thus, this validation study covers all the aspects of the software including heat transfer, natural convection and solidification.
As a first step, validation was attempted without incorporating the effects of uncertainty. Figure 3(a) plots the experimental measurements of temperature from the 25 thermocouples as a function of time beginning from solidification. Figure 3(b) plots the temperatures predicted by OpenCast simulations with a constant thermal conductivity. Subsequently, we included the temperature varying thermal conductivity which changes the temperature time curves to Fig. 3(c). This observation shows that the temperature dependency of thermal conductivity has an appreciable effect on the temperature evolution.
Local values of the thermophysical properties such as thermal conductivity and density depend on grain structure. At the length scale of current simulation, it is not possible to predict the grain structure and hence the properties from first principles. Thus, there can be some uncertainty in the input properties. Also, when the alloy solidifies and cools, there is a thermal contraction. This creates a gap between the mold wall and the casting and adds thermal contact resistance, reducing the amount of heat extracted. Further, the experimental measurements are also subject to sensor noise. For example, the thermocouples used by Quillet et al. (2007)
have a reported accuracy of 1K, thereby adding another uncertainty in the predictions. Thus, the validation must account for these uncertainties. We have therefore included the stochastic variation in the boundary conditions and material properties and propagated them using the stochastic collocation method. A confidence interval is estimated about the mean output parameters and the experimental and numerical predictions are compared within the confidence interval.
Uncertainty is therefore added to wall temperature, latent heat, density and thermal conductivity. Since it is difficult to estimate the thermal contact resistance, the wall temperature is specified as a Dirichlet boundary condition by adding an offset to the cooling rate of 5 K/min. The offset is estimated from the experimental temperature plot (Fig. 3(a)). In order to take into account the errors in the temperature measurement, an uncertainty is added on top of the offset. All the four input stochastic parameters are assumed to follow a normal distribution with mean (
) and standard deviation (
) as follows:
Wall temperature offset: C, C

Latent heat: J/kg,

Density: kg/m,

Temperature dependent thermal conductivity: W/mK, W/mK at temperature K
All the thermophysical properties are estimated as a weighted average of individual properties of Sn and Bi taken from the online version of Kaye and Laby (1921). These uncertainties are propagated using the stochastic collocation method described in section 2.7.
In order to make a comparison with the experimental data, the experimental temperaturetime data plot (Fig. 3(a)) is digitized. Since it is difficult to distinguish between the 25 temperatures, the thermocouple with the highest temperature is used for validation. Stochastic collocation is done with three accuracy levels of sample points in order to study its convergence. For estimating the interpolation error, 60 Latin Hypercube samples are used as test points. Deterministic simulations and polynomial chaos give two independent estimates of the same output parameter at the test points. The nondimensional error is defined as the root mean square of difference between these two estimates divided by the maximum value of the parameter. First column of Table 1 denotes the accuracy level of sample points. Accuracy level integrates polynomials upto order exactly Heiss and Winschel (2008). Second column is the number of sample points in two dimensional Smolyak sparse grid i.e., the number of deterministic simulations required ( in Eq. (28)). The last column lists the nondimensional RMS error in computation of the temperature. It can be seen that the error is of order for all the accuracy levels and it decreases with increasing level thus, showing convergence. Hence, level 6 is used for validation.
Accuracy Level  # Sample Pts.  Max Temp. RMS Error 
4  137  7.27E4 
5  385  4.11E4 
6  953  2.60E4 
The PolynomialChaosKriging (PCK) module of UQLab Marelli and Sudret (2014) is used to estimate the output (temperature–time history) as a function of the four input stochastic parameters. Kriging is a Gaussian process modeling tool based on stochastic interpolation algorithm to relate the inputs to outputs. PCK combines the polynomial chaos and kriging methods in a way that is more efficient than the individual methods.
Figure 5 plots the maximum temperatures (with error bands) from OpenCast simulations and experiments Quillet et al. (2007) together with cooling, solidus and liquidus lines. Before the solidification begins, the temperature drops at a rate similar to the cooling line (5 K/min). Thus, the temperature curve and cooling line are parallel. At the liquidus line, a plateau is observed in the temperature curve. This is due to the latent heat release inside the plate. Near the solidus temperature, it can be seen that the numerical simulation shows a small kink in the temperature curve whereas, the experimental curve drops smoothly. This can be attributed to the numerical smearing of the solid fractiontemperature relation (Eq. (21)). There is a sudden release of latent heat due to sharp change in solid fraction near the solidus temperature (Fig. 3). It is found that the kink reduces if the width of linear smear ( in Eq. (21)) is increased. The confidence interval about the mean is plotted as the simulation error band. Since the accuracy of experimental temperature measurement is 1 K Quillet et al. (2007), an experimental error band is plotted as K about the mean value. It can be seen that the experimental and numerical estimates are in good agreement thus validating the computational tool. This validation study also emphasizes the utility of uncertainty quantification as previously there was a mismatch between experiments and computations when no uncertainties were included.
4 Deterministic Solidification Results of Realistic Geometries
In order to demonstrate the utility of the developed software, some complex die casting geometries were simulated. Many die casting geometries have thin cross sections and solidify quickly in a few seconds. Thus, natural convection velocities are negligible and may have little effect on solidification. Thus, simulating without natural convection may be acceptable as it saves significant computational effort. Hence, in all the simulations of this section, natural convection was neglected.
Figure 6 shows the hexahedral grids of two selected geometries. Both the clamp and the pulley have approximately 300,000 control volumes each. Initially, the mould is filled with a liquid alloy at 1000 K and the walls are held at 500 K. Aluminum alloys with around 812% Silicon content are popular in die casting. Thus, for this simulation, an Al10%Si binary alloy is chosen. The liquidus, solidus and freezing temperatures are taken from the phase diagram.
Figure 7 plots the temperature contours for the clamp geometry along the midplane in Zdirection at different timesteps during the solidification. Figure 8 shows the corresponding solid fraction contours. Figure 9 plots microstructure parameters SDAS, grain size and yield strength estimated using the previously described empirical models. As time progresses, the cooling rates and temperature gradients decrease. Hence, the regions with highest thickness or the core take longer time to solidify. The core regions typically have higher SDAS, grain size and lower yield strength as the grains have longer time to grow. Such trends are seen in Fig. 9. The pulley is axisymmetric in geometry and loading conditions. Thus, SDAS, grain size and yield strength are plotted only along the XY plane as shown in Fig. 10. Even for this case, the core region has higher SDAS and grain size and lower yield strength.
5 Uncertainty Quantification Results
The clamp geometry is subsequently used to study the effect of stochastic variation in the process parameters on the outputs. The process parameters which can be tuned in a die casting foundry are alloy composition (the solute percentage), initial melt temperature and the wall temperature. It is assumed that these three parameters follow normal distribution () with the following means () and standard deviations ():

Solute concentration wt %

Initial molten alloy temperature K

Wall temperature K
Thermal properties like thermal conductivity, density and specific heat are estimated as a weighted linear combination of the properties of aluminum and silicon with solute concentration as the weight. Since there is a significant variation with temperature Kaye and Laby (1921), temperature varying properties are used. Alloy properties are obtained from the AlSi phase diagram (Fig. 11). The freezing and solidus temperatures are 660C and 577C, respectively. The solidus and liquidus curves can be approximated by straight lines. The partition coefficient can be estimated as the ratio of liquidus to solidus line slope Dantzig and Tucker (2001). The liquidus temperature is calculated from the liquidus line as a function of the solute concentration.
Accuracy Level  # Sample Pts.  Output 1  Output 2  Output 3  Output 4 

3  27  1.20E2  3.82E2  2.70E3  1.61E2 
6  216  7.41E3  3.35E2  2.42E3  1.42E2 
The impact of the three input parameters is studied on the following four outputs:

Total solidification time

Maximum SDAS

Minimum yield strength

Maximum grain size
where, the maximum or minimum is taken over the entire domain. The polynomial chaos method with stochastic collocation described in section 2.7 is used for uncertainty propagation. The error due to stochastic interpolation is estimated using 60 Latin hypercube samples. Estimates of the same output are obtained from the complete simulation and the polynomial chaos expansion (PCE) independently. The difference between these quantities nondimensionalized by their mean value is defined as the error. Since the errors in estimating all the four outputs reduce by using higher accuracy levels (Table 2), it can be seen that the collocation method has converged. The results presented correspond to the accuracy level 6. In order to get an idea of the interpolation error visually, both the simulation and PCE estimates can be plotted on the same graph. In the hypothetical scenario of exact interpolation, all the points should lie on the line. However, there is a deviation from the line because of the interpolation error. Figure 12 shows plots for each of the four outputs. Each output is normalized by subtracting its mean and dividing by its standard deviation and thus, is nondimensional. Most of the points follow the trend of the
line except some outliers. The outliers correspond to those random samples which are too far from the means of the input parameters. Table
2 and Fig. 12 prove that the polynomial chaos has converged and is accurate enough for further use.Sensitivity of each output with respect to each input can be easily estimated once an accurate polynomial chaos expansion is obtained. The sensitivity analysis tool of the software UQLAB Marelli and Sudret (2014) is used to estimate the total Sobol indices (Fig. 13). It can be seen that each output is highly sensitive only to one input parameter. Sensitivity is practically important as it gives an idea as to which input parameter should be tightly controlled. Thus, other parameters can be loosely controlled saving cost but yielding desired product quality at the same time. The amount of heat extracted from the wall is proportional to the temperature gradient near the wall and hence, solidification time is highly sensitive to the wall temperature. On the other hand, the microstructure parameters like SDAS, yield strength and grain size depend on the solute concentration. Thus, these three parameters are more sensitive to solute concentration than the wall and initial temperatures.
In the response surface (Fig. 14), the most important input is plotted on the xaxis. For solidification time response surface, the wall temperature (xaxis) and initial temperature (yaxis) are chosen. For other three outputs, solute concentration (xaxis) and wall temperature (yaxis) are chosen. It can be seen that all the response surface contours are nearly vertical. This confirms that the output is most sensitive to the input parameter plotted on the xaxis. The maximum grain size contours are nonlinear. This implies that the sensitivity varies locally in the input parameter space. For other three outputs, the contours are almost linear and thus, it can be concluded that the local sensitivity is similar everywhere and independent of the input parameter value.
6 Conclusions
This paper describes a numerical software framework OpenCast for simulations of solidification problems, including natural convection. Microstructure parameters and mechanical properties are estimated using published empirical relations. The flow equations are solved only in the liquid zone by modifying the coefficients of the discrete equations. The algebraic multigrid solver together with a Krylov subspace solver is used to solve the pressure Poisson equation. Complex geometries are meshed with unstructured hexahedral elements. Parameter uncertainty quantification is used as a wrapper over the deterministic simulations in order to assess the effect of stochastic variations in the inputs on the outputs.
The software is validated against published experimental results of solidification. This validation study shows the significance of stochastic analysis since it is observed that validation without uncertainty is unsuccessful. The validated software is used to simulate two practical die casting geometries. Sensitivity and uncertainty analysis is also performed. Sensitivity analysis shows that the product quality given by grain size and yield strength is highly sensitive to the solute concentration whereas, the productivity given by the solidification time is sensitive to the mold wall temperature. These results are practically useful as it gives an idea about the important input process parameters. The response surfaces show the variation of the outputs with the important input parameters. They can be used to get quick estimates of the outputs without running full deterministic simulations and also to get local sensitivities.
Although this paper demonstrates the ability of OpenCast to simulate die casting problems, it is a general purpose software. OpenCast can also be used to simulate other manufacturing processes such as sand casting, additive manufacturing, welding etc. The use of unstructured elements together with algebraic multigrid method adds complete flexibility to simulate arbitrary geometries. The coupling of uncertainty and sensitivity analysis tools enhance the power of the deterministic numerical method.
Acknowledgment
Authors acknowledge the financial support from Digital Manufacturing and Design Innovation Institute (DMDII Grant No. 150706). The authors also thank Steve Udvardy and Beau Glim of North American Die Casting Association (NADCA) for providing industry contacts. Technical discussions with Alex Monroe of Mercury Castings and his suggestions were helpful in this work.
References
References
 Bennon and Incropera (1987a) W. Bennon, F. Incropera, A continuum model for momentum, heat and species transport in binary solidliquid phase change systems—I. Model formulation, International Journal of Heat and Mass Transfer 30 (10) (1987a) 2161–2170.
 Bennon and Incropera (1987b) W. Bennon, F. Incropera, A continuum model for momentum, heat and species transport in binary solidliquid phase change systems—II. Application to solidification in a rectangular cavity, International Journal of Heat and Mass Transfer 30 (10) (1987b) 2171–2187.
 Voller and Prakash (1987) V. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection–diffusion mushy region phase–change problems, International Journal of Heat and Mass Transfer 30 (8) (1987) 1709–1719.
 Vynnycky and Kimura (2007) M. Vynnycky, S. Kimura, An analytical and numerical study of coupled transient natural convection and solidification in a rectangular enclosure, International Journal of Heat and Mass Transfer 50 (2526) (2007) 5204–5214.
 Bennett (2007) B. Bennett, Adaptive numerical modeling of natural convection, conduction, and solidification within mold cavities, Numerical Heat Transfer, Part A: Applications 51 (4) (2007) 313–342.
 Wang et al. (2010) S. Wang, A. Faghri, T. Bergman, A comprehensive numerical model for melting with natural convection, International Journal of Heat and Mass Transfer 53 (910) (2010) 1986–2000.
 Plotkowski et al. (2015) A. Plotkowski, K. Fezi, M. Krane, Estimation of transient heat transfer and fluid flow for alloy solidification in a rectangular cavity with an isothermal sidewall, Journal of Fluid Mechanics 779 (2015) 53–86.
 Hu et al. (2017) Y. Hu, D. Li, S. Shu, X. Niu, Lattice Boltzmann simulation for threedimensional natural convection with solidliquid phase change, International Journal of Heat and Mass Transfer 113 (2017) 1168–1178.
 Cleary et al. (2010) P. Cleary, J. Ha, M. Prakash, T. Nguyen, Short shots and industrial case studies: understanding fluid flow and solidification in high pressure die casting, Applied Mathematical Modelling 34 (8) (2010) 2018–2033.
 Gau and Viskanta (1986) C. Gau, R. Viskanta, Melting and solidification of a pure metal on a vertical wall, Journal of Heat Transfer 108 (1) (1986) 174–181.
 Quillet et al. (2007) G. Quillet, A. Ciobanas, P. Lehmann, Y. Fautrelle, A benchmark solidification experiment on an Sn–10% wt Bi alloy, International Journal of Heat and Mass Transfer 50 (34) (2007) 654–666.
 Hachani et al. (2012) L. Hachani, B. Saadi, X. Wang, A. Nouri, K. Zaidat, A. BelgacemBouzida, L. AyouniDerouiche, G. Raimondi, Y. Fautrelle, Experimental analysis of the solidification of Sn–3 wt.% Pb alloy under natural convection, International Journal of Heat and Mass Transfer 55 (78) (2012) 1986–1996.
 Davidson (1996) L. Davidson, A pressure correction method for unstructured meshes with arbitrary control volumes, International Journal for Numerical Methods in Fluids 22 (4) (1996) 265–281.
 Haselbacher (1999) A. Haselbacher, A gridtransparent numerical method for compressible viscous flows on mixed unstructured grids, Ph.D. thesis, 1999.
 Kim and Choi (2000) D. Kim, H. Choi, A secondorder timeaccurate finite volume method for unsteady incompressible flow on hybrid unstructured grids, Journal of Computational Physics 162 (2) (2000) 411–428.
 Ferziger and Peric (2012) J. H. Ferziger, M. Peric, Computational methods for fluid dynamics, Springer Science & Business Media, 2012.
 Mathur and Murthy (1997) S. Mathur, J. Murthy, A pressurebased method for unstructured meshes, Numerical Heat Transfer 31 (2) (1997) 195–215.
 Muzaferija and Gosman (1997) S. Muzaferija, D. Gosman, Finitevolume CFD procedure and adaptive error control strategy for grids of arbitrary topology, Journal of Computational Physics 138 (2) (1997) 766–787.
 Yang and Henson (2002) U. Yang, V. Henson, BoomerAMG: a parallel algebraic multigrid solver and preconditioner, Applied Numerical Mathematics 41 (1) (2002) 155–177.
 Ruge and Stüben (1987) J. Ruge, K. Stüben, Algebraic multigrid, in: Multigrid Methods, SIAM, 73–130, 1987.
 HYP (2017) HYPRE: scalable linear solvers and multigrid methods, URL http://www.llnl.gov/CASC/hypre/, 2017.
 Provatas et al. (1999) N. Provatas, N. Goldenfeld, J. Dantzig, Adaptive mesh refinement computation of solidification microstructures using dynamic data structures, Journal of Computational Physics 148 (1) (1999) 265–290.
 Karma and Rappel (1998) A. Karma, W. Rappel, Quantitative phasefield modeling of dendritic growth in two and three dimensions, Physical Review E 57 (4) (1998) 4323.
 Ramirez et al. (2004) J. Ramirez, C. Beckermann, A. Karma, H.J. Diepers, Phasefield modeling of binary alloy solidification with coupled heat and solute diffusion, Physical Review E 69 (5) (2004) 051607.
 Backer and Wang (2007) G. Backer, Q. Wang, Microporosity simulation in aluminum castings using an integrated pore growth and interdendritic flow model, Metallurgical and Materials Transactions B 38 (4) (2007) 533–540.
 Maxwell and Hellawell (1975) I. Maxwell, A. Hellawell, A simple model for grain refinement during solidification, Acta Metallurgica 23 (2) (1975) 229–237.
 Greer et al. (2000) A. Greer, A. Bunn, A. Tronche, P. Evans, D. Bristow, Modelling of inoculation of metallic melts: application to grain refinement of aluminium by Al–Ti–B, Acta Materialia 48 (11) (2000) 2823–2835.
 Desnain et al. (1990) P. Desnain, Y. Fautrelle, J.L. Meyer, J.P. Riquet, F. Durand, Prediction of equiaxed grain density in multicomponent alloys, stirred electromagnetically, Acta Metallurgica et Materialia 38 (8) (1990) 1513–1523.
 Xiu and Karniadakis (2002a) D. Xiu, G. Karniadakis, Modeling uncertainty in steady state diffusion problems via generalized polynomial chaos, Computer Methods in Applied Mechanics and Engineering 191 (43) (2002a) 4927–4948.
 Xiu and Karniadakis (2003) D. Xiu, G. Karniadakis, Modeling uncertainty in flow simulations via generalized polynomial chaos, Journal of Computational Physics 187 (1) (2003) 137–167.
 Maitre et al. (2001) O. Maitre, O. Knio, H. Najm, R. Ghanem, A stochastic projection method for fluid flow: I. basic formulation, Journal of Computational Physics 173 (2) (2001) 481–511.
 Carnevale et al. (2013) M. Carnevale, F. Montomoli, A. D’Ammaro, S. Salvadori, F. Martelli, Uncertainty quantification: A stochastic method for heat transfer prediction using LES, Journal of Turbomachinery 135 (5) (2013) 051021.
 Marepalli et al. (2014) P. Marepalli, J. Murthy, B. Qiu, X. Ruan, Quantifying uncertainty in multiscale heat conduction calculations, Journal of Heat Transfer 136 (11) (2014) 111301.
 Ganapathysubramanian and Zabaras (2007) B. Ganapathysubramanian, N. Zabaras, Sparse grid collocation schemes for stochastic natural convection problems, Journal of Computational Physics 225 (1) (2007) 652–685.
 Fezi and Krane (2017) K. Fezi, M. Krane, Uncertainty quantification in modeling metal alloy solidification, Journal of Heat Transfer 139 (8) (2017) 082301.
 Hosder et al. (2006) S. Hosder, R. Walters, R. Perez, A nonintrusive polynomial chaos method for uncertainty propagation in CFD simulations, in: 44th AIAA Aerospace Sciences Meeting and Exhibit, 891, 2006.
 Kumar et al. (2016) D. Kumar, M. Raisee, C. Lacor, An efficient nonintrusive reduced basis model for high dimensional stochastic problems in CFD, Computers & Fluids 138 (2016) 67–82.
 Fajraoui et al. (2017) N. Fajraoui, M. Fahs, A. Younes, B. Sudret, Analyzing natural convection in porous enclosure with polynomial chaos expansions: Effect of thermal dispersion, anisotropic permeability and heterogeneity, International Journal of Heat and Mass Transfer 115 (2017) 205–224.
 Spiegel and Veronis (1960) E. Spiegel, G. Veronis, On the Boussinesq approximation for a compressible fluid., The Astrophysical Journal 131 (1960) 442.
 Dantzig and Tucker (2001) J. Dantzig, C. Tucker, Modeling in materials processing, Cambridge University Press, 2001.
 Harlow and Welch (1965) F. Harlow, J. Welch, Numerical calculation of timedependent viscous incompressible flow of fluid with free surface, The Physics of Fluids 8 (12) (1965) 2182–2189.
 Geuzaine and Remacle (2009) C. Geuzaine, J. Remacle, Gmsh: A 3D finite element mesh generator with builtin preand postprocessing facilities, International Journal for Numerical Methods in Engineering 79 (11) (2009) 1309–1331.
 Artemiev (2015) M. Artemiev, TETHEX: Conversion from tetrahedra to hexahedra, URL https://github.com/martemyev/tethex, 2015.
 Patankar (1980) S. Patankar, Numerical heat transfer and fluid flow, CRC Press, 1980.
 Okayasu et al. (2015) M. Okayasu, S. Takeuchi, M. Yamamoto, H. Ohfuji, T. Ochi, Precise analysis of microstructural effects on mechanical properties of cast ADC12 aluminum alloy, Metallurgical and Materials Transactions A 46 (4) (2015) 1597–1609.
 Wiener (1938) N. Wiener, The homogeneous chaos, American Journal of Mathematics 60 (4) (1938) 897–936.
 Xiu and Karniadakis (2002b) D. Xiu, G. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM Journal on Scientific Computing 24 (2) (2002b) 619–644.
 Smith (2013) R. Smith, Uncertainty quantification: theory, implementation, and applications, vol. 12, SIAM, 2013.
 Smolyak (1963) S. Smolyak, Quadrature and interpolation formulas for tensor products of certain classes of functions, in: Doklady Akademii Nauk, vol. 148, Russian Academy of Sciences, 1042–1045, 1963.
 Heiss and Winschel (2008) F. Heiss, V. Winschel, Likelihood approximation by numerical integration on sparse grids, Journal of Econometrics 144 (1) (2008) 62–80.
 Marelli and Sudret (2014) S. Marelli, B. Sudret, UQLab: A framework for uncertainty quantification in MATLAB, in: Vulnerability, Uncertainty, and Risk: Quantification, Mitigation, and Management, 2554–2563, 2014.
 Kaye and Laby (1921) G. Kaye, T. Laby, Tables of physical and chemical constants: and some mathematical functions, Longmans, Green and Company, 1921.
 Murray and McAlister (1984) J. Murray, A. McAlister, The Al–Si (Aluminum–Silicon) system, Bulletin of Alloy Phase Diagrams 5 (1) (1984) 74.
Comments
There are no comments yet.