Stefan Problem Research Articles

Influence of Ground Ice Content on Thawing Depth Road Base

One of the important parameters that determine technical solutions in the design of roads in the permafrost zone is the depth of thawing of the road foundation soils. The aim of the work was to quantify the degree of influence of the ice content in the soil foundations of roads in the permafrost zone on the depth of their seasonal thawing. For the analysis, we used the classical formula for calculating the thawing depth for bodies of plane symmetry, obtained by solving the single-phase Stefan problem. The results of numerical calculations are presented in the form of 2D and 3D graphs, which make it possible to visually assess the effect of ground ice content and the degree of its change during the road operation period on the depth of thawing of the road foundation. It has been established, in particular, that the degree of change in the depth of thawing at the same value of the increase in ice content in different ranges (for example, from 10 to 20% and from 30 to 40%) for the considered typical case of the ground base of the permafrost zone, decreases by almost 1.3 times. It is shown that the greater the initial ice content of the soil, the lower the degree of decrease in the thawing depth when the ice content changes by a constant value.

Numerical simulation of selective laser melting by the SPH method

Abstract. Currently, additive manufacturing technologies develop actively. This requires creation of computational methods to describe physical processes occurring at the time of manufacturing. One of the methods used for the production of metal powder parts is the method of selective laser melting. This paper presents an SPH-based numerical technique for modeling the process of powder sintering under the influence of a laser beam. The flow of liquid formed as a result of melting is described by the Navier-Stokes equations. Pressure forces, viscous effects and surface forces at the interface are included in the force balance. The thermal state is determined from the energy conservation law, which takes into account thermal processes, volumetric absorption of laser radiation energy, convective heat exchange with the external environment and radiation. Phase transitions between solid and liquid phases are described in the framework of the generalized formulation of the Stefan problem. The calculation method is verified on tests specific to the class of problems under consideration. A comparison is made with the analytical solution, as well as with solutions obtained by other modifications of the SPH method, and with experimental data.

Adjoint-based optimization of two-dimensional Stefan problems

A range of optimization cases of two-dimensional Stefan problems, solved using a tracking-type cost-functional, is presented. A level-set method is used to capture the interface between the liquid and solid phases and an immersed boundary (cut cell) method coupled with an implicit time-advancement scheme is employed to solve the heat equation. A conservative implicit-explicit scheme is then used for solving the level-set transport equation. The resulting numerical framework is validated with respect to existing analytical solutions of the forward Stefan problem. An adjoint-based algorithm is then employed to efficiently compute the gradient used in the optimization algorithm (L-BFGS). The algorithm follows a continuous adjoint framework, where adjoint equations are formally derived using shape calculus and transport theorems. A wide range of control objectives are presented, and the results show that using parameterized boundary actuation leads to effective control strategies in order to suppress interfacial instabilities or to maintain a desired crystal shape.

Analysis of a fully discrete approximation to a moving-boundary problem describing rubber exposed to diffusants

We present a fully discrete scheme for the numerical approximation of a moving-boundary problem describing diffusants penetration into rubber. Our scheme utilizes the Galerkin finite element method for the space discretization combined with the backward Euler method for the time discretization. Besides dealing with the existence and uniqueness of solution to the fully discrete problem, we assume sufficient regularity for the solution to the target moving boundary problem and derive a a priori error estimates for the mass concentration of the diffusants, and respectively, for the position of the moving boundary. Our numerical results illustrate the obtained theoretical order of convergence in physical parameter regimes.

Convectively cooled solidification in phase change materials in different configurations subject to internal heat generation: Quasi-steady analysis

Phase change phenomena involving internal heat generation are of central interest in various applications, including freeze drying, preservation of biological tissues, microwave/Joule heating processes, materials processing using lasers, and nuclear power generation systems. In this work, mathematical models are developed to study the solidification processes driven by convective cooling effects on the exterior surface of phase change materials (PCM) in different geometric configurations, where both the solid and the liquid phases are subject to volumetric heating. By invoking a quasi-steady assumption and performing a simplified unified analysis, new semi-analytical approximate solutions for this unique Stefan problem in planar, cylindrical, spherical, and semi-infinite geometries of the PCM are obtained and studied. The interface evolution depends on the characteristic dimensionless parameters, such as the Biot number, Bi, internal heat generation parameter, Q, Stefan number, St, and the thermal conductivity ratio, κ. The effects of all the dimensionless parameters, viz., Bi, Q, St, and κ are systematically investigated for each geometrical configuration of the PCM. In particular, the results show the significant role of the Biot number on the transient evolution as well as the steady-state thickness of the solidifying front in the presence of an inhibiting internal heat generation. In addition, a simple and unified condition relating the dimensionless internal heating rate (Q) and the external convective cooling rate (Bi) for the onset of remelting for the special case of κ = 1 applicable for all the above geometries of the PCM is deduced.

