Applied Acoustics 67 (2006) 689–699

[Link]/locate/apacoust

Temperature gradient integration in

thermoacoustic stacks
a,*
Carl Jensen , Richard Raspet a, William Slaton b

a
National Center for Physical Acoustics, University of Mississippi, University, MS 38677, United States
b
Department of Physics and Astronomy, University of Central Arkansas, Conway, AR 72035, United States

Received 29 April 2005; received in revised form 8 November 2005; accepted 9 November 2005
Available online 27 December 2005

Abstract

Different approaches to numerical modeling of thermoacoustic effects are being investigated as a

step towards developing a model for inert gas–vapor mixtures in thermoacoustics. This paper
describes these approaches to numerical calculations, which differ primarily in the way that the tem-
perature profile is computed for the stack. Different aspects of the constant property of the thermo-
dynamic enthalpy in the stack are used by each approach: one maintains a constant value for the
enthalpy while calculating the temperature and the other keeps the spatial derivative of the enthalpy
equal to zero. The runtime and accuracy of the two approaches are compared by applying them to
both a thermoacoustic refrigerator and engine model, and it is found that the second methodÕs intro-
duction of an additional integration variable leads to more error than the first methodÕs use of large
numbers in the calculation.
Ó 2005 Elsevier Ltd. All rights reserved.

Keyword: Thermoacoustics

1. Introduction

Our research group is currently studying two different aspects of thermoacoustics which
includes inert gas–condensing vapor [1,2] thermoacoustics and devices of small scale with
high operating frequencies. A significant focus of this research is in developing numerical

*
Corresponding author. Tel.: +1 662 915 5843; fax: +1 662 915 7494.
E-mail address: crjensen@[Link] (C. Jensen).

0003-682X/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2005.11.004
690 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699

models from theory that can be used to accurately simulate devices of both types. The tem-
perature profile in the stack is of particular importance since the transport properties of an
inert gas–condensing vapor mixture vary rapidly with temperature fluctuations and
because the short stack in high-frequency devices leads to a high temperature gradient
where any errors are amplified by the small dimensions of the device.
Computer code that successfully modeled large, dry thermoacoustic devices developed
convergence problems when used to model small (high-frequency) dry thermoacoustic
devices. This difficulty along with the need for accurate temperature calculations moti-
vated an investigation into the sources of error in the method used to calculate the tem-
perature in the stack. The key concern with the original method was the mix of large
and small value acoustic variables that are added and subtracted to calculate the temper-
ature gradient. This raised concerns that taking the difference of large floating point num-
bers might generate errors. Taking the derivative of this expression and solving for the
derivative of the temperature gradient eliminates the large constants and reduces the large
values to smaller derivatives so that all of the terms are on the same scale. This produces a
new approach that is mathematically equivalent to the previous approach while eliminat-
ing the large numbers from the calculation.
This paper compares the two different approaches to temperature calculations in ther-
moacoustics. The first is the original technique used by DeltaE [3] and Arnott [4] which
employs the fact that the enthalpy is conserved in the steady state of an insulated stack
to evaluate the temperature gradient while using complex pressure and complex velocity
as integration variables. The second makes use of the constant enthalpy by setting the spa-
tial derivative of the enthalpy equal to zero in order to evaluate the second derivative of
the temperature in the integration. In this paper, we describe the two approaches, apply
them to inert gas thermoacoustics, and compare the results of integration.

2. Calculations

The thermoacoustic devices in this paper were modeled by solving the thermoacoustic
equations for a one-dimensional sound field along its long axis. The long and narrow
geometry of the designs ensures that the one-dimensional sound field will provide a suffi-
cient solution. The one-dimensional acoustic field presents a problem in which differential
equations for the variables that characterize the field must be solved with boundary con-
ditions that describe steady state operation. Computationally, the differential equations
are solved using a Runge–Kutta integration routine while the boundary conditions are
met by a Newton–Raphson zero search.
The Runge–Kutta [5] routine is a method of numerically integrating ordinary differen-
tial equations. The routine is used to evaluate the differential equations for the complex
acoustic pressure and velocity amplitudes that characterize the acoustic field being mod-
eled – these differential equations are discussed later in this section. Beginning on the left
end of the device (considering the device to be two ends of pipe separated by a stack as in
Fig. 1), an initial guess is made for the pressure amplitude as well as the temperature and
frequency. The acoustic impedance at the beginning (or end) of the device may be calcu-
lated from these guessed values because it is either a normal reflection at a solid surface [6]
or a transducer for the devices in this paper. This impedance is used to calculate the initial
velocity amplitude and a Runge–Kutta routine is used to integrate both the pressure and
velocity across the device to the beginning of the stack where it meets the heat exchanger.
 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699 691

