A cell-by-cell thermally driven mushy cell tracking algorithm for phase-change problems

A thermally driven mushy cell tracking algorithm for phase-change problems with a moving boundary is presented. The equation used to track the moving boundary is based on energy balance over the mushy cell and is applied to advance a moving front in a cell-by-cell manner. The efficacy of the tracking algorithm is demonstrated on specific problems solved using the finite volume method. An implicit scheme is adopted to ensure that the numerical solution is unconditionally stable in time. A preconditioned conjugated gradient (P-CG) solver is implemented to ensure that solutions converge in a finite number of iterations. Four benchmark cases are used to validate the algorithm including solidification in one dimensional space (two-region problem), melting of pure aluminum in two-dimensional (2D) space, solidification with periodic boundary conditions, and solidification of one-region problem. The results obtained show that the current algorithm is capable of converging to accurate solutions for moving fronts and the numerical predications are in excellent agreement with corresponding analytical solutions.     

A numerical simulation for two-dimensional moving boundary problems with a mushy zone

For materials such as alloy, organic phase-change materials and many others, the change of phases may take place over a temperature range. This leads to phase-change problems with the mushy zone in which the solid and liquid phases coexist. The present study introduces a numerical method combining the Laplace transform technique and the control volume method to solve two-dimensional phase-change problems with the mushy zone. The hybrid numerical method involves the control volume formulation for the space domain and the Laplace transform technique for the time domain. The Taylor's series approximation is applied to linearize nonlinear terms in the governing equation. The transfinite mapping method is used to generate control-volume meshes in each region. The growth of the mushy zone is unknown a priori and is predicted by using the least-square iteration scheme. It will be found that the present hybrid numerical method can be efficiently applied to solve two-dimensional phase-change problems with a mushy zone.     

A least-squares front-tracking finite element method analysis of phase change with natural convection

This paper focuses on the numerical modelling of phase-change processes with natural convection. In particular, two-dimensional solidification and melting problems are studied for pure metals using an energy preserving deforming finite element model. The transient Navier–Stokes equations for incompressible fluid flow are solved simultaneously with the transient heat flow equations and the Stefan condition. A least-squares variational finite element method formulation is implemented for both the heat flow and fluid flow equations. The Boussinesq approximation is used to generate the bulk fluid motion in the melt. The mesh motion and mesh generation schemes are performed dynamically using a transfinite mapping. The consistent penalty method is used for modelling incompressibility. The effect of natural convection on the solid/liquid interface motion, the solidification rate and the temperature gradients is found to be important. The proposed method does not possess some of the false diffusion problems associated with the standard Galerkin formulations and it is shown to produce accurate numerical solutions for convection dominated phase-change problems.     

Modified enthalpy method for the simulation of melting and solidification

Enthalpy method is commonly used in the simulation of melting and solidification owing to its ease of implementation. It however has a few shortcomings. When it is used to simulate melting/solidification on a coarse grid, the temperature time history of a point close to the interface shows waviness. While simulating melting with natural convection, in order to impose no-slip and impermeability boundary conditions, momentum sink terms are used with some arbitrary constants called mushy zone constants. The values of these are very large and have no physical basis. Further, the chosen values affect the predictions and hence have to be tuned for satisfactory comparison with experimental data. To overcome these deficiencies, a new cell splitting method under the framework of the enthalpy method has been proposed. This method does not produce waviness nor requires mushy zone constants for simulating melting with natural convection. The method is then demonstrated for a simple one-dimensional melting problem and the results are compared with analytical solutions. The method is then demonstrated to work in two-dimensions and comparisons are shown with analytical solutions for problems with planar and curvilinear interfaces. To further benchmark the present method, simulations are performed for melting in a rectangular cavity with natural convection in the liquid melt. The solid–liquid interface obtained is compared satisfactorily with the experimental results available in literature.     

Study on Heat Transfer Characteristics of Phase-Change Energy Storage Unit for Thermal Management

