The Local Radial Basis Function Collocation Method for Evaluating the Non-Fourier Heat Transfer in Thermal Metamaterials

1  Introduction

Recent advancements in micro/nanotechnology and ultrafast thermal processing have exposed the fundamental limitations of traditional Fourier’s law, particularly regarding its diffusion-based assumptions. Under these extreme spatial and temporal conditions, heat transfer behaviors diverge significantly from classical predictions. To address this, non-Fourier heat conduction models—most notably the Cattaneo-Vernotte (C-V) model—incorporate a thermal relaxation time. This modification characterizes heat transfer as a wave-like phenomenon with a finite propagation speed, providing a theoretical framework for understanding the wave-like thermal behavior in one-dimensional periodic structures, such as thermal wave crystals [1,2]. By relaxing the assumption of infinite propagation speed inherent in Fourier’s law, the C-V model offers a more accurate depiction of transient heat transfer in micro/nano-scale systems [3–5]. The versatility of this model is well-documented; it has been employed to analyze thermoelastic fields under thermal shock, which offers vital reference data for aerospace thermal protection systems [6]. Furthermore, the model has been extended to complex coupling problems, such as radiation-conduction heat transfer within unified lattice Boltzmann frameworks [7]. Moreover, the C-V model contributes to the thermal analysis of porous media, functionally graded materials, and coating structures, aiding in the optimization of thermomechanical properties and energy conversion efficiency [8,9]. Consequently, the precise manipulation of non-Fourier heat transfer represents not only a significant engineering challenge but also a fundamental scientific inquiry.

Strategies for controlling non-Fourier heat conduction have advanced significantly. Research demonstrates that one-dimensional periodical structures can induce phonon band gaps, effectively suppressing phonon propagation at specific frequencies and drastically reducing thermal conductivity [10–12]. Furthermore, the periodic structure has been widely used to manipulate acoustic or elastic waves [13,14]. In previous studies, the extended plane wave expansion method has facilitated the exploration of thermal rectification in nonlinear lattices and precise heat control in thermal metamaterials [15,16]. Despite these advancements, critical limitations persist. As a simplified hyperbolic equation, the C-V model does not fully account for microscale non-equilibrium effects; in ultra-fast or nanoscale regimes, its accuracy may trail behind more complex frameworks like the dual-phase lag (DPL) model [17,18]. Experimental validation also remains challenging, as direct observation and measurement of micro/nano-scale thermal waves present significant challenges, often resulting in discrepancies between theoretical predictions and empirical data [3,4]. Furthermore, numerical simulation becomes computationally prohibitive when addressing complex heterogeneities, defects, or multi-physics couplings (e.g., thermo-mechanical and thermo-electrical interactions) [7,8]. Consequently, there is a critical need for high-efficiency and accurate numerical methods to simulate non-Fourier heat conduction within these periodic artificial materials.

Research on metamaterials often involves complex geometric structures, material interfaces, and the coupling effects between heat and other physical fields such as electromagnetic, mechanical, acoustical fields, and so on [19–22]. Thanks to these pioneering metamaterials, the concept of thermal metamaterial emerges and paves new avenues for heat transfer manipulations [23]. Firstly, the thermal metamaterial is used for enhancing heat transfer [24]. Furthermore, the thermal metamaterial is used to break the reciprocity of heat transfer [25]. Recently, the original parabolic diffusion equation can be modified into the hyperbolic convection diffusion equation by introducing convection. By this method, the heat transfer can be used to realize the non-Hermitian phenomena [26–28]. From the perspective of exploration, thermal metamaterials not only developed novel functions absent in traditional conduction but also discovered new physical phenomena in diffusive systems. However, the calculation of the hyperbolic convection diffusion equation with wave-like heat conduction behaviors remains challenging.

In recent years, the local radial basis function collocation method (LRBFCM) has provided a new and effective approach for solving such high-dimensional problems. As a strong-form meshless method, the core advantage of LRBFCM lies in avoiding the complex mesh generation process required by traditional finite element methods (FEM) or finite difference methods (FDM) [29,30]. The LRBFCM offers exceptional flexibility in applying boundary conditions, enabling convenient handling of various complex boundary types encountered in numerical simulation. The LRBFCM is adopted to analyze thermal stress and deformation of engineering structures under temperature loads. For fluid dynamics, it serves as an effective numerical tool for simulating internal and external flow characteristics. Moreover, it also performs well in the numerical analysis of electromagnetic fields, covering related design and performance evaluation of electromagnetic equipment [20,31–34]. Compared with the traditional transfer matrix method (TMM) and plane wave expansion method, the LRBFCM is much easier to deal with high-dimensional and complex geometry problems; however, the TMM or PWE are not very accurate in dealing with irregular geometry.

For thermal wave crystals, heat conduction is non-instantaneous, propagating at a finite speed with distinct wave-like characteristics. The meshless nature of LRBFCM is particularly advantageous here, effectively capturing these hyperbolic heat conduction traits even in the presence of irregular geometries or material inhomogeneities [35–38]. Unlike global Radial Basis Function (RBF) approaches, the traditional LRBFCM approximates field quantities by locally collocating RBFs within overlapping influence domains, thereby circumventing common ill-conditioning issues [39,40]. This localization results in a sparse matrix while preserving super-convergence, simultaneously boosting computational efficiency and algorithmic stability [41,42]. Consequently, due to its meshless architecture, high precision, stability, and adaptability to complex geometries, LRBFCM demonstrates significant superiority in solving hyperbolic models of one-dimensional heat transfer in thermal wave crystals.

