Numerical Simulation of Time Domain Thermoacoustic Wave Equation Using the Finite Difference Time Domain Method with Comparison to the K-Space Pseudospectral Approach

1. Introduction

TA wave propagation has gained significant attention due to its wide range of applications in areas such as biomedical imaging, non-destructive testing, and acoustic sensing [1,2,3]. In these applications, localized energy deposition leads to the generation of pressure waves, which propagate through the medium and carry valuable information about the underlying physical properties [4,5,6]. Accurate numerical modeling of such wave phenomena is essential for understanding propagation characteristics and improving system performance [7,8,9,10,11,12].

Numerical methods such as the FDTD method have been extensively used for solving acoustic wave equations due to their simplicity and ease of implementation [11,13,14]. However, FDTD schemes are known to suffer from numerical dispersion, particularly at higher frequencies and coarse grid resolutions [12,15]. To address these limitations, spectral and pseudospectral methods, including k-Space pseudospectral approach, have been developed [9,10,16,17]. These methods compute spatial derivatives in the Fourier domain and significantly reduce dispersion errors while allowing larger time steps [7,8]. Previous studies have demonstrated the effectiveness of k-Space methods in simulating acoustic and TA wave propagation with improved accuracy [8,18].

Despite the advancements in numerical modeling techniques, most existing studies focus on single type of source geometry, typically Gaussian or simple symmetric distributions. Limited attention has been given to the systematic investigation of different source shapes and their influence on wave propagation characteristics. Furthermore a comprehensive comparison between FDTD and k-Space pseudospectral methods under varying CFL conditions and grid resolutions remains insufficiently explored. The impact of asymmetric and complex source geometries on both time domain and frequency domain responses also requires further investigation.

The present study addresses these gaps by performing a systematic analysis of TA wave propagation using multiple source geometries, including Gaussian, circular disk, Chebyshev-based, and asymmetric (rock-like) distributions. A direct comparision between FDTD and k-Space pseudospectral methods is conducted under consistent simulation conditions. Additionally, the combined effects of CFL number and spatial grid resolution on wave propagation, stability, and spectral behavior are investigated. This comprehensive approach provides new insights into the influence of source geometry and numerical parameters on acoustic wave simulations.

The main objectives of this study are as follows:

• To model thermoacoustic wave propagation using both FDTD and k-space pseudospectral methods.
• To investigate the effect of different source geometries on wave propagation behavior.
• To analyze the influence of the CFL number and grid resolution on stability and accuracy.
• To compare time-domain and frequency-domain responses obtained from both methods.
• To evaluate the performance of the numerical schemes in terms of accuracy and consistency.

2. Theory

2.1. Generation of Thermoacoustic Wave Equation

Under the assumptions of lossless, homogeneous, and isotropic medium, the propagation of thermoacoustic pressure field is governed by the classical acoustic wave equation. For a compressible fluid with density $\rho$, and the particle velocity $\mathbf{v}$, the conservation of mass (continuity equation) can be written as,

$${\partial }_t\rho +{\rho }_0\nabla .\mathbf{v}=0,$$

(1)

To analyze small acoustic perturbations $\rho ={\rho }_0+{\rho }_1$ and substituting in Equation (1) and neglecting higher order nonlinear terms,

$${\partial }_t{\rho }_1+{\rho }_0\nabla .{\mathbf{v}}_{\mathbf{1}}=0,$$

(2)

Using the linear acoustic relation between pressure ($p_1$) and density (${\rho }_1$) $p=c^2\rho$ and can rewrite the equation as,

$${\partial }_t{\rho }_1+{\rho }_0{c}^2\nabla .{\mathbf{v}}_{\mathbf{1}}=0,$$

(3)

From conservation of momentum and taking divergence on both sides,

$${\rho }_0{\partial }_t(\nabla .{\mathbf{v}}_{\mathbf{1}})=-{\nabla }^{2}p,$$

(4)

Substituting all terms and from Equation (1) the classical wave equation becomes,