T hot T cold
pˆ in pˆ out
Ẑ 1 vˆin vˆ out Ẑ 2
H H

Fig. 1. Simple diagram of a thermoacoustic resonator with hot and cold heat exchangers, a stack, and impedance
^ at both ends.
conditions, Z,

At the face of the stack, the heat flow from the heat exchanger and the intensity of the
sound wave are used to calculate the enthalpy entering the stack. The enthalpy is used to
calculate the temperature gradient and integrate the temperature as well as the pressure
and velocity through the stack as discussed below. The sound intensity exiting the stack
and the enthalpy are then used to calculate the heat flowing to the heat exchanger on
the right side. Integration of the pressure and velocity continues until reaching the end
of the device where the ratio of the two is used to calculate the acoustic impedance.
The routine above illustrates how initial guesses at the values of the frequency and the
acoustic pressure amplitude and temperature on the left side of the device are used to cal-
culate a variety of resulting properties including the heat flow, acoustic impedance, and
temperature on the right side of the device. These guessed values are matched to those
of the device being modeled by using a Newton–Raphson search to meet the boundary
conditions. A Newton–Raphson search [7] is a root-finding algorithm that finds the zeros
of a function from an initial guess near the root. In this case, an initial guess is made and
values of the temperature, pressure and frequency that meet conditions on the heat flow
and acoustic impedance on the right side are found. As mentioned before, the acoustic
impedance is known for the devices in this paper as well as the flow of heat, so these quan-
tities are known and can be targeted in the search.
The modeling procedure discussed above applies directly to the engine model presented
in the results section. The refrigerator model is very similar except that the values guessed
are the pressure amplitude, heat rejected on the hot side, and driver phase and the targets
are the complex end impedance and the heat pumped from the cold side.
The thermoacoustic equations are used to find the differential equation for each vari-
able used in the Runge–Kutta integration. First, the Navier–Stokes fluid field equations
provide differential equations for the area-averaged complex axial particle velocity, ^v,
and the complex acoustic pressure amplitude, ^ p. Inside the stack, the acoustic pressure
is constant across the cross-section and is a function only of position along the stack.
These considerations provide the following equations from Ref. [8] (a more thorough der-
ivation is presented in Ref. [4]):
d^p ixqðzÞ^vðzÞ
¼ ; ð1Þ
dz F ðkÞ
 
d^v 1 dT F ðkT Þ=F ðkÞ  1 F ðkÞ 2
¼ ^v  k ^p; ð2Þ
dz T dz 1r ixqðzÞ
where
x2 ½c  ðc  1ÞF ðkT Þ
k2 ¼ ; ð3Þ
c2 F ðkÞ
and q(z) is the density, x is the angular frequency, c is the ratio of specific heats, and r is
the Prandtl number. The thermoviscous function for parallel plates is given by
692 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699

 pffiffiffiffiffiffi
2 k i
F ðkÞ ¼ 1  pffiffiffiffiffiffi tanh ; ð4Þ
k i 2

where k = R(qx/g)1/2 is the shear wave number and kT = R(qxcP/j)1/2 is the thermal wave
number for a gas of density q, viscosity g, thermal conductivity j, and specific heat per
unit mass at constant pressure cP. R is the characteristic pore radius defined by twice
the transverse pore area divided by the pore perimeter and is equal to the spacing between
the plates for a parallel plate stack. The thermoviscous functions for other stack geome-
tries are available in Ref. [4].
The temperature profile in a thermoacoustic stack may be significantly non-linear for a
device running with high power output; [9] as a result, the temperature must also be used
as an integration variable in order to build a robust simulation. As before, a differential
equation for the temperature within the stack is needed and can be found by considering
the flow of energy down an insulated stack in the steady state. In such a state, the flow of
energy into one section of the stack must be the same as the flow of energy out, thus the
flow of energy must be a constant at every point in the stack. This acoustic energy flow is
given by the thermodynamic enthalpy, and thus the enthalpy must be constant throughout
the stack. Examining the equation for the enthalpy from Ref. [8] and combining the terms
for heat, work, and thermal conduction loss yields
     2
Agas 1  F ðkT Þ Agas dT cP  ^v 
H¼ Re ^v^p þr 
2 1þr F ðkÞ 2 dz x F ðkÞ
q dT
 ImðF  ðkT Þ þ rF ðkÞÞ  ðk gas Agas þ k stack Astack Þ . ð5Þ