This paper investigates the hyperbolic equation, non-Fourier heat conduction, within one-dimensional thermal metamaterials utilizing the LRBFCM. By integrating weighted techniques and novel derivative calculation schemes designed to enhance boundary stability, a robust and efficient numerical framework is constructed. This approach effectively mitigates the convergence challenges inherent in traditional methods, thereby providing a reliable computational tool for analyzing non-Fourier heat conduction in thermal metamaterials. Notably, the LRBFCM is capable of handling complex geometries and novel heat conduction phenomena, and can achieve computational efficiency comparable to that of analytical methods.

2  Algorithm of LRBFCM

This section provides a brief introduction to the relevant details of LRBFCM. A node with coordinate x = (x1, x2, x3) is selected as the center node within the computational domain, and Ns − 1 nearest neighboring nodes are identified around it to form a local subdomain containing Ns nodes. On each subdomain, the physical variable u at the ith node can be approximated as follows:

u(xi)=∑j=1Nsφ(‖xj−xi‖)αj,i=1,2,…,N(1)

where N denotes the total number of nodes within the computational domain, αj stands for the local unknown coefficient. In this paper, the multi-quadric (MQ) RBF is selected, where φ(‖xj−xi‖)=‖xj−xi‖2+c2, and c denotes the shape parameter. The shape parameter is a significant factor affecting the results of the LRBFCM, which is related to the node distance, interpolation accuracy, and other factors. There are various methods to search for the optimal shape parameter, such as the particle swarm optimization algorithm, cross-validation, and so on [43].

By performing interpolation calculations on nodes within local subdomains, the unknown coefficients can be obtained:

α=Φ−1u,(2)

where Φ=[φ(‖xj−xl‖)]1≤j,l≤Ns is a constant matrix if the positions of the nodes in the subdomain have already been determined, α=[α1,α2,…,αNs]T denotes the vector of unknown coefficients, and u=[u1,u2,…,uNs]T represents the vector containing physical variables at nodes within the local domain.

When solving partial differential equations, the mth-order derivative of the physical variable u in the n direction at any point in the computational domain can be represented by a linear combination of u at each node in the subdomain, that is:

∂mu(x)∂nm=∂mΘ(x)∂nmα=∂mΘ(x)∂nmΦ−1u=∂mw(x)∂nmu,(3)

where Θ=[φ(‖x1−x‖),φ(‖x2−x‖),…,φ(‖xNs−x‖)] and w=ΘΦ−1 is a weight vector.

When utilizing the LRBFCM to compute the variable derivatives at the boundaries, the algorithm often exhibits instability. Therefore, a direct method is employed to evaluate derivatives at the boundaries. For specific details, please refer to reference [44]. The direct method works by selecting local collocation points along the derivative direction, which effectively mitigates the numerical instability caused by first-order derivative evaluation in this approach.

3  The 1D Thermal Metamaterial

Consider a one-dimensional (1D) thermal metamaterial with unit cells as shown in Fig. 1a. The unit cell comprises a layer (sub-cell) A with thickness l1 and a layer (sub-cell) B with thickness l2, where the total thickness of the unit cell is l = l1 + l2. In 1D, the calculation domain of a unit cell with different boundaries and domains is shown in Fig. 1b.

Figure 1: Illustration of the 1D thermal metamaterial (a). The calculation domain of the unit cell (b). And the non-Fourier heat conduction in the finite thermal wave crystal with 10 unit cells (c).

3.1 Governing Equations

The ultra-fast laser heating process of nano-films is characterized by an ultra-short duration and ultra-small space size, in which the classical Fourier law based on the hypothesis of local equilibrium is no longer applicable [45,46]. When the thermal conduction process time scale is comparable to the relaxation time, the Cattaneo-Vernotte (CV) heat-conduction model [47,48], q+τqq˙=−κ∂T/∂x can be applied. By combining the energy conservation equation, the CV non-Fourier heat conduction equation can be written as

1τqiωT−ω2T=κρcpτq∂2T∂x2(4)

where q is the heat flux, i denotes the imaginary unit, ω is the angular frequency, T is temperature in the frequency domain, τq is the relaxation time for the phonon collision, κ is the thermal conductivity, ρ is the mass density, and cp is the specific heat. The heat propagates in the medium with a finite speed CCV=κ(ρcpτq)−1 [49]. Eq. (4) is the governing equation of the heat conduction in the frequency domain.

3.2 Interface Continuity Conditions

At the interface between two materials, the temperature continuity and heat conduction equilibrium conditions can be expressed as follows:

T1(x)=T2(x),x∈Γ0(5)

κ1∂T1(x)∂x−κ2∂T2(x)∂x=0,x∈Γ0(6)

where subscripts 1 and 2 represent two adjacent domains Ω1 and Ω2, respectively, as shown in Fig. 1b.

3.3 Periodic Boundary Conditions

According to Bloch’s theorem, a wave function in periodic structures can be expressed by a product of a plane wave and a periodic function. Specifically, the temperature and its surface traction components along the boundaries of the unit cell can be formulated by a lattice vector as follows:

T(x+a)=eiklT(x),x∈Γ1(7)

∂T(x+l)∂x−eikl∂T(x)∂x=0,x∈Γ1(8)