$$\{{\nabla }^2-(1/c^2){\partial }_t^2)p(\mathbf{r},t)=0,$$

(5)

In the thermo-acoustic (TA) effect, optical absorption produces a localized temperature rise in the medium. This temperature increase causes the medium to undergo thermal expansion. However, the surrounding medium restricts this expansion, resulting in the generation of a thermoelastic pressure rise. This pressure rise subsequently propagates through the medium as an acoustic wave can be expressed as,

$$p(\mathbf{r},t)=\beta K▵T^0,$$

(6)

where $p(\mathbf{r},t)$ denotes the time-dependent thermoelastic pressure rise, $\beta$ is the volumetric thermal expansion coefficient of the medium. K represents the bulk modulus, and $▵T^0$ is the temperature increase caused by optical energy absorption. When the conditions of thermal and stress confinement hold, the resulting TA pressure in a homogeneous and non-viscous acoustic medium can be formulated with thermal diffusion time by the following equation,

$$\{{\nabla }^2-(1/c^2){\partial }_t^2)p(\mathbf{r},t)=-(\beta /C_p){\partial }_t\mathcal{H}(\mathbf{r},t),$$

(7)

where $C_p$ is the specific heat capacity at constant pressure, c is the speed of sound of the ambient medium, and heat deposition function $\mathcal{H}(\mathbf{r},t)$ per unit time and volume. Equation (7) is called TA wave equation.

2.2. Discretization Using FDTD

To numerically solve the TA wave equation, the FDTD method is employed. The continuous spatial and temporal derivatives are approximated using second-order central difference schemes. The second order temporal derivative is approximated as in 2-D space and time,

$$(1/c^2){\partial }_t^{2}p(x,y,t)=({\partial }_x^2+{\partial }_y^2)p(x,y,t),(x,y)\in \Omega \subset {\mathbb{R}}^2,t\in [0,T]$$

(8)

where the domain specifies $\Omega =[0,N_x]\times [0,N_y]$. $N_x$ and $N_y$ are the rectangular sides of the domain. Including an external source term,

$$(1/c^2){\partial }_t^{2}p(x,y,t)=({\partial }_x^2+{\partial }_y^2)p(x,y,t)+(\beta /C_p){\partial }_t\mathcal{H}(x,y,t),$$

(9)

The temporal discretization is carried out using a second order central difference scheme derived from the Taylor series expansion,

$${\partial }_t^{2}p\approx {\displaystyle \frac{(p_{i,j}^{n+1}-2p_{i,j}^n+p_{i,j}^{n-1})}{d{t}^2}},$$

(10)

Similarly, for spatial discretization,

$$\begin{array}{cc}\hfill {\partial }_x^{2}p& \approx {\displaystyle \frac{p_{i+1,j}^n-2p_{i,j}^n+p_{i-1,j}^n}{d{x}^2}},\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill {\partial }_y^{2}p& \approx {\displaystyle \frac{p_{i,j+1}^n-2p_{i,j}^n+p_{i,j-1}^n}{d{y}^2}},\hfill \end{array}$$

(11b)

And the Laplacian,

$$({\partial }_x^2+{\partial }_y^2)\approx {\displaystyle \frac{(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n+p_{i,j-1}^n-4p_{i,j}^n)}{d{x}^2}},$$

(12)

From Figure 1 each grid point $(i,j)$ depends on its four neighbors $(i+1,j),(i-1,j),(i,j+1),(i,j-1)$ and time levels $(n+1),n,(n-1)$, respectively. Substituting all values in Equation (9),

$$(1/c^2){\displaystyle \frac{(p_{i,j}^{n+1}-2p_{i,j}^n+p_{i,j}^{n-1})}{d{t}^2}}={\displaystyle \frac{p_{i+1,j}^n-2p_{i,j}^n+p_{i-1,j}^n}{d{x}^2}}+{\displaystyle \frac{p_{i,j+1}^n-2p_{i,j}^n+p_{i,j-1}^n}{d{y}^2}}+(\beta /C_p){\partial }_t{\mathcal{H}}_{i,j}^n,$$