Simulation of Phase Change Material Absorbers for Passive Cooling of Solar Systems

One of the main challenges that face the reliable use of photovoltaic solar systems in hot arid regions is the prevailing high temperatures during the day. To overcome this issue, Phase Change Materials (PCM) are used for passive cooling providing different options to attain sufficient thermal management solutions for different applications. Passive cooling can be achieved by adjusting a heat sink to the solar PV module. This can be realized by attaching a PCM layer or sensible heat storage to the backside of PV panels. Few studies have reported on simplified modeling and numerical procedures using the apparent heat capacity formulation and volume averaging technique as a robust approach to solving such sophisticated problems with minimal computational efforts, high accuracy, and in a short period of time. However, there is still a need to bridge the large-scale gap between the macroscale within the PCM layer, with a moving melting front, and the length scale of PV modules. Hence, this work focuses on modeling and simulating PCM-Matrix Absorbers (PCM-MA) that consist of fibrous aluminum cellular structure filled with PCM for passive thermal management of PV panels using apparent heat capacity formulation and homogenization based on volume averaging technique. COMSOL Multiphysics FEM software was used for the numerical simulation of the phase change problem by using a Smoothed Heaviside step function to overcome the singularity of PCM challenge that arises with sharp melting temperatures. To validate the proposed model, it has been compared with a benchmark analytical solution for an ice melting problem, i.e., the Stefan problem, in a semi-finite slab, i.e., Neumann’s solution under the same assumptions and boundary conditions. The specific characteristics of phase change and evolution of melting front with time, heat capacity change with the temperature at different times, and with locations along the slab height are presented. As the phase change is modeled to take place over a mushy region, i.e., a narrow temperature interval rather than a sharp melting point, the results show a good coincidence of the heat capacity profile and its peak at different times and locations. The validated model can be used for the optimization of PCM-MA for any specific geographical location and other applications such as the passive cooling of buildings with PCM integrated with the outer envelope. To this end, the results of the simulation in this work are shown to be in agreement with those obtained from the analytical solution.

A hybrid level-set / embedded boundary method applied to solidification-melt problems

In this paper, we introduce a novel way to represent the interface for two-phase flows with phase change. We combine a level-set method with a Cartesian embedded boundary method and take advantage of both. This is part of an effort to obtain a numerical strategy relying on Cartesian grids allowing the simulation of complex boundaries with possible change of topology while retaining a high-order representation of the gradients on the interface and the capability of properly applying boundary conditions on the interface. This leads to a two-fluid conservative second-order numerical method. The ability of the method to correctly solve Stefan problems, onset dendrite growth with and without anisotropy is demonstrated through a variety of test cases. Finally, we take advantage of the two-fluid representation to model a Rayleigh–Bénard instability with a melting boundary.

Modeling of dendritic solidification and numerical analysis of the phase-field approach to model complex morphologies in alloys

