Thermal Simulation of PCM (Phase Change Materials) for Building Applications
Theory and Physics
What is PCM?
Professor, I heard PCM is used in construction. Does the wall melt? It sounds a bit scary...
The wall itself doesn't melt. PCM (Phase Change Material) is a material that absorbs a large amount of heat when it changes from solid to liquid. In construction, a type of paraffin wax is microencapsulated into particles 5–50 μm in diameter and mixed into plasterboard or mortar.
Microcapsules? So the paraffin is enclosed in a shell, and it melts and solidifies only inside that?
Exactly. For example, BASF's Micronal DS series encapsulates paraffin with a melting point around 23°C in a melamine resin shell. During the day, when the room temperature exceeds 23°C, the PCM in the wall begins to melt, storing solar heat as latent heat. At night, when the temperature drops, it solidifies and releases heat. This "storage → release" cycle can suppress room temperature fluctuations by 2–4°C.
Huh, so the wall becomes a "thermal buffer." But regular concrete also stores heat, right? What's the advantage of PCM?
Good question. Concrete stores heat as sensible heat—by increasing temperature. This results in low storage density. Water's sensible heat is 4.18 kJ/kg per 1°C rise, but the latent heat of fusion for paraffin C18 is about 244 kJ/kg. Comparing the same mass, PCM can store heat comparable to what concrete stores after a 50°C temperature rise, within a temperature change of just 1–2°C. Moreover, it operates near the "comfort zone" of room temperature, directly reducing heating and cooling loads.
The total enthalpy change during PCM's heat storage/release is expressed as the sensible heat of the solid phase + latent heat + sensible heat of the liquid phase.
Here, $m$ is mass, $c_s, c_l$ are the specific heats of the solid and liquid phases, $T_m$ is the melting point, $L_f$ is the latent heat of fusion, and $T_1, T_2$ are the initial and final temperatures. Representative property values for paraffin-based PCMs used in construction are shown below.
| Property | Value | Unit |
|---|---|---|
| Melting point $T_m$ | 21–26 | °C |
| Latent heat of fusion $L_f$ | 170–250 | kJ/kg |
| Density $\rho$ (solid) | 850–950 | kg/m³ |
| Thermal conductivity $k$ | 0.15–0.25 | W/(m·K) |
| Specific heat $c_p$ (solid/liquid) | 1.8–2.4 | kJ/(kg·K) |
Stefan Problem and Moving Boundary
When PCM melts or solidifies, the boundary between the solid and liquid phases moves, right? How is that formulated?
That's the Stefan problem. It's a moving boundary problem involving phase change, formulated in 1891 by the Austrian physicist Jozef Stefan. It solves the heat conduction equation in the solid and liquid regions separately, satisfying the energy conservation condition (Stefan condition) at the interface $x = s(t)$. In one dimension, it looks like this.
Solid region ($0 < x < s(t)$):
Liquid region ($s(t) < x < L$):
Stefan condition (energy balance at interface $x = s(t)$):
The interface position $s(t)$ is an unknown and moves over time. Can this be solved analytically?
Under special conditions like one-sided heating of a semi-infinite solid, there is an analytical solution like the Neumann solution $s(t) = 2\lambda\sqrt{\alpha t}$. But for problems with finite thickness, temperature variations on both sides, and multiple material layers like building envelopes, an analytical solution is impossible. Numerical methods are essential, and that's where the enthalpy method comes in.
Governing Equations of the Enthalpy Method
The enthalpy method is a convenient way to solve the Stefan problem? What's the idea behind it?
The tricky part of the Stefan problem is "tracking the interface position." The enthalpy method is a brilliant idea to avoid this by using volumetric enthalpy $H$ as the dependent variable instead of temperature $T$. The $H(T)$ function jumps by the latent heat amount across the phase change, allowing the solid, liquid, and mixed phases (mushy zone) to be treated uniformly with a single equation. There's no need to explicitly track the interface.
The governing equation for the enthalpy method can be written as follows.
Here, the volumetric enthalpy $H(T)$ is defined as a function of temperature.
For building PCMs, phase change does not occur at a single temperature but over a temperature range from $T_s$ to $T_l$ (typically 2–4°C wide). This range is the mushy zone, where solid and liquid phases coexist.
I see! Writing the equation in terms of temperature becomes discontinuous at the interface, which is problematic, but writing it in terms of enthalpy allows for smooth handling.
Exactly. As the mushy zone width $\Delta T = T_l - T_s$ approaches zero, it reduces to the ideal "sharp interface" Stefan problem. Actual building PCMs have a finite width of $\Delta T = 2\sim4$°C, so the enthalpy method is actually a more natural formulation.
Heat Transfer Mechanism of PCM Walls
Actual building walls aren't just a PCM layer; they have insulation, concrete, and various other layers, right? How do you model the whole thing?
Let me explain with an example of a typical PCM wall composition. From the outside, it's a multilayer structure like "exterior finish (10mm) → insulation (50mm) → PCM-containing plasterboard (12.5mm) → air gap → interior finish." We solve the coupled one-dimensional unsteady heat conduction for each layer.
For each layer $i$ of a multilayer wall:
Connection conditions between layers (temperature continuity and heat flux conservation):
Boundary conditions on the exterior and interior surfaces include convection ($h_\mathrm{ext}$, $h_\mathrm{int}$), solar absorption ($\alpha_s I_\mathrm{sol}$), and longwave radiation ($\varepsilon \sigma (T^4 - T_\mathrm{sky}^4)$).
So only the PCM layer has temperature-dependent $c_p(T)$, making it nonlinear. The other layers can remain as ordinary heat conduction equations.
Right. So the computational focus boils down to "how finely to discretize the PCM layer." Insulation and concrete are fine with coarse meshing, but the PCM layer requires fine discretization because enthalpy changes rapidly within the phase change temperature range $\Delta T$.
History of the Stefan Problem
The name "Stefan problem" comes from the Austrian physicist Jozef Stefan (1835-1893), who was also Boltzmann's teacher. He tackled the problem of "how thick Arctic sea ice grows." Given the surface temperature of the ice and the seawater temperature, the Neumann solution stating that ice thickness $s(t)$ grows proportionally to $\sqrt{t}$ is still used today as a benchmark problem for PCM simulation verification. Application to the building field became mainstream in the 1980s, spurred by research by Doris Boehm and others at the U.S. Department of Energy (DOE).
Physical Meaning of Each Term
- Enthalpy accumulation term $\partial H / \partial t$: Rate of enthalpy change per unit volume. Includes both sensible heat (temperature change) and latent heat (phase change). Near the PCM melting point, $H$ increases sharply for a small temperature change—this is the "heat absorption" effect. In everyday terms, ice water stays at 0°C for a long time because the ice continuously absorbs the latent heat of fusion.
- Heat conduction term $\nabla \cdot (k \nabla T)$: Heat diffusion based on Fourier's law. PCM's thermal conductivity is low at 0.15–0.25 W/(m·K), about 1/6 that of concrete (1.4 W/(m·K)). Therefore, PCM layers have the characteristic of being "good at storing heat but slow at transferring it." In wall design, it's important to balance utilizing PCM's storage capacity while ensuring a heat release path to the interior side.
- Heat source term $Q$: Internal heat generation (generally $Q=0$ for PCM walls). However, for active PCM panels with built-in electric heaters, Joule heating $Q = I^2R/V$ may be considered.
Assumptions and Applicability Limits
- Assumption of 1D heat transfer: When the thickness is sufficiently small compared to the wall area, heat conduction in the plane direction can be ignored, and it can be treated as a 1D problem in the thickness direction. Around window frames and corners, heat bridges occur, requiring 2D/3D analysis.
- Neglect of natural convection: In microencapsulated PCMs, the liquid phase is confined within tiny capsules, so macroscopic natural convection does not occur. However, for bulk PCMs (e.g., tank storage), convection in the liquid phase becomes dominant, requiring coupling with the Navier-Stokes equations.
- Neglect of volume change: Paraffin undergoes about 10% volume change during phase change, but this is absorbed by the elastic deformation of the microcapsules, so its structural impact is usually negligible.
- Neglect of hysteresis: Actual PCMs exhibit different temperature-enthalpy curves (hysteresis) for melting and solidification, but this is often ignored in many BES tools. Precise evaluation requires a hysteresis model.
Dimensional Analysis and Dimensionless Numbers
| Dimensionless Number | Definition | Physical Meaning | Typical Value for PCM Walls |
|---|---|---|---|
| Stefan number $\mathrm{Ste}$ | $c_p \Delta T / L_f$ | Ratio of sensible to latent heat | 0.02–0.1 (latent heat dominant) |
| Biot number $\mathrm{Bi}$ | $h L / k$ | Ratio of surface convection to internal conduction | 0.5–5 |
| Fourier number $\mathrm{Fo}$ | $\alpha t / L^2$ | Dimensionless time (degree of diffusion progress) | $\gg 1$ for annual cycles |
Numerical Methods and Implementation
Discretization of the Enthalpy Method
To actually solve the enthalpy method on a computer, how do you discretize it?
In building BES (Building Energy Simulation) tools, discretizing the wall in the thickness direction using the Finite Difference Method (FDM) is mainstream. Taking EnergyPlus's CondFD algorithm as an example, the PCM layer is divided into $N$ nodes, and the following implicit difference equation is solved for each node $j$.
Since both $H$ and $T$ are unknowns, they are connected by the $H(T)$ relation. But this is nonlinear, right?
Sharp observation. At each time step, $T^{n+1}$ must be found from $H^{n+1}$ using the inverse function $T = H^{-1}(H)$. Specifically, an iterative calculation is performed, starting with the temperature from the previous time step as an initial estimate and updating the temperature on the $H(T)$ curve. Within the mushy zone, $T = T_s + (H - H_s)\Delta T / (\rho L_f)$ can be solved explicitly, so iteration often converges in 1–3 steps.
Comparison with the Apparent Heat Capacity Method
I've also heard of the "Apparent Heat Capacity Method." How is it different from the enthalpy method?
The Apparent Heat Capacity Method expresses latent heat as a temperature-dependent equivalent specific heat $c_\mathrm{eff}(T)$.
Using the apparent heat capacity method, it can be solved as the usual heat conduction equation $\rho c_\mathrm{eff}(T) \partial T / \partial t = \nabla \cdot (k \nabla T)$, making it easy to implement in existing solvers. However, there's a big pitfall—if the time step is too coarse, it can "skip over" the phase change temperature range, causing the latent heat to be completely missed in the calculation.
Wow, that's scary! So if you jump from $T_s$ to $T_l$ in one step, the latent heat energy disappears?
Exactly. That's the biggest weakness of the apparent heat capacity method. On the other hand, the enthalpy method directly tracks $H$, so energy conservation of latent heat is guaranteed even with somewhat coarse time steps. In practice, they are used as follows.
| Method | Advantages | Disadvantages | Recommended Use |
|---|---|---|---|
| Enthalpy Method | Strict energy conservation, robust | Requires $H \to T$ inverse transformation | BES tools (EnergyPlus, etc.) |
| Apparent Heat Capacity Method | Easy to implement in existing solvers | Time step dependent, risk of missing latent heat | CFD/FEM (COMSOL, Fluent, etc.) |
Time Step and Mesh Discretization
What are appropriate time steps and mesh fineness for analyzing PCM walls?
This is a fairly critical point. PCM's thermal diffusivity is very small: $\alpha = k/(\rho c_p) \approx 8 \times 10^{-8}$ m²/s. The characteristic time for a 12.5mm thick PCM plasterboard is $\tau = L^2/\alpha \approx 2400$ seconds (40 minutes). Using this as a reference:
- Time step: 1–3 minutes (3 minutes is practical for annual simulations, 1 minute for precise evaluation)
- Spatial discretization: Minimum 4 divisions for PCM layer ($\Delta x \le 3$ mm), recommended 6–8 divisions
- Fourier number constraint: $\mathrm{Fo} = \alpha \Delta t / \Delta x^2 \le 0.5$ (for explicit methods. No constraint for implicit methods)
For an annual simulation, that's 365 days $\times$ 1440 minutes = 525,600 steps... Even with 3-minute steps, that's 175,200 steps. What about computational cost?
Since it's a 1D wall finite difference, the computational cost per step is tiny. For an annual simulation of a single-zone model in EnergyPlus, even with PCM walls, it finishes in tens of seconds to a few minutes. For a whole multi-zone building, about 10 minutes. Detailed 2D/3D analysis in COMSOL becomes orders of magnitude heavier, but 1D-BES is sufficient unless you need to see details like natural convection inside the wall or near the interface.
Handling Nonlinear Convergence
Does the nonlinearity of the PCM layer cause convergence problems?
Related Topics
なった
詳しく
報告