(13)

To compute, the pressure at the next time step, equal spatial step sizes are considered in both the x and y directions, i.e., $dx=dy$. The update equation as follows,

$$\begin{array}{c}\hfill p_{i,j}^{n+1}-2p_{i,j}^n+p_{i,j}^{n-1}=(c^{2}d{t}^2/d{x}^2)(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n+p_{i,j-1}^n-4p_{i,j}^n)+(\beta /C_p){\partial }_t{\mathcal{H}}_{i,j}^n,\end{array}$$

(14a)

$$\begin{array}{c}\hfill p_{i,j}^{n+1}=2p_{i,j}^n-p_{i,j}^{n-1}+CF{L}^2(p_{i+1,j}^n+p_{i-1,j}^n+p_{i,j+1}^n+p_{i,j-1}^n-4p_{i,j}^n)+(\beta /C_p){\partial }_t{\mathcal{H}}_{i,j}^n,\end{array}$$

(14b)

where CFL=$cdt/dx$ is the Courant-Friedrichs-Lewy [19] number.

2.3. First-Order Mur Absorbing Boundary Condition

In FDTD simulations, the computational domain is artificially truncated, where physical wave propagation occurs in an unbounded medium. As a result, outgoing waves reaching the domain boundaries are reflected back into the domain, leading to non-physical artifacts. To mitigate this issue, absorbing boundary conditions are employed. One of the simplest and most wisely used ABCs is the first-order Mur absorbing boundary conditions, which approximates the behavior of an outgoing wave at the boundary. The first-order Mur ABCs can be written as,

$$\begin{array}{c}\hfill p_{1,j}^{n+1}=p_{2,j}^n-R\left(p_{2,j}^{n+1}-p_{1,j}^n\right),i=1;\end{array}$$

(15a)

$$\begin{array}{c}\hfill p_{N,j}^{n+1}=p_{N-1,j}^n-R\left(p_{N-1,j}^{n+1}-p_{N,j}^n\right),i=1,2,3,...,N;\end{array}$$

(15b)

$$\begin{array}{c}\hfill p_{i,1}^{n+1}=p_{i,2}^n-R\left(p_{i,2}^{n+1}-p_{i,1}^n\right),j=1;\end{array}$$

(15c)

$$\begin{array}{c}\hfill p_{i,N}^{n+1}=p_{i,N-1}^n-R\left(p_{i,N-1}^{n+1}-p_{i,N}^n\right),i=1,2,3,...,N;\end{array}$$

(15d)

Using the same spatial and temporal discretization as in the FDTD formulation, the boundary nodes are updated based on the adjacent interior grid points and previous time-step values. The reflection coefficient is defined as $R=(dx-cdt)/(dx+cdt)$, which can also be expressed in terms of the CFL number as $R=(1-CFL)/(1+CFL)$. Accordingly, the Mur boundary condition is applied at all four edges of the computational domain, namely the left, right, top, and bottom boundaries. This approach ensures stable wave propagation and significantly reduces non-physical reflections at the boundaries, thereby improving the overall accuracy of the simulation.

2.4. Stability Parameter Analysis

The stability of the numerical scheme is governed by the CFL condition, which is defined as the ratio between the temporal and spatial discretization, given by CFL$=cdt/dx$. The selection of spatial grid spacing ($dx,dy$) and time step $dt$ must satisfy certain constraints to ensure stable numerical computaion. If these constraints are violated, the numerical solution becomes unstable and loses physical reliability. The stability criterion can be derived from the finite difference fomulation, where the amplification factor depends on the term (1-2CFL²). For stability, this term must remain bounded, which leads to the condition that the CFL number must not exceed a critical value. consequentl, the stability requirement can be expressed as CFL$\le 1/\sqrt{2}$ for the two dimensional case.

It is important to note that the CFL number must be positive in practical implementations, as it presents a physical ratio of wave propagation speed to numerical discretization. The influence of the CFL condition on wave propagation is illustrated through numerical experiments. When the CFL number exceeds the allowable limit (e,g., CFL=0.8), the numerical solution becomes unstable, and the wave propagation is no longer physically meaningful. These observations confirm that maintaining the CFL condition is essential for obtaining stable and accurate results.

