Actual impedance of nonreflecting boundary conditions:

Implications for computation of resonators

Laurent Selle, Franck Nicoud, Thierry Poinsot

To cite this version:

Laurent Selle, Franck Nicoud, Thierry Poinsot. Actual impedance of nonreflecting boundary con-
ditions: Implications for computation of resonators. AIAA Journal, 2004, 42 (5), pp.958-964.
�10.2514/1.1883�. �hal-00910165�

HAL Id: hal-00910165

https://hal.science/hal-00910165
Submitted on 27 Nov 2013

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
AIAA JOURNAL
Vol. 42, No. 5, May 2004

Actual Impedance of Nonreflecting Boundary Conditions:

Implications for Computation of Resonators
Laurent Selle∗
Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique, 31057 Toulouse, France
Franck Nicoud†
University Montpellier II, 34095 Montpellier Cedex, France
and
Thierry Poinsot‡
Institut de Mécanique des Fluides, Toulouse, 31400 Toulouse, France

Nonreflecting boundary conditions are essential elements in the computation of many compressible flows. Such
simulations are very sensitive to the treatment of acoustic waves at boundaries. Nonreflecting conditions allow
acoustic waves to propagate through boundaries with zero or small levels of reflection into the domain. However,
perfectly nonreflecting conditions must be avoided because they can lead to ill-posed problems for the mean flow.
Various methods have been proposed to construct boundary conditions that can be sufficiently nonreflecting for
the acoustic field while still making the mean flow problem well-posed. A widely used technique for nonreflecting
outlets is analyzed (Poinsot, T., and Lele, S., “Boundary Conditions for Direct Simulations of Compressible Viscous
Flows,” Journal of Computational Physics, Vol. 101, No. 1, 1992, pp. 104–129; Rudy, D. H., and Strikwerda, J. C., “A
Non-Reflecting Outflow Boundary Condition for Subsonic Navier–Stokes Calculations,” Journal of Computational
Physics, Vol. 36, 1980, pp. 55–70). It shows that the correction introduced by these authors can lead to large
reflection levels and resonant behavior that cannot be observed in the experiment. A simple scaling is proposed to
evaluate the relaxation coefficient used in these methods for a nonreflecting outlet. The proposed scaling is tested
for simple cases (ducts) both theoretically and numerically.

I. Introduction characterizing these instabilities is to study the forced response of

