Thermal Simulation of Thermal Storage Materials (PCM)

Good question. PCM (Phase Change Material) is a material that absorbs a large amount of heat when it changes from solid to liquid—that is, when it melts. The essence of thermal storage is using this "latent heat" to control temperature.

For example, when ice melts at 0°C, it absorbs 334 kJ/kg of heat without its temperature rising, right? PCM utilizes this principle. Paraffin-based PCMs have a latent heat of about 200 kJ/kg near a melting point of 28°C, and salt hydrates have about 260 kJ/kg near 32°C. Designs that fill the spaces between cells with PCM to keep battery pack temperatures below 40°C have been rapidly increasing in recent years.

Predicting how the solid-liquid boundary—the melting front—moves, when and how far it melts, and how the temperature rises after it's completely melted. This is a classic moving boundary problem called the "Stefan problem." In CAE, the mainstream approach is to solve it using the enthalpy method without explicitly tracking the boundary.

That's right. Let's consider a one-dimensional Stefan problem. Suppose there is a PCM slab heated to a high temperature on one side. If the position of the solid-liquid interface is $s(t)$, then at the interface, the condition "the interface advances while absorbing latent heat" holds:

Exactly. That's why in practice, the enthalpy method, which avoids interface tracking, is used most of the time. The idea is "we don't need to know where the interface is; as long as the energy balance is correct, the temperature field will come out right."

In the enthalpy method, we solve the energy conservation equation with enthalpy per unit volume $H(T)$ as the variable:

Correct! For pure substances, it jumps vertically at the melting point. For materials with a melting temperature range like paraffins or salt hydrates, it becomes a sloped step. The height of this jump corresponds to the latent heat $L$. In simulation, this sharp change causes nonlinearity, so care is needed with time steps and mesh.

Another well-known method is the Apparent Heat Capacity Method. It models latent heat as an apparent specific heat $c_{p,eff}$ over the melting temperature range:

Practically, the enthalpy method is overwhelmingly advantageous. With the apparent heat capacity method, if the time step is coarse, it can skip over the melting temperature range entirely and miss the latent heat—this is the so-called "latent heat skip". The phase change modules in Ansys Fluent and COMSOL internally use the enthalpy method; the apparent heat capacity method is a simplified technique when reusing a standard FEM thermal analysis solver.

The Stefan number is a dimensionless number representing the ratio of sensible heat to latent heat. It serves as an indicator to quantify the "difficulty" of PCM simulation:

Roughly speaking:

• $\mathrm{Ste} \ll 1$ (Latent heat dominant): Melting progresses slowly. PCM cooling for battery packs is in this region. If $\Delta T = 12$°C, $c_p = 2000$ J/(kg·K), $L = 200{,}000$ J/kg, then $\mathrm{Ste} = 0.12$. Time steps can be relatively large.
• $\mathrm{Ste} \gg 1$ (Sensible heat dominant): Almost like a normal heat conduction problem. Solidification in casting (hundreds of °C $\Delta T$) is an example. Time steps need to be fine to avoid numerical oscillation.
• $\mathrm{Ste} \sim 0.1$: The region where Neumann's approximate solution can be used. The melting front position can be estimated as $s(t) \approx 2\lambda\sqrt{\alpha t}$, which is convenient for verifying analysis results.

They are broadly classified into three types. Let's organize them along with the property values needed for simulation:

Yes, this is the biggest practical challenge in PCM simulation. Low thermal conductivity makes the melting rate extremely slow. So in actual products, aluminum fins are added or carbon nanotubes are dispersed to increase the effective thermal conductivity. In simulation, you either treat it as a composite of "PCM + fins" or use homogenized effective property values—this is where the designer's skill comes in.

With FVM (Finite Volume Method), you can directly treat the enthalpy per cell as a conserved quantity, so it's a good match. Ansys Fluent and STAR-CCM+ use this approach. With FEM (Finite Element Method), it's discretized after converting to the weak form.

For PCM, implicit methods (backward Euler) are the iron rule. There are two reasons:

• Paraffin has low thermal conductivity, so to satisfy the stability condition $\mathrm{Fo} = \alpha \Delta t / \Delta x^2 \leq 0.5$, a very small time step is required.
• The sharp change in enthalpy (latent heat jump) creates strong nonlinearity, but with implicit methods, you can iterate to convergence at each time step.

It's the intermediate state—neither solid nor liquid—where solid and liquid coexist. For paraffin mixtures, the melting point spreads over a range like 28~32°C, so that entire temperature band becomes the mushy zone.

Actually, selecting $C_{mush}$ is the "artisan" part of PCM simulation. For pure paraffin, $10^5$ to $10^6$ is common, adjusted to match experimental melting front positions. If $C_{mush}$ is too large, the mushy zone becomes too stiff and melting slows down; if it's too small, unphysical flow occurs in the solid region. Papers always perform a "$C_{mush}$ sensitivity analysis," so I recommend trying at least three levels.

It moves a lot. This is the tricky part of PCM simulation. Melted paraffin develops density differences due to temperature differences, causing buoyancy-driven natural convection. When the Rayleigh number exceeds $10^5$, convective effects become dominant, and melting can progress 2~5 times faster than predicted by pure conduction.