Dendrites are one of the most widely observed patterns in nature, and occur across a wide spectrum of physical phenomena—from snow flakes to river basins; from bacterial colonies to lungs and vascular systems; and in solidification and growth patterns in metals and crystals. The ubiquitous occurrence of these “tree-like” structures can be attributed to their excellent space-filling properties, and at times, dendritic structures also spatially manifest fractal-like distributions. As is the case with many fractal-like geometries, the complex multi-level branching structures in dendrites pose a modeling challenge, and a full resolution of dendritic structures is computationally very demanding. In the literature, extensive theoretical models of dendritic formation and evolution, essentially as extensions of the classical moving boundary Stefan problem exist. Much of this understanding is from the analysis of dendrites occurring during the solidification of metallic alloys, as this is critical for understanding microstructure evolution during metal manufacturing processes that involve solidification of a liquid melt. Motivated by the problem of modeling microstructure evolution from liquid melts of pure metals and metallic alloys during metal additive manufacturing, we developed a comprehensive numerical framework for modeling a large variety of dendritic structures that are relevant to metal solidification. In this work, we present a numerical framework encompassing the modeling of Stefan problem formulations relevant to dendritic evolution using a phase-field approach and a finite element method implementation. Using this framework, we model numerous complex dendritic morphologies that are physically relevant to the solidification of pure melts and binary alloys. The distinguishing aspects of this work are—a unified treatment of both pure metals and alloys; novel numerical error estimates of dendritic tip velocity; and the study of error convergence of the primal fields of temperature and the order parameter with respect to numerical discretization. To the best of our knowledge, this is a first of its kind study of numerical convergence of the phase-field equations of dendritic growth in a finite element method setting. Further, using this numerical framework, various types of physically relevant dendritic solidification patterns like single equiaxed, multi-equiaxed, single columnar and multi-columnar dendrites are modeled in two-dimensional and three-dimensional computational domains.

Butterflies

Butterflies

An approximate analytical solution for the movement of the phase change front in latent thermal energy storage heat exchangers

Latent thermal energy storage (LTES) heat exchangers are being applied in a wide range of energy systems. However, there is no analytical method to determine the outlet temperature of LTES heat exchangers from its operational conditions. This hinders experimental data reduction and sizing of these components for specific energy systems. The present paper develops an analytical model for the overall phase fraction and outlet heat transfer fluid temperature for planar LTES heat exchangers. The method is developed with solutions of the Stefan problem. The analytical results are validated based on a numerical model using the same assumptions as the analytical model. The analytical solution can be reduced to a limited set of parameters which are a dimensionless time and heat exchanger surface area. This limited set of parameters can form the basis for experimental data reduction of LTES heat exchangers.

Pattern formation and front stability for a moving-boundary model of biological invasion and recession

We investigate pattern formation in a two-dimensional (2D) Fisher–Stefan model, which involves solving the Fisher–KPP equation on a compactly-supported region with a moving boundary. The boundary evolves analogously to the classical one-phase Stefan problem, such that boundary speed is proportional to the local population density gradient. By combining the Fisher–KPP and classical Stefan theory, the Fisher–Stefan model alleviates two limitations of the Fisher–KPP equation for biological populations. Unlike the Fisher–KPP equation, solutions to the Fisher–Stefan model have compact support, explicitly defining the region occupied by the population. Furthermore, the Fisher–Stefan model admits travelling wave solutions with non-negative density for all wave speeds c∈(−∞,∞), and can thus model both population invasion and population recession. In this work, we investigate whether the 2D Fisher–Stefan model predicts pattern formation, by analysing the linear stability of planar travelling wave solutions to sinusoidal transverse perturbations. Planar fronts of the Fisher–KPP equation are linearly stable. Similarly, we demonstrate that invading planar fronts (c>0) of the Fisher–Stefan model are linearly stable to perturbations of all wave numbers. However, our analysis demonstrates that receding planar fronts (c<0) of the Fisher–Stefan model are linearly unstable for all wave numbers. This is analogous to unstable solutions for planar solidification in the classical Stefan problem. Introducing a surface tension regularisation stabilises receding fronts for short-wavelength perturbations, giving rise to a range of unstable modes and a most unstable wave number. We supplement linear stability analysis with level-set numerical solutions that corroborate theoretical results. Overall, front instability in the Fisher–Stefan model suggests a new mechanism for pattern formation in receding biological populations.

A generalized thermoelastic Stefan problem with spatio-temporal nonlocal effect in picosecond pulse laser material processing

The Stefan conditions are suitable for describing the solid-liquid phase change processes with moving interface. However, the existing literature only focuses on the Stefan problem of heat conduction and does not involve the thermo-mechanical coupling effects. In this work, the generalized thermoelastic Stefan problem with spatio-temporal nonlocal effect is investigated. The solid-liquid phase change of silicon with picosecond pulse laser machining is depicted by the Stefan interface conditions. We find the solutions making use of Laplace transformation and Laplace numerical inverse transformation. The temperature distribution at the moment of surface melting is seen as the initial state of the second phase. The distributions of transient thermoelastic responses under different laser intensities, pulse widths, fractional order parameters, spatial nonlocal parameters and time are obtained and expressed graphically.