the flame in the channel while disconnecting the coupling with the
D ERIVING nonreflecting boundary conditions for hyperbolic
or incompletely hyperbolic problems1 is a key problem in
multiple fields such as classical engineering fluid dynamics,2−8
acoustics of the channel. This can be done experimentally by chang-
ing the channel geometry upstream or downstream of the burner.
aeroacoustics,9−12 astrophysics,13−15 vibrations in solids,16 and Numerically, the same result can be achieved by diminishing the
electromagnetism.17 acoustic feedback of the channel to the burner, i.e., by using non-
In the field of fluid mechanics, acoustic phenomena can affect reflecting boundary conditions on inlets and outlets (Fig. 1b). This
the flow in a drastic manner in potentially unstable cases.18−20 The allows outgoing waves to leave the domain, but cancels the ampli-
numerical simulation of such flows in compressible codes is a ma- tude of waves entering the domain, thereby diminishing the possible
jor issue. Numerical boundary conditions must be treated with great effects of coupling mechanisms. In certain cases, making only the
caution to predict both the acoustic waves and the mean flow. Even outlet nonreflecting (Fig. 1c) may be sufficient to stabilize the flow,
though techniques have been developed to predict the mean flow as tested in Sec. V. There is no general rule indicating whether mak-
in steady compressible codes (when acoustic phenomena are sup- ing boundaries nonreflecting in a given configuration will indeed be
pressed by numerical or physical viscosity) or to control the acous- sufficient to damp the instability.
tic waves in linearized codes (where the mean flow is imposed), no The decomposition of the Euler equations into characteristic
method can handle both the mean flow and the acoustic waves in a waves is widely used to derive nonreflecting boundary conditions.
perfect way: boundary conditions that allow perfect control of the A specific point that is common to all characteristic techniques is to
mean flow generally reflect acoustic waves, whereas purely non- prescribe the amplitude of the incoming waves.3,5,13,14 The natural
reflecting conditions used, for example, for linearized simulations choice for nonreflecting boundaries is to set this amplitude to zero.
allow mean flow values to drift. This method is well suited to controlling the acoustic field (when
One example where such issues are critical is the prediction acoustic waves reach the boundary at normal incidence) but not the
of combustion instabilities in reacting flows.20−24 Flames exhibit mean flow: various authors2,25 indicate that this choice may lead to
strong combustion instabilities when they are placed in channels. ill-posed problems. Indeed, setting, for example, a zero amplitude
These oscillations are due to coupling between the flame and the for the incoming wave (L1 ) at a subsonic outlet (Fig. 2) may not al-
channel acoustics (Fig. 1a). Acoustic waves are reflected at inlets low the flow to retain a constant mean pressure. In the “real” world
and outlets into the domain where they can couple with the hydro- the mean pressure is imposed by the state at infinity downstream
dynamics and the unsteady heat release. A convenient method for of the burner, where the pressure is P∞ , and this information is fed
back into the domain through L1 . Setting this value to zero is usually
Received 16 April 2003; revision received 17 November 2003; accepted not a good choice because the information on P∞ is simply ignored
for publication 10 December 2003. Copyright  c 2004 by the American by the code, leading to possible drifts on the mean pressure. This
Institute of Aeronautics and Astronautics, Inc. All rights reserved. Copies problem is well known and various solutions have been proposed.
of this paper may be made for personal or internal use, on condition that the One of them is to provide a value for L1 , either using an analytical
copier pay the $10.00 per-copy fee to the Copyright Clearance Center, Inc.,
solution in region 2 or solving a linear problem between the outlet
222 Rosewood Drive, Danvers, MA 01923; include the code 0001-1452/04
$10.00 in correspondence with the CCC. section S and infinity (see recent review edited by Tourrette and
∗ Ph.D. Student, 42 avenue G. Coriolis; [email protected]. Halpern26 ). Such solutions can only be used in certain academic
† Professor, Math Department, ACSIOM Laboratory, CC51. cases for which an analytical solution may be derived between the
‡ Research Director, UMR 5502, Allee du Professeur Camille Soula. outlet of the computational domain and the “infinity” condition.16
Member AIAA. More practical (but less accurate) solutions have been proposed.
958
 SELLE, NICOUD, AND POINSOT 959

a)

b) c)
Fig. 1 Waves entering and leaving the computational domain.

and numerically (Sec. V); and 4) finally, to propose a method to

scale the K coefficient in practical computations (Sec. VI).
It will be shown that the values of K control the damping rate of
the modes but may also change the eigenvalues significantly. In other
words, a nonoptimal choice of K may lead to erroneous resonant
frequencies.
II. Wave Amplitudes
The developments in this section are based on the method derived
by Poinsot and Lele,2 called Navier–Stokes characteristic boundary
conditions (NSCBC). In all characteristic approaches the main is-
sue is the determination of the amplitudes of waves entering the
computational domain. In the NSCBC method the determination is
based on the assumption that the amplitudes can be obtained as if
the flow were laminar, one-dimensional, and inviscid (LODI). Note
Fig. 2 Influence of far-field conditions on the amplitude of the waves
that for one-dimensional nonviscous flows, the NSCBC method is
entering the computational domain.
equivalent26,29 to many other boundary treatments,3,13,14,30 so that
the results presented in Sec. III apply to all characteristic techniques
The simplest one (linear relaxation method; LRM) is to set the am- using LRM.
plitude of L1 as proportional to the static pressure difference,2,25 LODI equations link the wave amplitudes (Li ) and the temporal
L1 = K (P − P∞ ), where P is the predicted pressure at the outlet evolution of primitive Navier–Stokes variables (ρ, u, v, w, P). The
section S, P∞ is the far field pressure, and K has the dimension wave amplitudes L1 , L2 , L3 , L4 , L5 correspond, respectively, to the
of a frequency. “Appropriate” values of K are expected to provide left traveling acoustic wave (speed u − c), the entropy wave (speed
a quasi-nonreflecting boundary condition while avoiding pressure u), the first vorticity wave (speed u), the second vorticity wave
drifts; if P drifts away from P∞ , L1 acts as a spring force to relax (speed u), and the right-traveling acoustic wave (speed u + c). Their
P toward P∞ .27 expression is obtained through characteristic analysis2,27 :
Even though this relaxation approach is used today in numerous