1  r2 dz
We find that we can solve this equation for dT/dz, with Agas and Astack being the cross-
sectional area in the stack of the gas and the solid respectively, kgas and kstack being their
thermal conductivity, and * denoting complex conjugation. The resulting equation is
then a function of the acoustic pressure and velocity amplitudes, temperature, and
frequency.
The alternative approach takes the derivative of the enthalpy with respect to the long
axis, z, and solves for d2T/dz2 in order to eliminate the large, constant value of the
enthalpy from the calculations. Since kgas, kstack, q, k, and kT all vary with temperature
they also vary with position and must be considered variables in the differentiation as well
p, ^v, and T. This process yields an equation for d2T/dz2 (see appendix) that is a function
as ^
of the pressure, velocity and temperature as well as their derivatives with respect to z. Illus-
trating the two approaches symbolically:
Method 1:
 
dT dT dT
H ¼ H p ; v; T ;
^ ^ ¼ const. ! ¼ ð^
p; ^v; T ; x; H Þ.
dz dz dz
Method 2: Introduce a new variable G:
   
dH dH d^
p d^v dT d2 T dG d2 T d^p d^v dT
¼ p; ; ^v; ; T ;
^ ; 2 ¼0! ¼ 2 ^p; ; ^v; ; T ; ;x
dz dz dz dz dz dz dz dz dz dz dz
dT
¼ G.
dz
 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699 693

The Runge–Kutta routine operates in steps by using the derivative to trace the value of
the variable across the dimensions of the integration. A second-order routine evaluates the
derivative at the beginning of each step and then reevaluates it at the middle to increase
accuracy. So, for the first method, the derivative of the velocity, pressure and temperature
is evaluated at the beginning and then used to extrapolate the variables to the middle of
the step. In the middle, the derivatives are evaluated and the variables extrapolated again
to produce their value at the end of the step. This process is repeated step by step until the
end of the integration is reached. For the second method, the new variable G is integrated
using d2T/dz2 as dG/dz; then G is used as dT/dz in determining T. In order to integrate G,
an initial value entering the stack is calculated using the first method where dT/dz is found
using the enthalpy, pressure, velocity and temperature at the beginning of the stack. Now
G is an integration variable and is simply the derivative of the temperature, so the temper-
ature gradient has become an integration variable in the new method.
In this second approach, the enthalpy itself becomes an integration variable indirectly
through the second-order equation in the temperature. This approach is being examined
because eliminating the value of the enthalpy from the derivative calculation may alleviate
problems in taking the difference of several large numbers. The more complicated expres-
sion for the second-order temperature derivative contains many small terms that are each
scaled more closely to one another than the terms found in the first-order derivative, and
this may reduce the errors in the numerical model. This approach also has the advantage
that it provides a convenient method for estimating accuracy by checking the behavior of
the enthalpy in the stack – if the enthalpy is nearly constant then the integration is accu-
rate. However, these advantages may not outweigh the greater cost of computation.

3. Results

Two different example systems were modeled in order to examine any differences in
behavior between the two methods described above. The first is a solar-heated thermo-
acoustic engine with a STAR linear alternator designed by the Clever Fellows Inc. (CFIC)
[10]. The engine is comprised of a hemispherically shaped quartz heater head that supports
the high operating pressure of the device while focusing sunlight onto the hot end of the
stack. The stack is then cooled on the other side by a gas–water heat exchanger, and a little
further down from the cooler is the STAR alternator used for power conversion. The sec-
ond model used was a refrigerator designed as part of the NavyÕs Shipboard Electronic
Thermoacoustic Cooler project (or SETAC) [11]. The refrigerator uses a U-shaped half-
wavelength resonator with a stack and loudspeaker placed at each end as well as four iden-
tical gas–fluid heat exchangers. The relevant dimensions used for both devices are listed in
Table 1.
Fig. 2 presents the calculated temperature profiles in the stack of both devices during
operation. Fig. 2(a) is the temperature versus position in the stack for the engine model,
and Fig. 2(b) is the same plot for the refrigerator model. The position in both plots is nor-
malized to the stack length. These plots demonstrate that the temperature behavior in the
stack is significantly non-linear for a real operating device.
The plots in Figs. 3 and 4 demonstrate the behavior of the error in the stack integration
versus the percentage of the stack length covered per step. For these figures, the integra-
tion error is determined from the Newton–Raphson root-finding algorithm mentioned in
the calculations section which is used to find values of the temperature, pressure and
694 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699