2.5. Initial Pressure Rise

2.5.1. Gaussian Source

The initial pressure rise in the computational domain is modeled using a two dimensional Gaussian distribution to represent a localized TA excitation. This approach simulates a spatially confined energy deposition, where the pressure attains its maximum at the source center and gradually decays in all directions.

$$p(x,y,0)=e^{-{\displaystyle \frac{{(x-N_x/2)}^2+{(y-N_y/2)}^2}{{\sigma }^2}}},$$

(16)

where ($N_x\times N_y=200\times 200$, and standard deviation: $\sigma =10$ controls the spatial spread of the source. The Gaussian profile ensures a smooth and continuous initial condition, minimizing numerical artifacts and high frequency noise in the simulation. Mathematically, the pressure distribution is defined by a standard deviation parameter. A moderate value of the $\sigma$ is selected to achieve abalance between spatial localization and spectral bandwidth. The initial pressure perturbation generates radially propagating acoustic waves with symmetric wavefronts, providing a stable and physically consistent starting condition for the FDTD based TA analysis.

2.5.2. First Degree Chebyshev Polynomial

An asymmetric source distribution is modeled using a first oreder Chebyshev polynomial based shape to represent a non-uniform TA excitation [20]. The geometry of the source is defined in polar coordinates, where the radial boundary is modulated using a cosine function derived from a Chebyshev type formulation. This allows controlled deformation of a circular profile into a multi-lobed structure. The resulting shape mimics realistic irregular energy deposition, enabling the investigation of directional wave propagation and asymmetryeffects in the acoustic field. To ensure numerical stability and avoid sharp discontinuities, the generated shape is further smoothed using a Gaussian filter before being used as the initial condition. Distance from center with angular coordinate $d=\sqrt{(}{(x-N_x/2)}^2+{(y-N_y/2)}^2)$ and $\theta =co{s}^{-1}((x-N_x/2)/d))$.

$$\mathcal{R}\left(\theta \right)={\mathcal{R}}_0[1+ϵ{\mathcal{O}}_n\left(cos\theta \right)],$$

(17)

where ${\mathcal{O}}_n\left(cos\theta \right)=cosn\Theta$ is the Chebyshev polynomial of degree n sometimes called waviness parameter, ${\mathcal{R}}_0$ is the radius of unperturbed sphere and $ϵ$ being a deformation parameter.

3. Numerical Experiment

In this study, a 2D TA wave propagation model is implemented using the FDTD method. The computational domain consists of uniform square grid of size $200\times 200$ with a spatial resolution of $dx=1$ mm. The TA wave speed is taken as $c=1500$ m/s. The time step $dt$ is determined based on the CFL condition, with a chosen CFL number of $0.5$ to ensure numerical stability. The TA wave field is initialized using a localized source distribution. In the primary configuration, a Gaussian source centered at the midpoint of the domain is employed, defined by standard deviation of $\sigma =10$ grid points. Alternative source geometries, including circular disk, first order Chebyshev polynomial based, and asymmetric (rock-like) sources, are also considered to investigate the influence of source shape on wave propagation characteristics.

To simulate an unbounded domain, first order Mur ABCs are applied at all boundaries, effectively reducing artificial reflections. The wavefield is updated iteratively using a second order explicit FDTD scheme. At each time step, the field values are computed using spatial finite differences, followed by the application of boundary conditions and update of the field variables. Point sensors are placed from toward the center at 60 grid points to record the temporal evolution of the acoustic pressure. In particular, sensors are positioned at specified radial distances and angular orientations relative to the source to capture directional propagation effects. The simulation is executed for a total 500 time steps. During the simulation, snapshots of the wavefield are recorded at selected time intervals to visualize the propagation dynamics. Additionally, the time domain signals at the sensor locations are collected and analyzed. Post processing is performed to extract both temporal and spectral characteristics of the acoustic response. This recorded sensor signals are used to compute amplitude-time histories, and their frequency spectra are obtained from fast Fourier transform (FFT). The analysis focuses on the effects of grid resolution, CFL number, and source geometry on wave propagation behavior, numerical dispersion, and spectral charateristics. The detailed workflow of the numerical procedure is presented in Figure 2.