studies,28 two questions arise: ∂ρ 1 1
1) Are there optimal choices for K ? + 2 L2 + (L5 + L1 ) = 0 (1)
∂t c 2
2) What is the effect of this relaxation on the global acoustic
behavior of the boundary? ∂u 1
Obviously, low levels of K may be inefficient in terms of control- + (L5 − L1 ) = 0 (2)
∂t 2ρc
ling the mean pressure. Conversely, high levels of K ensure that P re-
mains very close to P∞ and then make the boundary partially or fully ∂v
+ L3 = 0 (3)
reflecting. Therefore a direct link exists between the magnitude of K ∂t
and the reflection coefficient R of the boundary. Actually K controls ∂w
not only the magnitude of R but also its phase. Choosing K without + L4 = 0 (4)
∂t
care may lead to numerical results in which the outlet boundary
condition introduces important biases in terms of acoustic behavior. ∂P 1
+ (L5 + L1 ) = 0 (5)
The objectives of this paper are 1) to derive a simple analytical so- ∂t 2
lution relating the magnitude of the relaxation coefficient K and the
complex reflection coefficient R of a boundary condition modeled where the wave amplitudes L1 , . . . , L5 are defined by
with LRM (Sec. III); 2) to verify this relation using corresponding  
numerical simulations in which the outlet of a simple tube is sub- ∂P ∂u 1
L1 = (u 1 − c) − ρc (6)
mitted to a harmonic propagating wave (Sec. IV); 3) to demonstrate ∂ x1 ∂ x1
the effect of nonreflecting boundary conditions on the global res-  
onance of a given duct using LRM, a simple case with fixed inlet ∂ρ ∂P
L2 = u 1 c2 − (7)
velocity and nonreflecting LRM outlet is studied both analytically ∂ x1 ∂ x1
960 SELLE, NICOUD, AND POINSOT

∂u 2 The transient term A0 e−K t/2 of Eq. (14) (A0 being a constant
L3 = u 1 (8) fixed by initial conditions) always vanishes with time since K > 0.
∂ x1
In further developments it will be assumed that a steady state has
∂u 3 been reached and this term will be omitted.
L4 = u 1 (9) Equations (14) (without transient part) and (12) make it possible
∂ x1
to reconstruct the incoming wave L1 and the complex reflection
 
∂P ∂u 1 coefficient of the boundary. The analytical value of the reflection
L5 = (u 1 + c) + ρc (10) coefficient R is
∂ x1 ∂ x1
Rout = L1 /L5 = −1/[1 − i(2ω/K )] (15)
Equations (1–5) provide a simple method to choose the incoming
wave amplitudes to be imposed at a boundary. For example, from The magnitude R and phase φ of a nonreflecting outlet modeled
Eq. (2), a fixed velocity inlet condition will require the incoming with LRM are derived from Eq. (15):
wave amplitude L5 to be equal to the outgoing wave L1 . Fixing
a constant pressure at an outlet will be achieved [from Eq. (5)] 
by setting L1 = −L5 . A more critical situation arises in mimicking R = 1 1 + (2ω/K )2 (16)
nonreflecting conditions; this would require setting the incoming
φ = −π − arctan(2ω/K ) (17)
waves to zero. As indicated earlier, such a perfectly nonreflecting
condition is not adequate because it may lead to a drift of the mean
The asymptotic behavior of R and φ is summarized in Table 1.
flow quantities. The next section illustrates the behavior of a duct
As expected, for a given pulsation ω, R goes to 0 when K is
outlet for which a nonreflecting condition is sought.
small, showing that the boundary condition is indeed nonreflecting
(R  0) when K is limited to small values. However, large values
III. Reflection Coefficient of a Linear Relaxation of K destroy the nonreflecting character of the boundary condition:
Method Boundary Condition when K goes to infinity R goes to 1, making the boundary fully
Let us consider the propagation of one-dimensional acoustic reflecting. Actually, as shown by Eq. (16), the control parameter for
waves in a semi-infinite tube of constant cross section (Fig. 3). The R is 2ω/K , so that one can define a cutoff pulsation ωc = K /2.
tube is infinite in the x < 0 direction and ends at x = x B in the other For a fixed value of K (which is the case in any computation),
direction, where a nonreflecting boundary condition must be imple- all frequencies will not be reflected with the same strength. High
mented. An harmonic wave propagating in increasing x direction frequencies will easily leave the computational domain (R → 0),
is imposed. This wave amplitude is chosen so that in the absence whereas very low frequencies will be strongly reflected (R → 1). In
of reflected wave, the inlet velocity signal would be u(t) = U0 e−iωt . practice a cutoff frequency f c separates waves that will be reflected
Note that the phase is set to zero at x = x B to simplify the algebra. ( f < f c ) from the ones that will leave the domain ( f > f c ). f c is
This has no influence on the result. The expression of this complex defined from Eq. (15) by
wave at x = x B is taken as
f c = ωc /2π = K /4π (18)
L5 = 2ρcU0 iωe−iωt (11)
√
This definition implies that R( f c ) = 1/ 2. In terms of energy this
The objective of the condition at x = x B is to be nonreflecting. In means that at the frequency f c half of the acoustic energy is fed back
practice, to avoid a drift of the mean pressure, the incoming wave into the computational domain.
amplitude L1 is not set to zero but to Figure 4 is a plot of Eqs. (6) and (17). The cutoff frequency f c
is represented by the vertical line. Equation (18) suggests that a
L1 = K (P − P∞ ) (12) proper interpretation of K is to view f c = K /(4π) as a frequency
below which the boundary condition will not let the waves leave the
Equations (11) and (12), together with LODI relations (2) and domain.
(5), lead to the system of equations
∂u 1 Table 1 Asymptotic behavior of a
+ (2ρcU0 iωe−iωt − K (P − P∞ )) = 0 nonreflecting outlet with LRM
∂t 2ρc
Pulsation R φ
∂P 1
+ (2ρcU0 iωe−iωt + K (P − P∞ )) = 0 (13) ω=0 1√ −π
∂t 2
ω = K /2 1/ 2 −π − π/4
The second equation of system (13) involves only P and can ω=∞ 0 −3π /2
easily be solved. (This system was derived at the outlet boundary;
consequently, coordinates are fixed and both u and P are functions
of time only.) The solution for P is
ρcU0 iω −iωt
P(t) = P∞ + A0 e−K t/2 − e (14)
K /2 − iω

