Quantitative and Qualitative Analysis of Main Effects and Interaction Parameters in Thermoacoustic Refrigerators Performance

1. Introduction

Thermoacoustic refrigerators are devices, that use sound waves to generate cooling without moving parts. Due to the complex relationship between the multiple variables that determine its performance, optimization studies are required from a practical point of view. Many proposals are currently being explored to optimize thermoacoustic refrigerators, not only the performance response to environmental concerns related to many refrigerants but also to obtain fast and simple engineering estimates. The performance response depends primarily on the system’s design parameters, including operating conditions, the structural geometry of stack/regenerator, thermo-physical properties of the stack, working fluid, and acoustic driver. These design parameters can be studied and optimized from different perspectives and study groups/areas such as:

• Energy efficiency and thermodynamics: Improving the energy efficiency of the thermoacoustic cycle involves studying heat transfers, energy losses and conversion efficiencies in the cooling process [1]–[11,25].
• Materials and structural design: Optimizing the materials of the acoustic resonator and internal components can enhance the durability and efficiency of the thermoacoustic cooler. This involves studies on composite materials, aerogels [17], thermal conductivity, mechanical resistance, and durability under operating conditions [12]–[14].
• Acoustics and vibrations: Since thermoacoustic refrigerators use sound waves to generate pressure changes and thus cooling, they can be studied from the perspective of acoustics and vibrations. This includes optimizing acoustic resonator geometry and system design to maximize acoustic-to-thermal energy transfer efficiency [9,10,15,16,20].
• Control and electronics: Electronic control of the thermoacoustic cycle can be crucial to optimize the stability and performance of the refrigerator. Control algorithms, temperature and pressure sensors, and the integration of electronic components can be studied to improve system accuracy and efficiency. [18]

In the literature, the efficiency of a thermoacoustic device is given by the Coefficient of Performance ($COP$) [1,17] defined as the ratio between heat flow ($Q_c$) and the acoustic work (W) on the stack. From an energy efficiency optimization perspective, the $COP$ generated in the stack must be maximized. Regarding geometrical optimization parameters, several studies have been carried out, including intelligent optimization techniques such as Genetic Algorithms [1,7] Response Surface methodology (RSM) [15,16], Nonlinear programming models with discontinuous derivatives (DNLPs) [6], among others. Most of these works have been based on the linear thermoacoustic theory [9]-[10]. Using dimensionless analysis and short stack and boundary layer approximations, Tijani [9] obtained dimensionless expressions for $Q_c$ and W; which depend on eight dimensionless parameters. Despite these simplifications, to the authors’ knowledge, no study has yet been carried out to identify and quantify the contribution of each of these parameters, as well as the cross-effect of their interactions. When the working gas is chosen, the $COP$ depends on five parameters, namely D, B, $\Delta T_m$, $X_s$, $L_s$, see Figure 1.

In this work we present results from the application of the design of experiments (DOE) [21] using MINITAB [22], to quantitatively and qualitatively estimate the individual contribution of the five parameters described above to the $COP$ and their interactions. This study builds upon the methodology outlined by Tijani [9], readers interested in a detailed development of the theory are encouraged to refer to this work.

2. Methods

2.1. COP Estimation for Thermoacoustic Devices

The $COP$ is a measure of the efficiency of a refrigeration system, heat pump, or heat engine. In the context of the efficiency of transferring acoustic energy from the speaker to the environment, the $COP$ and acoustic power are related [9]-[10].

The performance of a thermoacoustic system depends strongly on the characteristics and properties of the porous medium (referred to in the Literature as the stack) used to exchange heat with the working fluid, see Figure 1.

For the development of a thermoacoustic refrigerator, the basic theory and equations for linear standing-wave thermoacoustic systems have been extensively developed by Tijani [9]-[10]. The term “regenerator” is used when the pore size of the porous element is small compared to the depth of thermal penetration. The thermal penetration depth refers to the layer around the regenerator where the thermoacoustic phenomenon occurs.

The thermal ${\delta }_{\kappa }$ and viscous ${\delta }_{\nu }$ penetration depths are defined as follows:

$${\delta }_{\kappa }=\sqrt{{\displaystyle \frac{2K}{\rho c_p\omega }}}$$

(1)

and