The objective of the study was to investigate the heat transfer characteristics of a phase-change energy storage unit for thermal management. Considering the conduction in the solid and natural convection in the liquid, a physical and mathematical model for heat transfer was formulated. The governing conservation equations were solved using the finite-volume method on fixed grids. An enthalpy-porosity method was used for modeling the melting phenomenon of a phase-change energy storage unit. The time and space movement of the phase front, the temperature distribution, and the heat dissipation rate have been analyzed based on the model. The influence of the unit geometry, heat source location, and types of phase-change materials on the thermal performance of the energy storage unit were investigated. The model and numerical method were evaluated by comparing the numerical predictions with the experimental results. There was found to be excellent agreement between the calculation and experiment, indicating that the numerical method for heat transfer simulation of a phase-change energy storage unit is accurate. The results from the analysis elucidate the thermal performance of the phase-change energy storage unit and will provide the basis for the design and optimization of thermal management systems.     

A fixed-grid numerical modelling of transient liquid phase bonding and other diffusion-controlled phase changes

A transient liquid phase (TLP), in which a liquid layer is formed and subsequently solidifies, and other diffusion-controlled phase changes are generally associated with moving phase-change interfaces. Both fixed and variable grid discretization models have been formulated to investigate these diffusion-controlled problems. However, all numerical efforts to date have employed one of the approaches explicitly to track the moving interfaces across which there exist step changes in concentrations. In this article, the fixed-grid source-based method originally developed to simulate the temperature fields for melting-solidification phase change processes has been adopted to simulate diffusion-controlled dissolution and solidification. This method solves a unique diffusion equation for the different phases and the moving interfaces using implicit time integration. Compared with previously developed models, it is not only simpler in numerical formulation and procedure, but also more convenient to extend to many phases and high-dimensional problems. We report here the detailed formulation of the relevant equations, and compare and validate the model using experimental data and previous modelling predictions for several systems available from the existing literature.     

An approximate method for evaporation problem with mushy zone

In the present paper, a simple mushy zone model is used to track the moving boundaries in an evaporation problem in which the vapor is removed upon formation. Two main parameters for the mushy zone model are analyzed as well as their effect on the movement of the moving boundaries and the thickness of the mushy zone. A new approximate method is developed for analysis and tracking the moving boundaries appears throughout the process. The proposed method mainly based on applying the boundary integral equation corresponding to each phase in such a way that the associated boundary and initial conditions as well as energy equations at the moving boundaries achieved with minimum error and low number of iterations. The results of the present paper seem to be good because there are neither analytical or numerical solutions available.     

A mixed-variable continuously deforming finite element method for parabolic evolution problems. Part II: The coupled problem of phase-change in porous media

The conceptual framework of a least squares rate variational approach to the formulation of continuously deforming mixed-variable finite element computational scheme for a single evolution equation was presented in Part I.¹ In this paper (Part II), we extend these concepts and present an adaptively deforming mixed variable finite element method for solving general two-dimensional transport problems governed by a system of coupled non-linear partial differential evolution equations. In particular, we consider porous media problems that involve coupled heat and mass transport processes that yield steep continuous moving fronts, and abrupt, discontinuous, moving phase-change interfaces. In this method, the potentials, such as the temperature, pressure and species concentration, and the corresponding fluxes, are permitted to jump in value across the phase-change interfaces. The equations, and the jump conditions, governing the physical phenomena, which were specialized from a general multiphase, multiconstituent mixture theory, provided the basis for the development and implementation of a two-dimensional numerical simulator. This simulator can effectively resolve steep continuous fronts (i.e. shock capturing) without oscillations or numerical dispersion, and can accurately represent and track discontinuous fronts (i.e. shock fitting) through adaptive grid deformation and redistribution. The numerical implementation of this simulator and numerical examples that demonstrate the performance of the computational method are presented in Part III² of this paper.     

Numerical modeling of convection-diffusion phase change problems