Fig. 3 One-dimensional harmonic wave impacting on the outlet Fig. 4 Modulus and phase of reflection coefficient vs frequency at
boundary of a tube. K = 2000 s−1 .
 SELLE, NICOUD, AND POINSOT 961

Table 2 Summary of numerical Consider a simple duct filled with a homogeneous gas (in which
simulation parameters the speed of sound c is constant), as represented in Fig. 6. The
inlet speed is imposed: u(x = 0, t) = U0 . The steady flow solution
Parameter Value
is of course u(x, t) = U0 . This system is acoustically defined by its
Size of the domain L, m 0.5 reflection coefficients at the inlet and the outlet. Assuming that the
Number of cells 400 duct inlet is the phase reference, one can write the relations between
Mean inlet velocity, m · s−1 10 the amplitude of the acoustic waves in the duct and the reflection
Sound speed c, m · s−1 348 coefficients:
Forcing frequency f , Hz 500
Forcing amplitude, m · s−1 0.1 L5 L1 e−iω(L/c)
Rin = , Rout = (19)
L1 L5 eiω(L/c)

which can be written as

Rin Rout eiω(2L/c) − 1 = 0 (20)

Calculating the eigenfrequencies of this duct consists in giving

the values of ω that are solutions of Eq. (20).

B. Analytical Solutions
For this test, the inlet corresponds to fixed velocity conditions,
which implies that Rin = 1. The outlet is modeled using a nonreflect-
ing condition with LRM. Therefore Eq. (15) shows that the effective
impedance of this outlet is Rout = −1/[1 − i(2ω/K )]. In this case
Eq. (20) degenerates to

eiω(2L/c) + [1 − i(2ω/K )] = 0 (21)

The values of ω that satisfy Eq. (21) are complex. The real part
Fig. 5 Comparison of numerical and theoretical reflection coefficients
at a “partially reflecting” outlet.
of ω (ωr = (ω) = 2π fr ) is the eigenpulsation and the imaginary
part (ωi = (ω) = 2π f i ) is the damping (or amplification) rate. The
temporal evolution of pressure and velocity at an eigenfrequency is
At constant K , the phase φ of the reflected signal is that of a fixed- proportional to e−iωt = e−iωr t eωi t (in the linear regime):
pressure outlet for low values of ω. But as shown in Fig. 4, when the 1) If ωi > 0, the mode is amplified: ωr is an unstable pulsation.
frequency increases, the boundary treatment induces a delay that 2) If ωi < 0, the mode is damped: ωr is a stable pulsation.
moves the phase from −π to −3π/2. This phase shift can induce Figure 7 is a plot of the eigenfrequencies of Eq. (21) vs K . At
dramatic changes in the acoustic properties of the domain, as shown very high values of K the system responds as if the pressure were
in Sec. V. fixed at the outlet. Combining Eqs. (11), (12), and (14) one can

