Comparison of Different Numerical Models for Solidification of Paraffin-Based PCMs and Evaluation of the Effect of Variable Mushy Zone Constant

1. Introduction

Among renewable energy technologies, converting solar energy into thermal energy is currently one of the most widely adopted methods [1,2,3,4]. To maximize the efficient use of renewable energy in daily life and industrial applications, thermal energy storage technology is often necessary to stabilize the energy output and ensure its effective utilization [5]. Phase change materials (PCMs) possess a robust capacity for energy storage and exhibit a noteworthy ability to maintain a consistent temperature while absorbing or releasing heat.

Numerous studies have investigated the solidification process of PCMs within a range of encapsulation geometries, including spheres, rectangular enclosures, and horizontal cylinders [6,7,8,9,10,11,12]. The melting and solidification behavior of PCM with oscillations initiated at different times was investigated by Liu et al. [13]. They examined the effects of oscillation frequency and amplitude on melting performance and found that oscillation accelerates the phase transition process by enhancing convective heat transfer. Sparrow and Broadbent studied the solidification of a PCM within a vertical cylinder [14]. They experimentally investigated the freezing process of initially superheated PCMs. The research concludes that the impact of initially superheating the liquid on energy extraction and frozen mass is moderately effective, particularly for short freezing durations. An experimental setup was designed to monitor temperatures within a cylindrical vessel containing PCM by Shmueli et al. [15]. The obtained experimental data were then compared with results from a computational model utilizing the enthalpy‒porosity method. The influence of container geometry on energy storage was investigated experimentally by Menon et al. [16]. Additionally, they proposed a numerical model designed to forecast heat fluxes throughout the processes of solidification and melting of the PCM within a vertical cylinder. Mahamudur et al. [17] carried out an experimental study to explore the impact of the cylinder size, thermal loading, and inclusion of a graphene additive on the melting and solidification dynamics of PCMs within a millimeter-scale vertical cylinder. Kapilow et al. [18] conducted an experimental study on the melting and solidification of PCMs within a small-diameter vertical cylinder operating under convective boundary conditions. They employed air as the heat transfer fluid and studied the impacts of different air velocities and temperatures in their investigation. the impact of natural convection on heat transfer during the solidification and melting of PCMs within a vertical cylinder was investigated by Kim et al. [19]. They presented that natural convection in the liquid phase of a PCM is significant primarily during the melting process, whereas its impact is negligible during the solidification phase. A three-dimensional model for evaluating various aspects of the discharge of phase change material (PCM) enclosed within a vertical cylindrical enclosure was developed by Izgi and Arslan [20]. The computational model considers the influence of natural convection. Their findings showed that natural convection played a substantial role in the initial stages of the freezing process. Conduction subsequently became the dominant mechanism throughout the entire process. Furthermore, they reported that the diameter of the cavity significantly affects the freezing process, whereas changes in the cavity height do not have a noticeable effect.

Computational fluid dynamics (CFD) is an influential tool particularly adept at addressing complex, nonlinear problems involving the concurrent presence of mass, heat, and momentum transfer [21]. Consequently, it proves invaluable for simulating the thermal behavior of PCMs. The leading commercial CFD software packages include Ansys Fluent and COMSOL Multiphysics. Ye et al. [22] applied Ansys Fluent software to research the fluid flow and heat transfer within a plate-fin unit. Ansys Fluent’s solidification and melting model is based on the enthalpy‒porosity method, which is a widely employed tool for investigating phase change phenomena. Nowadays, substantial research efforts have been dedicated to utilizing this model, yielding valuable agreements with experimental findings [23,24,25,26].

Numerical phase change models in the literature are generally categorized into two main groups: fixed grid (or continuum) models and variable grid (or two-phase) models [27]. Fixed grid models apply a single computational domain, where continuity, energy conservation, and momentum equations are used uniformly across the entire domain. In contrast, variable grid models use separate computational domains for the solid and liquid phases, coupled through transfer terms between the phases [28]. As a result, fixed grid models are simpler to implement than variable grid models, while still providing accurate results [29]. Nazzi Ehms et al. [30] provide a clear classification of fixed grid models for the solidification and melting of PCM, organized by latent heat modeling and velocity transition modeling. Latent heat models, applied within the energy equation, can be based on a source term, enthalpy, or temperature. Temperature-based models are further categorized into three types: apparent heat capacity, effective heat capacity, and heat integration. Similarly, velocity transition models are classified into three main groups: the switch-off method, the source term method, and the variable viscosity method. Moreover Li et al. [31] present axisymmetric, enthalpy-based Lattice Boltzmann models for simulating solid–liquid phase change. Their findings indicate that these models effectively simulate axisymmetric solid–liquid phase transitions.

