Empirical Modeling of Subcritical Hopf Bifurcation of Thermoacoustic Stirling Engine

Abstract

This study models the subcritical Hopf bifurcation in thermoacoustic Stirling engines using the Stuart–Landau model, highlighting the role of nonlinear dynamics. By inducing self-sustained oscillations and measuring pressure fluctuations across different temperature gradients imposed on the regenerator, we reveal the engine’s transition to a nonlinear domain, characterized by heightened oscillation amplitudes and unique periodic patterns. Interpreted Landau constants and growth rates illuminate the stabilizing effects of nonlinear dynamics, demonstrating the Stuart–Landau model’s applicability in thermoacoustic engine analysis. Our research confirms that this empirically refined model reliably describes oscillation amplitudes and transient phenomena, contributing valuable perspectives for advancing thermoacoustic engine design and operational understanding.

1. Introduction

In dissipative fluid systems, spontaneous thermoacoustic oscillations undergo transitions between dynamical states through bifurcations, leading to significant amplitude oscillations indicative of a self-sustained state, known as thermoacoustic instability. This phenomenon arises from thermally nonlinear interactions between sound waves and solid walls in confined flow channels, driven by steep axial temperature gradients, and precipitates through Hopf bifurcations. It is possible to apply these bifurcations of thermoacoustic oscillations as a natural heat engine [1].

The evolution of thermoacoustic systems into qualitatively distinct behaviors is termed bifurcation [2]. Changes in system parameters, such as the axial temperature gradient, can dramatically shift their qualitative nature. These shifts are crucial as they can create or destroy fixed points, such as stable fixed points or limit cycle oscillations, fundamentally altering the system’s dynamics. The certain values of parameters at which these changes occur are known as bifurcation points. Understanding bifurcations is essential for gaining insights into the dynamic behavior of thermoacoustic engines and crucial for predicting and managing thermoacoustic instability

Within thermoacoustic engines, both supercritical and subcritical Hopf bifurcations are observed [3]. An increase in the temperature gradient results in a supercritical Hopf bifurcation, leading to a gradual increase in oscillation amplitude, whereas a subcritical Hopf bifurcation causes an abrupt amplitude jump at the Hopf point, where the system’s stable steady-state becomes unstable and a stable limit cycle emerges. Subcritical cases have been extensively documented across various thermoacoustic engine configurations, including standing wave, traveling wave, phase change, and free-piston Stirling types [3,4,5,6,7,8,9,10,11]. Furthermore, a decrease in the temperature gradient reveals hysteresis in oscillation amplitude with subcritical Hopf bifurcations, eventually ceasing through a saddle-node bifurcation at a distinct fold point [2], which signifies the minimum temperature gradient necessary for engine operation.

The Bifurcation of the thermoacoustic engine system, in terms of triggering and ceasing thermoacoustic instability, can be described by Rott’s linear theory [12]. Based on the linearized hydrodynamic equations, Rott’s theory provides a theoretical framework for analyzing thermoacoustic engine performance using wave equations in the frequency domain. This framework enables us to determine the growth/attenuation rate of oscillation amplitudes, crucial for estimating non-combustion-driven thermoacoustic instability.

While Rott’s linear theory lays a solid foundation for analyzing thermoacoustic phenomena, predicting frequencies and growth rates consistent with supercritical Hopf bifurcations [13], it inadequately addresses the nonlinear dynamics crucial for understanding subcritical bifurcations. This oversight underscores the importance of delving into the nonlinear dynamics that govern the complex behaviors observed in thermoacoustic engines, particularly those leading to hysteresis effects.

Further, thermoacoustic systems exhibit rich nonlinear dynamics [14], such as limit cycles, quasi-periodic oscillations, harmonics, shock waves, synchronization, and more. While computational fluid dynamics (CFD) methods effectively capture these nonlinear phenomena [15], they do come with significant computational costs. Recently, we have refined Rott’s model by integrating empirically derived nonlinear flow resistance [3]. This addition temporarily characterizes the reduction of the linear growth rate, enhancing Rott’s frequency theory to better describe the evolution of amplitude saturation from onset to limit cycles. Still, integrating every nonlinear effect emerging in thermoacoustic systems into the linear theoretical framework presents significant challenges. This difficulty is perceived due to the complex and interdependent nature of these nonlinear effects. It highlights the necessity of developing a simplified nonlinear model to depict their collective impact on the thermoacoustic system.