IV. Numerical Evaluation of the Reflection Coefficient

In this section Eq. (15) is compared to numerical simulations.
This is done by simply running a full Euler code on the configura-
tion of Fig. 3 and using Eq. (12) at the outlet. The code AVBP was
used for this validation. (The AVBP home page for detailed informa-
tion is http://www.cerfacs.fr/cfd/CFDWeb.html.) AVBP is a three-
dimensional fully compressible Navier–Stokes equation solver us-
ing characteristic boundary conditions.2
The outlet condition is nonreflecting (with LRM) as described in
Sec. III and the inlet is a nonreflecting pulsed inlet.31 The numerical
scheme is second order in space and uses a three-step Runge–Kutta Fig. 6 Configuration.
method (third order) for time integration. Calculation parameters
are summarized in Table 2.
For technical reasons, it is more convenient to vary K at the outlet
than to change the forcing frequency f at the inlet. Therefore the
following results present R(K ) at fixed f instead of R( f ) at fixed
K . The reflection coefficient R is obtained in AVBP by R = L1 /L5 .
The amplitudes L1 and L5 are measured at the outlet using Eqs. (6)
and (10) and one-sided spatial derivatives.
The agreement between numerical simulations and the theory
[Eq. (15)] is extremely good (Fig. 5). Both modulus and phase of the
reflection coefficient at the outlet are superimposed on the theoretical
curves.

V. Computing a Closed Domain with Nonreflecting

Boundary Condition
A. Description of the Configuration
The preceding sections have shown that nonreflecting boundary
conditions using LRM can have impedances that make them par-
tially reflecting. This section shows the implications of these results
for the computation of flows where one or more sections are modeled Fig. 7 Eigenfrequencies of the first three modes vs K for the duct of
using such nonreflecting conditions. Fig. 6.
962 SELLE, NICOUD, AND POINSOT

Fig. 9 Summary of the influence of σ on the acoustics and mean flow

quantities.

is lower than all acoustic frequencies expected in the computational

Fig. 8 Imaginary part of the first three modes vs K for the duct of domain.
Fig. 6. As an example, this strategy is now applied to the duct of Sec. V.
Since velocity is imposed at the inlet and pressure at the outlet, the
lowest acoustic frequency is that of a quarterwave mode. At a given
show that L5 + L1 goes to zero when K goes to infinity, and the Mach number M the quarter wave mode frequency is
eigenfrequencies
f 0 = (1 − M2 )(c/4L) (24)
f n = (2n + 1)(c/4L) (22)
From Eqs. (18) and (24) the highest value of K for this calculation
correspond directly to the 14 wave (n = 0), 34 wave (n = 1), . . . modes that will significantly damp the 1/4 wave mode (as well as other
of a duct with fixed inlet velocity and outlet pressure (represented higher order modes) is
by horizontal lines in Fig. 7): the boundary is fully reflecting.
When K decreases the eigenfrequencies of the duct also decrease K max = π(1 − M2 )(c/L) (25)
due to the phase shift induced by the boundary treatment shown by
Eq. (17). At very low K values the 1/4 wave mode even disappears. This scaling of K is very similar to the theory of Rudy and
At very high values of K the imaginary part f i of all eigenfre- Strikwerda25 given in Eq. (23) but with a σ coefficient of π in-
quencies is zero (Fig. 8), indicating that they are not damped. This stead of 0.58. This approach suggests that the highest admissible
is consistent with the fact that there is no source term in the duct and value of σ to prevent acoustic feedback is
boundary conditions are acoustically closed (no fluctuation of P or
u allowed). Obviously at such large K values the outlet boundary σ max = π (26)
condition fails to evacuate acoustic waves. For lower values of K
the modulus of Rout is lower than 1 and thus all modes are damped On the other hand, K (or σ ) has to be chosen large enough to
( (ω) < 0). prevent a drift in the mean values (pressure, mass flow, etc.). Fig-
Figures 7 and 8 show that the value of K must remain small to ure 9 offers a simple summary of these findings. (Note that Fig. 9
provide damped modes but also that the modes that appear are not is actually a one-dimensional plot.) The minimum admissible value
physical: their frequencies are not the eigenfrequencies of the duct. σmin is not fixed by acoustics. It strongly depends on the computa-
tional parameters (Reynolds number, three-dimensional effects such
VI. Scaling Strategy for the Relaxation Coefficient K as swirl and geometry, etc.). Consequently, the boundary σmin is not
clearly defined in Fig. 9. Other tests (not reported here) suggest that
As shown in the preceding sections, the value of K has a drastic choosing σ lower than 0.1 often increases convergence times and
influence on the results of a numerical simulation. This has already sometimes does not allow mass fluxes and pressure to reach a steady
been pointed out by different authors for either steady25 or unsteady2 state.
calculations. The main problem for several authors is to know how
to choose K in practical cases. Rudy and Strikwerda25 suggested VII. Numerical Computations
the following scaling of K for optimal convergence of steady cal-
culations: Section V.B has provided a theoretical stability analysis of the duct
of Fig. 6. In the present section, these results are tested numerically.
K = σ (1 − M2 )(c/L) (23) The duct inlet is a fully reflecting characteristic inlet with imposed
velocity, achieved through LODI relation 2 by imposing L5 = L1 .
where M is the Mach number of the mean flow, c the sound speed, A nonreflecting outlet with LRM is applied at the outlet; the value
and L the domain size. of σ is varied to illustrate the results of Secs. V and VI. A Gaussian
Numerical simulations2,25,27 showed that an optimum is reached perturbation is superimposed on the initial pressure field to excite
optim
for σNum = 0.58, whereas the theory of Rudy and Strikwerda25 sug- all acoustic modes.
optim
gested an optimum value of σRS = 0.27. Rudy and Strikwenda’s Two values for σ are investigated to illustrate the influence on
definition of the optimum value for σ is based on the convergence both the damping of acoustics and the shift of the eigenfrequencies
of steady-state calculations, whereas in Refs. 2 and 27 the definition of the duct.
is based on both the convergence of mean values and the evacuation
of acoustic waves in unsteady calculations. A. Case 1: σ = 10π
π
The present work gives a new interpretation of K , which is now In this section, the relaxation factor at the outlet σ is fixed to
linked to the cutoff frequency of the boundary by Eq. (18). Thus σ = 10π, which is too high to evacuate the first eigenmodes of the
“good” values of K are those that allow all duct acoustic modes to duct. The pressure perturbations recorded at the inlet show exponen-
leave the domain. It was shown in Sec. III that frequencies lower tial decay (Fig. 10) before t = 0.1s, corresponding to the evacuation
than f c = K /(4π ) are reflected, and frequencies higher than f c leave of high-frequency modes. The signal then remains almost constant,
the computational domain. It is then relevant to choose K so that f c because low-frequency modes are reflected by the outlet condition.
 SELLE, NICOUD, AND POINSOT 963

Fig. 10 Time evolution of pressure perturbation at the inlet (σ = 10π). Fig. 12 Time evolution of pressure perturbation at the inlet (σ = π).

Fig. 11 Spectrum of pressure signal (σ = 10π) compared to analytical

values of Sec. V and Eq. (22). Fig. 13 Spectrum of pressure signal (σ = π) compared to analytical
values of Sec. V and Eq. (22).

Figure 11 is the Fourier transform of the inlet pressure signal. The

spectrum is compared to the expected eigenfrequencies analytically is no source term in the domain and that the simulation is quite
derived in Sec. V. The first three modes predicted by Eq. (21) and long compared to f c (35 cycles), so that even a reflection coefficient
displayed in Fig. 11 remain in the computational domain, as sug- slightly√lower than 1 (the definition of the cutoff frequency f c is that
gested from Fig. 9 for σ = 10π. Moreover, the frequencies of these R = 1/ 2) attenuates the mode.
modes are not the true eigenfrequencies as predicted by Eq. (22), Figure 13 is the Fourier transform of the pressure signal of Fig. 12
which gives the eigenfrequencies of an acoustically closed duct (im- between t = 0 and 0.05s. The spectrum is compared to the expected
posed velocity at the inlet and imposed pressure at the outlet). This eigenfrequencies analytically derived in Sec. V and to the true eigen-
confirms the biasing effects induced by the outlet condition when frequencies of the duct given by Eq. (22). The spectrum is not very
large values of σ are used: sharp, due to the exponential decay of the signal, and yet the agree-
1) The outlet is strongly reflecting, so that certain modes are not ment is fairly good. Again the eigenfrequencies of the duct are
damped. significantly shifted compared to those of an acoustically closed
2) The frequencies of the modes differ from the true eigenfre- duct.
quencies of the duct.
Obviously, this would not be a good choice for σ for practical VIII. Conclusions
computations. A simple acoustic theory was proposed to characterize the ac-
The agreement between numerical and analytical results is good; tual reflection coefficient of numerical “nonreflecting” boundary
the shift between the eigenfrequencies of the closed duct and the conditions using LRM (linear relaxation method) as proposed by
duct with LRM is correctly predicted. Rudy and Strikwerda25 or Poinsot and Lele.2 It has been shown
that large values of the relaxation coefficient K used in these meth-
B. Case 2: σ = π ods to link the ingoing wave amplitude L1 to the pressure differ-
In this section the relaxation factor σ is set to π, which is the ence (P − P∞ ) make the boundary condition fully reflecting and
highest value that enables the evacuation of all acoustic modes of that K can be linked to a cutoff frequency f c by K = 4π f c . At a
the duct. given value of K , all modes such that f < f c are not damped. A
Figure 12 shows the exponential decay of the pressure perturba- proper strategy to minimize acoustic coupling is to choose K such
tions. At t = 0.1s all acoustic modes have vanished. Note that there that f c is lower than the first duct acoustic mode. In the case of a
964 SELLE, NICOUD, AND POINSOT