where k = kreal + ikimag is the Bloch wave number vector, kreal is the real part of the Bloch wave vector, and kimag is the imaginary part of the Bloch wave vector. By giving different frequencies, the wave vector can be obtained by solving a generalized eigenvalue problem, which is presented in the next section. For the complex wave vector, the real part kreal represents the propagation constant of the thermal wave and determines the phase change of the non-Fourier thermal wave. The imaginary part kimag represents the attenuation coefficient, determining the amplitude attenuation rate of the non-Fourier thermal wave. Therefore, the attenuation caused by the band gap can be described by the e−kimagl.

3.4 Response of Periodic Structure

In order to validate the results of the band structure, a 10-unit cell periodic structure is proposed as shown in Fig. 1c, where extra material with the length of l is used to alleviate the influence of the boundary conditions. The governing equation in Eq. (4), the interface conditions in Eqs. (7) and (8) are considered. While the boundary conditions on Γs and Γend are considered as

T(x)=1,x∈Γs(9)

T(x)=0,x∈Γe(10)

The Tout at an inner node close to the Γo is obtained to show the response of the finite periodic structure in Fig. 1c.

4  Numerical Discretization

In this section, the details of the numerical discretization of the LRBFCM are presented.

4.1 Governing Equations

Considering the discretization in Eqs. (3) to (4) yields,

1τqiωw(xm)Tm−ω2w(xm)Tm−κρcpτq∂2w(xm)Tm∂x2=0,x∈Ωm,m=1,2(11)

Tm is the temperature vector related to the computational domain Ωm,m=1,2, as shown in Fig. 1b.

4.2 Continuity Conditions

Substituting Eq. (3) into Eq. (6), the following expression is obtained

κ1∂w(x)∂xT1−κ2∂w(x)∂xT2=0,x∈Γ0(12)

where Γo is the interface in Fig. 1b. By assuming the gradient of temperature is continuous, the heat flux including relaxation time is also continuous.

4.3 Periodic Boundary Conditions

Taking Eqs. (3)–(8), the following can be obtained

∂w(x+l)∂xT1=eikl∂w(x)∂xT1,x∈Γg(13)

Using Eqs. (11)–(13) the following form is obtained

AT=eiklHT,(14)

where A and H are the RBF matrices given as

A=[1τqiωw(x1)−ω2w(x1)−κρcpτq∂2w(x1)∂x2001τqiωw(x2)−ω2w(x2)−κρcpτq∂2w(x2)∂x2κ1∂w(xΓ0)∂x−κ2∂w(xΓ0)∂xκ1∂w(xΓ2)∂x0]

H=[000000κ1∂w(xΓ1)∂x0]

T=[T1,T2]T

If the discrete temperature unknowns at boundary and interior nodes are ordered, then the following can be defined:

T=[T1,T2,TΓ1,TΓ2,TΓ01,TΓ02]T(15)

A=[A1,A2,AΓ1,AΓ2,AΓ01,AΓ02](16)

H=[H1,H2,HΓ1,HΓ2,HΓ01,HΓ02](17)

where A1,AΓi(i=0,1,2) and AΓ01 are respectively related to the nodes at xn located on Ω1, Γi(i=0,1,2), the left- and right-hand sides of Eq. (14) can be written explicitly as

AT=A1T1+A2T2+AΓ1TΓ1+AΓ2TΓ2+AΓ01TΓ01+AΓ02TΓ02,(18)

HT=H1T1+H2T2+HΓ1TΓ1+HΓ2TΓ2+HΓ01TΓ01+HΓ02TΓ02,(19)

Now the temperature periodicity condition (7) of the unit cell and the temperature continuity condition (5) on the interface have to be taken into consideration, which can be rewritten into the following forms

A~T~=eiklH~T~(20)

where

A~=[A1,A2,AΓ1,AΓ01+AΓ02]

H~=[H1,H2,HΓ2+HΓ1−AΓ1,HΓ01+HΓ02]

T~=[T1,T2,TΓ1,TΓ01]T

For the band structures, the eigenfrequencies are calculated by given specific wave vector combinations [50,51]. However, the wave vector is calculated by giving different frequencies by solving Eq. (20), which is different from previous work in [44]. Furthermore, to avoid ill-conditioning and improve high-frequency stabilities, a similar strategy is applied in our calculation process [52,53].

5  Numerical Results and Discussions

As a numerical example, a thermal metamaterial with l1 = l2 = 0.01 mm is considered. The layer A is stratum-like material, and the layer B is dermis-like material. For the material constants, Ref. [54] is cited; the values are given in Table 1. For the porous and biomaterials, the relaxation time varies from 1 to 103 s. The corresponding experimental data can also be found in this literature [55].

5.1 Convergence Rate Analysis

All the problems are tested using MATLAB R2026a on a Windows 10 operating system with an Intel(R) Xeon(R) Platinum 8272CL CPU running at 2.60 GHz and 192 GB of memory. To test the convergence rate, the stability is analyzed by considering different node distributions with a fixed shape parameter c2 = 0.01. The nodes 5, 10, 15, and 20 are used in the unit cell to obtain the numerical results in Fig. 2. To directly illustrate the convergence, the kreal is plotted with respect to the frequency, while kimag is plotted as the color map of the kreal. It can be seen from Fig. 2 that when only 5 nodes are used to calculate the complex dispersion curves, the eigenvalues in the low-frequency part seem acceptable. The colors are different from those of the 20 nodes. That is because the limited number of nodes cannot capture the high-frequency plots. As the number of nodes increases from 5 to 20, the situation changed at high-frequency parts. When the number of nodes is increased from 15 to 20, the band structures are very similar, which indicates that the results are convergent. More than 20 nodes are suggested in one unit cell to evaluate the band structures.