The Stuart–Landau equation (SLE) [16,17,18] is instrumental in describing the bifurcation dynamics within fluid dynamic systems, offering insights into the nonlinear stabilizing effects within the system. It is instrumental in studying a variety of physical systems undergoing bifurcation, especially in systems that exhibit oscillatory behavior, such as jet and shear flow instabilities and bluff body wake instabilities [19,20,21,22,23,24]. Additionally, the SLE plays a crucial role in modeling significant nonlinear phenomena like frequency lock-in [25] and resonance in flow systems [26]. In studies of combustion-driven thermoacoustic systems [27,28], the SLE captures the weakly nonlinear dynamics near bifurcation points. The criticality of the bifurcation is determined by the sign of the coefficient of the leading nonlinearities, known as a Landau constant, which underlines the equation’s utility in assessing the stabilizing effects of nonlinearities within the system.

For the case of the supercritical Hopf bifurcation observed in a standing wave thermoacoustic system, Biwa et al. [29] experimentally derived the evolution equation by measuring energy and energy flow under external oscillatory perturbations at the limit cycle frequency. This study seeks to expand the SLE model’s application, also referred to as the evolution equation, for the empirical determination of Landau constants to explore the critical behavior of Hopf bifurcations leading to self-sustained oscillations in a thermoacoustic Stirling engine. By focusing on the nonlinear dynamics, especially those governing subcritical Hopf bifurcations, we aim to shed light on their implications for engine operation. In the following, we begin by detailing the theoretical background and methodology of our study. Section 3 outlines the experimental setup designed for our thermoacoustic Stirling engine. Section 4 then presents our findings and their analysis. Finally, we conclude by summarizing our key insights in Section 5.

2. Stuart–Landau Model

The Landau model [18] offers a framework for analyzing the nonlinear dynamics around phase transitions. It is extensively used to describe and categorize the bifurcation behaviors in fluid systems. Additionally, there have been several studies on using this model for understanding combustion-driven thermoacoustic instabilities in Rijke tubes. A detailed explanation of the Stuart–Landau model, which is derived from the Navier–Stokes equations, is provided in Stuart [16,17]. The model is read as

$$\begin{array}{c}\hfill \frac{dA}{dt}=A(\sigma +i\omega )+A\sum _{k=1}^{\infty }{(-1)}^k{l}_k{\left|A\right|}^{2k},\end{array}$$

(1)

where $A\left(t\right)$ represents the complex amplitude of the oscillation mode, $\sigma$ is the linear temporal growth rate, $\omega$ is the linear angular frequency, and the complex coefficients $l_k$ are called the kth Landau constants that adjust the oscillation frequency at saturation. Landau constants $l_k$ describe nonlinear interactions in systems near instability. Derived by linearizing and then expanding system equations around a base state, they show how a system’s behavior changes near critical points. Note that due to the inherent symmetries of the oscillating system’s wavefront, only odd terms of the complex amplitude A are permitted on the right-hand side of the SL model [18]. The Landau constants, as key nonlinear coefficients, moderate system responses significantly by counteracting amplitude growth, thus ensuring steady oscillation states. These coefficients are critical in quantifying nonlinear effects and are typically determined through experimental data [28] or theoretical calculations [27,30] based on physical models.

Presenting the complex amplitude $A\left(t\right)$ in modulus and argument form as

$$\begin{array}{c}\hfill A\left(t\right)=\left|A\right(t\left)\right|exp\left[i\varphi \right(t\left)\right],\end{array}$$

(2)