3.1. k-Space Pseudospectral Method

In the k-Space pseudospectral method, acoustic wave propagation is modeled by decomposing the second order wave equation into a set of coupled first order equations. This formulation enables efficient computation of spatial gradients in the spectral domain while maintaing accurate time evolution. The acouustic field is described by pressure p and particle velocity components $u_x$ and $u_y$. The governing equations are,

$${\displaystyle \frac{\partial u_x}{\partial t}}=-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\partial p}{\partial x}},\&{\displaystyle \frac{\partial u_y}{\partial t}}=-{\displaystyle \frac{1}{{\rho }_0}}{\displaystyle \frac{\partial p}{\partial y}},$$

(18)

If equilibrium density ${\rho }_0$ and pressure density relation $p={\rho }_0{c}^2$ the continuity equation yeild as,

$${\displaystyle \frac{\partial p}{\partial t}}=-{\rho }_0{c}^2({\displaystyle \frac{\partial u_x}{\partial x}}+{\displaystyle \frac{\partial u_y}{\partial y}}),$$

(19)

In this framework, spatial derivatives are evaluated in the Fourier domain. The pressure field is transformed using the Fourier transform,

$$\tilde{p}(\mathbf{k},t)=\mathcal{F}\left[p(\mathcal{r},t)\right],$$

(20)

where $\mathbf{k}$ represent the wave number vector. The Laplacian operator in Fourier space becomes ${\nabla }^{2}p\to -k^2\tilde{p}$. To improve temporal accuracy, a correction factor is introduced. The pressure update equation in the spectral domain can be expressed as,

$${\tilde{p}}^{n+1}=2{\tilde{p}}^n-{\tilde{p}}^{n-1}+c^{2}d{t}^2(-k^2{\tilde{p}}^n)sin{c}^2(ckdt/2),$$

(21)

To prevent artificial reflections at the boundaries, a perfectly matched layer (PML) is implemented. Within the PML region, additional damping terms are introduced into the governing equations to absorb outgoing waves. The momentum equation with damping,

$${\displaystyle \frac{\partial u}{\partial t}}=-{\displaystyle \frac{1}{{\rho }_0}}\nabla p-\Lambda .\mathbf{u},$$

(22)

where $\Lambda ={\Lambda }_x+{\Lambda }_y$. ${\Lambda }_x$ and ${\Lambda }_y$ are the absorption coefficients at the PML region.

4. Results

The finite dimain responses obtained using both the FDTD method and k-Space pseudospectral method are presented in Figure 3, illustrating the variation of acoustic pressure with respect to time for different CFL numbers. For this analysis, a single point sensor is considered and positionsed at a distance of 60 grid points from the source center to capture the propagating wave signal in a simple and controlled configuration. The total propagation time is defined as the product of the time step ($dt$) and the total number of time iterations ($n_t=500$), i.e., $t=n_{t}dt$. This allows consistent comparison of temporal signals across different simulation setup. Figure 3a shows the comparison of pressure signals obtained using two different spatial grid resolutions, namely $dx=1$ mm and 50 mm, while maintaining a constant CFL number that satisfies the stability condition. Under these conditions, it is observed that the numerical scheme preserves the propagation characteristics when the CFL number is properly maintained.

This demonstrates the consistency and reliability of the numerical implementation across different spatial discretization. In Figure 3b presents the results when different CFL numbers are used for the two grid resolutions. Since the CFL condition directly depends on both the time step and spatial grid size, varying the CFL number alters the stability and accuracy of the numerical scheme.As a result the coressponding pressure signals exhibit noticeable differences, including variations in amplitude and phase.