Figure 2: Results of the complex dispersion curves with 5, 10, 15, and 20 nodes.

The summation relative error is defined as

Errsum=∑|u1−u2|∑|u2|(21)

where the u2 is the eigenvalues obtained by using 100 nodes, and u1 is the eigenvalues obtained with 5, 10, 15, 20 nodes, and so on. While the relative error for every frequency can be derived by

Errevery=|u1−u2||u2|(22)

The Errsum and Errevery are calculated and shown in Fig. 3. Fig. 3a,b illustrates the Errsum with respect to the node number, where both error curves follow an oscillatory convergence pattern when the node number is larger than 20. This means the error magnitude declines rapidly with the increase in nodes, indicating that finer discretization consistently reduces the overall Errsum. Then Fig. 3c,d presents the Errevery of the real and the imaginary parts of the wave vector, respectively. For both real and imaginary wavenumber components, error magnitude uniformly decreases with increasing node count across all frequency modes: the 5-node configuration yields the largest deviations, while the over 15 nodes configuration delivers enough accuracy, and errors generally become more substantial at higher frequencies. In a word, Fig. 3 confirms that increasing the number of discretization nodes effectively improves numerical accuracy.

Figure 3: The Errsum of the real part (a) and imaginary part (b) of the wave vector. The Errevery of the real part (c) and imaginary part (d) of the wave vector with respect to frequency. The Errevery of different nodes 5 nodes (black), 20 nodes (red), and 60 nodes (blue).

At last, the mean absolute error,

MAE=1mtotal∑1mtotal|u1−u2|(23)

where mtotal is the total calculation time of the eigenfrequencies in Eq. (20). The calculation time is summarized in Table 2. From Table 2, it can be seen that a configuration of over 15 nodes delivers enough accuracy.

5.2 Numerical Validations

To validate the numerical results, the LRBFCM results are compared with those from the transfer matrix method [35], as shown in Fig. 4a,b. In these figures, circles denote the LRBFCM results, and the solid line denotes the transfer matrix results. The temperature response is also shown in Fig. 4, and the results are given in Fig. 4c. Totally, 100 nodes are distributed in the unit cell, and 1000 nodes are used in the whole computational domain, as shown in Fig. 4; the shape parameter c2 = 0.01 according to the test function and experience [43,44]. The LRBFCM with 100 nodes takes 20.622 s to complete the computation, whereas the transfer matrix method reported in Ref. [35] requires only 0.877 s. This notable gap in computational cost stems from their fundamental methodological differences. The transfer matrix method is an analytical approach with closed-form formulations, while the LRBFCM incurs additional overhead from local radial basis function approximation and collocation-based solving as a meshfree numerical scheme. Nevertheless, the LRBFCM achieves computational efficiency on par with the analytical method when the number of collocation nodes is reduced to 15. This result indicates that the LRBFCM can match the speed of analytical solutions while retaining the numerical flexibility to handle irregular geometries and complex boundary conditions.

Figure 4: Complex dispersion curves with (a) for the real part, (b) for the imaginary part, (c) the temperature responses. The black curves are from previous results [35], while the red dots are calculated by the LRBFCM.

The complex dispersion curves of the real and imaginary parts are illustrated in Fig. 4a,b, respectively. However, both results can obtain similar conclusions as previous in Ref. [35], where the pass bands are located at the frequency intervals of 0–1.8, 3.5–4.5, and 6.7–7.9 Hz. While the band gaps are in the frequency intervals of 1.8–3.5, 4.5–6.7 and 7.9–9.3 Hz with larger kimag, as shown in Fig. 4b. The temperature at the right boundary of the finite thermal wave crystal is shown in Fig. 4c. The temperature in Fig. 4c at the right end is quite small at the band gap areas, and very large at the pass band parts, which validates the results of the proposed LRBFCM.

Similarly, we can define the mid-gap frequency, fmid of a band gap from the Bragg condition,

fmid=m2(lAκA/ρAcpAτqA+lBκB/ρBcpBτqB)(24)

in which m = 1, 2, 3… is an integer representing the order of the band gap, and the subscripts A and B represent the material parameters from materials A and B, respectively. Based on Eq. (24), one can derive fmid = 2.8 Hz for m = 1, fmid = 5.6 Hz for m = 2, and so on. By comparing the fmid to the results in Fig. 4, one can find the fmid is very close to the center of the band gap, which testifies that the band gaps originate from the Bragg scattering.

At last, the temperature at Γo in Fig. 1c is calculated and shown in Fig. 4c as Tout. The Tout displays a clear alternating pattern of pass bands and band gaps, consistent with the dispersion characteristics of the thermal wave crystal. Three pass bands appear in the ranges of approximately 0–1.8, 3.5–4.5, and 6.7–7.9 Hz. Within these pass bands, Tout is larger than that of the band gaps, meaning the non-Fourier heat conduction can propagate through the thermal wave crystal. The significant drop of Tout demonstrates that non-Fourier heat conduction is strongly attenuated, with the strongest attenuation occurring near the center of each band gap.