where $\left|A\right(t\left)\right|$ represents the instantaneous amplitude of the mode and $\varphi \left(t\right)$ its phase, both being real variables. Note that $\left|A\right|$ and $\varphi$ are experimentally determined in the following Section 3 for estimating the Landau constants. Substituting $A\left(t\right)$ into Equation (1) and focusing on the real part yields the Landau amplitude equation:

$$\begin{array}{c}\hfill \frac{1}{\left|A\right|}\frac{d\left|A\right|}{dt}=\frac{d\left(ln\left|A\right|\right)}{dt}=\sigma -{l\left|A\right|}^2+q{\left|A\right|}^4-\cdots ,\end{array}$$

(3)

where l and q are the real part of Landau constants for the cubic and quintic terms, respectively. The linear growth rate $\sigma$ is emphasized to depend on the bifurcation control parameter [22]: positive values indicate instability and amplitude growth, while negative values suggest stability and decay of oscillations. The equation is typically truncated after the cubic term to represent supercritical Hopf bifurcations effectively. This approach stems from the cubic term, with l being positive, being generally sufficient to moderate the initial exponential growth, leading to saturation, thereby ensuring the model accurately captures the nonlinear dynamics near transitions, especially for small saturation amplitudes [21,22]. If l is negative, however, higher-order terms become essential for depicting saturation due to enhanced perturbation growth by the cubic term q. The sign of l determines the transition type: positive l indicates a supercritical (smooth) transition, while negative l signifies a subcritical (hysteretic) transition [31]. The linear growth rate, $\sigma$, marks the stability change at the transition point by shifting the sign. In non-combustion-driven thermoacoustic systems, this control parameter is the temperature difference $\Delta T$ across the regenerator. Clearly, Equation (3) describes the dynamics of the temporal growth rate, which are tempered by high-order terms corresponding to oscillation amplitudes. It should be noticed that, in the scenario of a subcritical bifurcation, the temporal growth of the amplitude is saturated by the higher-order nonlinear terms, necessitating the inclusion of at least quintic terms.

3. Experimental Setup and Mehod

Figure 1 presents a schematic diagram of the experimental setup for the present study. This experimental setup was aimed to investigate subcritical bifurcation behaviors for a thermoacoustic Stirling engine, featuring a looped tube and a branch resonator made from a 0.03-m-inner diameter stainless-steel cylindrical pipe. The loop has an average length of $L=1.68$ m, with the branch resonator of $L_b=$1.875 m opened to the atmosphere, operating with air at atmospheric pressure. The looped section includes various components: resonator tube 1, a cold heat exchanger, a stacked-screen regenerator, a hot heat exchanger, a thermal buffer tube, and resonator tube 2, arranged along the axial coordinate x from the T-junction. The 20 mm long regenerator comprises stainless-steel woven mesh screens of # 30 mesh number, with a wire diameter of 0.22 mm, hydraulic diameter of 0.82 mm, and a volume porosity of 0.79.

To impose temperature gradients on the regenerator, electrical cartridge heaters and cooling water are used in the hot and cold heat exchangers, respectively, both fitted with brass fins 1 mm apart and a porosity of 0.667. Type-K thermocouples monitor temperatures at both ends of the regenerator to determine $T_H$ and $T_C$ and hence the temperature difference $\Delta T=T_H-T_C$, recorded by a temperature recorder (LR8450, HIOKI, Ueno, Japan), ensuring precise temperature control. Pressure oscillations within the gas column are measured at $x=0.24$ m using a pressure transducer (PMS-5M-2-1M, JTEKT, Nagoya, Japan). The signals are then amplified by a DC Amplifier (AA6210, JTEKT, Nagoya, Japan) and analyzed with a fast Fourier transform analyzer (OR35, OROS, Grenoble, France) to track the transition from initial to steady oscillation states. For data collection, we utilize the FFT spectrum analyzer to capture the time series of the pressure signal $p^{\prime }\left(t\right)$ at a sampling rate of $1/\tau =3276.8$ Hz. It is important to note that we preprocess the measured data by applying band filtering to $p^{\prime }\left(t\right)$ to reduce noise. We isolated measured $p^{\prime }\left(t\right)$ within the 22 to 26 Hz range, determined by the engine’s spontaneous oscillation frequency, roughly estimated at $24.1$ Hz using $a/(L+L_b)\times 1/4$ [3], where a is the adiabatic sound speed at room temperature.