Numerical simulation of a PCM requires inputting material properties and parameters as input values. The mushy zone is where the melting or solidification process starts and progresses, during which latent heat is absorbed or released. Therefore, accurate representation of the solid‒liquid phase change and the associated latent heat storage performance is crucial. This parameter is part of a source term in the Carman–Kozeny equation, which is used in the enthalpy‒porosity technique for the momentum equation to simulate the flow resistance through the mushy zone. This approach addresses the main drawback of the fixed grid method, which struggles with handling the zero-velocity condition as the liquid transitions to a solid within the mushy zone. The main questions are how the mushy zone constant $A_{mushy}$ affects the predicted solidification performance and how to determine $A_{mushy}$. While the Carman–Kozeny equation is commonly used in most models, the effects of the $A_{mushy}$ parameter continue to be a topic of frequent investigation [32]. Values for $A_{mushy}$ ranging from ${10}^3$ to ${10}^{10}$ have been recommended in commercial software guidelines and by various researchers [33]. The melting of lauric acid was numerically studied by Abdi et al. [34] and they proposed an $A_{mushy}$ value of ${10}^{6}kg/\left(m^3.s\right)$ because it was closer alignment with experimental data than did the values of ${10}^5$ and ${10}^{7}kg/\left(m^3.s\right)$. Similar research has been conducted by Hong et al. [35]. Fadl et al. [36] reported that the optimal $A_{mushy}$ values are ${5\times 10}^5$ and ${2\times 10}^5$for the vertical and horizontal orientations, respectively. The most common value for $A_{mushy}$ is the recommended ${10}^{5}kg/\left(m^3.s\right)$ [37,38,39,40]. Kheirabadi et al. [33] suggested that for lauric acid, the optimal value for $A_{mushy}$ is ${10}^6$ when ANSYS FLUENT and COMSOL Multiphysics solvers are used.$A_{mushy}$ significantly impacts the accuracy and correctness of the predicted melting performance. However, its role in predicting solidification performance remains poorly reported to date. It is quite interesting to examine how this flow resistance parameter affects solidification and to understand the extent of its influence on the related heat recovery performance.

This research aims to investigate the impact of the numerical model and the mushy zone parameter on phase change solidification simulations, using paraffin as a representative phase change material. The cooling process was simulated using three methods implemented in the CFD code ANSYS Fluent 2020 R2 [41]: the enthalpy–porosity method, the AHC method, and a new model proposed by the authors, in which heat capacity is directly introduced into ANSYS Fluent. Furthermore, the impact of the〖 A〗_mush on the accuracy and correctness of the predicted solidification performance has been investigated by testing different values of 〖 A〗_mush (constant and variable). Then the results such as temperature profiles, solidification time, and solidification profiles analyze scientifically. This paper offers a detailed discussion of the observed differences and how they change with varying model parameters, providing explanations in physical terms. This research aims to provide a more accurate approach to study the effects of simulation methods and the mushy zone constant on the solidification process in a vertical cylinder under real boundary conditions.

2. Methodology

2.1. Physical Model

The physical model in this study is based on the authors' experimental work, which examined the solidification of paraffin-based PCMs under T-history conditions. [42,43]. The schematic diagram of the 2-D simulation geometry is presented in Figure 1. A rectangular computational domain was defined, taking into account two distinct fluid regions: one containing air and the other comprising the working substance, either water or PCM. The tube has a diameter of 30 mm and a height of 170 mm, with 140 mm filled with liquid PCM. In the upper part, a rubber plug was modeled via a thin-wall thickness approach, incorporating combined heat transfer around it. Moreover, the model was numerically simulated, assuming axisymmetric conditions, incompressible laminar flow, and the Boussinesq approximation for buoyancy effects. For solidification, the external wall was designed to account for heat losses via conduction, convection, and radiation to the surrounding environment, as depicted in the equivalent thermal circuit in Fig. 1. To calculate surface emissivity and radiation heat transfer, the authors [42] proposed a numerical method that reduces computational costs and simulation time, improving the efficiency of thermal analysis while maintaining accuracy.

2.2. Governing Equations