A numerical methodology is presented for the modeling of convection-diffusion controlled mushy region change problems. An efficient and accurate non-staggered control volume method, based on the momentum interpolation practice and on a high-order convection differencing scheme, is proposed for the solution of the continuum model equation. Suitable numerical techniques are implemented to overcome the numerical instability problems resulting from the strong coupling between the equations of the model. Special attention is given on the efficient treatment of the latent head evolution in the energy equation. A new numerical technique is developed which accounts for the dependence of the latent heat on the variation of temperature and concentration fields. The proposed method is applied on two phase change problems. Satisfactory agreement with previously published results is observed.     

On finite element modelling of surface tension: Variational formulation and applications – Part II: Dynamic problems

In Part I, a finite element model of surface tension has been discussed and used to solve some quasi-static problems. The quasi-static analysis is often required to find not only the initial shape of the liquid but also the static equilibrium state of a liquid body before a dynamic analysis can be carried out. In general, natural and industrial processes in which surface tension force is dominant are of dynamic nature. In this second part of this work, the dynamic effects will be included in the finite element model described in Part I.A fully Lagrangian finite element method is used to solve the free surface flow problem and Newtonian constitutive equations describing the fluid behaviour are approximated over a finite time interval. As a result the momentum equations are function of nodal position instead of velocities. The resulting ordinary differential equation is integrated using Newmark algorithm. To avoid overly distorted elements an adaptive remeshing strategy is adopted. The adaptive strategy employs a remeshing indicator based on viscous dissipation functional and incorporates an appropriate transfer operator.The validation of the model is performed by comparing the finite element solutions to available analytical solutions of a droplet oscillations and experimental results pertaining to stretching of a liquid bridge.     

Numerical simulation of cryogenic cavitating flow by an extended transport-based cavitation model with thermal effects

Thermodynamic effects on cryogenic cavitating flow is important to the accuracy of numerical simulations mainly because cryogenic fluids are thermo-sensitive, and the vapour saturation pressure is strongly dependent on the local temperature. The present study analyses the thermal cavitating flows in liquid nitrogen around a 2D hydrofoil. Thermal effects were considered using the RNG k-ε turbulence model with a modified turbulent eddy viscosity and the mass transfer homogenous cavitation model coupled with energy equation. In the cavitation model process, the saturated vapour pressure is modified based on the Clausius-Clapron equation. The convection heat transfer approach is also considered to extend the Zwart-Gerber-Belamri model. The predicted pressure and temperature inside the cavity under cryogenic conditions show that the modified Zwart-Gerber-Belamri model is in agreement with the experimental data of Hord et al. in NASA, especially in the thermal field. The thermal effect significantly affects the cavitation dynamics during phase-change process, which could delay or suppress the occurrence and development of cavitation behaviour. Based on the modified Zwart-Gerber-Belamri model proposed in this paper, better prediction of the cryogenic cavitation is attainable.     

A mixed-variable continuously deforming finite element method for parabolic evolution problems. Part I: The variational formulation for a single evolution equation

The objective of the research presented here was to develop a generic adaptive computational method for porous media evolution problems that involve coupled heat flow, fluid flow and species transport processes with sharply defined phase-change interfaces. In this paper we examine the general least squares variational approach and develop the conceptual framework for a rate least squares variational formulation of a continuously deforming mixed variable finite element method for solving highly non-linear time-dependent partial differential equations. In Part II of this paper¹ we extend the formulation given here for a single evolution equation to a system of coupled evolution equations. In Part III² we discuss in detail the numerical procedures that were implemented in a computer program and present several numerical examples that demonstrate the performance of this computational method.     

Low-velocity impact-induced damage of continuous fiber-reinforced composite laminates. Part I. An FEM numerical model