$${\delta }_{\nu }=\sqrt{{\displaystyle \frac{2\mu }{\rho \omega }}}$$

(2)

Where:

• $\omega =2\pi f_0$ is the angular frequency of the sound wave, $f_0$ is the design frequency.
• $\mu ,\rho$, K and $c_p$ are the viscosity, density, thermal conductivity, and isobaric specific heat of the gas respectively.

Thus, the performance of the thermoacoustic regenerator depends on many parameters. Olson and Swift [19] proposed the application of Similitude in thermoacoustics without using the acoustic approximation.

The equations which are important to thermoacoustics (continuity, motion, and heat transfer) are rewritten in dimensionless form, verifying that the list of dimensionless variables obtained from similitude is complete.

The dimensionless cooling power obtained from similitude $Q_{cn}=Q_c/P_{m}aA$ and the acoustic power $W_n=W/P_{m}aA$, where $P_m$ is the mean pressure; a is the sound velocity, and A is the cross-sectional area of the stack proposed by [19] are mentioned by Tijani et al. [10] that rewrites $Q_{cn}$ and $W_n$ in dimensionless form by using dimensionless parameters as:

$$\begin{array}{c}{Q}_{cn}=-\left[{\displaystyle \frac{{\delta }_{\kappa n}{D}^{2}sin\left(2x_n\right)}{8\gamma (1+\sigma )\Lambda }}\right]\hfill \\ \times \left[{\displaystyle \frac{\Delta T_{mn}tan\left(x_n\right)}{\left(\gamma -1\right)B{L}_{sn}}}{\displaystyle \frac{1+\sqrt{\sigma }+\sigma }{1+\sqrt{\sigma }}}-\left(1+\sqrt{\sigma }-\sqrt{\sigma }{\delta }_{\kappa n}\right)\right]\end{array}$$

(3)

$$\begin{array}{c}{W}_n=\left[{\displaystyle \frac{{\delta }_{\kappa n}{D}^2{L}_{sn}\left(\gamma -1\right)B{cos}^2\left(x_n\right)}{4\gamma }}\right]\hfill \\ \times \left[{\displaystyle \frac{\Delta T_{mn}tan\left(x_n\right)}{B{L}_{sn}\left(\gamma -1\right)\left(1+\sqrt{\sigma }\right)\Lambda }}-1\right]\\ \hfill -\left[{\displaystyle \frac{{\delta }_{\kappa n}{L}_{sn}{D}^2}{4\gamma }}\times {\displaystyle \frac{\sqrt{\sigma }{sin}^2\left(x_n\right)}{B\Lambda }}\right]\end{array}$$

(4)

Where $\Lambda =\left(1-\sqrt{\sigma }{\delta }_{\kappa n}+{\displaystyle \frac{1}{2}}\sigma {\delta }_{\kappa n}^2\right)$, ${\delta }_{\kappa n}={\delta }_{\kappa }/y_0$ is the dimensionless thermal penetration depth, $B=y_0/(y_0+l)$ is the regenerator blocking ratio, $y_0$ is half the distance between the stack plates and l is half the thickness of the regenerator planes, $\Delta T_{mn}=\Delta T_m/T_m$ is the dimensionless temperature difference, k is the wave number, $x_n$ is the dimensionless stack position, $L_{sn}$ is the dimensionless stack length, $\gamma$ is the ratio of isobaric to isochoric specific heats, ${\delta }_{\kappa \nu }={\delta }_{\kappa }/y_0$ is the dimensionless viscous penetration depth, and D is the drive ratio.

Thus, the performance of the regenerator depends explicit on dimensionless variables, as shown in Equations (3) and (4). It is important to note that the Prandtl number $\sigma$ is proportional to the ratio of the thermal ${\delta }_{\kappa }$ and viscous ${\delta }_{\nu }$ penetration depths, see Equations (1) and (2) respectively. By other hand, the dimensionless viscous penetration depth ${\delta }_{\nu n}$ is implicit in the $COP$ calculation even though, it is not present in Equations (3) and (4). The Prandtl number is an essential parameter for understanding and analyzing the convective heat transfer characteristics of fluids. If gas working is chosen, then the $COP$ depends on five dimensionless parameters. For the analysis in this study, five of these dimensionless variables are selected from Table 1: drive ratio D, the blocking ratio B, the dimensionless temperature difference $\Delta T_{mn}$, the dimensionelss stack position $x_n$, the dimensionless stack length $L_{sn}$. The analysis focuses on these variables through the $COP$ defined as $COP=Q_{cn}/W_n$. Our analysis considers a quarter-wavelength resonator using air as a working gas.