2.2.1. Enthalpy‒Porosity Method

The enthalpy‒porosity technique has been extensively applied in simulating the melting and solidification phenomena within PCM enclosures, particularly when natural convection plays a significant role. When dealing with fusion and solidification within a fixed domain, adjustments to the momentum equation are necessary to accommodate phase transitions, ensuring zero velocities in the solid phase. A common approach to address this involves incorporating a high-magnitude source term ($S_u$) into the momentum equation.

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho u_i\right)+{\displaystyle \frac{\partial }{\partial X_j}}\left(\rho u_i{u}_j\right)=-{\displaystyle \frac{\partial p}{\partial X_i}}+\mu {\displaystyle \frac{{\partial }^2{u}_i}{\partial X_j\partial X_j}}+\rho g_i+S_u$$

(1)

where ρ is the density and where μ is the dynamic viscosity. The term S_u modifies the moment balance depending on the phase. The term varies, ranging from a high value that enforces complete immobility in the solid region to eventually reaching a limit of zero as the material transitions into a fully liquid state. The source term of Eq. (1) is defined as

$$S_u=-A\left(\gamma \right)u_i$$

(2)

Here, A (γ) is the porosity function, which allows the moment equation to mimic the behavior of the flow in the mushy region of the PCM, defined by the Carman–Kozeny equation derived from Darcy's law for a porous medium [44]:

$$A\left(\gamma \right)={\displaystyle \frac{C{(1-\gamma )}^2}{{\gamma }^3+\epsilon }}$$

(3)

where $\gamma$ is the liquid volume fraction and where $C$ is the constant of the mushy zone, which varies according to its morphology. Higher values of $C$ result in stronger velocity damping, leading to a thinner mushy region and approaching a situation resembling that of a pure substance. The constant $\epsilon$ is introduced to prevent division by zero.

Voller [44] initially introduced the enthalpy‒porosity technique to address phase change problems related to convective diffusion-controlled heat transfer. The energy equation can be written as:

$${\displaystyle \frac{\partial \rho H}{\partial t}}+\overrightarrow{\nabla }.\left(\rho \overrightarrow{\vartheta }H\right)=\overrightarrow{\nabla }.\left(\lambda \overrightarrow{\nabla }T\right)$$

(4)

The governing energy equation thereby contains one enthalpy H that includes latent and specific heat

$$H\left(T\right)={\int }_{T_m}^T{C}_{p}dT+\rho \gamma L$$

(5)

where $\rho$ corresponds to the density, $C_p$ corresponds to the specific heat and $L$ corresponds to the latent heat of the PCM. $L$ is associated with the liquid fraction, $\gamma$, which allows the computation of the change in enthalpy from the energy in the material during the phase change [43]. The general form of $\gamma$ can be written as