The experimental procedure begins with the initiation of spontaneous gas column oscillations by establishing steep temperature gradients across the regenerator, ensuring thermal steadiness by maintaining a constant heat power supplied to the cartridge heaters. Throughout the experiment, we closely monitor the stability limits and pressure oscillation behaviors. To experimentally observe the subcritical bifurcation behaviors, we systematically reduce the supplied heat power once steady pressure oscillations are achieved, recording the smaller pressure oscillations steadily in the reverse path. This approach facilitates the observation of hysteresis loops of the thermoacoustic Stirling engine.

To experimentally obtain Landau coefficients, we examine the transient behavior under a subcritical bifurcation scenario. This involves setting up constant temperature differences across the regenerator within bistable regions, located between the fold and the Hopf points, as well as above the Hopf bifurcation point. Once bistable regions are reached, we introduce external airflow disturbances at the branch resonator’s open end to induce pressure oscillations $p^{\prime }\left(t\right)$. Note that the airflow uses a manually operated air blow gun. Air is applied at a pressure of 5 kg-f/cm² for approximately 0.2 seconds to initiate the engine. Depending on the magnitude of external disturbances and their distance from subcritical Hopf bifurcation points, oscillation amplitudes can either grow over time and saturate finite amplitudes or decrease and eventually fade away. In contrast, we maintain steady temperature differences above the Hopf point while sealing the end of the branch resonator with a plate to prevent the excitation of spontaneous oscillations. Subsequently, we remove the plate to observe the initial growth of pressure oscillations, leading to finite amplitudes.

With the measured pressure fluctuations $p^{\prime }\left(t\right)=[p_j^{\prime },p_{j+\tau }^{\prime },p_{j+2\tau }^{\prime },\cdots ,p_{j+(d-1)\tau }^{\prime }]$, where j signifies the initial time index and d represents the number of recorded data points, we obtain its normalization $p^{*}\left(t\right)$ using the mean value $\overline{p}$ and the standard deviation $sd$, which can be expressed as follows:

$$\begin{array}{c}\hfill p^{*}=\frac{p^{\prime }-\overline{p}}{sd}.\end{array}$$

(4)

$p^{*}\left(t\right)$ are characterized by a complex angular frequency, ${\omega }_n={\omega }_{\mathrm{R}}+i{\omega }_{\mathrm{I}}$, with ${\omega }_{\mathrm{R}}$ representing the oscillation’s angular frequency and ${\omega }_{\mathrm{I}}$ its growth rate. We perform a Hilbert transform (HT) on the normalized pressure fluctuations $p^{*}\left(t\right)$ using MATLAB’s built-in function (version R2021b, Mathworks, Inc., Natick, MA, USA). This converts the normalized pressure signal $p^{*}\left(t\right)$ of the real-axis to an imaginary-axis signal $p_{\mathrm{H}}\left(t\right)$, enabling the calculation of the normalized oscillation magnitude $\left|A\left(t\right)\right|=\sqrt{{p^{*}}^2\left(t\right)+p_{\mathrm{H}}^2\left(t\right)}$ from the experimental data. The obtained time series of the instantaneous amplitude $\left|A\left(t\right)\right|=\left[\left|A_j\right|,\left|A_{j+\tau }\right|,\left|A_{j+2\tau }\right|,\cdots ,\left|A_{j+(d-1)\tau }\right|\right]$ enables us to numerically calculate the time derivative $d\left(ln\left|A\right|\right)/dt$. The time derivative $d\left(ln\left|A\right|\right)/dt$ is approximated from its finite difference $\Delta (ln\left|A\right|)/\Delta t$ in the time series as

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d\left(ln\left|A\right|\right)}{dt}}\approx {\displaystyle \frac{\Delta (ln\left|A\right|)}{\Delta t}}\hfill \\ \hfill =& \left[\frac{ln\left|A_{j+\tau }\right|-ln\left|A_j\right|}{\tau },\frac{ln\left|A_{j+2\tau }\right|-ln\left|A_{j+\tau }\right|}{\tau },\cdots ,\frac{ln\left|A_{j+(d-1)\tau }\right|-ln\left|A_{j+(d-2)\tau }\right|}{\tau }\right]\hfill \end{array}.$$