Table 1
Dimensions for thermoacoustic engine and refrigerator models
CFIC engine SETAC refrigerator
Geometry
Device length (cm) 24.89 39.06
Device radius (cm) 4.37 2.23
Stack length (cm) 8.29 4.03
Stack material Stainless Steel Mylar
Stack porosity 0.81 0.80
Stack plate spacing (mm) 0.3 0.28
Stack plate thickness (mm) 0.025 0.052
Hot heat exchanger
Length (cm) 7.5 0.25
Porosity 0.20 0.72
Cold heat exchanger
Length (cm) 3.5 0.25
Porosity 0.35 0.72
Operating conditions
Mean pressure (atm) 40 20.7
Temperature
Hot side (K) 752.74 300
Cold side (K) 342.08 286.26
Frequency (Hz) 178.2 320

frequency (in the case of the engine model) that match boundary conditions on the acous-
tic impedance and heat flow. An accurate solution for the temperature, pressure, and fre-
quency is found by using a large number of steps in the integration and tight tolerances in
the root-finding routine. These values are then used in the integration with a varying num-
ber of step sizes and the deviation of the heat flow and acoustic impedance from the
boundary conditions of the system are recorded as a straightforward measure of the error
generated by the entire calculation. Lastly, in these figures, the error is recorded as the root
squared sum of the relative error in the acoustic impedance and the heat flow.
Figs. 3 and 4 quickly reveal both that the engine model generates more error than the
refrigerator and that method 2 is at a significant disadvantage in terms of accuracy as its
error grows much faster than that of method 1. The larger error in the engine model cal-
culation is caused by the more significantly non-linear temperature profile of the device. As
can be seen in Figs. 2a and b, both devices exhibit a great deal of curvature but that of the
engine model is greater in magnitude due to the larger temperature gradient and much
greater heat flow. On the other hand, the most likely cause of the large errors shown by
method 2 is the additional integration variable used by the new method where the temper-
ature gradient is integrated as well as the temperature, pressure and velocity. The added
error in the fourth integration variable overcomes any possible improvements in other
sources of numerical error and puts method 2 decidedly behind method 1.
However, the new approach has an advantage in being able to use the enthalpy to esti-
mate the total error in the integration. The Runge–Kutta routine that was used has an
inherent error thatÕs second order in the step size, as can be seen in Figs. 3 and 4 where
the error varies quadratically from the accepted value as a function of the step size. The
error in the stack integration variables also behave this way and, with the second method,
the enthalpy does as well; though, the error in the enthalpy is determined by its deviation
 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699 695

800

Temperature (K)
700

600

500

400

300
0.00 0.20 0.40 0.60 0.80 1.00
(a) Position (z/L)

305

300
Tempe ra ture (K)

295

290

285

280
0.00 0.20 0.40 0.60 0.80 1.00
(b) Position (z/L)

Fig. 2. Temperature profile in the stack for the CFIC solar-heated engine (a) and the SETAC refrigerator (b).
The position in the stack is recorded as the fraction of the stack length, L. The plots demonstrate that the
temperature profile in these devices is significantly non-linear.

0.080

0.070
2
0.060
Total error (% )

y = 0 .03x - 0.00x + 0.00

2
0.050 R = 1.00

0.040
0.030

0.020
2
0.010 y = 0.000 4x - 0.0000x + 0.0000

0.000
0.2 0 0.4 0 0.60 0.80 1 .00 1 .20 1.40 1.60
Step size (%)

Fig. 3. The total error for the CFIC thermoacoustic engine is calculated as the magnitude of the Newton–
Raphson target vector with an accepted solution for the guess vector found using small integration steps and tight
tolerances in the root-finding routine. The step size is recorded as the percentage of the stack length and the two
are plotted against each other demonstrating the second-order behavior of the Runge–Kutta integration.
s – fixed enthalpy approach – method 1; h – zero derivative approach – method 2.
696 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699

8.00

7.00 2
y = 5.13x - 0.09 x + 0 .01
6.00 2

Total error (%)

R = 1.0 0
5.00

4.00

3.00

2.00
2
1.00 y = 0 .06x + 0.00 x - 0.00

0.00
0.20 0.4 0 0.6 0 0.80 1.00 1 .20
Step size (%)