Overall, the LRBFCM results exhibit excellent agreement with the reference previous results across the full frequency range in Fig. 4. The minor observable discrepancies originate from inherent differences between the two numerical methods, rather than any physical difference in the thermal wave crystal system. The previous results are calculated by an analytical solution. By contrast, LRBFCM is a meshfree collocation method that constructs the solution via local radial basis function approximation on scattered nodes. This fundamental difference in spatial discretization and numerical formulation naturally introduces small deviations in the computed eigenfrequencies [16].

To further analyze the non-Fourier heat conduction process in the thermal metamaterial, the temperature distribution of the periodic structure in Fig. 1c is presented for several frequencies, as shown in Fig. 5. The amplitudes of T(x) are obtained by considering the boundary conditions in Eqs. (9) and (10), the T(x) is changing with the position of x. In Fig. 5a, the heat can propagate through the periodic structure because the excitation frequency, f = 1.44 Hz, is in the pass band shown in Fig. 4. However, the heat fast decays when propagating through the thermal wave crystal, as shown in Fig. 5b–d. This is because these excitation frequencies, f = 2.59, 4.90, and 8.29 Hz, are in the band gaps in Fig. 4. In a word, good agreement is shown between the present results and those of previous work [35].

Figure 5: Temperature distribution in 1D thermal metamaterial with different excitation frequencies, f = 1.44 Hz (a), f = 2.59 Hz (b), f = 4.90 Hz (c) and 8.29 Hz (d).

Furthermore, one can find that the non-Fourier heat conduction decays faster at 4.90 Hz than at 2.59 and 8.29 Hz. That is because the imaginary part of the wave vector in the band gap, 4.5–6.7 Hz, is larger than that of the 1.8–3.5 and 7.9–9.3 Hz. However, the non-Fourier heat transfer still meets the laws of thermodynamics. Therefore, the thermal waves cannot be localized to certain positions, i.e., kreal = 0, in the thermal wave crystal. That is the main difference from that of the photonic and phononic crystals.

Later, the influence of different numbers of unit cells on the non-Fourier heat conduction is studied. Then, the thermal metamaterial is composed of different numbers of unit cells and an extra material A with the length of l. Similarly, we define the Tout as the temperature at Γo. The governing equations of Eq. (4), interface continuity conditions of Eqs. (5) and (6), and temperature conditions of Eqs. (9) and (10) are considered. The Tout of different unit cells is evaluated as shown in Fig. 6. The axis-x in Fig. 6 is the range of the frequency. From Fig. 6, it is not easy to find that, as the number of unit cells increases, the Tout decreases when the frequency is located at the band gap. For those results in the pass band, the Tout decreases as frequency increases. At last, it can be seen from Fig. 6 that there appear to be strong boundary reflections around 5 Hz, which is different from previous studies. With fewer periodicities, there appears to be a resonance peak at 5.6 Hz.

Figure 6: The Tout of different periodicities.

To analyze the boundary reflection of the non-Fourier heat conduction. The temperature distribution of different periodicities is calculated and shown in Fig. 7. Fig. 7 compares the temperature distribution in the thermal metamaterial at 4.9 Hz (blue) against the higher-frequency 5.6 Hz case (red) for 3 to 20 periodicities. It can be seen that the temperature at boundary Γo is largely affected by the periodicity. This frequency is located in the band gap, and the attenuation coefficient of thermal waves, kimag, is large. However, reflection occurs when thermal waves propagate to the end of the structure. Then the incident wave and reflected wave interfere with each other, forming standing non-Fourier thermal waves. Thus, a resonance peak appears, leading to an increase in temperature amplitude, when the structure length satisfies the resonance condition. But, the attenuation in the band gap is enhanced as the number of periods increases, the influence of boundary reflection gradually weakens, and the resonance peak disappears. Beyond numerical validation, the rapid decay of amplitude near the boundary across all stabilized configurations serves as tangible evidence of the thermal attenuation capabilities of the band gap. And it confirms the successful blocking of heat propagation.

Figure 7: Temperature distribution of different periodicities around 4.9 Hz (blue) and 5.6 Hz case (red).

6  Conclusions

In this paper, the LRBFCM is applied for the first time to evaluate the band gaps by using eigenvalues rather than the wave vector within the LRBFCM framework. The detailed formulation of the LRBFCM is presented, and the complex dispersion curves are calculated and validated against results from the transfer matrix method. The transient temperature response in finite systems is further simulated using the LRBFCM. The results reveal that Bragg scattering gives rise to distinct band gaps that significantly suppress heat conduction in non-Fourier thermal transport processes. From a numerical perspective, increasing the node density markedly enhances stability in the high-frequency regime. The LRBFCM is shown to be a meshless and geometrically flexible approach. The convergence analysis indicates that MAE decreases steadily as the node density increases, demonstrating reliable accuracy for the problems investigated. These features render it suitable for studying complex geometries and unconventional heat conduction phenomena.

Looking ahead, we plan to extend the LRBFCM to higher-dimensional thermal wave crystals and thermal metamaterials, as well as to structures with complex geometries that are beyond the capabilities of transfer matrix and plane wave expansion methods. Additionally, the LRBFCM holds potential for designing novel thermal devices such as thermal concentrators, cloaks, Luneburg lenses, and thermal mirages, which could find applications in laser hardening, laser cladding, and metal additive manufacturing.