Development of a modified mass transfer model based on height function method to capture the pressure oscillation occurred in direct contact condensation process

For direct contact condensation (DCC) simulation, it is a challenge to capture the numerical characteristic of pressure oscillation precisely, especially the frequency and amplitude. It is the most crucial key issue to calculate the mass transfer across the vapor-liquid interface. Up to now, the pressure characteristic is generally obtained by experiment and the existing mass transfer models always lead to unacceptable deviation from the experiment result. In order to estimate the accurate pressure characteristic, a modified mass transfer model is originally proposed based on the height function method in this paper. An estimation method for the interfacial area density in this model is derived from the algebraic spherical model with the mean bubble diameter. This method correlates the interfacial area concentration with the measurable interfacial curvature and the interfacial curvature is calculated by the height function method. The proposed model is then employed based on the VOF method for the verification of Stefan's analytical problem and steam bubble condensation problem. The corresponding numerical results are both in great agreement with the analytical and experimental results. In addition, the proposed model achieves its advantage for the steam flow condensation to obtain more accurate results of pressure oscillation than other mass transfer models. The numerical result of pressure oscillation frequency is 8.9 Hz, close to the experimental data of 10 Hz with a deviation of 11%. However, the energy balance model and Lee model predict that the frequencies are 31.5 Hz and 19 Hz relatively with the corresponding deviation of 215% and 90%. It is validated to be convincing with excellent accuracy and shows promise for extensive numerical simulation for direct contact condensation applications.

A study of dual phage lag model of heat transfer for solidification of binary eutectic alloy

In this paper, we developed a generalized non-Fourier model of energy transfer for solidification of a binary eutectic system (Al-Cu alloy). The mathematical model is divided into three stages: stage I (liquid region), II (liquid and mushy region), and III (liquid, mushy and solid region). The analysis of the model has been done in dimensionless form. This model is a boundary value problem of the system of second-order hyperbolic partial differential equations for Stage I, the moving boundary problem of the System of second-order hyperbolic partial differential equations for Stage II and III. For the solution, we applied Laplace transform technique in Stage I and Legendre wavelet spectral Galerkin method in Stage II and III. The effect of unit-less parameters: Relaxation time, Predvoditelev numbers, dimensionless thermal conductivity, Fourier numbers, Stefan number, solid fraction on temperature, and position of moving interface are discussed in detail.

Exact solution for non-classical one-phase Stefan problem with variable thermal coefficients and two different heat source terms

Exact solution for non-classical one-phase Stefan problem with variable thermal coefficients and two different heat source terms

Numerical analysis of fully non-linear sloshing waves in an arbitrary shape tank by meshless method

This paper is devoted to develop a truly meshless numerical procedure for the fully non-linear analysis of sloshing phenomenon in an arbitrary-shape tank. For this purpose, the potential theory is considered as a liquid sloshing model, and the Lagrangian form of kinematic and dynamic boundary conditions in a moving coordinate system fixed to the tank is used to accurately capture the free surface. The multi-quadric radial basis function (MQ-RBF) augmented with the polynomial terms is employed to determine RBF-FD weights on neighboring nodes surrounding the center point for the spatial derivatives. These weight coefficients are used to solve both boundary value problems for velocity and acceleration potential. Compared with mash-based methods, the present method has the advantage of being easy to construct and flexible in dealing with the moving boundary problems. The free surface elevations, hydrodynamic pressure, and first-order natural sloshing frequencies have been reported for various geometries tanks under horizontal excitations at resonance and non-resonance cases. Fairly satisfactory agreement is observed in the numerical results and exiting results, especially for free surface elevation and first-order natural sloshing frequency.

Fractional Stefan Problem Solving by the Alternating Phase Truncation Method

The aim of this paper is the adaptation of the alternating phase truncation (APT) method for solving the two-phase time-fractional Stefan problem. The aim was to determine the approximate temperature distribution in the domain with the moving boundary between the solid and the liquid phase. The adaptation of the APT method is a kind of method that allows us to consider the enthalpy distribution instead of the temperature distribution in the domain. The method consists of reducing the whole considered domain to liquid phase by adding sufficient heat at each point of the solid and then, after solving the heat equation transformed to the enthalpy form in the obtained region, subtracting the heat that has been added. Next the whole domain is reduced to the solid phase by subtracting the sufficient heat from each point of the liquid. The heat equation is solved in the obtained region and, after that, the heat that had been subtracted is added at the proper points. The steps of the APT method were adapted to solve the equations with the fractional derivatives. The paper includes numerical examples illustrating the application of the described method.