Figure 4a illustrates the initial two dimensional Gaussian source distribution at time $t=0$ sec, showing the localized pressure excitation within the computational domain. In this case the CFL number is fixed at 0.5 to ensure numerical stability. The sensor position is kept consistent with the previous configuration, located at a distance of 60 grid points from the source center. Figure 4b presents the wave propagation map at time $t=6.67\times {10}^{-5}$ s, demonstrating the outward radial expansion of the wavefront and its interaction with the domain boundaries. The propagation remains smooth and symmetric, indicating stable numerical behavior. The corresponding time domain signals, shown in Figure 4c, reveal excellent agreement between the FDTD and k-Space pseudospectral methods, confirming the accuracy of both approaches under stable conditions. Furthermore, the frequency domain representation in Figure 4d highlights the spectral characteristics of the propagated signal. A distinct minimum is observed around $f=1.0\times {10}^5$ Hz, indicating consistent frequency behavior between the two methods. Overall, the results demonstrate strong agreement in both time and frequency domains. A similar analysis is performed for the circular disk source with a radius of $r=20dx$, as shown in Figure 5. In this case, both numerical methods again exhibit excellent agreement in the time domain signals. However, slight discrepancies are observed at higher frequencies in the spectral domain. The frequncy spectrum shows multiple minima that become more prominent with increasing frequency, with the first minimum occuring around $f=0.5\times {10}^5$ Hz. Figure 6 and Figure 7 extend the analysis to more complex source geometries, including a first order Chebyshev polynomial based source and an asymmetric (rock-like) source. In both cases, the k-Space pseudospectral method shows very good agreement with the FDTD results.

Figure 8a–f are shown at different time instances for a fixed CFL number of the pressure distributions. Figure 9 shown variation of wave propagation distance with time under (a) satisfied CFL condition and (b) violated CFL condition, illustrating the impact of stability on wave propagation behavior. Additionally, the Figure 10 is the time domain responses of the TA pressure recorded at two distinct sensor locations for a Gaussian source under a fixed CFL condition.

5. Discussion and Conclusions

The present study demonstrates that both the FDTD and k-Space pseudospectral methods provide consistent and reliable solutions for TA wave propagation under stable numerical conditions. For a constant CFL number within the stability limit the time domain signals and wavefield distributions show excellent agreement across different grid resolutions and source configurations. The Gaussian source produces symmetric wave propagation, while more complex geometries such as circular, Chebyshev-based, and asymmetric (rock-like) sources introduce variations in wavefront structure and directional propagation. The observed agreement between the two numerical methods indicates that both approaches accurately capture the underlying physics of wave propagation. The radial symmetry of the Gaussian source leads to uniform wavefront expansion, whereas deviations from symmetry in other source types result in directional energy distribution and asymmetric wavefields. The presence of multiple minima in the frequency spectra reflects the inherent spectral characteristics of the source geometry. The k-Space pseudospectral method exhibits slightly improved performance in the frequency domain, particularly at higher frequencies due to its reduced numerical dispersion. The FDTD method simpler and computationally efficient shows minor discrepancies at higher frequencies especially for complex source geometries. Both methods show strong agreement in the time domain under stable CFL conditions. When the CFL condition is violated significant deviations are observed in the wave propagation behavior, including distorted signals and unstable numerical solutions. The local approximation used in FDTD introduces phase errors that become more pronounced at higher frequencies and for coarse grid resolutions. Both methods are suitable for modeling TA wave propagation.

The simulations are performed on a two dimensional domain which may not fully capture three dimensional wave behavior. The use of first order Mur ABCs may introduce minor reflections, affecting long time simulations. The assumption of homogeneous medium also limits the applicability of the results to more complex real world scenarios. Future work may focus on extending the model to large scale and three dimensional domains and incorporating heterogeneous media to better represent realistic environments. The implementation of more advanced ABCs such as PML could further improve boundary performances.

The study demonstrate that careful selection of numerical parameters, particularly the CFL condition, enables accurate and stable simulation of TA wave propagation across a range of source configurations and numerical methods.