Acknowledgement: None.

Funding Statement: This work is supported by the National Natural Science Foundation of China (No. 12402099), the Interdisciplinary Research Project for Young Teachers of USTB (Fundamental Research Funds for the Central Universities) (No. FRF-IDRY-GD21-01), the State Key Laboratory of Nuclear Power Safety Technology and Equipment Open Fund (No. SKL-2024-TS-04).

Author Contributions: Study conception and design: Anyu Hong, Zheng-Yang Li; Data collection: Anyu Hong; Analysis and interpretation of results: Anyu Hong, Zheng-Yang Li; Draft manuscript preparation: Anyu Hong, Zheng-Yang Li. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Nie BD, Cao BY. Reflection and refraction of a thermal wave at an ideal interface. Int J Heat Mass Transf. 2018;116(1):314–28. doi:10.1016/j.ijheatmasstransfer.2017.09.043. [Google Scholar] [CrossRef]

2. Manzanares-Martinez J, Esquivel-Sirvent R. Non-Fourier omnidirectional thermal mirror in finite biological multilayer. Results Phys. 2026;80(22):108563. doi:10.1016/j.rinp.2025.108563. [Google Scholar] [CrossRef]

3. Chen G. Non-Fourier phonon heat conduction at the microscale and nanoscale. Nat Rev Phys. 2021;3(8):555–69. doi:10.1038/s42254-021-00334-1. [Google Scholar] [CrossRef]

4. Wang FF, Wang B. Current research progress in non-classical Fourier heat conduction. Appl Mech Mater. 2013;442:187–96. doi:10.4028/www.scientific.net/amm.442.187. [Google Scholar] [CrossRef]

5. Vedavarz A, Kumar S, Moallemi MK. Significance of non-Fourier heat waves in conduction. J Heat Transf. 1994;116(1):221–4. doi:10.1115/1.2910859. [Google Scholar] [CrossRef]

6. Guo SL, Zhang YX, Wang KF, Wang BL, Zhang CW. Effects of non-Fourier heat conduction and surface heating rate on thermoelastic waves in semi-infinite ceramics subject to thermal shock. Ceram Int. 2021;47(12):17494–501. doi:10.1016/j.ceramint.2021.03.067. [Google Scholar] [CrossRef]

7. Yan W, Wei Y, Liu X, Zhu K, Huang Y. Unified lattice Boltzmann framework for coupled non-Fourier conduction and thermal radiation based on the C-V model. Appl Therm Eng. 2024;252(2):123609. doi:10.1016/j.applthermaleng.2024.123609. [Google Scholar] [CrossRef]

8. Cui YJ, Li WJ, Wang KF, Wang BL, Guo SL. Thermal shock fracture of honeycomb-based porous thermoelectric materials with non-Fourier heat conduction. Ceram Int. 2024;50(1):2151–61. doi:10.1016/j.ceramint.2023.10.328. [Google Scholar] [CrossRef]

9. Yang Y, Dai HL, Ye C, Xu WL, Luo AH. Investigation of the one-dimensional transient heat conduction problem of a coated high strength steel plate. Math Mech Solids. 2019;24(11):3472–84. doi:10.1177/1081286519847709. [Google Scholar] [CrossRef]

10. Zen N, Puurtinen TA, Isotalo TJ, Chaudhuri S, Maasilta IJ. Engineering thermal conductance using a two-dimensional phononic crystal. Nat Commun. 2014;5(1):3435. doi:10.1038/ncomms4435. [Google Scholar] [PubMed] [CrossRef]

11. Anufriev R, Maire J, Nomura M. Review of coherent phonon and heat transport control in one-dimensional phononic crystals at nanoscale. APL Mater. 2021;9(7):070701. doi:10.1063/5.0052230. [Google Scholar] [CrossRef]

12. Maire J, Anufriev R, Nomura M. Thermal conduction control by phononic crystal nanostructures. JSAP Annu Meet Ext Abstr. 2015:2804. [Google Scholar]

13. Wang YF, Wang YZ, Wu B, Chen W, Wang YS. Tunable and active phononic crystals and metamaterials. Appl Mech Rev. 2020;72(4):040801. doi:10.1115/1.4046222. [Google Scholar] [CrossRef]

14. Fu C, Tang XL, Liu YD, Ma TX, Wang YS. Simultaneous manipulation of elastic and acoustic waves in acousto-elastic metamaterial beams. Extreme Mech Lett. 2025;75(12):102286. doi:10.1016/j.eml.2024.102286. [Google Scholar] [CrossRef]

15. Dhar A. Heat transport in low-dimensional systems. Adv Phys. 2008;57(5):457–537. doi:10.1080/00018730802538522. [Google Scholar] [CrossRef]

16. Li ZY, Mellmann M, Wang Y, Ma TX, Yan D, Golub MV, et al. Non-Fourier heat conduction in 2D thermal metamaterials. Mater Today Commun. 2024;38:107828. doi:10.1016/j.mtcomm.2023.107828. [Google Scholar] [CrossRef]

17. Yang W, Chen Z. Nonlocal dual-phase-lag heat conduction and the associated nonlocal thermal-viscoelastic analysis. Int J Heat Mass Transf. 2020;156(8):119752. doi:10.1016/j.ijheatmasstransfer.2020.119752. [Google Scholar] [CrossRef]