(5)

Note that Equation (5) can also be determined using MATLAB’s built-in function “gradient”. This experimental data analysis provides insight into the oscillation’s growth rate, ${\omega }_{\mathrm{I}}$, realizing from the Landau amplitude Equation (3). By plotting $d\left(ln\left|A\left(t\right)\right|\right)/dt$ against ${\left|A\left(t\right)\right|}^2$, as shown in the latter, we can ascertain the Landau constants from Equation (3) through least-squares fitting of polynomial curves. In this study, we utilize MATLAB’s built-in function “lsqcurvefit” to perform this fitting.

4. Result and Discussion

Figure 2 depicts the hysteresis loop of saturated pressure oscillations measured across the regenerator under various steady temperature differences $\Delta T$. Two distinct paths, the forward (indicated by red markers △) and reverse (indicated by blue markers ▽), illustrate the steady oscillations of acoustic pressures. Notably, upon steadily increasing temperature differences above the engine onset, the forward path reveals significant finite amplitudes of acoustic pressure, culminating at the Hopf bifurcation point. Conversely, during the reverse process, characterized by decreasing temperature differences after engine initiation, thermoacoustic oscillations persist until approaching the lower critical temperature difference, denoted as the saddle-node bifurcation point. Figure 2 demonstrates nearly identical steady amplitudes of pressure oscillations in both forward and reverse paths above the Hopf point, indicative of subcritical Hopf bifurcation behavior in the current setup of the thermoacoustic Stirling engine.

Within the region delimited by the saddle-node and Hopf points, commonly referred to as the “bistability region”, we set steady $\Delta T$ before thermoacoustic oscillations to explore Landau constants, aiming to investigate the potential for empirical modeling. This exploration includes acknowledging the variations in temperature differences once oscillations are initiated, which arise due to the axial heat transfer caused by acoustic pumping [32], and can impact the attainment of saturated steady-state pressure amplitudes and steady temperature differences in the system. Within this bistability region, Figure 3 illustrates the triggering of original self-sustained pressure oscillations by external airflow disturbances and bandpass-filtered $p^{\prime }$ mainly oscillating with 23 Hz. Notably, the filtered $p^{\prime }$ excludes the pulse resulting from external disturbances, and its clear wavefront, as depicted in the inlet of Figure 3, enables us to pursue further investigations. Capitalizing on this, we utilize least-squares fits with the Stuart–Landau model to determine Landau constants.

We examine the evolution of pressure oscillations at specific temperature differences: within the bistability region at 339 K and 341 K, and above the Hopf bifurcation point at 343 K, 345 K, and 347 K. These conditions are set to explore the system’s response in thermally steady states before the initiation of oscillations. The experiments were replicated up to three times to ensure reliability. The analysis employs Hilbert transform techniques to reveal the nonlinear transition characteristics of the normalized pressure oscillations $p^{*}$. These transitions are categorized into subfigures of Figure 4 for clear comparison: Figure 4a,b represent the conditions within the bistability region at 339 K and 341 K, respectively; while Figure 4c–e correspond to conditions above the Hopf bifurcation point at 343 K, 345 K, and 347 K, respectively. Each subfigure consists of two parts: (i) traces the evolution of the normalized instantaneous amplitudes $\left|A\right(t\left)\right|$ from their initiation to the maximum amplitude phase, and (ii) displays the corresponding growth rates versus ${\left|A\right|}^2$.

Due to the necessity of external perturbations for the excitation of self-sustained oscillations for the bistability region (339 and 341 K), the analysis indicates that the growth rates exhibit positive values at specific amplitudes ${\left|A\right|}^2$, decreasing to zero as the amplitude reaches saturation. In contrast, for the conditions above the Hopf bifurcation point (343 K, 345 K, and 347 K), growth rates are positive even at lower amplitudes and increase to a peak before declining to zero with further amplitude growth.