A numerical model for simulating the process of low-velocity impact damage in composite laminates using the finite element method is presented in this paper, i.e. Part I of this two part series on the study of impact. In this model, the 9-node Lagrangian element of the Mindlin plate with consideration of large deformation analysis is employed. To analyze the transient response of the laminated plates, a modified Newmark time integration algorithm previously proposed by the authors is adopted here. We also proved that the impact process between a rigid ball and laminated plates is a stiff system, therefore a kind of A(α) stable method has been advocated here to solve the motion equation of the rigid ball. Furthermore, various types of damages including delamination, matrix cracking and fiber breakage, etc. and their mutual influences are modeled and investigated in detail. To overcome the difficulty of numerical oscillation or instability in the analysis of the dynamic contact problem between delaminated layers using the traditional penalty methods, we have employed dynamic spring constraints to simulate the contact effect, which are added to the numerical model by a kind of continuous penalty function. Moreover, an effective technique to calculate the strain energy release rate based on the Mindlin plate model is proposed, which can attain high precision. Finally, some techniques of adaptive analyses have been realized for improving the computational efficiency. Based on this model, a program has been developed for numerically simulating the damage process of cross-ply fiber-reinforced carbon/epoxy composite laminates under low-velocity impact load. In Part II, this numerical model will be verified by comparing with the experimental results. Also the impact damage will be investigated in detail using this numerical approach.     

Flow instabilities of liquid and mushy regions during alloy solidification and under high gravity environment induced by rotation

Linear flow instabilities of the liquid and mushy regions during directional solidification of a binary alloy are studied under a high gravity environment where the rotation axis is inclined with respect to the high gravity vector. Stability analysis and numerical computation are carried out to determine the results for the stationary disturbances at several values of the rotation rates and for given values of the other parameters. The results provide information about the effects of Coriolis force on various flow features in the liquid and mushy layers. The preferred structure of the mush-liquid interface is found to be that of longitudinal rolls. The main mode of convection is found to be able to generate double-cell structure in the vertical direction and is strengthened in the mushy layer by the Coriolis force. The Coriolis force appears to be generally stabilizing in the sense that the motion in the liquid zone is significantly weakened, the tendency for the chimney formation in the mushy zone is reduced and the critical values of the liquid and mush Rayleigh numbers and the wave number increase with increasing the rotation rate.     

Measurements of the Thermal Conductivity and Thermal Diffusivity of Liquids. Part II: “Convective and Radiative Effects”

The second part (Part II) of this work is concerned with coupling in the transient regime of conduction with convection and radiation in the experimental bench developed and presented in Part I for the measurement of the thermal conductivity and thermal diffusivity of fluids by an impulse technique. The first section will analyze heat transfer in the liquid by conduction and convection. This will help to define the optimal extension of the measuring cell to reduce the influence of natural convection for the case of impulse heat flux stimulation. The second section is about coupled conductive–radiative heat transfer and will show how to deal with radiative effects in the problem of parameters estimation.Paper presented at the Seventh Asian Thermophysical Properties Conference, August 23–28, 2004, Hefei and Huangshan, Anhui, P. R. China.     

Meshless simulations of the two-dimensional fractional-time convection–diffusion–reaction equations

The aim of this work is to propose a numerical approach based on the local weak formulations and finite difference scheme to solve the two-dimensional fractional-time convection–diffusion–reaction equations. The numerical studies on sensitivity analysis to parameter and convergence analysis show that our approach is stable. Moreover, numerical demonstrations are given to show that the weak-form approach is applicable to a wide range of problems; in particular, a forced-subdiffusion–convection equation previously solved by a strong-form approach with weak convection is considered. It is shown that our approach can obtain comparable simulations not only in weak convection but also in convection dominant cases. The simulations to a subdiffusion–convection–reaction equation are also presented.     

Thermal convection in the presence of a vertical magnetic field

Summary The interaction between thermal convection and an external uniform magnetic field in the vertical is numerically simulated within a computational domain of a horizontally periodic convective box between upper and lower rigid plates. The numerical technique is based on a spectral element method developed earlier to simulate natural thermal convection. In this work, it is extended to a magnetoconvection problem. Its main features are the use of rescaled Legendre-Lagrangian polynomial interpolants in expanding the flow variables except the pressure for which a modal expansion in terms of lower order polynomials is used to avoid the complicated staggered grid approach. The technique is validated in the steady roll and oscillatory convective regimes where various experimental and numerical results are available in the literature. The effect of a vertical magnetic field in such a way to inhibit the convective motions has been demonstrated.     