Fig. 4. The total error for the SETAC thermoacoustic refrigerator is calculated as the magnitude of the Newton–
Raphson target vector with an accepted solution for the guess vector found using small integration steps and tight
tolerances in the root-finding routine. The step size is recorded as the percentage of the stack length and the two
are plotted against each other demonstrating the second-order behavior of the Runge–Kutta integration.
s – fixed enthalpy approach – method 1; h – zero derivative approach – method 2.

from constant and itÕs the standard deviation that decreases with the step size. This makes
the error in the enthalpy integration simple to determine and it yields a direct estimate for
the accuracy of all the integration variables if their errors are correlated.
In Fig. 5, the relationship in the engine model between the error in the enthalpy and the
total error is examined by plotting them against one another. The linear correspondence
shown implies that the error in the stack integration is directly related to the error in
the enthalpy. Thus the error in the stack variables as well as the total error can be deter-
mined by analyzing the behavior of the enthalpy through the stack without an actual solu-
tion being known.
Lastly, in order to compare the computation cost of the two methods, they each were
run several times and the average runtime is recorded in Table 2 at different step sizes. The

0.080
0.070
y = 0 .85x + 0.0 0
2
0.060 R = 1.00
Total error (%)

0.050
0.040
0.030
0.020
0.010
0.000
0.0 00 0.02 0 0 .040 0.060 0.08 0 0.100
Enthalpy error (%)

Fig. 5. A graph illustrating the correspondence in the engine model between error in the enthalpy, which is
determined as the standard deviation of the enthalpy calculated at each point in the stack, and the total error
determined against the accepted solution. The linear correspondence implies that the different integration errors
are directly related.
 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699 697

Table 2
Runtime comparison of both calculation methods applied to both the refrigerator and engine models
Number of steps Engine Refrigerator
Method 1 (s) Method 2 (s) Method 1 (s) Method 2 (s)
300 1.02 1.11 1.05 1.10
600 1.96 2.20 2.02 2.19
900 2.95 3.31 3.04 3.27

differences in runtime shown in Table 1 amount to the first method having roughly a ten
percent performance advantage.

4. Conclusion

The new approach to thermoacoustic temperature calculations described in this paper

was found to have no advantage in accuracy for the situations examined. The two cases
compared demonstrated that while the two methods converged to comparable solutions
for a reasonable step size, the second method generates more error for a given step size
and has about a ten percent performance penalty. So, in the end, the new method is a via-
ble solution for calculating the temperature profile with the added benefit of offering a sim-
ple estimate of numerical error, but with a notable impact on numerical error as well. For
these reasons, both methods will be retained as research progresses but the constant
enthalpy method will be used predominantly.

Acknowledgements

The authors thank Steve Garrett of Pennsylvania State University and John Corey of
Clever Fellows Innovation Consortium, Inc. for providing guidance and dimensions for
modeling the thermoacoustic refrigerator and engine. In addition, the authors thank Hank
Bass and Elliot Hutchcraft of the University of Mississippi for helpful suggestions and
discussions.

Appendix A

In order to find the differential equations for the temperature used in the two calcula-
tion methods discussed, we start first with the enthalpy given in Eq. 5:
     2
Agas 1  F ðkT Þ Agas dT cP  ^v 
H¼ Re ^v^p þr 
2 1þr F ðkÞ 2 dz x F ðkÞ
q dT
 ImðF  ðkT Þ þ rF ðkÞÞ  ðk gas Agas þ k stack Astack Þ . ðA1Þ
1  r2 dz
If we solve the enthalpy for dT/dz we find the differential equation for the temperature
used in method 1:
  
Agas 1
dT H  2 1þr
Re p FF ðk
^v^ TÞ
ðkÞ
þr
¼   2 . ðA2Þ
dz A cP  ^v  q 
 gas2 x  F ðkÞ  1r2 ImðF ðk T Þ þ rF ðkÞÞ  ðk gas A gas þ k stack A stack Þ
698 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699

For method two, the equations for d2T/dz2 are simplified if we rename the repeated term
in the enthalpy:
1
D¼ ðrF ðkÞ þ F  ðkT ÞÞ; ðA3Þ
1þr
which produces
    2
Agas  D Agas dT cP  ^v  q
H¼ Re ^v^ p    ImðDÞ
2 F ðkÞ 2 dz x F ðkÞ 1r
dT
 ðk gas Agas þ k stack Astack Þ . ðA4Þ
dz
To find the derivative of the enthalpy we recall that ^v; ^p, and T vary directly with position
whereas kgas, kstack, r, k, and kT vary with temperature. Also, we assume a simple power
law for the sheer and thermal wave numbers and the Prandtl number:

kT / T a ; ðA5aÞ
b
k/T ; ðA5bÞ
d
r/T . ðA5cÞ

Taking the derivative of H:

(    !
dH Agas d^v  p
d^ D  0 dT D
¼ Re p þ ^v
^  Re ^v^p F ðkÞbk
dz 2 dz dz F ðkÞ dz T ðF ðkÞÞ2
  0  )  2
D Agas dT cP  ^v  q
þRe ^v^ p 
F ðkÞ 2 dz x F ðkÞ 1  r
   
2 d^v k dr 1 dT
 Re  2F 0 ðkÞb þ 1 ImðDÞ þ ImðD0 Þ
^v dz F ðkÞ ð1  rÞ T dz
 2  2  2
dk gas dT dk stack dT Agas d2 T cP  ^v 
 Agas  Astack 
dT dz dT dz 2 dz2 x F ðkÞ
q d2 T d2 T
 ImðDÞ  k gas Agas 2  k stack Astack 2 ; ðA6Þ
1r dz dz
where
 
dD dT 1 rk dF  ðkT Þ
D0 ¼ ¼ ðF ðkÞ  F  ðkT ÞÞ þ rbF 0 ðkÞk þ a kT ðA7Þ
dz dz ð1 þ rÞT 1þr dz
and
dF ðkÞ
F 0 ðkÞ ¼ . ðA8Þ
dk
d2 T
Grouping the terms with dz2
we get
2
dH dT
¼AB 2 ðA9Þ
dz dz
where
 C. Jensen et al. / Applied Acoustics 67 (2006) 689–699 699
(    !
Agas d^v  p
d^ D  0 dT D
A¼ Re p þ ^v
^  Re ^v^p F ðkÞbk
2 dz dz F ðkÞ dz T ðF ðkÞÞ2
  0  )  2
 D Agas dT cP  ^v  q
þRe ^v^p 
F ðkÞ 2 dz x F ðkÞ 1  r
   
2 d^v k dr 1 dT
 Re  2F 0 ðkÞb þ 1 ImðDÞ þ ImðD0 Þ
^v dz F ðkÞ ð1  rÞ T dz
 2  2
dk gas dT dk stack dT
 Agas  Astack ðA10Þ
dT dz dT dz
and
  2
Agas cP  ^v  q
B¼ ImðDÞ  k gas Agas  k stack Astack . ðA11Þ
2 x F ðkÞ 1  r
And finally
d2 T A
¼ . ðA12Þ
dz2 B

References

[1] Raspet R, Slaton W, Hickey C, Hiller R. Theory of inert gas–condensing vapor thermoacoustics:
propagation equation. J Acoust Soc Am 2002;112:1414–22.
[2] Slaton W, Raspet R, Hickey C, Hiller R. Theory of inert gas–condensing vapor thermoacoustics: transport
equations. J Acoust Soc Am 2002;112:1423–30.
[3] Ward WC, Swift GW. Design environment for low amplitude thermoacoustic engines (DeltaE). J Acoust Soc
Am 1994;95:3671–2. Software and userÕs guide available either from the Los Alamos thermoacoustics web
site at [Link]/thermoacoustics/ or from the Energy Science and Technology Software Center, US
Department of Energy, Oak Ridge, Tennessee.
[4] Arnott WP, Bass HE, Raspet R. General formulation of thermoacoustics for stacks having arbitrarily
shaped pore cross sections. J Acoust Soc Am 1991;90:3228–37.
[5] Flowers BH. An introduction to numerical methods in C++. New York: Oxford University Press Inc.; 2000.
[6] Pierce AD. Acoustics: an introduction to its physical principles and applications. New York: American
Institute of Physics; 1989.
[7] Bhatti MA. Practical optimization methods with mathematicaÒ applications. New York: Springer; 2000.
[8] Raspet R, Brewseter J, Bass HE. A new approximation method for thermoacoustic calculations. J Acoust
Soc Am 1998;103:2395–402.
[9] Hebert D, Atchley A. Measurements of the evolution of the temperature profile in a parallel plate stack. J
Acoust Soc Am 1996;100:2846.
[10] Information and specifications provided by John Corey of Clever Fellows Innovation Consortium Inc.
(CFIC). More information available in US patents no. 6,578,364 and 6,487,859.
[11] Ballister SC, McKelvey DJ. Shipboard electronics thermoacoustic cooler. Master of Science Thesis, Physics
Department, Naval Postgraduate School; 1995.