one-dimensional duct with imposed velocity inlet, this is obtained in the High Sonar Corona,” Astronomy and Astrophysics, Vol. 362, 2000,
by writing K = σ (1 − M2 )(c/L) and choosing 0.2 < σ < π. For pp. 342–358.
16 Barry, A., Bielak, J., and MacCamy, R. C., “On Absorbing Bound-
more complex cases (three dimensions, complex geometries, com-
bustions, etc.), an extension of this strategy is to solve the general ary Conditions for Wave Propagation,” Journal of Computational Physics,
acoustic equation in the domain to find its eigenfrequencies f i . The Vol. 79, 1988, pp. 449–468.
17 Mur, G., “Absorbing Boundary Conditions for the Finite-Difference of
maximum value of σ is then given by choosing f c = min( f i ). the Time-Domain Electromagnetic-Field Equations,” IEEE Transactions on
Electromagnetic Compatibility, Vol. 23, 1981, pp. 377–382.
References 18 Ho, C. M., and Huerre, P., “Perturbed Free Shear Layers,” Journal of
1 Engquist, B., and Majda, A., “Absorbing Boundary Conditions for the Fluid Mechanics, Vol. 16, 1984, pp. 365–424.
19 Lucas, L., and Rockwell, D., “Self-Excited Jet: Upstream Modulation
Numerical Simulation of Waves,” Mathematics of Computations, Vol. 31,
1977, pp. 629–651. and Multiple Frequencies,” Journal of Fluid Mechanics, Vol. 147, 1984,
2 Poinsot, T., and Lele, S., “Boundary Conditions for Direct Simulations of pp. 333–352.
20 Crighton, D. G., Dowling, A., Heckl, M., Leppington, F., and Williams,
Compressible Viscous Flows,” Journal of Computational Physics, Vol. 101,
No. 1, 1992, pp. 104–129. J., Modern Methods in Analytical Acoustics, Springer-Verlag, Berlin, 1992,
3 Hirsch, C., Numerical Computation of Internal and External Flows, Vol. Chap. 13.
21 Yang, V., and Culick, F. E. C., “Analysis of Low-Frequency Combustion
2, Wiley, New York, 1988, Chap. 19.
4 Hagstrom, T., and Hariharan, S. I., “Accurate Boundary Conditions for Instabilities in a Laboratory Ramjet Combustor,” Combustion Science and
Exterior Problems in Gas Dynamic,” Mathematics of Computation, Vol. 51, Technology, Vol. 45, 1986, pp. 1–25.
22 Poinsot, T., Trouvé, A., Veynante, D., Candel, S., and Esposito, E.,
1988, pp. 581–597.
5 Tourrette, L., “Artificial Boundary Conditions for the Linearized Com- “Vortex Driven Acoustically Coupled Combustion Instabilities,” Journal of
pressible Navier–Stokes Equations,” Journal of Computational Physics, Fluid Mechanics, Vol. 177, 1987, pp. 265–292.
23 McManus, K., Poinsot, T., and Candel, S., “A Review of Active Control
Vol. 137, 1997, pp. 1–37.
6 Hayder, M. E., and Turkel, E., “Nonreflecting Boundary Conditions for of Combustion Instabilities,” Progress in Energy and Combustion Science,
Jet Flow Computations,” AIAA Journal, Vol. 33, 1995, pp. 2264–2270. Vol. 19, 1993, pp. 1–29.
7 Tsynkov, S. V., “Numerical Solution of Problems on Unbounded 24 Angelberger, C., Egolfopoulos, F., and Veynante, D., “Large Eddy Sim-

