Gradient based porosity calculation in casting simulation

Treball realitzat per: Daniel Díez Rodríguez Dirigit per: Dr. Riccardo Rossi and Jordi Rubio Màster en: Métodos Numéricos en la Ingeniería Barcelona, 15/06/2017 Departamento de Ingeniería Civil y Ambiental Gradient Based Porosity Calculation in Casting Simulation UNIVERSITAT POLIT` ECNICA DE CATALUNYA Escola T`ecnica Superior d’Enginyers de Camins Canals i Ports de Barcelona MASTER THESIS: GRADIENT BASED POROSITY CALCULATION IN CASTING SIMULATION A Master Thesis submitted by Daniel D´ıez for the Master in Numerical Methods in Engineering Supervised by: Dr. Riccardo Rossi and Jordi Rubio Abstract Porosity is one of the most common defects in casting production. When some parts of the casting begin to solidify they contract due to the difference of densities between the liquid and solid phases. The areas that are still liquid feed the solidifying regions until at some point feeding is not longer possible due to the growth of the solid phase and porosity begins to form. If the feeding cut is due to the growth of dendrites, pores will have a very small size and will be usually mentioned as microporosity. When a pool of liquid is isolated, surrounded by solid material, a big void will appear known as macroporosity. If the liquid pool is connected to the exterior it will form pipe shrinkage. The aim of this work is the implementation of a microporosity model inside the commercial casting software Click2Cast, based on the general finite elements framework Kratos Multiphysics. The model that has been implemented is based on the Dimensionless Niyama method developed by Carlson and Beckermann in 2008. This approach requires the recover of the gradient of the temperature. In order to do this, several methods to compute the gradients were developed and implemented inside the code. Finally, another model to calculate macroporosity and pipe shrinkage is being currently implemented. This method is based on the calculation of the pressure during the solidification. No results are provided for this model as still further work is needed to expect a correct under different circumstaces .The objective is that after validation, this model will be also included in Click2Cast. 2 CONTENTS Contents 1 Introduction 14 1.1 CastingDefects............................... 14 1.2 Objectives.................................. 16 1.3 Structureofthework............................ 16 2 Porosity Literature Review 18 2.1 Classification and Formation of Porosity . . . . . . . . . . . . . . . . . 18 2.2 Shrinkage Porosity Simulation . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Microporosity Models : Criterion Functions . . . . . . . . . . . . . . . 22 2.4 Macroporosity Models : Thermal / Volume Calculation Models . . . . 24 2.5 Coupled Microporosity and Macroporosity Models : Thermal/ Fluid FlowModels................................. 25 3 Gradient Reconstruction 30 3.1 Introduction................................. 30 3.2 L2 Projection Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2.1 Results................................ 35 3.3 Newton-Raphson iterative scheme . . . . . . . . . . . . . . . . . . . . . 37 3.3.1 Results................................ 39 3.4 ZhangMethod................................ 41 3.4.1 Description of the method . . . . . . . . . . . . . . . . . . . . . 42 3.4.2 Results................................ 44 4 CONTENTS 3.5 Pouliot Gradient Recover Method . . . . . . . . . . . . . . . . . . . . . 46 3.5.1 Description of the Method in 1-D . . . . . . . . . . . . . . . . . 47 3.5.2 The Multidimensional Case . . . . . . . . . . . . . . . . . . . . 49 3.5.3 Results................................ 50 3.6 Comparison of results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4 Microporosity Model 58 4.1 Introduction................................. 58 4.2 Dimensionless Niyama Criterion . . . . . . . . . . . . . . . . . . . . . . 61 4.3 StudyCases................................. 66 4.3.1 Lever................................. 67 4.3.2 BrakeCaliper............................ 68 4.3.3 OilPump .............................. 70 4.3.4 Comments on the results . . . . . . . . . . . . . . . . . . . . . . 72 4.4 GradientRecovery ............................. 73 4.5 Sensitivity to the element size . . . . . . . . . . . . . . . . . . . . . . . 74 4.6 Sensitivity to different evaluation temperatures . . . . . . . . . . . . . 76 4.7 Casting parameters that have influence on microporosity . . . . . . . . 77 4.8 Conclusions ................................. 80 5 Macroporosity and Pipe Shrinkage Model 84 5.1 Introduction................................. 84 5.2 GoverningEquations ............................ 84 5 CONTENTS 5.3 ShrinkagePorosity ............................. 88 5.4 OpenShrinkage............................... 89 5.5 WeakForm ................................. 89 5.6 Resolution Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6 Conclusions and Future Work 94 6 1 INTRODUCTION 1 Introduction 1.1 Casting Defects Metal casting is one of the oldest manufacturing processes where liquid metal is poured into a mold which has the shape of the desired final part that will be obtained once the liquid metal solidifies. The main advantage of this processes with respect to other manufacturing processes is that very complex geometries can be easily produced. Only limitation in this aspect is that we have to be able to open the mold and get the final casting. Besides the mold itself, other components can be used in order the obtain the final shape such as the cores. Using casting, a part made of practically any metal and with any size can be produced with relatively no effort compared to others manufacturing processes. For these reasons, casting processes have a huge influence in the manufacturing market and its volume of production has been regularly increasing. However, the internal processes that occur inside the casting are very complex and the appearance of defects is very hard to eliminate. According to American Iron and Steel Institute, it is estimated that between 3% and 7% of the total produced castings have to be scraped due to the presence of defects. These defects can be classified as follows: •Filling related defect: Filling defects include misruns, cold shuts and inclusions. A misrun occurs when a front of liquid metal solidifies before the casting is completely filled leaving an unfilled portion. Cold shuts are a result of two fronts of liquid metal which do not fuse properly inside the cavity of the mold due to their differences, specially in temperature, that will lead to weak areas where the material is not regular and cracks can form. Both misruns and cold shuts are caused by a lack of fluidity in the molten metal. An inclusion is any impurity in the pour metal like oxides which will result in worse material properties of the final casting. •Gas Porosity: Gas porosity is the formation of bubbles within the casting after it has cooled. Liquid metals usually carry dissolved gas in the melt, but when the metal starts solidifying the amount of gas that the solid phase can contain is much lower than the amount of gas that was present in the liquid. The consequence of this is that gas will try to create its own phase nucleating and 14 1 INTRODUCTION forming gas porosity (see figure 1.1). Besides, gasses can be present inside the mold if as the liquid flow fills the mold cavity, the air that was inside the mold has not enough time to exit ending up as entrapped air. •Shrinkage Defect: Shrinkage defects due to the differences of densities between the liquid and solid phases. As the metal solidifies it might happen that there is no available feed metal to compensate for the shrinkage and voids will appear. These empty areas are classified as pipes or caved surfaces if they are present on the surface of the casting or shrinkage porosity if it appears inside the casting. •Metallurgical defects: When the solid metal is still hot, its mechanical properties are reduced and residual stresses due to the solidification might end up cracking the material and forming the so-called hot tears. Figure 1.1: Gas Porosity [10] In order to eliminate these defects and achieve optimum operating conditions industrial experiments are usually conducted. This kind of experiments are complicated due to the nature of the problem and can be very expensive. Though there is steel a need of these experimental results, the casting producers understand more and more that they can improve the product quality through computer simulation. This approach facilitates the optimization of the process by increasing energy efficiency and avoiding defects with less cost. 15 1 INTRODUCTION In this work we will focus on the prediction of porosity as one of the most persistent and common complains of casting users. Figure 1.2: Shrinkage Porosity [10] 1.2 Objectives The main objective of this work is to implement a gradient based model to compute microporosity in the finite elements commercial code Click2Cast which will be included in next versions of the code. The code has been implemented using the multidisciplinary framework Kratos Multiphysics. As the main key of this calculation is the gradient of the temperature, the implementation and comparison of different gradient recover methods has been also part of the work. At the current time, another model to compute macroporosity is being implemented with the objective of being included as well in future versions of Click2Cast. 1.3 Structure of the work In section 2 a review of the available literature about porosity is exposed, explaining the different approaches to the problem. One of the models for microporosity presented in this section has been implemented in Click2Cast and explained in detail in section 4. For the implementation of this method, it is previously needed to compute the gradient of the temperature. Different ways to compute it are shown and compared in section 3 . Finally, in section 5, it is presented the current and future work in a macroporosity and pipe shrinkage model which is being tested and implemented in these moments. 16 1 INTRODUCTION 17 2 POROSITY LITERATURE REVIEW 2 Porosity Literature Review 2.1 Classification and Formation of Porosity Porosity is one of the major defects in metal castings. Its presence results in a deterioration of the mechanical properties of the final part. Particularly tensile strength and fatigue resistance decrease significantly when porosity forms, resulting in a failure in the performance of the material for lower loads or in a lower number of cycles that we could expect. Despite all the efforts and studies that have been carried out in order to fully understand how porosity forms and grows, there is no clear agreement on which are the mechanisms that rule this process. One of the biggest difficulties in this problem is the considerable number of variables that play a role in the solidification. We can classify them in one of the following groups: •Material Properties •Composition of the Alloy •Casting Conditions •Mold Properties •Geometry of the part All these variables will somehow affect the solution and all of them will have to be taken into account. However, it is very hard to know the relative effect that each of them will have compared to the rest. For this reason, despite casting has been part of the human activities for thousands of years, its most complex aspects are not yet well known and it is almost impossible to consider all of them in one single model. These difficulties should not discourage us to investigate in this field. More the opposite, this means that this problem requires more investigation and efforts as porosity is known to be one of the most important aspects in casting simulation and a major cause of casting rejections. Porosity can be classified according to different criteria, usually origin and size. The soundness of a casting depends on the flow of liquid metal to the solidifying regions 18 2 POROSITY LITERATURE REVIEW Figure 2.1: Types of Porosity [30] that will contract as they become solid. The primary cause of porosity is therefore the difference in the densities of liquid and solid phases and the obstruction of the fluid flow as the solidification advances. When this feeding is cut from a region, shrinkage porosity will form. Besides the contraction of the metal there are other factors that can affect to the amount of porosity present. The presence of gas dissolved in the metal will help the pore nucleation to begin and grow. The solubility of gas in the solid metal is lower than the solubility in the melt and therefore, as the solidification advances more gas will be present in the liquid favoring the nucleation of the pores. By the size of the pores we can distinguish between microporosity or macroporosity. Macroporosity As the solidification advances from the exterior to the center of the filling cavity it might happen that large pools of liquid end up surrounded by solid metal. As there is no possibility of feeding through the solid walls, when the whole pool becomes solid there will be a big region of void in the center called macroporosity 19 2 POROSITY LITERATURE REVIEW Figure 2.2: Microporosity in ductile iron [20] Microporosity On the other hand, microporosity forms in regions where there is a poor feeding due to the growth of dendrites. If the gas forms in the mushy zone after the dendrite coherency point has been reached, it will be entrapped in the dendrite network, in the so called interdendritic spaces and will nucleate forming small pores called microporosity or microshrinkage. The size of this pores ranges between 1-100 µm . In this case, these pores will grow following the dendritic shape instead of being round. Open Shrinkage All the previous forms of porosity are classified as closed shrinkage defect. However, if the casting part is open somewhere to the exterior it will form open shrinkage defects or pipe shrinkage. In these regions the free surface is capable of moving, lowering its level as the rest of the metal solidifies. As long as the solid fraction is smaller than the coherency point it may be assumed that the solid crystal move along with the liquid. In this situation the air will enter in the cavity forming characteristic pipe forms. Once solidification is complete, only limited solid feeding through elastic and plastic deformation is possible. The latter is responsible in part of the formation of caved surfaces on the exterior of the part. 20 2 POROSITY LITERATURE REVIEW 2.2 Shrinkage Porosity Simulation There are many commercial software packages that give results for shrinkage porosity. Each of them have to select the physics that will be modeled in order to have an accurate model. However, users usually have a limited computational capacity and a trade-off between accuracy and computational time must be satisfied. For this reason one must be very careful deciding which aspects of the solidification and porosity modeling will be taken into account. In general, we can number the following aspects that will have influence in the nucleation and growth of the pores. •Thermal variables such as the temperature and solid fraction. •Solid phases growth •Flow variables such as the pressure or the velocity. •Interfacial characteristics between the gas and liquid or solid phases. We will review later porosity models that will have in account some of the phenomena mentioned above. However, none of them will have all of them into consideration as such a model would be too costly and inefficient for commercial purposes. Besides, we have to be aware that many of the parameters that are required to perform a simulation are unknown and have to be estimated. Even if we do experimental measures, the related error of the own measurement procedure might be higher than any other inaccuracy that we might introduce in the model due to simplifications. This fact only reinforces the idea that it was stated before that the final model that will be implemented must ensure a not only precise but also efficient calculation. In general, the existing models to compute the porosity can be classified into different groups according to Lee et.al [20] and Stefanescu [30]. •Criterion Function Models •Thermal/Volume Calculation Models •Thermal/Fluid Flow Models We will further classify them according to type of porosity that they aim to model. 21 2 POROSITY LITERATURE REVIEW •Microporosity Models : Criterion Functions •Macroporosity Models : Thermal/Volume Calculation Models •Coupled Microporosity and Macroporosity Models : Thermal/Fluid Flow Models For each of the families above it will be shown different implementation methods in the following sections. 2.3 Microporosity Models : Criterion Functions Criterion functions are simple rules that relate the local conditions to the propensity to form micropores. The obtained criterion are functions of thermal parameters such as the cooling rate, the thermal gradient or the solidification time and material parameters as the shrinkage factor, the viscosity etc. These models take into account simplified formulations for the transport problems, however they usually ignored the contribution of gas rejected by the solidifying melt to microporosity. The application of these criterion started a long time ago with Pellini in 1953 [26]. Pellini was one of the first persons to realize the correlation between the presence of porosity and the gradient of the temperature. He stated that the temperature gradient Gshould be greater than a critical value Gcr to avoid centerline shrinkage porosity. Niyama took this idea from Pellini and developed his own criterion function simply known as Niyama Criterion [24]. This model has become quite popular specially between Japanese foundry men and it has been proved to very useful for detecting microporosity in ferrous alloys. On the other hand, its implementation for non-ferrous alloys has been frequently questioned[30] . As Pellini did before, Niyama et al. took measures in a number of commercial castings and confirmed the previous experiments that had related the formation of shrinkage porosity to the gradient of the temperature. However, they found that the critical value under which microporosity could be expected depended largely in other factors as well such as the shape of the casting or its size so it could not be predicted in advanced [23]. They proposed therefore that the critical value of the gradient of the temperature depended also on the solidification time ts. More specifically they stated that it should be proportional to 1/√tS. Approximating the solidification time as tS=4Tf/˙ T, where 4Tfis the difference 22 2 POROSITY LITERATURE REVIEW Carlson Carlson presented in 2002 [6] a multi-phase model that predicts melt pressure, feeding flow and porosity formation and growth in steel castings during solidification. Carlson derived a momentum equation by combining Darcy’s Law and Stokes flow that is valid everywhere in the mushy zone. ∇2v=fl Kv+fl µl∇P−fl µl ρg (2.16) The permeability Kis computed in the usual way using the Kozeny-Carmen equation. K=KrK0 (1 −fs−fp)3 (fs+fp)2(2.17) Where K0= 6 ×10−4λ2 1, being λ1the primary dendrite arm spacing, and Kris the relative permeability between air and liquid Kr=fl/(fl+fa). When the solid fraction is low, the permeability tends to infinite and the momentum equation reduces to Stokes’ flow equation. When the solid fraction is high at the end of the solidification the laplacian term in the left hand side is very small and the equation is equivalent to Darcy’s Law. Finally the species equation accounts for macro segregation of gas species due to the flow. Once the total gas pressure is high enough to cause pore nucleation, the amount of porosity that forms is determined from the continuity equation. This multi-phase model has been successfully implemented in a general-purpose casting simulation code. This model has been adapted and will be describe in detail in section 5. 29 3 GRADIENT RECONSTRUCTION 3 Gradient Reconstruction As stated in the introduction and in section 2 several methods for computing porosity are present in the literature about casting. In this section and the following we will focus on microporosity. Due to the nature of this phenomena, where different processes affect each other in a very small scale, criterion function models are usually chosen, using only local values of the variables of interest. Between these models, the Dimensionless Niyama model proposed by Carlson and Beckermann [5] has gained notoriety during the last years and it has been implemented in some commercial software packages specialized in casting simulation. Unlike other criterion function models like [24, 26] this method accounts for both thermal and material properties. The key aspect of this criterion is its dependency on the gradient of the temperature. Therefore, an accurate and efficient method has to be implemented to recover the gradients from the known temperature field that we have already calculated. With this aim, several methods will be explained in the following sections and will be latter compared in terms of precision and computational cost. 3.1 Introduction In many problems in the field of engineering the recovery of gradients is needed for different reasons. For example in the field of the fluid dynamics one might be interested in calculating the drag coefficient or the vorticity of the field. This quantities can be determined as a post-process of the previously obtained solution for the pressure and the velocity. In our case, we are interested in calculating the microporosity of a casting part during solidification. It has been experimentally proved that the presence of microporosity is closely related to areas where the gradient of the temperature is low. The most simple way to obtain the gradient of a known field is just to take derivatives inside each element. If our finite element solution was of degree k, then the gradient field that we obtain in such a way is of degree k-1.This means that if we use linear finite elements then our solution will be discontinuous on the nodes. We want of course our solution to be continuous so we will focus on methods to obtain a gradient field such that it belongs to the same space of interpolation as the original 30 3 GRADIENT RECONSTRUCTION solution. There are many authors that have been interested in this problem in the field of error estimation. In finite elements problems it can be proved that the error of the solution is a function of the element size h, therefore as happroaches zero so does the error. Evidently, this will imply also a higher computational cost. The approach taken in the called adaptive mesh refining methods is just to refine those regions where the errors are higher. Following the lines firstly proposed by Zienkiewicz and Zhu [25], the usual way to compute these a posteriori errors is to evaluate them as the difference between the piecewise discontinuous gradients (if we use linear elements) and the so called superconvergent post-processed recovered gradients. The reason for this name is that according to the authors the rate of convergence of these methods is at least quadratic, that is, order 2 with respect to the element size. With this purpose, some authors proposed their own method to recover the gradients [32, 28, 25]. In this chapter we will explore some of these proposed methods. After each models is explained in detail, we will find out which one is the best for our needs. We will compare the following methods in terms of the L2 errors with respect to analytical functions and the computational time. •L2 Projection with Consistent Mass Matrix •L2 Projection with Iterative Newton Raphson •L2 Projection with Lumped Mass Matrix •Zhang & Naga Method •Pouliot Recover Method These methods will be tested for a structured regular mesh with different element sizes. Element Size Number of Nodes Number of Elements 0.100 1331 6000 0.050 9261 48000 0.025 68921 384000 0.020 132651 750000 Table 3.1: Meshes used for testing gradient recover methods 31 3 GRADIENT RECONSTRUCTION Two different analytical functions will be used in order to asses the accuracy of each method. Case 1 f1=sin(πx)sin(πy)sin(πz) (3.1) Case 2 f2=e−25x+e−25y+e−25z(3.2) With the corresponding analytical gradients: ∇f1=   πcos(πx)sin(πy)sin(πz) πsin(πx)cos(πy)sin(πz) πsin(πx)sin(πy)cos(πz)   (3.3) ∇f2=   −25e−25x −25e−25y −25e−25z   (3.4) The error ewill be measured using the L2 norm, defined as: e=ZΩ (∇Tana −φh)2dV (3.5) This integral is going to be approximated using the obtained nodal values of the gradient and the nodal volumes Vi. e=X nodes ∇Ti ana −φi2Vi(3.6) Where ∇Tana is the analytical gradient obtained either from equation 3.3 or 3.4 as it corresponds and φhis the solution of a specific gradient recovery method. In Kratos, the nodal volume is defined as the sum of the volumes of the elements to which that specific node belongs divided by the number of nodes of that element, four for our case since only tetrahedron are used in Click2Cast. 32 3 GRADIENT RECONSTRUCTION 3.2 L2 Projection Method The fist and most simple approach that we can think of is to project the discontinuous gradient of the temperature field that we have already calculated into the L2 space. The method and the subsequent implementation follows the same principles explained in [15]. If we have previously calculated a temperature field Th, we can obtain a function for the gradient ∇Tby simply taking the spatial derivative. If we are using linear elements, our function Thwill belong to the linear finite element space Vh∈ L2and ∇Twill be a piecewise continuous function. Th=XNiTi(3.7) ∇T=XNi,jTi(3.8) The most simple approach is to project ∇Tinto the finite element space using an L2 Projection. That is, to find a solution φh∈ Wh, being Wha finite element space Wh∈(L2)3, such that minimizes the error ea in a Least Squares sense. e(x) = ZΩ (ϕh(x)−∇T(x))2dΩ (3.9) This minimization problem can be solved in the following equivalent way. We first take derivative with respect to ϕand multiply by the test function N, arriving to the weak form of the problem. ZΩ N·(ϕh(x)−∇T(x))dΩ = 0 (3.10) The arbitrary test functions Nare the shape functions of our finite element space Wh. 33 3 GRADIENT RECONSTRUCTION ZΩ N·ϕh(x)dΩ = ZΩ N·∇T(x)dΩ (3.11) Knowing that ϕh=XNjϕj(3.12) We get ZΩ NiNjϕjdΩ = ZΩ Ni∇T(x)dΩ (3.13) Using the following compact notation MC=ZΩ NiNjdΩ (3.14) f=ZΩ Ni∇T(x)dΩ (3.15) Rewriting the system we obtain Mij Cϕj=fi(3.16) We can now build the whole system and solve it. This global system is going to be expensive in terms of computational cost and must be avoided when possible. The problem can be simplified using a lumped mass matrix MLinstead of the consistent mass matrix MC. For a tetrahedral element: MC=V olume 20       2 1 1 1 1 2 1 1 1 1 2 1 1 1 1 2       (3.17) 34 3 GRADIENT RECONSTRUCTION ML=V olume 4       1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1       (3.18) As we see, the lumped mass matrix is obtained adding the values of the rows in diagonal terms. This lumped mass matrix is diagonal and therefore the global system will be diagonal as well. We can avoid then to solve the global system as each node can be solved independently. The diagonal terms of the left hand side will coincide with the Nodal Volume Vi, so we can simply calculate the nodal values as the sum of the element contributions divided by the Nodal Volume. ϕi=Pcont elems fi Vi (3.19) 3.2.1 Results Figure 3.1: L2 Projection - L2 Error for Case 1 35 3 GRADIENT RECONSTRUCTION Figure 3.2: L2 Projection - L2 Error for Case 2 Figure 3.3: L2 Projection - Computational Time for Case 1 36 3 GRADIENT RECONSTRUCTION Figure 3.4: L2 Projection - Computational Time for Case 2 Method Case 1 Case 2 L2 Lumped 1.52 1.35 L2 Consistent 1.51 1.45 Table 3.2: L2 Projection - Convergence Rates 3.3 Newton-Raphson iterative scheme When we substitute the mass matrix Mcby the lumped mass matrix Mlwe are doing an approximation, and therefore we will obtain a different solution from the one we would have using the consistent mass matrix. This issue can be overcome using a iterative Newton-Raphson scheme. From the equation 3.16 we can define the following residual R(ϕ) = fi−Mij Cϕj(3.20) Our objective is to find the values of ϕthat make R(ϕ) = 0. If we start from an arbitrary solution φi, we can approximate our solution using a 1st order Fourier 37 3 GRADIENT RECONSTRUCTION expansion R(ϕi+1) = R(ϕi) + ∂R(ϕi) ∂ϕi(4ϕi+1) = 0 (3.21) Rearranging terms we get ∂R(ϕi) ∂ϕi4ϕi+1 =−R(ϕi) (3.22) From the equation 3.20 we calculate ∂R(ϕi) ∂ϕi=−MC(3.23) Including this result into the previous equation we get the following iterative scheme MC4ϕi+1 =f−MCϕ(3.24) Now, we can approximate the mass matrix of the left hand side with the lumped mass matrix to obtain the following iterative scheme 4ϕi+1 =M−1 L(f−MCϕ) (3.25) The solution of this system will converge to ϕ=MC −1fin a number of steps. Other methods have been proposed following similar lines, like D.M Hawken et al. proposed in [13]. In this paper a scheme with optimal convergence was given starting from the same equation that we formulated in 3.16. Mij Cϕj=fi(3.26) In this occasion the iterative scheme that Hakwen proposes has an iterative convergence acceleration parameter wsuch that 38 3 GRADIENT RECONSTRUCTION Figure 3.10: Zhang - L2 Error for Case 2 Figure 3.11: Zhang - Computational Time for Case 1 45 3 GRADIENT RECONSTRUCTION Figure 3.12: Zhang - Computational Time for Case 2 Method Case 1 Case 2 Zhang 2 1.88 Table 3.4: Zhang - Convergence Rates 3.5 Pouliot Gradient Recover Method As well as Zhang and some others before, Pouliot et al. became interested in the gradient recovery in the field of a posteriori error estimators. This error is simply measured as the difference between the recovered gradient and the piecewise discontinuous solution gradients that we can directly obtain from our finite elements solution. Therefore, a number a methods have been developed with the aim of approximate the gradient of the solution called super-convergent methods. These error estimators are then used to drive isotropic or even anisotropic mesh adaptation. For a linear FEM solution the error is directly linked to the Hessian matrix, so not only the gradient but the Hessian matrix (the gradient of the gradient) needs to be recovered. With a more suitable method to recover the gradients, this can be achieved more efficiently, with a very practical application in almost all engineering 46 3 GRADIENT RECONSTRUCTION problems that involve finite elements. In the paper published by Pouliot et al. [28], a review of the methods that were beings used up to that date. As it was pointed out by Zhang and Naga in [32], the most popular of the superconvergent methods, the Zienkiewicz and Zhu patch method, lost its superconvergence property when non regular meshes were used. In the same paper, Pouliot mentioned the Zhang method from [32], recognizing its capability to maintain its superconvergence in a wide number of cases. Despite the fact that this method can provide very good results inside the domain, Pouliot et al. pointed out that as in most of the gradient recovery methods, annoying oscillations were observed near the boundaries of the domain. Finally, Pouliot also mentioned the L2 Projection method with a consistent mass matrix that we explained in section 3.2 though this method was disregarded by Pouliot due its linear convergence rate and its relatively high computational cost. In the Pouliot Recover Method we have to build a global system that must be solved using an iterative method. The method constructs linear equations using superconvergent points on each edge of the mesh. This method, as well as the others presented in this section is superconvergent but besides it should be particularly good at reducing the oscillations near the boundaries. 3.5.1 Description of the Method in 1-D To expose the method we will begin with a 1D case where the domain is an interval [a, b] divided in n+ 1 elements and nnodes such that a=x0< x1< ... < xi< ... < xN=b(3.40) We start from a piecewise linear solution uh, which nodal values are already known. The derivative of this function is a piecewise constant function. The objective is to reconstruct the function using the nodal values to obtain dh, a piecewise linear and continuous approximation of the derivative of our function uh. Pouliot achieves this by introducing a bubble function ubidefined by ubi(x) = Ci(x−xi)(x−xi+1) in the interval of the element [xi, xi+1] and vanishing elsewhere and where Ciare unknown DOFs. This function is added to our previous solution obtaining a new solution 47 3 GRADIENT RECONSTRUCTION up=uh+ubIf the coefficients Ciare properly chosen we can expect our solution to go from second-order to a third-order approximation of uin L2. The derivative our function is now simply u0 p(x) = dh(x) = u0 h(x) + u0 b(x) (3.41) Notice that the bubble function is unknown so we cannot directly compute dh. However if we want to evaluate the derivative in the middle point of the element xs= (xi+xi+1)/2, the derivative of the bubble function is always equal to zero. These points are called superconvergent points. dh(xs) = u0 h(xs) (3.42) Knowing that dhand uhare linear functions, we can write this equation just in term of the nodal values di+di+1 2=ui+1 −ui xi+1 −xi (3.43) Which can be rewritten in an elementary form as "1 1 1 1#" di di+1#=−2 hi"ui−ui+1 ui−ui+1#(3.44) Notice that we are imposing one equation for each edge of the mesh. It usually happens that meshes have more nodes than edges and therefore our system will be singular and we will not be able to solve it. To overcome this problem, Pouliot et. al. propose to modify the system adding small terms in order ti stabilize the formulation. As an example, the propose to use a Laplacian operator to the derive multiplied by a stabilization parameter h. ("1 1 1 1#+h"1−1 −1 1 #)" di di+1#=−2 hi"ui−ui+1 ui−ui+1#(3.45) The parameter his usually chosen small (h= 10−6h) 48 3 GRADIENT RECONSTRUCTION 3.5.2 The Multidimensional Case The generalization of the method to the 2 and 3D cases is quite straight from the 1D case. In this cases, the superconvergent points are located in the middle of the edges of the mesh. We have to replace now the derivatives at nodes by the directional derivatives along the edges. We will show the formulation of the 2D case, being the 3D case very similar. For each edge we calculate: Le x=xe 2−xe 1(3.46) Le y=ye 2−ye 1(3.47) he=q(Le x)2+ (Le y)2(3.48) le= (le x, le y) = 1 he(Le x.Le y) (3.49) Figure 3.13: Element Edge [28] Now we write the value of the directional derivative in the center of the edge ass a function of the nodal values, arriving to: de 1+de 2 2·le=ue 2−ue 1 he(3.50) 49 3 GRADIENT RECONSTRUCTION This can be written as:       (le x)2le xle y(le x)2le xle y le xle y(le y)2le xle y(le y)2 (le x)2le xle y(le x)2le xle y le xle y(le y)2le xle y(le y)2             de x,1 de y,1 de x,2 de y,2       =−2 he       le x(ue 1−ue 2) le y(ue 1−ue 2) le x(ue 1−ue 2) le y(ue 1−ue 2)       (3.51) This elementary system is computed for each element of the mesh and then assembled in the global system. As a remark, in general, for 2D and 3D meshes, there are more nodes than edges so the problem will be non-singular. However, this is not the case for regular meshes, for which stabilization has to be added just as it was explained in the previous section. 3.5.3 Results Figure 3.14: Pouliot - L2 Error for Case 1 50 3 GRADIENT RECONSTRUCTION Figure 3.15: Pouliot - L2 Error for Case 2 Figure 3.16: Pouliot - Computational Time for Case 1 51 3 GRADIENT RECONSTRUCTION Figure 3.17: Pouliot - Computational Time for Case 2 Method Case 1 Case 2 Pouliot 2 1.82 Table 3.5: Pouliot - Convergence Rates 52 3 GRADIENT RECONSTRUCTION 3.6 Comparison of results Figure 3.18: Comparison - L2 Error for Case 1 Figure 3.19: Comparison - L2 Error for Case 2 53 3 GRADIENT RECONSTRUCTION Figure 3.20: Comparison - Computational Time for Case 1 Figure 3.21: Comparison - Computational Time for Case 2 54 4 MICROPOROSITY MODEL 4.2 Dimensionless Niyama Criterion The dimensionless Niyama criterion is derived from a solidifying 1-D system. The flow of a liquid inside a porous media is described with Darcy’s Law: flvl=−K µl ∂P ∂x (4.2) where lis the liquid fraction, vlis the velocity of the liquid phase, µlis the dynamic viscosity and Pis the melt pressure. The permeability of the mushy zone Kis calculated from the Kozeny-Carman relation: K=K0 f3 l (1 −fl)2(4.3) K0=λ2 2/180 (4.4) being λ2is the secondary dendrite arm spacing (SDAS). Using the shrinkage factor β= (ρs−ρl)/ρlto simplify the mass conservation equation and then integrating the result, we can arrive to the conclusion that the velocity of the liquid phase or shrinkage velocity is constant and can be expressed as vl=−βR, where Ris the isotherm velocity which can be calculated at the same time as the relation between the gradient of the temperature ∇Tand the cooling rate ˙ T. The whole formula will be thus: vl=β∇T ˙ T(4.5) Inserting this equation into equation 4.2: flβ∇T ˙ T=−K µl ∂P ∂x(4.6) reordering terms 61 4 MICROPOROSITY MODEL ∂P ∂x=−flβ∇T ˙ T µl K(4.7) As the solidification advances, the melt pressure starts falling from an initial value Pliq down to a critical value Pcr at which porosity begins to form. The initial value of the pressure Pliq is simply determined by the hydro-static pressure due to the head of liquid metal h. In gravity casting: Pliq =Patm +ρgh (4.8) In this term we can include the packing pressure for high pressure simulations. This pressure is made by the piston when the casting is already solidifying in order to make the void areas collapse and reduce the porosity. We can simply include it by using: Pliq =Ppacking +ρgh (4.9) Where ρis the density of the liquid, his the head of liquid and gis the gravity. This pressure can be determined using a mechanical equilibrium. In the inside we have the pore pressure Ppand in the outside the critical pressure Pp. We can also take into account the capillary pressure Pσwhich also tries to close the pore. Therefore, the equilibrium can be written as: Pcr +Pσ=Pp(4.10) The capillary pressure is usually described as: Pσ=2σ r0 (4.11) Where σis the surface tension between the liquid and the pore and r0is the initial radius of curvature of the pore. If the presence of gasses in the melt is negligible the pressure inside the pore will be null and the critical pressure will simply be: 62 4 MICROPOROSITY MODEL Pcr =−2σ r0 (4.12) Notice that the value that we assign to the capillary pressure is rather arbitrary as the initial radius of the pore will depend on many factors and cannot be known in general. It is sometimes evaluated as r0=λ2/4. In more rough estimations, the whole term is neglected and they simply use Pcr = 0 With this idea in mind we can now integrate the equation 4.7 from the beginning of the solidification when fl= 0 and P=Pliq until porosity forms when fl=fl,cr and P=Pcr. We define as well 4Pcr =Pliq −Pcr ZPliq Pcr dP =4Pcr =Z0 xcr µlβ˙ Tfl K∇Tdx =ZTliq Tcr µlβ˙ Tfl K∇TdT =µlβ˙ T (∇T)2Z1 fl,cr fl K d˙ T dfl dfl (4.13) In order to have a non-dimensional formula, a dimensionless temperature is introduced θ= (T−Tsol)/4Tf, where Tsol is the solidus temperature and 4Tfis the alloy freezing range. We can define the integral value now as I(fl, cr) = Z1 fl,cr 180(1 −fl)2 f2 l dθ dfl dfl(4.14) 63 4 MICROPOROSITY MODEL Figure 4.3: Mushy Zone 1D Analysis from [5] Putting altogether we arrive to 4Pcr =µlβ4Tf λ2 2 ˙ T (∇T)2I(fl, cr) (4.15) We are going to define here the dimensionless Niyama Ny∗=∇Tλ2√4Pcr qµlβ4Tf˙ T =pI(fl, cr) (4.16) The dimensionless Niyama criterion given by 4.16 accounts not only for the local thermal conditions, but also for the properties and solidification characteristics of the alloy and the critical pressure drop across the mushy zone. If we determine the SDAS as a function of the cooling rate from the relation λ2=Cλ˙ T−1/3the dimensionless Niyama can be expressed in the following alternate form 64 4 MICROPOROSITY MODEL N∗ y=Cλ∇T ˙ T5/6s4Pcr µlβ4Tf (4.17) Notice that in this equation, unlike the usual Niyama Criterion which is proportional to ∇T/ ˙ T1/2, the dimensionless Niyama criterion is proportional to ∇T/ ˙ T5/6 which according to the authors will result in a more generally applicable criterion for varying section thicknesses and mold materials. The approach here is to determine the liquid fraction that is still present in a point when the pressure reaches the critical value. From this point, feeding is assumed to be impossible and the amount of porosity will be calculated applying the conservation of mass. ∂ρ ∂t +∇·(ρv) = 0 (4.18) Where ρis the mixture density, defined as ρ=flρl+fpρp+fsρsand each subscript s,land paccounts for the solid, liquid and porosity contribution respectively. Making use of the fact that 1 = fl+fs+fpand neglecting the density of the porosity we can define the density as: ρ=ρs−fpρs+fl(ρl−ρs) (4.19) If there is no feeding ∇·(ρv) = 0 and therefore the equation 4.18 is simplified to: ∂(ρs−fpρs+fl(ρl−ρs)) ∂t = 0 (4.20) ∂fp ∂t =(ρl−ρs) ρs ∂fl ∂t (4.21) Integrating from the point where porosity begins to form, that is when the pressure reaches the critical value Pcr ,fl=fl,cr and fp= 0 until the end of the solidification when fp=fend pand fl= 0. 65 4 MICROPOROSITY MODEL fend p=(ρs−ρl) ρs fl,cr (4.22) Recalling that β= (ρs−ρl)/ρlwe end with: fend p=β β+ 1fl,cr (4.23) The algorithm to solve the problem is the following for each point: 1. Compute the local thermal conditions and obtain the properties of the material. 2. Calculate the value of the Dimensionless Niyama using equation 4.16 or 4.17 3. Compute the critical fluid fraction from equation 4.14 knowing that from equation 4.16. N∗ y=pI(fl, cr) 4. Compute the porosity using the mass continuity equation. There is another aspect that we have not mentioned yet and might have a great relevance. We refer to the moment of the solidification at which we compute the gradient of the temperature and the cooling rate. In the paper, Beckermann and Carlson suggest to compute these values when the temperature is a ten per cent of the freezing range above the solidus temperature, this is TNy =Tsol + 0.14Tf(4.24) They based this decision on the assumption that porosity forms late on the solidification, however there are other authors that claim that porosity could nucleate and grow for values of the temperature much higher than this one. In the following sections we will check the importance of this factor in different simulations. 4.3 Study Cases In this section we are going to show some results for microporosity using the method previously exposed using different parts and materials. All the simulations were made 66 4 MICROPOROSITY MODEL using Click2Cast. The results are shown as void fraction volume using Hyperview, a post-processing software from Altair Engineering. For all the results, the variables were evaluated at TNy =Tsol + 0.34Tf, the packing pressure was set to Ppacking = 0 and the secondary dendrite arm spacing was SDAS = 40µm. To recover the gradients, the iterative L2 Projection method proposed in section 3.3 was used with 5 iterations. The filling was not computed for any simulation. 4.3.1 Lever In the first case, the simulation corresponds to the lever of a motorbike. The material used for the casting was Aluminum AlSi7Mg with a temperature at the beginning of the simulation Tinitial = 700oC. The mold material was Steel-H13 with an initial temperature of 150 C. The solidification was simulated using 1,185,840 elements. Figure 4.4: Lever Part 67 4 MICROPOROSITY MODEL Figure 4.5: Lever - Regions over 2% of Microporosity Figure 4.6: Levers - Midplane As we can see in the figures above most of the porosity is located at the center of the casting part with a maximum value of 3.29% of porosity. It is also present with lower intensity in some areas that have small thicknesses in both the extremes of the casting. 4.3.2 Brake Caliper The second study case correspond to a brake caliper made of Magnesium AE42 with a temperature at the beginning of the simulation Tinitial = 727oC. The mold material was Steel-H13 with an initial temperature of 150 C. 68 4 MICROPOROSITY MODEL The solidification was simulated using 1,425,004 elements. Figure 4.7: Brake Caliper Part Figure 4.8: Brake Caliper - Regions with over 0.9% Microporosity 69 4 MICROPOROSITY MODEL Figure 4.9: Brake Caliper - Section Cut In this casting the microporosity is distributed on the more massive areas at the top and the bottom as well as in the thin walls of the center cylinder. In this case, the maximum porosity value is 3.33%. 4.3.3 Oil Pump The last case that we will show in this section corresponds to a oil pump made of Carbon-Steel ASTM-A216. The initial temperature of the liquid was set to Tinitial = 1592C. Same as in the other cases, the mold material was Steel-H13 with an initial temperature of 150 C. For this casting, a mesh of 1,183,087 elements was created. 70 4 MICROPOROSITY MODEL Figure 4.16: Evolution of the temperature gradient - Red: Exterior node - Blue: Interior Node In this graphic we see the evolution of the gradient of the temperature for an interior and an exterior node. It is clear that though the values in the interior are always lower, the point at which we take the value will have a big influence. Peaks are present due to the shape of the temperature-solid fraction curves, which are usually not linear, and the effect of the latent heat during the solidification. Although Beckermann et al. suggested in [5] to use Teval =Tsol + 0.14Tfas evaluation point, we have obtained more realistic results using Teval =Tsol + 0.34Tf. This is no surprise as each commercial code uses different formulations and parameters and therefore different results will be obtained. 4.7 Casting parameters that have influence on microporosity In this subsection we are going to see how different casting parameters affect the final amount of microporosity. This value is going to be calculated using the nodal values in the following way: Microporosity(%) = Pnodes Vifp,i VT (4.25) Where Viis the nodal volume, fp,i is the nodal value of the microporosity and VT 77 4 MICROPOROSITY MODEL is the total volume of the part. Heat Transfer Coefficient (HTC) Figure 4.17: Microporosity - HTC Influence The heat transfer coefficient, usually known as HTC is the relation between the heat flux and the temperature difference between two points: HTC W/m2oC=q 4T(4.26) It basically gives us an idea of how fast heat is extracted. It depends on the material of both the casting and the mold and the gap that might form as solidification advances. This graphic clearly shows us than low HTC values will lead to a higher amount of microporosity. This is caused by lower temperature gradients that will make the flow of liquid metal through the mushy zone more difficult. For this reason, we can expect in general to have lower levels of microporosity in steel molds than in sand molds. 78 4 MICROPOROSITY MODEL Packing Pressure Figure 4.18: Microporosity - Packing Pressure Influence In many occasions, the metal is forced to be introduced in the cavity of the mold by a piston instead if letting it flow simply by gravity. When the filling has concluded, the piston keeps pressuring the metal, in a process known as packing pressure, in order to ensure that metal does not flow back outside, but also to reduce the level of porosity as the pressured metal will try to close the empty spaces. The effect of this pressure can be introduced in our model and as we expected, it helps to reduce the amount of microporosity. Feeding Temperature Figure 4.19: Microporosity - Feeding Temperature Influence 79 4 MICROPOROSITY MODEL The feeding temperature is an important factor in the casting since it must be ensure that the metal does not solidify anywhere during the filling, which would provoke defects like misruns as explained in section 1. In counterpart, as we see in the graphic, high feeding temperatures will lead to slightly higher levels of microporosty, though its effect is very low. Mold Temperature Figure 4.20: Microporosity - Mold Temperature Influence In this last graphic we can see how higher mold temperatures have a considerable negative effect on microporosity. Again, this is due to the presence of lower temperature gradients that as we stated is the key aspect in microporosity. All these results agree with [22], where it was experimentally shown an increase of the mold temperature would be followed by an increment in the total amount of porosity in the casting part. Moreover, in the same work, Y.M. Lee proved the little influence of the temperature of the liquid feeding. This conclusion also matches our results. 4.8 Conclusions The Dimensionless Niyama and the original Niyama Criterion have been object of many debates within the academic world. It cannot be denied its importance as practically all casting simulation packages offer at least the Niyama values and some 80 4 MICROPOROSITY MODEL of them based their microporosity prediction on it. However one must be aware of the important limitations that this approach has and the need of some interpretation by the user. In this last section we will try to put together some of the conclusions that can be found along the literature about Niyama and Dimensionless Niyama together with our own conclusions from the results that we have obtained. First thing that we have to take into account is that the nucleation and growth of pores is a very complex problem and many factors have to be taken into account. As stated in section 2.2 and [30] the following aspects should be included in an ideal model. •Thermal Field •Solid phases growth •Flow Field •Interfacial characteristics between the gas and liquid or solid phases. The Niyama Criterion and in general most of the criterion functions only computes the thermal field and neglect the importance of the rest of the fields. The Darcy’s Law based assumption for the flow does not cover the whole physics behind the defects formation. For example in clean melts, porosity only appear because of the segregation of gases which is completely neglected in Niyama. In spite of being true that the velocities in the solidification phase are small, they are going to have a key role since the liquid will flow from the interior areas to those where the solidification has already begun. Besides, though some researchers have approved its application for simple steel casting tests [7, 3] in general its use is not advised if working with complex geometries as it does not take into account the geometry of the part or other types of alloys which might have very different freezing ranges [31, 30]. One could wander how this criterion has acquired so much notoriety with the limitations mentioned above in sight. We have to realize first that it is very difficult to validate microporosity results as most of the radio graphic techniques that are usually used to measure porosity fail to detect microporosity as it belongs to a very small size scale. Secondary, as its simplicity might make it questionable for general applications, it is also a reason why it became so popular since results can be obtained very quickly 81 4 MICROPOROSITY MODEL and with little implementation effort. This was specially important at the moment when it was developed by Niyama et al. in the early 80’s where computers were much less powerful and complex simulations could not be run as easily as nowadays. Beckermann and Carlson created the Dimensionless Niyama Criterion with the aim to address some of the doubts that have been formulated by other authors. They tried to come with a solution where other factors could be taken into account. Besides the local thermal variables (the cooling rate and the gradient of the temperature), their formulation took into account also material parameters in order to be able to use it with materials with different freezing ranges. Other of the factors that allows us to have a better control of the simulation is the critical pressure, which defines a threshold under which porosity will appear. Although this pressure is sometimes simply taken as the atmospheric pressure, a more precise approximation can be taken including the effect of the gravity and the capillary pressure. As we have shown in subsection 4.7 this model is able to deal with different parameters that will affect to the problem in a realistic way. We have seen in 4.3 that results require some degree of interpretation. Dimensionless Niyama and general Niyama criterion also mark areas where macroporosity will be present. These results do not have to be considered as stated in [18] and for that reason, another model to compute macroporosity is required. In this study we have shown some calculations for some example castings were microporosity seemed to give reasonable qualitative results. Regarding to the quantitative aspect, many variables such as SDAS or the evaluation temperature might require some tuning due to its high influence in the final result. Beckermann and Carlson show some experimentally agreement in their paper, but only for very simple geometries. It is actually very hard to quantify even experimentally the actual values of microporosity, so it can be quite difficult to validate the model in this aspect. Future work must be focused on this direction. 82 4 MICROPOROSITY MODEL 83 5 MACROPOROSITY AND PIPE SHRINKAGE MODEL 5 Macroporosity and Pipe Shrinkage Model 5.1 Introduction As some parts of the casting begin to solidify, a mass deficit is created as a consequence of the difference densities between the solid and liquid metal. When this mass feeding is cut, shrinkage porosity begins to form. It has been previously mentioned that when the cause of this cut is the growth of dendrites, this kind of porosity is called microporosity. When this effect happens in a larger scale, with isolated pools of liquid, porosity can be seen the naked eye and is referred as macroporosity. When the isolated pool of liquid is open to the air, porosity is usually called open shrinkage or pipe shrinkage because of its characteristic shape. Due to its impact on mechanical properties such as ductility, dynamic properties or fatigue life, macroporosity is considered a major cause of casting rejection. In many occasions, if the macro-shrinkage is not limited to the riser the casting is not accepted. In this section, a pressure-based model will be presented capable to compute macroporosity and pipe shrinkage. Following other approaches as the ones introduced in [8, 6, 27, 2], pressure is calculated using Darcy’s Law model along the whole casting and shrinkage is assumed to form in those points where pressure falls below a certain value. As this model is still being implemented, no numerical results or simulations are shown, which will be reserved for future work. 5.2 Governing Equations The present model states that each element is composed of a combination of air(a), porosity(p), liquid metal (l) and solid metal(s), such that the sum of the volume fractions is equal to 1. The air fraction acorresponds to the air that enters through the top of the riser and forms the open shrinkage. a+p+l+s= 1 (5.1) The mixture properties are obtained as a function of the properties of each phase 84 5 MACROPOROSITY AND PIPE SHRINKAGE MODEL multiplied by the corresponding volume fraction. ρ=aρa+pρp+lρl+sρs(5.2) Energy Conservation The energy conservation equation is written as: ρc −ρLds dT ∂T ∂t =∇·(λ∇T) (5.3) Where ρis the mixture density, cis the specific heat capacity, Tis the temperature and λis the conductivity. With this equation we can compute the temperature and the solid fraction s, as we suppose that the solid fraction is a known function of the temperature. Mass Conservation For the mass conservation we are only going to calculate the liquid phase, assuming that the porosity and the solid phase are stationary. Subtracting the air phase equation (∂(aρa)/∂t +∇·(aρava) = 0) to the continuity equation we get ∂ ∂t(pρp+lρl+sρs) + ∇·(ρlv) = 0 (5.4) where vis the superficial velocity v=lvl The last term can be rewritten as ∇·(ρlvs) = ρl∇·v+v·∇ρl(5.5) where the second term is very small and can be neglected. Introducing this modification into the previous equation and reordering terms 85 5 MACROPOROSITY AND PIPE SHRINKAGE MODEL ∇·v=1 ρl∂ ∂t [p(ρp−ρl) + s(ρs−ρl)−aρl+ρl](5.6) Momentum Conservation The momentum equation corresponds to the Stokes’ equation where an additional Darcy term has been added. As the velocity of the fluid and therefore the Re number are very small, the inertial term can be neglected as well as the convective term. ∇2v=l Kv+l µl∇P−l µl ρref g(5.7) Where gis the gravity, ρref is the density of the liquid at the liquidus temperature, µlis the dynamic viscosity and Kis the permeability given by: K=KrKo 3 l (1 −l)2(5.8) with Kr=l l+a (5.9) Ko= 6.10−4λ2 1(5.10) being λ1the primary dendrite arm. Notice that the momentum equation reduces to the Stokes equation when the permeability tends to infinite, that is, when liquid is the only present phase. Likewise, when the liquid fraction is very small, the left hand side term becomes very small and equation reduces to Darcy’s Law. 86 5 MACROPOROSITY AND PIPE SHRINKAGE MODEL 93 6 CONCLUSIONS AND FUTURE WORK 6 Conclusions and Future Work The objective of this work was the implementation of a model to compute porosity in order to implement it in the commercial software Click2Cast. In this last section we are going to summarize the different tasks that were carried out and the conclusions that we can extract from each of them as well as the current and future work that is being done in each specific area. In section 2 we did a review of the literature available on porosity and particularly, the different models that exist to compute each specific type of porosity. We classified them as follows •Microporosity Models : Criterion Functions •Macroporosity Models : Thermal/Volume Calculation Models •Coupled microporosity and macroporosity Models : Thermal/Fluid Flow Models and Stochastic Models Following this classification, we developed a Microporosity model based on Beckermann and Carlson’s paper from 2008 [5]. In order to do this we first had to implement a method to recover the gradient of the temperature in an accurate and efficient way. In section 3 we explained and tested 5 different methods to recover gradients: •L2 Lumped Projection •L2 Consistent Projection •Iterative L2 Projection •Zhang Method •Pouliot Method For each method, we obtained results for the error and the calculation time in order to be able to compare them. Although Pouliot method was proved to be the most accurate, we finally decided to go along with the iterative method which demonstrated its capacity to calculate gradients much faster than Pouliot with a reasonable accuracy. 94 6 CONCLUSIONS AND FUTURE WORK The microporosity model was explained in detail in section 4 along with results for different parts and materials. By doing a sensibility analysis of the mesh size, we showed that in order to have precise results a fine mesh is needed, in general finer than the mesh that we would require if only the solidification simulation was taken into account. Another sensibility analysis of the evaluation temperature also warned us about the importance of this parameter in the simulation. The tuning of this and other parameters such as the Secondary Dendrite Arm Spacing will be an important part of the future work. It would be desirable to have the chance to validate this model with some experimental results, although we know that due to the small size of these defects they cannot be detected by usual radio-graphic procedures and especial techniques are required. Once this validation work is finished, the current implementation will be included in future versions of Click2Cast. Finally, in the last part we studied how different parameters such as the mold temperature or the heat transfer coefficient affected the final results with satisfactory results that matched experimental findings. [22] In section 5, we presented a macroporosity model based on a simplified flow model, following the same approach as other methods like [6, 8, 19]. In this model, only the pressure was calculated throughout the solidification. When pressure falls under a certain value, pressure is no longer calculated and porosity is assumed to grow. This model is currently implemented but some work is still required to correct some wrong behavior under certain circumstances. We expect to solve this problems and validate the results in the next months in order to implement it along with the microporosity model in Click2Cast. 95 REFERENCES References [1] J. Beech, M. Barkhudarov, K. Chang, and s. B. Chin. Modeling of casting welding and advanced solidification processes. B.G Thomas and C. Beckermann, eds, The Minerals, Metals and Materials Soc. p.1071, 1998. [2] S. Bounds, G. Moran, K. Pericleous, M. Cross, and T.N. Croft. A computational model for defect prediction in shape castings based on the interaction of free surface flow, heat transfer, and solidification phenomena. Metallurgical and Materials Transactions B January 2000, 2000. [3] Kent D. Carlson and Christoph Beckermann. Use of the niyama criterion to predict shrinkage-related leaks in high-nickel steel and nickel-based alloy castings. Proceedings of the 62nd SFSA Technical and Operating Conference Paper No. 5.6 SFSA, Chicago, IL, 2008. [4] Kent D Carlson and Christoph Beckermann. Authors reply to discussion of prediction of shrinkage pore volume fraction using a dimensionless niyama criterion. Metallurgical and Materials TransactionsA, 40(13):3054 3055, 2009. [5] Kent D. Carlson and Christoph Beckermann. Prediction of shrinkage pore volume fraction using a dimensionless niyama criterion. Metall. Mater. Trans. A, Vol. 40A, pp. 163 175, 2009. [6] Kent D. Carlson, Zhiping Lin, Richard A. Hardin, and Christoph Beckermann. Modelling of porosity formation and feeding flow in steel casting. Proceedings of the 56th SFSA Technical and Operating Conference Paper No. 4.4, Chicago IL, 2002. [7] Kent D. Carlson, Shouzhu Ou, Richard Hardin, and Christoph Beckermann. Development of a methodology to predict and prevent leaks caused by microporosity in steel castings. Proceeding of the 55th technical and operating conference, SFSA, Chicago, 2001. [8] A.V. Catalina and C.A. Monroe. Simplified pressure model for quantitative shrinkage porosity prediction in steel castings. Materials Science and Engineering 33 (2012) 012067, 2012. 96 REFERENCES [9] John. P. Cornthwaite. Pressure Poisson Method for the Incompressible NavierStokes Equations Using Galerkin Finite Elements. PhD thesis, Georgia Southern University, 2003. [10] Pengfei Du. Numerical Modeling of Porosity and Macrosegregation in Continuous Casting of Steel. PhD thesis, University of Iowa, 2013. [11] P.N. Hansen and P.R. Sahm. How to model and simulate the feeding process in casting to predict shrinkage and porosity formation. Modeling of Casting and Welding Process IV, pages 33 42. TMS AIME, 1988. [12] R.A. Hardin, S. Ou, K. Carlson, and C. Beckermann. Relationship between casting simulation and radiographic testing: Results from the sfsa plate casting trials. 1999 SFSA Technical and Operating Conference, 1999. [13] D. M. Hawken, P. Townsend, and M. F. Webster. A comparison of gradient recovery methods in finite-element calculations. Communications in Applied Numerical Methods, VOl. 7, 195 204, 1991. [14] I. Imafuku and K. Chijiiwa. A mathematical model for shrinkage cavity prediction in steel castings. AFS Transactions 91: 527 540, 1983. [15] Ute Israel. Implementation of an algorithm for data transfer on the fluid-structure interface between non-matching meshes in kratos and an algorithm for the generation of an interface between gid and kratos. Master’s thesis, Technische Universit¨at M¨unchen, 2006. [16] Neelesh Jain, Kent D. Carlson, and Christoph Beckermann. Round robin study to assess variations in casting simulation niyama criterion predictions. Proceedings of the 61st SFSA Technical and Operating Conference Paper No. 4.4 SFSA, Chicago, IL, 2007. [17] J Jakumeit, S Jana, B B¨otger, R Laqua, M Y Jouani, and A B˜ AŒhrig Polaczek. Simulation-based prediction of micro-shrinkage porosity in aluminum casting: Fully-coupled numerical calculation vs. criteria functions. IOP Conference Series Materials Science and Engineering 27(1):2066, 2012. [18] Maodong Kang, Haiyan Gao, Jun Wang, Lishibao Ling, and Baode Sun. Prediction of microporosity in complex thin-wall castings with the dimensionless niyama criterion. Proceedings of the 8th Pacific Rim International Congress on Advanced Materials and Processing., 2013. 97 REFERENCES [19] Kimio Kubo and Robert D.Pehlke. Mathematical modeling of porosity formation in solidification. Metallurgical Transactions B, June 1985, Volume 16, Issue 2, pp 359 366, 1985. [20] P. D. Lee, A. Chirazi, and D. See. Modeling microporosity in aluminum-silicon alloys: a review. Journal of Light Metals 1 (2001) 15 30, 2001. [21] Y.W. Lee, E. Chang, and C.F. Chieu. Met. trans. b,. 1990, 21b, 715, 1990. [22] Y.M. Li and R.D. Li. Effect of the casting process variables on microporosity and mechanical properties in an investment cast aluminium alloy. Science and Technology of Advanced Materials 2, 2001. [23] E. Niyama, T. Uchida, M. Morikawa, and S. Saito. Predicting shrinkage in large steel castings from temperature gradient calculations. AFS Int. Cast Metals Journal, 6(2):16 22,, 1981. [24] E. Niyama, T. Uchida, M. Morikawa, and S. Saito. A method of shrinkage prediction and its application to steel casting practice. AFS Cast Metals Research J., 1982, 7, 52, 1982. [25] Zienkiewicz OC and Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. International Journal for Numerical Methods in Engineering 1987; 24:337357., 1987. [26] W.S. Pellini. Factors which determine riser adequacy and feeding range. AFS Transactions, 61: 61 81, 1953. [27] Ch. Pequet, M.Gremaud, and M. Rappaz. Modeling of microporosity, macroporosity, and pipe-shrinkage formation during the solidification of alloys using a mushy-zone refinement method: Applications to aluminum alloys. Metallurgical and Materials Transactions A, July 2002, Volume 33, Issue 7, pp 2095 2106, 2002. [28] B. Pouliot, M. Fortin, A. Fortin, and ´ E. Chamberland. On a new edge-based gradient recovery technique. Int. J. Numer. Meth. Engng (2012), 2012. [29] G. K. Sigworth. Discussion of prediction of shrinkage pore volume fraction using a dimensionless niyama criterion. Metallurgical and Materials Transactions A, 40(13):3051 3053, 2009. 98 REFERENCES [30] Doru M. Stefanescu. Computer simulation of shrinkage-related defects in castings - a review. Int. J. Cast. Met. Res. 18 129 143, 2005. [31] R. Tavakoli. On the prediction of shrinkage defects by thermal criterion functions. nt J Adv Manuf Technol (2014) 74: 569, 2009. [32] Zhimin Zhang and Ahmed Naga. A new finite element gradient recovery method: Superconvergence property. SIAM Journal on Scientific Computing 26(4):1192 1213, 2005. 99