Furthermore, the growth rates derived from the time series pressure data analyzed via the Hilbert transform are fitted using a model that extends to quintic terms of the Landau amplitude Equation (3). This approach facilitates the determination of linear growth rate $\sigma$ and the Landau constants l and q, which are used to plot fitted dash-dotted curves within each subfigure (ii) in Figure 4. These curves, based on coefficients averaged at every experimental iteration, suggest that while the fit to quintic terms may not perfectly align with the experimental observations at 345 K and 343 K, the latter comparison of the calculated instantaneous amplitudes, based on the Stuart–Landau equation, with measured values justifies the model’s overall appropriateness.

The derived linear growth rate $\sigma$ and Landau constants l and q are presented in

Table 1. Linear growth rate and Landau constants determined from Equation (3), extended to incorporate quintic terms.

	$\mathbf{\Delta }\mathit{T}$ (K)	$\mathit{\sigma }$	l	q
Within the bistability region	339	−0.015326	−0.028968	−0.0070583
	341	0.010792	−0.061028	−0.02404
Above the Hopf bifurcation point	343	0.010044	−0.058972	−0.013558
	345	0.024856	−0.10616	−0.039699
	347	0.10994	−0.04836	−0.037495

. Observations indicate that $\sigma$ increases with rising $\Delta T$. The cubic term l consistently shows negative values, indicating an acceleration in growth as oscillation amplitudes increase. Consequently, the negative quintic term q plays a crucial role in stabilizing the system as amplitudes grow. While the growth rate $\sigma$ can potentially be described within Rott’s linear theory framework due to its linear nature [3], the constants l and q represent the system’s nonlinear aspects. Furthermore, we observe that the constants l and q are dependent on temperature differences across the regenerator. This dependence might be attributed to nonlinear viscous and heat transfer effects in the tiny flow channels of the regenerator, which are sensitive to local temperatures. The weakly nonlinear theoretical study of the thermoacoustic engine [33], though limited in describing the supercritical bifurcation, shows that the analytical evolution equation clearly exhibits a temperature dependence on perturbation parameters. This would explain the observed $\Delta T$ dependence of the Landau constants obtained in this study. Addressing the theoretical implications of l and q, like weakly nonlinear theoretical studies conducted by Subramanian et al. [27] and Orchini et al. [30], will be an important focus of future research. In contrast to theoretical studies, the empirical modeling results presented in this study would contribute to guide the theoretical analysis for developing the evolution equation using the Stuart–Landau model for the thermoacoustic Stirling engine.

Upon empirically deriving Landau constants from Figure 4 and listed in Table 1, we utilized MATLAB’s “ode89” numerical solver to tackle the Landau amplitude equation:

$$\begin{array}{c}\hfill \frac{d\left|A\right|}{dt}=\sigma \left|A\right|-{l\left|A\right|}^3+q{\left|A\right|}^5.\end{array}$$

(6)

This equation allowed us to model the evolution of acoustic pressure oscillation amplitudes, adjusting these by scaling with the standard deviations $sd$ and adding the mean pressure value $\overline{p}$ to account for normalization in the time series pressure analysis. Note that we use the default settings for the ode89 solver, with a relative error tolerance (RelTol) of ${10}^{-3}$ and an absolute tolerance (AbsTol) of ${10}^{-6}$. These adjusted computational results were then directly compared with the time series pressure measurements shown in Figure 5. Additionally, the figures include the evolution of $\Delta T$ for context, indicating how the initiation of oscillation pressures leads to acoustic pumping of the heat flux, which slightly reduces $\Delta T$ over time. It is also observed that the amplitude reaches finite saturation after achieving peak pressure values, highlighting the system’s dynamic stability. Notably, all changes in $\Delta T$ over time are kept within a 5 K range, emphasizing the system’s consistent response to induced oscillations.