Domains: A Review,” Applied Numerical Mathematics, Vol. 27, 1998, ulations of Chemical and Acoustic Effects on Combustion Instabilities,”
pp. 465–532. Flow Turbulence and Combustion, Vol. 65, No. 2, 2000, pp. 205–222.
8 Vasilyev, O. V., and Bowman, C., “Second-Generation Wavelet Collo- 25 Rudy, D. H., and Strikwerda, J. C., “A Non-Reflecting Outflow Bound-

cation Method for the Solution of Partial Differential Equations,” Journal of ary Condition for Subsonic Navier–Stokes Calculations,” Journal of Com-
Computational Physics, Vol. 165, No. 2, 2000, pp. 660–693. putational Physics, Vol. 36, 1980, pp. 55–70.
9 Kim, J. W., and Lee, D. J., “Generalized Characteristic Boundary Con- 26 Tourrette, L., and Halpern, L., Absorbing Boundaries and Layers, Do-

ditions for Computational Aeroacoustics,” AIAA Journal, Vol. 38, 2000, main Decomposition Methods: Applications to Large Scale Computation,
pp. 2040–2049. Nova Science, Huntington, NY, 2001.
10 Colonius, T., Lele, S., and Moin, P., “Boundary Conditions for Direct 27 Poinsot, T., and Veynante, D., Theoretical and Numerical Combustion,