$$\gamma =\left\{\begin{array}{c}0ifT<T_s\\ {\displaystyle \frac{T-T_s}{T_l-T_s}}if{T}_s<T<T_l\\ 1ifT>T_l\end{array}\right.$$

(6)

In this formulation, the porosity is equated to the liquid fraction, $f_l$, within a given element, as dictated by the enthalpy balance. The liquid fraction is recalculated for each iteration. The "mushy zone" refers to the region where the liquid fraction ranges from 0 to 1. Consequently, when a material is completely solidified within an element, the porosity decreases to zero. During the phase change from liquid to solid, the velocity decreases to zero, making it crucial to accurately incorporate this phenomenon in the numerical formulation. Researchers have employed various methods [47] to implement this aspect of physics.

2.2.2. Apparent Heat Capacity (AHC)

The Apparent heat Capacity Method (AHC) takes into account the influence of enthalpy and its temporal evolution by incorporating an apparent heat capacity during thermal phase transitions. The approach is based on the following relationship:

$${\displaystyle \frac{\partial H}{\partial t}}={\displaystyle \frac{\partial H}{\partial T}}{\displaystyle \frac{\partial T}{\partial t}}$$

(7)

where ${\displaystyle \frac{\partial H}{\partial T}}=\rho c_{app}\left(T\right)$ represents the temperature-dependent apparent heat capacity. The apparent specific heat capacity is calculated as the derivative of the specific enthalpy, taking into account a temperature-dependent phase change enthalpy [45].

$$c_{app}\left(T\right)=\left(1-\gamma \right)c_p^s+\gamma c_p^l+{\displaystyle \frac{d\gamma }{dT}}(h^l\left(T\right)-h^s\left(T\right)$$

(8)

where, $h^s\left(T\right)$ and $h^l\left(T\right)$ represent the extrapolated single-phase specific enthalpies for the solid and liquid states, respectively [46].

In the enthalpy porosity model, $\gamma \left(T\right)$ represents a linearly increasing liquid fraction between the solidus and liquidus temperatures, assuming constant values (zero and one) outside the phase change temperature range. Buruzs et al. [47] indicate that the non-differentiable nature of this piecewise linear liquid fraction function can lead to convergence problems, so they introduced an alternative relation for $\gamma \left(T\right)$ has been introduced based on the smooth step function and cubic Hermite spline functions [46].

2.2.2. New Model

In the proposed model, the heat capacity of the PCM is determined through Differential Scanning Calorimetry (DSC) and implemented in ANSYS FLUENT software using a User-Defined Function (UDF) file. It is important to note that DSC results are highly sensitive to varying heating and cooling rates; therefore, heat capacity is calculated at heating rates of 1, 5, 10, and 20 ◦C/min. The methodology is simple: the obtained specific heat capacity curve as function of temperature ($c_p-T$) from DSC is used as a variable, specific heat capacity of the material in ANSYS FLUENT software. One of the primary heat transfer phenomena in phase change is natural convection, which arises from temperature-induced density differences. Without accounting for variable density, the model reduces to a conduction problem, addressing only pure diffusion and neglecting the movement within the liquid phase. Since convection is incorporated in the enthalpy-porosity model, and it is likely to be important in our problem, we use the following generalization [48,49]. The technique of a source term $S_u$ in the momentum equation will relate linearly the enthalpy with the liquid fraction, and thus lead to numerical results similar to those of the enthalpy-porosity method. Instead, we model a liquid with a viscosity that varies with temperature according to our measurements on the liquid PCM above the solidus point, and below we fair a temperature function reaching a significantly high value. This assures that velocity of the medium decreases in the mushy region and finally stops when it is a solid.

$${\mu }_l=\left\{\begin{array}{c}-\left(2.273\times {10}^{-1}\right)T+33.83T\le T_s\\ \\ \left(2.283\times {10}^{-7}\right)T^3-\left(3.501\times {10}^{-5}\right)T^2+\left(1.73\times {10}^{-3}\right)T-0.02587T>T_s\end{array}\right.$$

(9)

2.3. Numerical Model Setup

The numerical model setup parameters were chosen on the basis of ANSYS FLUENT guidelines [41] and previous research papers [27,30,40,42,50] that have numerically studied the phase change process.

The ANSYS Fluent 2022 R2 code was employed to simulate the solidification process. The conservation equations for both momentum and mass were addressed via the semi-implicit method for pressure-linked equations (SIMPLE) algorithm. The momentum and energy equations were discretized via a second-order upwind scheme, and pressure correction was carried out via the PRESTO! (PRESure STaggering Option) method. Relaxation factors of 0.3, 0.7, and 1 were applied to pressure, momentum, and energy, respectively. At each time step, the convergence of the solution was verified throughout the simulation. These modeling choices align with the recommendations in [40 and 42].

2.4. Mesh Tests

To confirm that the results are independent of the chosen mesh (grid) resolution, various element sizes were tested and compared with experimental results over the entire process duration. Structured meshes with 65844, 23800, 6120, and 2332 quadrilateral elements were employed. Parameters such as the maximum aspect ratio, orthogonality, cell quality, and obliquity were meticulously adjusted to meet the recommended criteria. In Figure 2, a close alignment is observed between the experimental findings and the numerical model when water is used as the working fluid. Evidently, there are no substantial disparities in the simulated cooling curves across various mesh sizes compared with the experimental results.

3. Numerical Results

In this section, numerical results are presented to examine the influence of the thermal properties of the PCM and the mushy zone parameters on the solidification process. The numerical modeling was conducted under the assumptions that (1) the flow is laminar and incompressible, (2) the material properties in each simulation remain constant, as indicated in Table 1, (3) the Boussinesq approximation was applied, and (4) the domain is modeled as a fluid in which the solid phase has a high viscosity.

Figure 3 illustrates the PCM cooling process employing the enthalpy‒porosity method with a mushy zone constant (C) value of 105. The PCM properties were inputted on the basis of the manufacturer's provided values, as detailed in Table 1. The Boussinesq approximation was used for buoyancy due to variable density. The graph reveals two noteworthy observations. First, the numerical cooling curves exhibited a consistent pattern across the various mesh sizes, with no notable discrepancies. Second, there is a noticeable abrupt shift in the cooling curve at the start and end points of solidification (at 46 and 40 °C, respectively) in the numerical simulation. In contrast, under real experimental conditions, a much more gradual change in slope is observed, especially at the end.

This study is divided into two main sections. The first part focuses on investigating the influence of the numerical model on phase change solidification simulation using paraffin as a sample phase change material to determine which model is more accurate. The second part examines the effect of the mushy zone parameter on solidification, comparing constant and variable mushy zone parameters to identify the optimal values for use in numerical modeling of the phase change process for paraffin material.

3.1. Comparing the Three Numerical Methods

As outlined in section 2.2, three different numerical models were used to simulate the solidification of PCM-paraffin under T-history conditions. Figure 4 shows the temperature during the solidification process for all models. It is evident that the new model aligns most closely with the experimental data compared to the others, while the S&M and AHC models exhibit a slight discrepancy with the experimental results. Specifically, for the S&M model, there is a noticeable abrupt shift in the cooling curve at the start and end points of solidification (at 46 and 40 °C, respectively). In contrast, under real experimental conditions, a much more gradual change in slope is observed, especially at the end. Although the AHC model has a good agreement with experimental result, the new model demonstrates a more precise change in the solidification slope, which aligns more closely with the experimental results. To compare the results obtained with respect to the experimental data, the mean absolute percentage error (MAPE), which is also known as the mean absolute percentage deviation (MAPD), was used. The variations in the cooling curves of New Model and AHC model compared with the experimental curve showed percentage errors (EPAM) of 3.1% and 8.3%, respectively.

Figure 5 presents the comparison of liquid fraction for the three models. Although the liquid fraction for all models are similar at the beginning of the solidification process, a noticeable discrepancy between them emerges towards the end of solidification. Moreover, the solidification profile of the new model differs significantly from the others. Initially, the solidification rate of the new model is faster than the other models, but it slows down as the process progresses, ultimately resulting in a bit shorter total solidification time compared to the other two models. Moreover, the new model shows a change in the gradient of the liquid fraction, with a distinct maximum point for the solidification rate. In other words, the new model shows that the solidification rate initially rises sharply to a peak, then decreases through the middle of the process. Following this, a gradual increase occurs until a local maximum is reached, after which the rate declines to zero. The same pattern has been observed AHC but no for S&M model (Fig.5b). This phenomenon has been noted by Milad et al. [43], who also showed that the solidification rate reaches a maximum point during the process. Furthermore, an examination of the slope of the liquid fraction curve revealed that the rate of solid formation initially increased, peaked, and then steadily decreased to zero as the solidification process progressed. Reference [51] explained that this phenomenon occurs because, initially, numerous nuclei form, and the crystals grow larger due to liquid subcooling. It causes an enhancment in the solid formation rate. However, as solidification progresses, continuous nucleation events are effectively inhibited. This suppression occurs because the thermal gradient driving the phase change decreases as the nucleation site moves farther from the cooling wall, and the thermal resistance rises. As a result, the entire solidification process gradually slows.

The key characteristic of the T-history experiment is that the Biot number satisfies the lumped system analysis condition (Bi < 0.1). This implies that temperature varies with time but not with position, meaning the temperature of the medium changes uniformly over time.

Figure 6 illustrates the temperature difference between the centerline and the wall for the three points (Y = 3.5, 7, and 10.5 cm). The new model shows the smallest temperature difference, fluctuating between 1.25 and 0.5 K and the same pattern can be seen for the AHC model where temperature difference is put in a narrow range. In contrast, the S&M model causes more uneven temperature changes, with discrepancies ranging from 2.0 to 0.4 K. Another interesting observation is that, although a slight variation is seen in the S&M model when moving from the bottom to the top of the tube, both the AHC model and the new model show an almost unchanged temperature pattern along the tube’s length. In fact, with these two models, the temperature difference between the centerline and the wall remains consistent, regardless of the measurement point’s location. Moreover, during the solidification process, the cooling rate of the PCM decreases, likely due to the release of latent heat, increased thermal resistance from the growing solid phase, or a combination of both factors, as the temperature of the sample and ambient environment approach each other. This implies that the Biot number decreases during solidification, leading to a reduction in the temperature difference between the centerline and the wall of the tube. This phenomenon is clearly illustrated by the new model.

Figure 7 compares the solidification profiles for three models over a 150-200 minute period. The primary observation is that the solidification profiles differ based on the simulation model. As anticipated, the solidification rate in the S&M model is lower than in the other models during this timeframe. After 150 minutes, solidification begins to progress at the bottom of the tube in the S&M model, a pattern not observed in the other two models. Although solidification develops at the tube’s bottom in both the S&M and ACH models, in the new model, it slows in the bottom layers and instead concentrates near the cold wall. This discrepancy in the solidification profile behavior is directly related to the heat transfer characteristics in the mushy region.

3.2. The Effect of the Mushy Zone

In order to examine the impact of the mushy zone parameter, several simulations for different mushy zone parameters of ${10}^5$-${10}^8$ were conducted. Liquid fraction throughout the solidification process were investigated for each set of simulations, as they are directly influenced by the mushy zone parameter.

Numerous studies suggest a reliable method, grounded in the physical characteristics of the mushy zone, to determine the mushy zone constant, $A_{mush}$. For various PCMs, the $A_{mush}$ value can be determined as a function of dynamic viscosity and the characteristic length of solid microstructures within the mushy zone. In fact, $A_{mush}$ is parameter that related the between the microstructural characteristics of the mushy zone and the macroscopic heat transfer during phase change. Yang et al. [52] evaluated the $A_{mush}$ parameter on the basis of measured microstructural features and proposed the $A_{mush}$ value function of liquid fraction for paraffin.

For solidification:

$$A_{mush}={\displaystyle \frac{4\mu }{125{\gamma }^{0.3}}}\times {10}^{12}$$

(10)

For melting:

$$A_{mush}={\displaystyle \frac{4\mu }{125\gamma }}\times {10}^{12}$$

(11)

Figure 8 shows the liquid fraction of the PCM for different values of $A_{mush}$ and compare it related to simulation methods. The liquid fraction follows a similar trend across varying $A_{mush}$ values, with only slight differences toward the end of the solidification process. The solidification rate for $A_{mush}={10}^8$ is faster than that for $A_{mush}={10}^5$ in both the S&M and AHC methods, predicting a solidification time that aligns more closely with the new model. Furthermore, when $A_{mush}$ is variable, the liquid fraction profile remains the same as when $A_{mush}$ is constant at ${10}^8$. Notably, the effects of different $A_{mush}$ values on the melting rate by applied S&M method have been examined [40], and they reported that a significant variation in the mushy zone constant, from ${10}^5$ to${10}^8$, led to a greater temperature difference.

Figure 9. The effect of $A_{mush}$ on liquid fraction.

Figure 9. The effect of $A_{mush}$ on liquid fraction.

4. Conclusions

The main goal of the current study was to examine the numerical model of simulation and the mushy zone constant on the simulation of the solidification process inside a vertical cylinder for thermal energy storage applications. Three simulation models—enthalpy-porosity, AHC, and the new model—were applied to simulate solidification, and their results were compared. In the new numerical model, the specific heat capacity curve as a function of temperature, obtained from DSC, is directly introduced into ANSYS FLUENT. Additionally, a variable viscosity method is used to account for the increased viscosity associated with the solid phase. Moreover, the effect of the mushy zone parameter on solidification was examined across the simulation models, and the results were compared.

This study revealed that the proposed model shows excellent agreement with experimental data for cooling temperature and provides a more accurate prediction than the other models. The liquid fraction pattern for the new model and AHC is almost identical, whereas the S&M model displays a distinct pattern, with the solidification rate decreasing continuously until it reaches zero. Furthermore, the new model better satisfies the lumped system condition, showing less sensitivity to position in temperature variation compared to the other models. $A_{mushy}$ is a critical parameter for accurately modeling phase change phenomena, as it gradually influences the prediction of solidification processes Although the impact of $A_{mushy}$ is more pronounced in the melting process, where heat transfer is primarily dominated by natural convection, its influence is generally less prominent in the solification process, where conductive heat transfer is dominant. Increasing the mushy parameter value results in a reduction of solidification time and similar results observed for both a variable $A_{mushy}$ and $A_{mushy}={10}^8$. A mushy zone constant of ${10}^8$ has been proposed for the solidification process, which is in good agreement with the real results. The criteria established from this study can guide the accurate and effective modeling of paraffin’s PCM behaviors across various applications.