Numerical difference solution of moving boundary random Stefan problems

This paper deals with the construction of numerical solutions of moving boundary random problems where the uncertainty is limited to a finite degree of randomness in the mean square framework. Using a front fixing approach the problem is firstly transformed into a fixed boundary one. Then a random finite difference scheme for both the partial differential equation and the Stefan condition, allows the discretization. Since statistical moments of the approximate stochastic process solution are required, we combine the sample approach of the difference schemes together with Monte Carlo technique to perform manageable approximations of the expectation and variance of both the approximating stochastic process solution and the stochastic moving boundary solution. Qualitative and reliability properties such as positivity, monotonicity and stability in the mean square sense are treated. Feasibility of the proposed method is checked with illustrative examples of a melting problem and a binary metallic alloys problems.

Rationale Behind EU-DEMO Limiter’s Plasma-Facing Component Design Under Material Phase Change

The protection strategy adopted for the European DEMOnstration (EU-DEMO) fusion power reactor foresees the use of sacrificial components–referred to as limiters–dealing with plasma-wall contacts. Their aim is to protect the first wall (FW) against the huge amount of plasma energy (up to GigaJoules) released in a few milliseconds during disruptive events, which might lead to melting and/or vaporization of the foreseen plasma-facing tungsten armor. The current limiter design concepts rely on actively water-cooled plasma-facing components (PFCs) made of tungsten. As water is not allowed inside the main chamber, limiter’s PFCs must be designed to preserve the cooling system integrity under any scenario, therefore the estimate of the thickness of material undergoing any phase change is crucial. Given the initial assessment of plasma magnetic configurations during off-normal events, this article describes the procedure followed for designing the limiter’s PFC, which includes a novel approach for estimating the molten thickness of material under high heat flux. A simplified 1-D model has been implemented in MATLAB, which deals with multiphase moving boundary problems, hence its name Thermal Analysis foR Tracking InterFaces under meLting&vaporizaTion-induced plasma Transient Events (TARTIFL&TTE). This model takes its inspiration from the way phase change interface tracking problems are tackled in food industry (i.e., freeze-drying processes) and solute concentration diffusion-controlled problems. It overcomes the complexity of solving a strongly coupled non-linear system of partial and ordinary differential equations in moving spatial domains by adopting a change of coordinate system based on the Landau transformation. As a result, an equivalent and fixed spatial coordinate system is defined where the spatial domain boundaries of the different phases are constrained, and the systems of governing equations are easy to solve. Both the solid and liquid phases are modeled, while the vapor is assumed to be removed once formed. Its benchmark against computational results found in the literature has shown a very good agreement, which paves the way to further development of it. For phase change interface tracking problems in 2-D/3-D and more complex geometries, a commercial Multiphysics software adopting the Lagrangian–Eulerian approach in moving mesh frames will be used for tackling problems where material is removed following a phase change. Although the vapor domain is not simulated, a set of gas kinetics boundary conditions couples the interface between vapor and liquid phases, driving its position over time. This will be detailed described in a future companion article.

Interface-resolved simulation of solid-liquid phase change with thermal convection using discrete unified gas kinetic scheme

This work presents an interface-resolved model for simulating the solid-liquid phase transition problems. In the method, the natural convection flow is handled by the recently proposed discrete unified gas kinetic scheme (DUGKS). Compared with the conventional lattice Boltzmann equation, DUGKS has more advantages on numerical stability and multiscale computation. On the other hand, the ghost-cell (GC) immersed boundary scheme is applied to capturing the solid-liquid interface of pure material. In contrast to the phase-field or enthalpy-based method, the interface here with zero thickness is resolved and explicitly tracked on the fixed Cartesian grid. It is worth noting that due to the adsorption or release of latent heat, the phase front should be evolved according to the Stefan condition. Since compatible with the GC method, the position of the melting front is directly identified by the intersection point, and the complicated treatment of marker inserting and/or deleting in the previous studies is removed. Three numerical examples, including the one-dimension conduction melting, convection melting in a square cavity and melting in a circular cavity with natural convection are conducted to test the present model. The result demonstrates the accuracy of the method for solid-liquid phase change problems.