Computation of Aerodynamic Sound Generation,” AIAA Journal, Vol. 31, R. T. Edwards, Flourtown, PA, 2001.
28 Wall, C., Pierce, C. D., and Moin, P., “A Semi-Implicit Method for
No. 9, 1993, pp. 1574–1582.
11 Rowley, C. W., and Colonius, T., “Discretely Nonreflecting Bound- Resolution of Acoustic Waves in Low Mach Number Flows,” Journal of
ary Conditions for Linear Hyperbolic Systems,” Journal of Computational Computational Physics, Vol. 181, 2002, pp. 545–563.
29 Nicoud, F., “Defining Wave Amplitude in Characteristic Boundary
Physics, Vol. 157, 2000, pp. 500–538.
12 Freund, J. B., “Proposed Inflow/Outlow Boundary Condition for Di- Conditions,” Journal of Computational Physics, Vol. 149, No. 2, 1998,
rect Computation of Aerodynamic Sound,” AIAA Journal, Vol. 35, 1997, pp. 418–422.
30 Giles, M., “Non-Reflecting Boundary Conditions for Euler Equation
pp. 740–742.
13 Thompson, K. W., “Time Dependent Boundary Conditions for Hy- Calculations,” AIAA Journal, Vol. 28, No. 12, 1990, pp. 2050–2058.
31 Kaufmann, A., Nicoud, F., and Poinsot, T., “Flow Forcing Techniques
perbolic Systems,” Journal of Computational Physics, Vol. 68, 1987,
pp. 1–24. for Numerical Simulation of Combustion Instabilities,” Combustion and
14 Thompson, K. W., “Time Dependent Boundary Conditions for Hy- Flame, Vol. 131, 2002, pp. 371–385.
perbolic Systems,” Journal of Computational Physics, Vol. 89, 1990,
pp. 439–461. S. Mahalingam
15 Grappin, R., Lorat, J., and Buttighoffer, A., “Alfvén Wave Propagation Associate Editor