Within the bistability region at $\Delta T=339$ K and 341 K, we utilized the initial experimental pressure amplitude as the initial condition for solving the Landau model described by Equation (6). The close alignment between the numerical solutions and experimental data validates the Landau model’s capability to accurately capture the essence of self-sustained pressure oscillations. For scenarios above the Hopf bifurcation point, specifically at 343 K, 345 K, and 347 K, the model was solved using minimal initial amplitudes. It effectively captures the system’s transition into a state of increased nonlinearity, evidenced by the growth in oscillation amplitudes. This phase is marked by augmented amplitudes and clear periodic patterns. Despite slight variances in transient phases between the model’s predictions and experimental findings, the critical aspects, such as the trend and periodicity of the oscillations, are precisely depicted. The slight discrepancies observed are likely due to the system’s inherent nonlinearity, which suggests the need for incorporating terms beyond quintic in the Stuart–Landau Equation for accurately modeling the transient phase. However, the current model up to quintic terms suffices for depicting the saturated amplitude.

Also, to further validate our model’s accuracy in depicting bifurcation behavior, we conducted numerical calculations using Equation (6) under varying initial conditions within the bistable region at $\Delta T=339$ K. Figure 6 illustrates the system’s responses, calculated using the ode89 solver by varying sets of the relative error RelTol and AbsTol. The results demonstrate consistent finite amplitudes or stable states (no oscillation) across different initial conditions of pressure when setting smaller error tolerances, thereby reinforcing the robustness and reliability of our approach in modeling the critical dynamics of thermoacoustic systems. These findings clarified the potential impacts of initial conditions on the numerical integrations and highlighted the necessity of rigorous numerical methods to ensure accurate empirical modeling. This convinces us of the empirical modeling approach’s effectiveness in capturing the dynamics of the system near critical thresholds.

Besides SLE, studies of combustion-driven thermoacoustic instability offer alternative approaches for nonlinear modeling into thermoacoustic Stirling engines. Noiray and Schuermans [34], Noiray [35] explore deterministic quantities and linear growth rates within thermoacoustic systems, applying a noise-driven Van Der Pol oscillator and analyzing acoustic signal envelopes to elucidate the complex interactions in combustors. Seshadri et al. [36] further this exploration by modeling the effects of vortex shedding on oscillation amplitudes, providing insights into instability control through system design. Dutta et al. [37] employ the Kuramoto model to examine synchronization dynamics in combustors, highlighting transitions between unstable and stable states. Additionally, Bonciolini et al. [38] study how rate changes can delay transitions in thermoacoustic systems. Together, these studies deepen our understanding of nonlinear dynamics in thermoacoustics and could offer valuable insights for advancing research on thermoacoustic engines.

5. Conclusions

This study effectively utilized the Stuart–Landau equation for empirical modeling to investigate the dynamics of the thermoacoustic Stirling heat engine. Through the analysis of saturated pressure oscillations, from their initiation to finite amplitudes using filtered pressure signals, and by observing the oscillation evolution across different temperature gradients, we have verified the equation’s capability in capturing the intricate nonlinear behavior of the thermoacoustic system. The study shows that the Stuart–Landau equation, particularly when extended to quintic terms, effectively captures the system’s shift to a more nonlinear regime, marked by increased oscillation amplitudes and distinct periodic patterns. In summary, this research validates the Stuart–Landau equation’s applicability in thermoacoustic engine modeling, emphasizing the importance of the empirical modeling approach in understanding the underlying dynamics associated with subcritical Hopf bifurcations. This study highlights how the Stuart–Landau equation, with quintic term extensions, precisely captures the transition to a more nonlinear regime in thermoacoustic engines, evident from increased oscillation amplitudes at temperature differences exceeding the Hopf bifurcation point. It emphasizes the equation’s effectiveness in detailing the dynamics specific to subcritical Hopf bifurcations and underscores the significance of empirical modeling. These findings are crucial for enhancing the design and control of thermoacoustic systems, offering a path forward for advanced theoretical modeling of nonlinear behaviors under the Stuart–Landau framework.