Compressible flow analysis of filling and postfilling in injection molding with phase-change effect

In order to predict the shrinkage, warpage and mechanical properties of the injection molded parts, it is necessary to know the history of the flow field during injection-molding processes. In the present investigation a numerical simulation program was developed to predict the flow field in filling and post-filling stages of injection molding. To simulate the real molding conditions more accurately, a generalized Hele-Shaw model for a non-Newtonian fluid was assumed considering the effects of phase change and compressibility of the resin. A finite-element-finite-difference (FEM-FDM) hybrid scheme with control volume approach was employed as the solving technique. For modeling the viscosity of the resin, a modified Cross model was used with a double-domain Tait equation of state being employed in describing the compressibility of the resin during molding. The energy balance equation, including latent-heat dissipation for semicrystalline materials, was solved in order to predict the solidified layer and temperature profile in detail. For verification of the numerical results obtained from the developed program, the simulation results were compared with the experimental results obtained from the test mold set designed in the current study using commercial-grade PP and the data available in the literature. Based on a comparison between experiments and simulations, it was found that the currently developed program was useful in unified simulations of filling and post-filling in injection-molding processes when considering the phase-change effect.     

Oscillatory instabilities of the liquid and mushy layers during solidification of alloys under rotational constraint

Summary Linear flow instabilities due to oscillatory disturbances of the liquid and mushy regions during solidification of binary alloys are investigated under a rotational constraint where the rotation axis is inclined to gravity vector. Results of stability analyses and numerical computations for a preferred centrifugal mode of general oscillatory disturbances at zero and non-zero rotation rates are determined which provide information about the preference of oscillatory flow and its role on the solidification system as modified by the rotational effects. The main results are due to a preferred oscillatory mode of convection which is more significant for non-zero rotation case and is restricted mostly to the mushy region. The preferred oscillatory mode of convection is a traveling wave in the presence of rotation, but it is a standing wave in the absence of rotation. The results for different Prandtl numbers indicate that the freckles formation tendency for metallic alloys is less than that for aqueous solutions. Freckles are imperfections that reduce the quality of the solidified materials.     

Modelling convection in solidification processes using stabilized finite element techniques

Solidification of dendritic alloys is modelled using stabilized finite element techniques to study convection and macrosegregation driven by buoyancy and shrinkage. The adopted governing macroscopic conservation equations of momentum, energy and species transport are derived from their microscopic counterparts using the volume‐averaging method. A single domain model is considered with a fixed numerical grid and without boundary conditions applied explicitly on the freezing front. The mushy zone is modelled here as a porous medium with either an isotropic or an anisotropic permeability. The stabilized finite‐element scheme, previously developed by authors for modelling flows with phase change, is extended here to include effects of shrinkage, density changes and anisotropic permeability during solidification. The fluid flow scheme developed includes streamline‐upwind/Petrov–Galerkin (SUPG), pressure stabilizing/Petrov–Galerkin, Darcy stabilizing/Petrov–Galerkin and other stabilizing terms arising from changes in density in the mushy zone. For the energy and species equations a classical SUPG‐based finite element method is employed with minor modifications. The developed algorithms are first tested for a reference problem involving solidification of lead–tin alloy where the mushy zone is characterized by an isotropic permeability. Convergence studies are performed to validate the simulation results. Solidification of the same alloy in the absence of shrinkage is studied to observe differences in macrosegregation. Vertical solidification of a lead–tin alloy, where the mushy zone is characterized by an anisotropic permeability, is then simulated. The main aim here is to study convection and demonstrate formation of freckles and channels due to macrosegregation. The ability of stabilized finite element methods to model a wide variety of solidification problems with varying underlying phenomena in two and three dimensions is demonstrated through these examples. Copyright © 2005 John Wiley & Sons, Ltd.