2.2. Design of Experiments (DOE)

The Design of Experiments (DOE) or factorial design, is a regression analysis extensively employed in industry to identify and control critical parameters for enhancing product quality [21]. Experimental studies are normally time-consuming, and requiere examining the effects of two or more factors. Factorial design allows for the investigation of all possible combinations of the levels of specified factors (see Figure 2). Due to the time and expensive of the experiments, the selection of the factors to be analyzed must be made with careful consideration. One of the main reasons to apply a DOE is to minimize the number of experimental trials and to derive a transfer function. Another benefit of using DOE is the seamless integration of experimental data with numerical methods, such as the finite element method (FEM) to simulate alterations in physical parameters (factors).

A complete factorial design incorporates all possible "levels" of a given number of "factors". These levels are categorized as "low" (in red) and "high" (in blue) and are represented by the symbols "-" and "+", respectively (see Figure 2). The effect of a factor is defined as the change in the response cause by altering the level of that factor. This is referred to as the main effect1 because it represents the primary influence in the experiment. Furthermore, the difference in response in some experiments between levels of one factor may vary across the levels of the other factors. When this variation occurs, it indicates an interaction between the factors. Analyzing the main effects and interactions between factors in the results the archivement of one or more of the following objectives:

• Determine the optimal condition within a specified range.
• Assess the contribution of individual factors and their interactions.
• Predict the response under optimal conditions.

The estimation of the contribution of each factors and their interactions will be determined within a specified range or limit.

In this study, the parameter range is set to derive a response function ($\widehat{y}$) which is defined by the parameters or factors ($x^{\prime }s$), coefficients ($b_i$) and the error term ($e_{error}$), see Equations (5)–(7):

$$\widehat{y}=b_0+\sum _{i=0}^k{b}_1{x}_1+\sum _{i=0}^k{b}_2{x}_2+\sum _{i=0}^k{b}_{12}{x}_1{x}_2+\dots +e_{error}$$

(5)

where:

$$Y=Xb+e_{error}$$

(6)

$$\begin{array}{ccccccc}\hfill Y& =\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right],\hfill & \hfill X& =\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1k}\\ a_{21}& a_{22}& \cdots & a_{2k}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nk}\end{array}\right],b\hfill & \hfill =\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right],& e\hfill & \hfill =\left[\begin{array}{c}{e}_{erro{r}_1}\\ e_{erro{r}_2}\\ ⋮\\ e_{erro{r}_n}\end{array}\right]\end{array}$$

(7)

