A meshless approach to thermomechanics of DC casting of aluminium billets

The ability to model thermomechanics in DC casting is important due to the technological challenges caused by physical phenomena such as different ingot distortions, cracking, hot tearing and residual stress. Many thermomechanical models already exist and usually take into account three contributions: elastic, thermal expansion, and viscoplastic to model the mushy zone. These models are, in a vast majority, solved by the finite element method. In the present work the elastic model that accounts for linear thermal expansion is considered. The method used for solving the model is of a novel meshless type and extends our previous meshless attempts in solving fluid mechanics problems. The solution to the problem is constructed using collocation on the overlapping subdomains, which are composed of computational nodes. Multiquadric radial basis functions, augmented by monomials, are used for the displacement interpolation. The interpolation is constructed in such a manner that it readily satisfies the boundary conditions. The discretization results in construction of a global square sparse matrix representing the system of linear equations for the displacement field. The developed method has many advantages. The system of equations can be easily constructed and efficiently solved. There is no need to perform expensive meshing of the domain and the formulation of the method is similar in two and three dimensions. Since no meshing is required, the nodes can easily be added or removed, which allows for efficient adaption of the node arrangement density. The order of convergence, estimated through an analytically solvable test, can be adjusted through the number of interpolation nodes in the subdomain, with 6 nodes being enough for the second order convergence. Simulations of axisymmetric mechanical problems, associated with low frequency electromagnetic DC casting are presented.


Introduction
The thermomechanical phenomena that occurs during DC casting of aluminium alloys plays an important role in designing the casting device and determining the process parameters. Therefore, there exists a great interest in developing numerical models that are able to describe this kind of phenomena. The effort reflects not only in improvement of casting process itself but also in a substantial amount of related publications (see for example [1,2] and references therein).
The majority of existing numerical models are based on the established finite element method (FEM) [3], which allows for efficient solutions of involved physical models which take into account the thermal expansion, linear elasticity and viscoplastic phenomena on a transient growing domain of the DC cast ingot [4,5]. By the best knowledge of the present authors, the meshless numerical methods have not been used for the related purpose until now. It is the objective of the present work to fill this gap.
The meshless methods, as compared to the established mesh-based methods, have several advantages such as: the meshless formulation requires no additional topological structure to be established (the spatial positions of nodes suffice), allowing for efficient addition or removal of the computational nodes. The method allows for unstructured node arrangements to be used and for the solution in complicated geometries to be easily obtained. Last but not least, the method can be generalized to higher dimensional problems needing virtually no modifications in formulation [3].

Physical model
The primary focus of the present paper is in applying a new numerical method to the considered engineering problem. Therefore, the physical model is kept quite simple. A linear thermoelastic model, resulting from the stress equilibrium, Hooke's law for inhomogeneous isotropic body and linear isotropic thermal expansion is considered. The following governing equation, written in terms of the deformation field u , is obtained [3] Here G stands for shear modulus,  for Lamé parameter and f for the body force. The coupling with the temperature field is described by the coefficient  defined as T stands for reference temperature at which the thermal expansion is considered to be zero.
The displacement, traction and symmetry boundary conditions are considered. The boundary condition that limits the displacement in the contact with the rigid mold is solved by iteration.
Since the metal undergoes a phase change, the elastic properties substantially vary with temperature. The related dependencies for the considered alloy were obtained from the JMatPro database [3].

Solution procedure
The present thermomechanics model is used in a post-processing mode after the calculations with inhouse developed model [4] of heat and mass transfer are performed. The stationary state temperature field, obtained from the model, is used as a known input.

