Energy transport and relaxation of phonon polaritons (PhPs) are studied, based on a macroscopic phonon model, for atomic layers of hexagonal boron nitride (hBN) and transition metal dichalcogenides (TMDs). The velocity of the energy flow (energy velocity) is derived from the energy flow and density; it equals the group velocity, similar to the results of bulk and surface PhPs. In electrostatic approximation, valid once the frequency is slightly above ω0 (e.g., ω>1.002ω0∼224cm−1 for pentalayer MoTe2; ω0 is the zone-center optical-phonon frequency), simple formulas are obtained for the energy velocity and relaxation rate (ERR). While the energy velocity increases proportionally with the number of layers N, the ERR is independent of N. The ERR equals the phonon damping rate in freestanding layers, but it is slightly decreased in SiO2-supported layers and has a non-monotonic frequency dependence (the decrease is smaller than 5.4% for hBN layers and negligible for TMD layers). The energy velocity decreases significantly with frequency in both freestanding and supported layers. Near ω0, however, the PhP properties should be calculated rigorously, and they all depend on N as well as the dielectric environment. High-frequency screening should be included to study the energy transport. The energy velocity can be engineered by varying N and the dielectric environment; it also can be tuned together with the propagation quality factor by the incident light frequency. The MoTe2 layers should be exploited for a far-infrared PhP material (wavelengths 43–45 μm); this is just an example of application of the proposed model to the considered materials.