The response of the $COP$ will be represented by $\left[Y\right]$, where $\left[X\right]$ is an $n\times p$ matrix of independent variables (defined as A,B,C,D, and E in Table 2. The vector b contains the regression coefficients and $e_{error}$ is a vector of random errors. The system’s $COP$ will be calculated in the following section using the dimensionless parameters, expressed in Equations (3) and (4).

Additionally, once the parameters are established, the DOE is primarily analyzed in three stages, encompassing all the theorical calculation. These stages are:

• Planning stage
• Analysis stage
• Results stage

The most crucial part of the DOE is the planning stage. It is vital the development of an orthogonal matrix to account for the effects of various factors influencing the $COP$ response. According the literature, the relationship of these factors depends on the physical parameters of the thermoacoustic refrigerator, as detailed in Table 2. These factors will be analyzed quantitatively and qualitatively in the DOE.

Therefore, five parameters, each with two levels (low and high), are put forward as factors to create a full $2^5$ factorial. Table 3 illustrates the matrix to be followed based on the $2^5$ levels by factors of Table 2. The determination of the interaction signs is calculated by multiplying the signs of each interaction shown in Table 3 up to the third interaction for matrix size ratios. For example, the first value in column AB is ”+” derived by multiplying the signs $A-$ and $B-$. For the higher interactions, the remaining signs are calculated similarly. Using these interaction signs and the $COP$ response, the average $COP$ values (low and high) are computed to obtain the coefficients. These coefficients are calculated with MINITAB 19 [22] and used to generate the main interaction effects plots. It is crucial that the low and high levels remain orthogonal to each other, preserving the same proportion between the factors. These values depend on the nominal values chosen for the DOE. For intance, if the nominal value of a factor is 0.10, the low and high limits will be 0.05 and 0.15, respectively.

Several statistical decision approaches constitute the analysis phase: examining the significance of the distribution, performing analysis of variance (ANOVA), and evaluating the R-Sq coefficient (which represents the ratio of the variation explained by each variable). These methods enable the identification of significant coefficients and facilitate the contruction of the transfer function in two stages.

First, the significance of the $COP$ distribution is assessed to establish a confidence level of 95%, which corresponds to a significance level of $\alpha$=0.05 for a single proportion of population. The Alpha ($\alpha$) represents the maximum acceptable probability of incorrectly rejecting the null hypothesis when the alternative hypothesis is true. The $COP$ results are evaluated using average analysis2$\mu ={\sum }_{i=0}^N{\displaystyle \frac{Y_i}{N}}$, the standard deviation3$s={\sum }_{i=1}^N{\displaystyle \frac{{(Y_i-\mu )}^2}{N}}$ the sum of squares deviations$\left(SST\right)=\sum s^2$ and the variance $s^2={\displaystyle \frac{SST}{DF}}$ of the factors affecting the $COP$, where N represents the number of runs and "DF" stands for degrees of freedom$=n-1$.

Second, the ANOVA test evaluates the averages that affect the $COP$ to identify the variance among factors and their interactions. A P-value ≤ 0.05 indicates the significance of the main parameters and their interactions in the ANOVA analysis. Coefficients and errors are determined through on statistical decision approximations using the Anderson-Darling Test method with MINITAB 19, that assess the significance and normality of the data distribution [22]. As shown in Table 3, the ANOVA is conducted using the $COP$ responses derived from the orthogonal matrix,. Numerous visualizations, such as main effects plots, interaction effects plots, Pareto charts, cube plots, contour plots, and surface plots, are utilized to analyze the factors and interactions affecting the $COP$.

Finally, the results phase entails constructing the transfer function using the coefficients derived from the R-Sq value. The relationship between the variation of each variable and the total variation (variation of the Mean Sum of Squares (MSS)) is quantified with the coefficient of determination R-Sq. The R-Sq is also referred to as a measure of the goodness of fit. A highher R-Sq value, approaching 100 %, suggests a better representation of the total variation of the $COP$ response.

3. Results

The parameters intervals for the $COP$’s are defined, as shown in Table 4. The $COP$’s values for each trial are computed using a numerical rutine in Python 3.6 [23] based on the Equations (3 and 4 for each established interval, see Table 5. A representative plot of the $COP$ is observed in Figure 3 as a function of $x_n$ and $L_{sn}$ defined in Table 1. Typically, the $COP$ increases as the dimensionless stack position $x_n$ is decreased. Based on the maximum $COP$’s values obtained numerically, the DOE analysis is carried in three steps:

• Establish the $2^5$ orthogonal design to determine the $COP$ of a thermoacoustic cooling system, as outlined by Tijani [10].
• Perform Pearson analysis: Assess main and interactions effects, Pareto chart, cube plots, contours, and surfaces plots.
• Conduct ANOVA analysis and develop transfer function.

3.1. COP’s Orthogonal Array

The full factorial design, is created by generating a set of combinations based on the number of parameters and their specified values.

Table 4 provides the values for the lower and upper limits. Table 5 presents the matrix of $2^5$ = 32 runs plus the center point, along with the $COP$ values obtained for each combination. The order of runs is randomized to prevent aliasing4. The $COP$ response is calculated for each combination of "low" and "high" values of the parameters.

The DOE is conducted using each combination to determine the $COP$ response. The subsequent step involves analyzing the main effects and interactions of the parameters affecting the $COP$. Note that Table 5 includes 33 runs due to the inclusion of a central point to assess curvature. This central point marked in red, is shown in Figures, Figure 4 and Figure 5. If the central point deviates significantly from the line that joins the means of the vertex points, it indicates a potential curvature in the relationship, and its P-value will be assessed for statistical significance. If the contribution of curvature is found to be insignificant, it implies that the variation in $COP$ response with the chosen factors does not exhibit substantial changes.

3.2. Pearson Analysis

A main effects analysis can be conducted to identify the critical parameters and interactions using a Pearson correlation5, back when the $COP$ values are obtained based on the combinations in Table 5. Pearson correlation assesses the linear relationship between defined variables for a given significance level (in this case $\alpha$=0.05).

Moreover, interactions between parameters can also be assessed. It is crucial to understand that the range of values was established by multiplying each of the parameters at the central point of Table 5 (run 33) by a value of 0.5, and then adding and subtracting this result from the center point to define the upper and lower limits, respectively.

Multiplying the factors by this value to adjust for the lower and upper limits is referred to as orthogonality in the full factorial design, which helps to maintain the "balance" of the DOE. The effects of each chosen parameter on the $COP$ are depicted in Figure 4. Main effects are determined by the slope of the $COP$ values obtained between the defined limits for each parameter: a steeper slope indicates a stronger effect of the parameter. According to the plot results, the dimension parameters with the most significant impact on the $COP$ are: the temperature difference $\Delta T_m$, the length of the regenerator $L_s$, the Blocking ratio B as well as it’s position $X_s$. Counterintuitively, the parameter D does not have a notable effect on the $COP$. However a comprehensive quotient analysis of Equations (3) and (4) reveals that this outcome is anticipated. As shown in the main effect results for $COP$ in Figure 4 , parameter D does not affect the $COP$.

Additionally, Figure 5 illustrates the interactions between parameters. These interactions are identified based on the slope observed between the specified limits. The columns in the interaction graph represent the response for a specific interaction of a given factor (- or +), while the rows represent the average $COP$ response for a particular configuration (- or +) of factor interactions (columns), as shown in Figure 5. Parallel lines in the interactions suggest no interaction effect, whereas a larger difference in slope between the lines indicates a stronger interaction effect. Consequently, the ”low” values within the defined range show a positive effect on $COP$ results compared to the ”high” values. The interactions with the most significant influence, depicted in Figure 5, are between the following parametes pairs in the lower section: $X_s$*B, $\Delta T_m$*$L_s$, $L_s$*B. These interactions, especially those in the lower section of the graph, positively impact the average $COP$ value due to their interaction slopes. Only the pair $X_s$*$L_s$ in the upper section also shows a positive effect on improving the average $COP$ value.

The positive effect on both the upper and lower limits of the interaction between $X_s$*B factors is emphasized. However, the interaction graph alone does not reveal the statistical interaction significance of these interactions. To further understand the impact of main and interaction effects, a Pareto6 chart can be used. A Pareto chart helps identify which factors have the most significant influence on the $COP$ response. It shows the absolute values of the standardized effects from the largest to the smallest. Standardized effects are t-statistics used to test the null hypothesis that the effect is zero [22]. The reference or threshold line for statistical significance depends on the significance level (denoted by $\alpha$). Unless a stepwise selection method is used to define an alpha value, the significance level is calculated as 1 minus the confidence level for the analysis. The reference line indicates which effects are statistically significant. Factors such as $\Delta T_m$, $L_s$, the interaction $X_s$*B and the others with bars extending above the threshold red dashed line (value of 12.7), see Figure 6, are typically the main contributors to be focused on the subsequent Analysis of variance (ANOVA) subsection. With the current model terms, these factors are statistically significant at the 0.05 level.

Additionally, the cube plot, contour plots and surface plots will offer a visual representation of the trends of each interaction based on the established intervals. A cube plot (see Figure 7) illustrates the relationships among the five factors and their impact on the $COP$ response measure for 2-level factorial designs or Plackett-Burman designs. This plot visually represents the factorial design, showing how different factors affect the $COP$ values. By examining the cube plot, one can observe the most significant impacts on the $COP$ values at the lower limits of the factors.

Contour plots and surfaces are used to visualize and explore the potential interaction among all parameters relative to the $COP$ response. Contour plots provide a two-dimensional view of the lower and upper limits, with the x and y axes representing these limits, while the contours indicates the $COP$ response. Similar to a topographic map, where longitude, latitude, and elevation are plotted and connected by contour lines, the contour plot displays x, y, and z values connected by lines of constant $COP$ response. In the contour plots, a maximum $COP$ response of 1.76 can be observed in the products for the lower limits, such as $B*D$, $B*X_s$, $B*\Delta T_m$, and $L_s*B$, as shown in Figure 8 and Figure 9. Notably, at the center point in Figure 8, the maximum $COP$ value is 1.4, specifically for the product $\Delta T_m*D$.

The $COP$ surface plots, illustrated in Figs. Figure 10 and Figure 11, display the variable response surfaces for the center and lower limit hold values respectively. Surface plots are particularly useful for examining the trends or behavior of $COP$ concerning a specific factor between center and lower limit hold values. A notable maximum ridge pattern response of 1.76 in the $COP$ is observed for the products $\Delta T_m*B$ in Figure 11 at the lower limits hold values. The highest $COP$ results for the lower limit values are seen in Run numbers 1 and 2 in Table 5. At the center point hold values depicted in Figure. Figure 10, the peak $COP$ value is 1.4 for the product D*$\Delta T_m$. Surface plots, along with contour paths, provide a comprehensive visualization of the results. Adjusting the hold values will lead to changes in the $COP$ response surface.

3.3. COP’s ANOVA Analysis

The ANOVA is employed to assess the significance of the design parameters affecting the $COP$. This analysis is conducted with a 95% confidence level, corresponding to $\alpha$=0.05. The results obtained in Table 6, demostrate that dividing the sum of squares by main effects or interactions can reveal their impact on the $COP$ response. The main effects contribute significantly (55.69%) to the $COP$’s response, followed by the 2-Way Interactions (33.30%), 3-Way Interactions (7.36%), and the 4-Way Interactions (3.35%). A parameter or interaction is deemed statistically significant if the P-value is less than 0.05. According the ANOVA results in Table 6, the parameters, as well as the second, third and fourth interactions in the model, have a significant impact on the $COP$ response.

The influence of these interactions can be visualize using the coefficient of determination R-Sq, which represents the ratio of the variation explained by the model to the total variation. An R-Sq value close to one indicates a better-fitting model

[24]. For this model, the R-Sq value is 100.00%, accounting for up to 4-Way Interactions, reflecting a high degree of variability in the data based on the factors and their interactions as demostrated by the $COP$ values obtained in the previous subsection.

3.3.1. Transfer Function of the $COP$’s response

After performing the ANOVA, the transfer function can be determined using the estimated uncoded coefficients derived from MINITAB for $COP$ responses. The uncoded coefficients shown in Table 6 are based on the ANOVA analysis including up to 4-Way Interactions. It is important to emphasize the importance of including the appropriate number of interactions to minimize error in the transfer function. The estimated uncoded coefficients, which consider up to 4-interactions, help in reducing the error. As indicated in Table 6, the R-Sq, which measure the goodness of fit, decreases to 55.69% when only the main effects are considered. However, it improves to 88.99%, when 2-Way Interactions are included. The R-Sq value is calculated by dividing the Seq-SS of the main effects by the total sum of squares.

The Regression Equation in Uncoded Units with all terms is observed in Eq.8:

$$\begin{array}{c}COP=5.5922+0.00656D-0.41781x_n-0.33544\Delta T_{mn}-5.5487B\hfill \\ -0.25144Ls-0.000438D·x_n-0.000563D·\Delta T_{mn}-0.01625D·B\\ -0.000813D·L_{sn}+0.024912x_n·\Delta T_{mn}+0.69325X·B+0.021062x_n·L_{sn}\\ +0.44975\Delta T_{mn}·B+0.018988\Delta T_{mn}·L_{sn}+0.22775B·L_{sn}\\ +0.000038D·x_n·\Delta T_{mn}+0.000750D·x_n·B+0.000038D·x_n·L_{sn}\\ +0.001250D·\Delta T_{mn}·B+0.000062D·\Delta T_{mn}·L_{sn}+0.002250D·B·L_{sn}\\ -0.047450x_n·\Delta T_{mn}·B-0.001553x_n·\Delta T_{mn}·L_{sn}-0.032450x_n·B·L_{sn}\\ -0.024550\Delta T_{mn}·B·L_{sn}-0.000050D·x_n·\Delta T_{mn}·B\\ -0.000002D·x_n·\Delta T_{mn}·L_{sn}-0.000050D·x_n·B·L_{sn}\\ \hfill -0.000150D·\Delta T_{mn}·B·L_{sn}+0.002790x_n·\Delta T_{mn}·B·L_{sn}+0.14844Ct·Pt\end{array}$$

(8)

The Regression Equation in Uncoded Units with P-values ≤ 0.05 is defined in Eq.9:

$$\begin{array}{c}COP=5.5922-0.41781x_n-0.33544\Delta T_{mn}-5.5487B\hfill \\ -0.25144L_{sn}+0.024912x_n·\Delta T_{mn}+0.69325x_n·B+0.021062x_n·L_{sn}\\ +0.44975\Delta T_{mn}·B+0.018988\Delta T_{mn}·L_{sn}+0.22775B·L_{sn}\\ -0.047450x_n·\Delta T_{mn}·B-0.001553x_n·\Delta T_{mn}·L_{sn}-0.032450x_n·B·L_{sn}\\ \hfill -0.024550\Delta T_{mn}·B·L_{sn}+0.002790x_n·\Delta T_{mn}·B·L_{sn}=1.759313\end{array}$$

(9)

The transfer function in Eq. 8 includes all coefficients, while Eq. 9 features only those coefficients with a P-Value ≤ 0.05. Both functions are constructed using these uncoded units and highlight the most significant coefficients affecting the $COP$’s response. Specifically, Eq. 9 contains only the significant coefficients, including the 4-Way interactions identified through the ANOVA analysis of the $COP$’s results. Given the ANOVA results, the R-sq value from a response surface methodology (RSM) was not included, as it indicated a minimal contribution from curvature in the Full Factorial (less than 1%), see Table 6. Consequently, the optimal $COP$ value for the analyzed space can be calculated and validated using the "Solver" function in ’Excel’ or ’Libreoffice calc’, (see Figure 12) based on the transfer function derived from the $COP$’s parameters $x_n$, $\Delta T_{mn}$, B, $L_{sn}$ and their interaction with the uncoded coefficients. The residual error related in the ANOVA analysis for this case is zero (0.00%). Additionally, the transfer function indicates that the parameter, D and its interactions are not significant, based on the P-values presented in Table 6.

4. Conclusions

In this study, the $COP$ of a thermoacoustic device, dependent of five dimensionless parameters was analyzed using the equations from linear thermoacoustic theory. By applying DOE methodology and tools, the contribution of each parameter and their interactions at various levels were quantified and qualified. Although the case study is illustrative, the methodology can be extended to analyze other $COP$ parameter spaces under different operating conditions. A $COP$ value of 1.76 was archived for $X_s$=0.05m, $\Delta T_m$= 5 °C, B=0.25, $L_s$=0.05m and D=1 in this study. The Full factorial design allowed for a systematic variation of all factors at all levels, providing comprehensive quantitative and qualitative insights into the main effects and their interactions. This approach facilitates easy interpretation due to the clear separation of effects. The graphs generated in this DOE study, including Pareto, cube plot, contour, and surface plots are particularly valuable for visualizing the impact and influence of the $COP$s factors and their interactions. Contours and surface plots, in particular, are effective for optimizing responses involving multiple factors. They offer a flexible approach to model fitting, exploring the response surface, and identifying optimal process conditions, and understanding the curvature of the response surface, which in this case, contributes 0.29%.

The $COP$ results derived with the Pearson correlation and the ANOVA analysis are utilized to assess the sensibility of the stack parameters selected for a thermoacoustic refrigerator. Evaluating the sensibility of the $COP$ parameters suggests that opting for the lower values rather than the upper ones can increse the $COP$ to 1.76. Finally, with the transfer function allows for archieving the best possible $COP$ value within the limits established by the DOE using the “Solver” option in Excel to determine the optimal final values of the factors. In this case, the best $COP$ value achievable with the transfer function align with the $COP$ value obtained in run 1, as shown in Table 5. Finally, the performance relative to Carnot COPR vs the COP values obtained in Table 6 will be analyzed in future work to find out the parameter differences applying the proposed DOE methodology.