Local radial basis function collocation method
The local radial basis function collocation method (LRBFCM) is used to solve the boundary value problem, given by equation (1). A simplified method formulation is given in the following subsections, with a more elaborate formulation stated in [6].
where c is the shape parameter and s  the maximum distance between subdomain points in  coordinate. The purpose of adding the monomials to the set of basis functions is to allow for exact reproduction of constant and linear deformation fields. The interpolant is given by  determined by collocation and m denoting the number of augmentation monomials. One of the key points of the method is the way the boundary conditions are treated. Let us introduce the linear boundary conditions operator at point j specified by the equation . In case the subdomain of a certain point includes a point that lies on the boundary, the corresponding collocation equation is replaced by the equation specifying the boundary condition.
To obtain a non-degenerated system, the augmentation equations are added. The system of equations can be written in a matrix form as , , where indices ,j l count discretization nodes, i basis functions and ,,    the spatial directions. The matrix elements of interpolation matrix A and the components of the right hand side vector  are given by The differential operators can be evaluated, by using the described interpolation, in an arbitrary point of the computational domain. By inverting the collocation matrix A , the interpolant can be written in terms of unknown values in collocation nodes. This fact can be used to elegantly discretize the governing equation.
The governing equation is written in a short form as Du g     . To allow for a reasonably compact notation of the discretized system of equations, the domain In this manner a system of equations with dimensions 2 2 N N  is obtained. Since the interpolation is local, the obtained system of equations is sparse. This allows the use of highly efficient specialized solvers [5].

Treatment of the mold contact
An axisymmetric billet DC casting is considered in the present work, described in a related coordinate system

Convergence test
The performance of the method is tested on a two-dimensional domain specified by a unit square The temperature field is assumed to be ) sin 2 ( The convergence order is shown in figure 2. We can see that the order of convergence depends on the number of the nodes used in the subdomain. The method is of the second order if 6 or 10 nodes are used and of the fourth order if 15 or 21 nodes are used in each of the subdomains.
Another important issue in meshless methods based on radial basis functions is that they employ a shape parameter [6]. The dependence of the error on the shape parameter is shown in figure 3. We can see that for each number of subdomain points, there exists a certain optimal shape parameter. The optimal value gets smaller as the number of subdomain points is increased.

The DC casting example
The temperature field is assumed from our compatible meshless heat and mass transport model (HMM) for low frequency electro-magnetic DC casting, elaborated in [4]. The computational domain extends from 0 to 1 m height with the diameter of 0.15 m. The mechanics domain is smaller than the one used in the HMM, since only the area below the highest level of liquidus isotherm is considered, resulting in top 0.03 m of the domain being cut off. The mold extends from 0.905m to 1m. The procedure described in subsection 3.2 is used to treat the boundary conditions. The rest of boundary conditions are illustrated in figure 5. The material properties used in the model correspond to AA6082 alloy and are displayed in figure 4.
In figure 6 the deformation in radial direction is shown, along with all three diagonal components of the stress tensor. The shrinkage at the bottom of the computational domain agrees with the expected shrinkage due to the difference in density between the solid and the liquid phase the temperature field for this case is shown in the leftmost plot in figure 7. The fact that the stress state is very sensitive to the changes in temperature field is demonstrated in figure 7. In the figure we see the comparison between stress fields obtained for two different amplitudes of the driving current. The changes in temperature field are quite hard to notice. The sump depth is slightly larger and the mushy zone is a bit thinner for the larger driving current. The change in the stress state is however more dramatic. In the centre of the billet the nature of the stress changes from compression to tension, which may reduce the quality of the produced billet due to hottearing.

Conclusions
It is evident that the developed LRBFCM is a viable novel method for modelling the linear thermoelastic problems as used in a basic thermomechanics model of DC casting. As is well known from the literature, the order of the convergence depends on the number of the nodes in each subdomain and the overall convergence strongly depends on the value of the shape parameter. Solutions obtained by the method for DC casting in the presence of low frequency electro-magnetic field are shown. The method performs well for the simple material model. The described  The temperature fields were obtained using in-house HMM [4].
4th International Conference on Advances in Solidification Processes (ICASP-4) IOP Publishing IOP Conf. Series: Materials Science and Engineering 117 (2016) 012032 doi:10.1088/1757-899X/117/1/012032 6 developments thus successfully complement our previous work on meshless methods for fluid mechanics problems [8,9]. In future we plan to add a viscoplastic constitutive model to the framework in order to more realistically describe the alloy behaviour at elevated temperatures.