18. Fan QM, Lu WQ. A new numerical method to simulate the non-Fourier heat conduction in a single-phase medium. Int J Heat Mass Transf. 2002;45(13):2815–21. doi:10.1016/S0017-9310(01)00364-7. [Google Scholar] [CrossRef]

19. Dong HW, Shen C, Zhao SD, Qiu W, Zheng H, Zhang C, et al. Achromatic metasurfaces by dispersion customization for ultra-broadband acoustic beam engineering. Natl Sci Rev. 2022;9(12):nwac030. doi:10.1093/nsr/nwac030. [Google Scholar] [PubMed] [CrossRef]

20. Dong HW, Zhao SD, Wang YS, Zhang C. Topology optimization of anisotropic broadband double-negative elastic metamaterials. J Mech Phys Solids. 2017;105:54–80. doi:10.1016/j.jmps.2017.04.009. [Google Scholar] [CrossRef]

21. Wen Z, Gao P, Mao J, Dai S, Marti-Sabaté M, Jin Y, et al. Compact metaplate with bound state in the continuum: from quasisymmetry to symmetry. Phys Rev Lett. 2025;135(8):087201. doi:10.1103/sbk1-szl7. [Google Scholar] [PubMed] [CrossRef]

22. Fu M, Miao H, Kang G. Nondestructive subwavelength metablocks for manipulation of shear horizontal waves. Mech Syst Signal Process. 2026;247:113936. doi:10.1016/j.ymssp.2026.113936. [Google Scholar] [CrossRef]

23. Ju R, Xu G, Xu L, Qi M, Wang D, Cao PC, et al. Convective thermal metamaterials: exploring high-efficiency, directional, and wave-like heat transfer. Adv Mater. 2023;35(23):2209123. doi:10.1002/adma.202209123. [Google Scholar] [PubMed] [CrossRef]

24. Xu LJ, Huang JP. Transformation thermotics and extended theories. Berlin/Heidelberg, Germany: Springer; 2023. [Google Scholar]

25. Li Y, Li J, Qi M, Qiu CW, Chen H. Diffusive nonreciprocity and thermal diode. Phys Rev B. 2021;103(1):014307. doi:10.1103/physrevb.103.014307. [Google Scholar] [CrossRef]

26. Xu G, Li W, Zhou X, Li H, Li Y, Fan S, et al. Observation of Weyl exceptional rings in thermal diffusion. Proc Natl Acad Sci U S A. 2022;119(15):e2110018119. doi:10.1073/pnas.2110018119. [Google Scholar] [PubMed] [CrossRef]

27. Cao PC, Ju R, Wang D, Qi M, Liu YK, Peng YG, et al. Observation of parity-time symmetry in diffusive systems. Sci Adv. 2024;10(16):eadn1746. doi:10.1126/sciadv.adn1746. [Google Scholar] [PubMed] [CrossRef]

28. Li Y, Peng YG, Han L, Miri MA, Li W, Xiao M, et al. Anti–parity-time symmetry in diffusive systems. Science. 2019;364(6436):170–3. doi:10.1126/science.aaw6259. [Google Scholar] [CrossRef]

29. Fu Z, Tang Z, Xi Q, Liu Q, Gu Y, Wang F. Localized collocation schemes and their applications. Acta Mech Sin. 2022;38(7):422167. doi:10.1007/s10409-022-22167-x. [Google Scholar] [CrossRef]

30. Siraj-ul-Islam, Vertnik RSarler B,. Local radial basis function collocationmethod along with explicit time stepping for hyperbolic partial differential equations. Appl Numer Math. 2013;67(3):136–51. doi:10.1016/j.apnum.2011.08.009. [Google Scholar] [CrossRef]

31. Mavrič B, Šarler B. Application of the RBF collocation method to transient coupled thermoelasticity. Int J Numer Methods Heat Fluid Flow. 2017;27(5):1064–77. doi:10.1108/hff-03-2016-0110. [Google Scholar] [CrossRef]

32. Mramor K, Vertnik R, Šarler B. Simulation of laminar backward facing step flow under magnetic field with explicit local radial basis function collocation method. Eng Anal Bound Elem. 2014;49:37–47. doi:10.1016/j.enganabound.2014.04.013. [Google Scholar] [CrossRef]

33. Hon YC, Šarler B, Yun DF. Local radial basis function collocation method for solving thermo-driven fluid-flow problems with free surface. Eng Anal Bound Elem. 2015;57(1):2–8. doi:10.1016/j.enganabound.2014.11.006. [Google Scholar] [CrossRef]

34. Dobravec T, Mavrič B, Šarler B. Improved finite difference method for phase-field modelling of dendritic solidification. J Comput Phys. 2026;553(6):114716. doi:10.1016/j.jcp.2026.114716. [Google Scholar] [CrossRef]

35. Chen AL, Li ZY, Ma TX, Li XS, Wang YS. Heat reduction by thermal wave crystals. Int J Heat Mass Transf. 2018;121(4):215–22. doi:10.1016/j.ijheatmasstransfer.2017.12.136. [Google Scholar] [CrossRef]

36. Idesman A, Dey B. A high-order numerical approach with Cartesian meshes for modeling of wave propagation and heat transfer on irregular domains with inhomogeneous materials. Comput Meth Appl Mech Eng. 2020;370(2):113249. doi:10.1016/j.cma.2020.113249. [Google Scholar] [CrossRef]

