results in the system of non-dimensional equations
whereM = ¯u(x)/c0 is the local Mach-number.
Such an nondimensionalization is useful to reduce the different scalings in the system, and can be appropriate in many cases when a physical problem is dealt with numerically, since it reduces the potential scale disparities. For the purpose here, where the acoustics are simulated in time to obtain time-series, this would be benificial if a low-order numerical scheme was used, where the low scales of acoustic perturbations might drown in numerical noise (round-off errors and such). But such an nondimensionalization also scales the frequency (by the time scaling), and most importantly the resulting wavelength. A general rule-of-thumb in numerical analysis is that the domain discretization should allow for at least 5 grid-points per wavelength in order to resolve a wave, and thus the spatial discretization becomes very fine and hence computationally more expensive. So later when then governing equations are discretized, they will be in dimensional form, although the scaling presented above are used in some aspects of the method.
When the heat release is introduced in the equations, an appropriate scaling could be taken as energy density per time.
2.4 Heat Release
The acoustical theme of this project is thermoacoustics, a term implying that we consider the conversion of energy in form of heat into acoustic energy.
In general, the medium is considered a premixed gas of reactants, which is ignited at some point in the domain, so that the combustion will raise the tem-perature downstream of the combustion zone. The combustion process itself is not taken into account. It is assumed that the unburnt and burnt gas both behaves as perfect gases, and that no molecular weight change occurs in the com-bustion process. If the comcom-bustion process were to be included more carefully, a model describing the chemical reaction along with evolution equations for the gas components would be needed. Here, the combustion is treated as kind of black-box process, the only outcome being heat addition and its consequences.
A portion of the gas in thermodynamic equilibrium will not have its temperature affected be the passing through of a sound wave i.e. a compression wave. This means that the compression acts as an adiabatic process, so that the entropy
remains constant (no entropy fluctuation). Which is why the acoustic field gen-erally is considered isentropic.
Contrary, an unsteady heat release will cause temperature fluctuations unre-lated to the pressure wave, which are convected by the flow. Therefore down-stream of the flame there will also be an entropy waves present in the system, which is the mentioned convective wave with speed ¯uin the Euler equations in non-conservative form. These entropy fluctuations will appear in the density fluctuation, which can have different effects.
The energy equation would be the most natural place to consider the heat release, since it is an energy source.
However, the path towards inclusion of the heat release source term presented below is the same as in [1]. This derivation is chosen, because it is based on thermodynamic considerations rather than just adding an energy source, and more clearly presents the deviation from the usual assumption of isentropicity in acoustics.
The heat input will affect the thermodynamics of the gas, so if the density is considered a state function of both pressure and entropyρ=ρ(p, s), the chain rule will give
Since the gas is considered inviscid and non heat-conducting, the change in entropy is given by the second law of thermodynamics
ρTDs
Dt =q(x, t) (2.29)
whereq(x, t) is the heat input per unit volume.
The gas is considered perfect, hence ∂ρ/∂s|p =−ρ/cp =−ρT(γ−1)/c2, and the expression (2.28) results in
Dρ
Combined with the continuity equation (2.1), this results in a evolution operator on the pressure similar to the previous derived expression (2.15c), but now with the sought source term caused by the heat release
∂p
∂t +u∂p
∂x+γp∂u
∂x = (γ−1)q . (2.31)
This heat release is considered the cause of temperature changes in the medium.
So a change in the stationary temperature field is the consequence of a stationary
2.4 Heat Release 15
heat release, which changes the properties of the medium but does not act as a source of sound.
However, if the heat release has an unsteady component, this will in turn act as an acoustic source. So when the heat input enters the linearized equations, the source term will consist of the unsteady part, while the steady part is entered implicitly through the mean values.
When the axial extent dof the combustion zone is very small compared to the acoustic wavelength – i.e. kd1,k being the wavenumberω/¯c – it is said to be acoustical compact. Thus this energy source can be viewed as a source onto a plane, represented by a Diracδ-function. Although the stationary heat release could be distributed, the unsteady part is only considered non-zero at the flame position. If distributed steady heat release is part of the model, one would have to consider the conservation laws across this distribution. But here, the general expression used for the heat release is
q(x, t) = ¯q(x) +q0(x, t)
δ(x−β) (2.32)
where x = β is the position of the flame. The steady heat release ¯q hence is discontinuous at the flame position, reflecting a sudden temperature change. If isothermal walls are considered, the steady heat release represents the energy necessary to divide the domain up- and downstream ofx=β, into regions with different but uniform mean flow values. If losses due to cooling by the duct walls was included, this could be modelled by a smoothly varying steady heat release downstream of the combustion zone. A schematic representation of the model and the notations used is presented in Figure 2.1.
Note that the use of compactness is based on the acoustical wavelength, so in terms of the sound field it is well-founded. But the convective wave has a wavelength dependent on the mean flow velocity, and for low Mach number, the assumption of compactnessωd/¯u1 for the entropy wave is valid for very very low frequencies.
When the source enters the system, mathematical and numerical, it does so through balance of fluxes across the source, and as such on assumption of compactness. The entropy wave can have a strong influence on the acoustics under certain circumstances, and as such this lack of compactness is a factor that should be considered, since it might affect the numerical simulations.
The steady heat release changes the temperature, which in turn changes the properties of the gas regarding sound propagation. The characteristic impedance of the medium has a jump at the flame position, which causes an incoming sound wave to be partially transmitted and reflected. Considering the mass and momentum equation in a small volume around the impedance jump, it can be seen that the jump in pressure and velocity perturbations vanishes, as one could expect.
x= 0 x=β x=L
inlet outlet
flame
T¯1, u¯1, ρ¯1, c¯1 T¯2, u¯2, ρ¯2, ¯c2
Figure 2.1: Scenario of the model set-up.
Also, the mean flow velocity changes across the heat release (by conservation of mass), so that all in all, the acoustic modes are affected by such temperature inhomogeneity.
The term thermoacoustic oscillation is related to the general scenario when heat exchange occurs in physics, i.e. the inhomogeneity of the medium caused by (steady) heat exchange and an eventual unsteady heat release. A heat re-lease perturbation will locally cause volumetric expansion, acting as an acoustic monopole source.
Whenever there exists an coupling between the unsteady heat input and the acoustic field in the resonator, it is possible that this leads to the generation of self-sustained high amplitude oscillations in the chamber, the so-called ther-moacoustic instability.
A coupling implies that the unsteady heat release is affected by the acous-tic perturbations though the combustion zone somehow. If for example, the combustion process depends on the mass-flux fluctuation (ρ0u¯+ ¯ρu0) – which could represent a premixed gas where the ”supply” of reactants into the com-bustion zone is proportional to the acoustic field – then a higher sound level gives a higher (unsteady) heat release, giving more power to the acoustic field, increases the heat release, etc. This coupling is what can cause the acoustic instability, it acts as a positive feed-back mechanism. Or where it inverse pro-portional to an acoustic perturbation, it might lead to flame extinction, but this is not considered here.
A condition implying stability is the Rayleigh criterion. This states that when the pressure fluctuations and the unsteady heating are in phase, the acoustic mode will turn unstable. The criterion is given as
Z T 0
q0(t)p0(t) dt >0 (2.33) where T is the period of the oscillation. This criterion being satisfying might not be sufficient to guarantee stability, since also the interaction between the entropy wave and boundaries can lead to thermoacoustic instabilities ([1], [4]).
2.4 Heat Release 17
Modelling the combustion reaction can be a tricky affair, and many differ-ent suggestion can be found in the literature. This often involves combustion chemistry, with the medium composition and more complex combustion reac-tion expressions. In this project, the unsteady heat release model serves the purpose of creating the coupling and thus amplifying the resonant modes in the chamber, i.e. as the phase-locking mechanism so that the fundamental modes are the most dominant. The flame model primarily used here is that of the reference solution, since it is against this that the varying factors are compared and the method evaluated.
Balance conditions
Consider the source term the equation (2.31) and let this be represented by a delta function q(x) =q(x)(x−β), so that the solution is discontinuous across the source position. Hence the equation are not valid across the source, since the solution is not continuous nor differentiable across the source. But this is not of great concern, since a weak solution is sought. So instead, the physical principles are supposed to hold, by requiring the flux of basic conservation quantities, expressed in integral formulation (or equivalently, a weak formulation), to be conserved or balanced in the limit over the source position. These are the Rankine-Hugoniot conditions across this imposed, stationary shock. For the heat release, this reads
[ρu]ββ+−= 0 (2.34a)
[ρuu+p]ββ+−= 0 (2.34b)
u ρE+pβ+
β−=q (2.34c)
where the factor (γ−1) has been dropped with respect to (2.31). When dealing with gas dynamics, these conditions tells something about the state on either side of a shock and here they express the differences in fluxes caused by a source term.
The energy balance (2.34c) tells us something about the scale of the heat release.
The model has the temperature as a tunable parameter, so it would be natural to look at the total enthalpy density fluxρu(e+p/ρ+u2/2). For a perfect gas, this gives
u ρE+p
=ρu cPT+12u2
(2.35) so it is a bit intuitive to take ¯q∝ρ¯1u¯1 cP T¯2−T¯1
+12 u¯22−u¯21
considering in/out of the combustion zone.
The governing equations for the acoustics are linearized, so the Rankine-Hugoniot conditions above is not used as presented, but in a linearized form. The steady-unsteady splitting is introduced in these expressions (2.34), and the steady state form of the conditions is subtracted, again by the steady field implicitly having conservation of mass and momentum, and balanced energy, for the heat release.
These RH-conditions for acoustic perturbations are presented here, and their use is described in later section on the numerical scheme.
[ρ0u¯+ ¯ρu0]β
+
β− = 0 (2.36a) p0+ρ0u¯2+ 2¯ρ¯uu0β+
β− = 0 (2.36b) cPT¯+12u¯2
ρ0u¯+ ¯ρu0
+ ¯ρ¯u cPT0+ ¯uu0β+
β− =q0 (2.36c) The energy flux can also be describes in term of the used dependent variables.
Either by the reformulationu ρE+p
= γ−γ1pu+12ρu3
and linearizing, or equivalently by use of (2.36c) and the relation cPT0 =p0/¯ρ+ ¯T s0 and cPT¯= γp/¯¯ ρ(γ−1). This results in the equation
γ
γ−1 up¯ 0+ ¯pu0
+ ¯u2 12uρ¯ 0+32ρu¯ 0x+ x−
=q0 (2.37)
which also is the used form in the implementation. The case of a stagnant medium (i.e. no flow) is also considered. The RH-conditions on the conserved quantities are not directly applicable in this case. The jump conditions, when the medium is at rest reduces significantly (see [1] for derivation), such that for unsteady heat release, the conditions becomes
[p0]xx+−= 0 and [u0]xx+− =(γ−1)
¯
ρ1¯c21 q0 . (2.38) These jump conditions hence describes the coupling between the variables across a source, which is how they are used when a numerical model of the problem is constructed. As mentioned, the unsteady heat release is a sound source through volumetric expansion, and as such, the jump in primitive variables are expected to be largest (relatively) in the density fluctuations.
When a mean flow is present, the density fluctuation downstream of the point of application of these conditions will carry an ”excess density” (deviation from ρ0−p0/¯c2), which is the entropy wave.
Driving sources
As mentioned, the unsteady heat release acts as an acoustic source, but only in presence of an acoustic field by the coupling of these. So there is need for an
2.4 Heat Release 19
extra source in order to generate the sound field.
Maybe the most simple way to generate an acoustic perturbation in the chamber is to let an inflow boundary condition (at x = xIN) act as an abstraction of a piston, still allowing a mean flow to pass through i.e. an inlet boundary condition on the formu(xIN, t) = ¯u(xIN) +u0e−iωtwithω being the frequency andu0 the amplitude of this acoustic perturbation.
In this project, again, the specific type of boundary condition is important for the resonator and it seemed reasonable not to start mixing a driver signal and the more physical boundary conditions. It might even deteriote the simulation, since incoming waves would ”feel” the piston as well.
A modification of this approach, still in order to generate an acoustic field, is to use the hyperbolicity of the system, decomposing the system into its character-istics, and then imposing the driving source as the incoming wave and letting the outgoing passing through. If the model had non-reflecting boundaries as a requirement, this method would be more appropriate.
The inclusion of driving sources is done by adding source terms to either the mass or momentum conservation equation, which corresponds1 to either externally injection or force representation ([3, p. 28]). So for an external mass-injection, the inhomogeneous conservation law becomes
∂ρ
∂t + ∂
∂x(ρu) = ∂ρfρ
∂t (2.39)
where the right hand side is the rate of injected mass by volume fraction fρ. And for an external force source per unit volume
∂
∂t(ρu) + ∂
∂x(ρuu+p) =fu (2.40)
These two sources are fixed in space and considered compact (δ-sources), and are used to drive the system with a broadband signal. The exact expression of this signal is given later.
During the linearization process, the momentum source does not change (and hence is considered small-scale), while the mass-injection term results in ¯ρ dfρ/dt, again withfρbeing small.
The inclusion of a mass term in (2.39) affects the formulation of the momen-tum conservation in non-conservative form (2.15b) by adding a term−u∂(ρfρ)/∂t or in linearized form−u¯¯ρ∂fρ/∂t. However, when the source is included in the numerical scheme, it enters through considerations about the conservative form of the system. So the term does not matter as such, and is ignored here.
1these are simplifications, the actual physical process can be rather complicated
Either one of these source terms can be applied through the appropriate jump conditions (2.36a) or (2.36b), and be included into the numerical method in the same way as the unsteady heat release, as will be explained later.
By considering the conservation laws (2.39) and (2.40) in a integral formula-tion, one can arrive at more simple jump conditions/estimates for such external sources [3, p.84]. For the mass injection at instant time this reads
[u0]+−=fρ and [ρ0]+−= [p0]+−= 0 (2.41) with the superscripts±denoting either side of the source position. Very similar, the jumps describing an external force are
[p0]+−=fu and [u0]+−= [ρ0]+−= 0 . (2.42) The jump expression (2.41) is used in stagnant medium to represent a mounted loud-speaker with some volume velocityfρ, and when an isentropic model with mean flow is considered (presented in later section), these are the conditions used to enter the driver in the numerical scheme.
The signal used to drive the chamber is considered in greater detail in the following section
The placement of the acoustic drivers for the stagnant medium is not impor-tant, this just creates an acoustic disturbance. But when there is a flow, the placement matters. The driver is imposed via jump conditions, so the solution is discontinuous across the driver. This discontinuity creates entropy fluctua-tions since entropy cannot be conserved across a shock. This entropy fluctuation convects with the flow, so this has to taken into account if the model expects that the acoustic wave entering the combustion zone is isentropic. Placing the driver downstream of the heat release source removes this particular problem with respect to the heat source.
Some of the forms for including a driver source generates stronger entropy waves than others. The simple velocity jump condition (2.41) generates only weaker entropy waves. Hence, if outlet boundary conditions would have an effect by the impinging entropy wave, it would be more appropriate to use the simpler (2.41) in order not to add to this boundary effect unnecessarily..
Driving signal
The reason for solving this acoustic problem in the time-domain, relies on its non-restriction (in regard of formulation) to single frequency perturbations. So the idea is to excite the chamber with a signal which covers some frequency band,
2.4 Heat Release 21
so that the unknown, and possibly more, resonance frequencies are excited. A linear broadband signal is chosen here. In principle, a block signal (or binary signal) would have a broader spectrum, but the propagation of such a signal is more complex in regard of the numerical treatment. Besides, it is expected that the different conditions under consideration, will shift the fundamental, geometric eigenfrequency, but not with large jumps. So the excitation spectrum can be set thereafter, and a signal with a relatively narrow spectrum would still excite the sought instability frequencies, and give, in the end of the process, spectras just as good. So occam’s razor applies, and the simpler, linear driving signal is used.
The expression for the driver signal, also known as a sweep signal, is given as y(t) =Asin (∆ω t+ωmin)t
(2.43) where ∆ω = 2π(fmax−fmin) is the interval of angular frequency and ωmin the lower frequency. The amplitude is set to the scale of the acoustic perturbations.
This signal is used over intervals of 1s, so the simulations are usually run over this time-interval. In Figure 2.2 is shown an example of such a sweep signal.
0 0.2 0.4 0.6 0.8 1
Figure 2.2: Sweep signal example, f ∈[0,10].