37. Yan XB, Zheng H, Zhang C, Wen PH, Sladek J, Sladek V. A time-domain local radial basis function collocation method for the band structure analysis of 2D anti-plane phononic crystals. Eng Anal Bound Elem. 2024;162(30–31):203–19. doi:10.1016/j.enganabound.2024.01.034. [Google Scholar] [CrossRef]

38. Yun DF, Hon YC. Improved localized radial basis function collocation method for multi-dimensional convection-dominated problems. Eng Anal Bound Elem. 2016;67(2):63–80. doi:10.1016/j.enganabound.2016.03.003. [Google Scholar] [CrossRef]

39. Vertnik R, Šarler B. Meshless local radial basis function collocation method for convective-diffusive solid-liquid phase change problems. Int J Numer Meth Heat Fluid Flow. 2006;16(5):617–40. doi:10.1108/09615530610669148. [Google Scholar] [CrossRef]

40. Jiang P, Zheng H, Xiong J, Rabczuk T. The localized radial basis function collocation method for dendritic solidification, solid phase sintering and wetting phenomenon based on phase field. J Comput Phys. 2025;520(8–9):113515. doi:10.1016/j.jcp.2024.113515. [Google Scholar] [CrossRef]

41. Zheng H, Zhang C, Yang Z. A local radial basis function collocation method for band structure computation of 3D phononic crystals. Appl Math Model. 2020;77(13):1954–64. doi:10.1016/j.apm.2019.09.006. [Google Scholar] [CrossRef]

42. Jiang P, Zheng H, Xiong J, Zhang C. A stabilized local RBF collocation method for incompressible Navier–Stokes equations. Comput Fluids. 2023;265(3):105988. doi:10.1016/j.compfluid.2023.105988. [Google Scholar] [CrossRef]

43. Zheng H, Yao G, Kuo LH, Li X. On the selection of a good shape parameter of the localized method of approximated particular solutions. Adv Appl Math Mech. 2025;10(4):896–911. doi:10.4208/aamm.oa-2017-0167. [Google Scholar] [CrossRef]

44. Zheng H, Zhang C, Wang Y, Sladek J, Sladek V. A meshfree local RBF collocation method for anti-plane transverse elastic wave propagation analysis in 2D phononic crystals. J Comput Phys. 2016;305:997–1014. doi:10.1016/j.jcp.2015.10.020. [Google Scholar] [CrossRef]

45. Mao Y, Liu S, Liu J, Yu M, Li X, Yang K. Non-Fourier heat conduction of nano-films under ultra-fast laser. Materials. 2023;16(14):4988. doi:10.3390/ma16144988. [Google Scholar] [PubMed] [CrossRef]

46. Camacho de la Rosa A, Becerril D, Gómez-Farfán MG, Esquivel-Sirvent R. Bragg mirrors for thermal waves. Energies. 2021;14(22):7452. doi:10.3390/en14227452. [Google Scholar] [CrossRef]

47. Cattaneo C. Sur une forme de l’équation de la Chaleur éliminant le paradoxe d’une propagation instantanée. Comptes Rendu. 1958;247:431–3. [Google Scholar]

48. Vernotte P. Les paradoxes de la théorie continue de l’équation de la Chaleur. Comptes Rendu. 1958;246:3154–5. doi:10.1007/2-287-26785-9_4. [Google Scholar] [CrossRef]

49. Tzou DY. Macro- to microscale heat transfer: the lagging behavior. New York, NY, USA: John Wiley & Sons; 2014. doi:10.1002/9781118818275. [Google Scholar] [CrossRef]

50. Xu C, Wei P, Li Z, Guo X. Shear horizontal wave propagation on a piezoelectric substrate with periodic gradient local resonant metasurfaces. Appl Math Model. 2025;137:115733. doi:10.1016/j.apm.2024.115733. [Google Scholar] [CrossRef]

51. Yan X, Zheng H, Yan D. Analysis of the band structure of transient in-plane elastic waves based on the localized radial basis function collocation method. Appl Math Model. 2024;125:468–84. doi:10.1016/j.apm.2023.09.002. [Google Scholar] [CrossRef]

52. Jančič M, Kosec G. Stability analysis of RBF-FD and WLS based local strong form meshless methods on scattered nodes. In: Proceedings of the 2022 45th Jubilee International Convention on Information, Communication and Electronic Technology (MIPRO); 2022 May 23–27; Opatija, Croatia. p. 275–80. [Google Scholar]

53. Mishra PK, Fasshauer GE, Sen MK, Ling L. A stabilized radial basis-finite difference (RBF-FD) method with hybrid kernels. Comput Math Appl. 2019;77(9):2354–68. doi:10.1016/j.camwa.2018.12.027. [Google Scholar] [CrossRef]

54. Xu F, Seffen KA, Lu TJ. Non-Fourier analysis of skin biothermomechanics. Int J Heat Mass Transf. 2008;51(9–10):2237–59. doi:10.1016/j.ijheatmasstransfer.2007.10.024. [Google Scholar] [CrossRef]

55. Massaguer A, Teixidor M, Leroy M, Goeminne J, Suñol JJ, Massaguer E. Method of obtaining dual-phase-lag times in aggregate porous materials with a nonhomogeneous inner structure. J Therm Anal Calorim. 2025;150(19):15251–60. doi:10.1007/s10973-025-14605-x. [Google Scholar] [CrossRef]