Modeling the Generation of Entropy Waves of a Technically Premixed Swirl Burner Setup 

S

Modeling the Generation of Entropy Waves of a Technically Premixed Swirl Burner Setup

Autor:

Anne Berger Matrikel-No:

03637603 Betreuer:

Dipl.-Ing. Thomas Steinbacher Prof. Wolfgang Polifke, Ph. D.

September 29, 2017

Professur für Thermofluiddynamik

Hiermit versichere ich, die vorliegende Arbeit selbstständig verfasst zu haben. Ich habe keine anderen Quellen und Hilfsmittel als die angegebenen verwendet.

Ort, Datum Anne Berger

I would like to express my gratitude to every person who offered their help, which was crucial for the accomplishments during this work.

I would like to especially thank

Thomas Steinbacherfor patiently teaching me the applied CFD and his supervision through- out the process,

Alexander Avdoninfor helpful discussions facilitating the start with the Fluent software, Prof. Wolfgang Polifkefor the opportunity to conduct this semester thesis at his chair

A more detailed insight into the transport mechanisms of entropy waves in a turbulent tech- nically premixed swirl combustor is established by a detailed analysis of the entropy transfer function with an underlying basic approach known from system identification. This novel method, so far only applied on laminar flames, is for the first time implemented on a turbu- lent burner using the Fluent software package. A computationally simple state space model was applied to account for the advection of fluctuations while the combustion chemistry was included via the static gain. To set a basis for the analysis, the simulation and validation of the turbulent flow field for the burner setup was conducted by applying the probability density function theory. The entropy transfer function was determined by analyzing the step response of the system. Using this approach, the influence of a diabatic boundary condition as well as the impact of diffusion could be evaluated. It was shown that the realistic diabatic boundary condition for the fluctuations leads to a dampening of the response of the system, but does not change its general properties. After an analysis of laminar, turbulent, and overall diffusion, it was found that laminar diffusion has only a minor impact in contrast to turbulent diffusion.

The overall diffusion was found to be a major contributor to the penetration of fluctuations into the inner recirculation zone. In an approach to model the impulse response of the system, the data shows a superposition of multiple phenomena: most of the transport of the pertu- bation takes place similar to a 1D duct over the major streamline and the response can be described with a Gaussian distribution. The phenomenon of advective dispersion was found to be the main contributor to the distribution of the step response front. The tailing behavior impulse response can be explained by the occurance of backmixing in the inner recirculation zone. The influence of diffusion was narrowed down to a higher degree of distribution of the signal in the 1D duct and more backmixing in the inner recirculation zone.

Nomenclature viii

1 Introduction 1

2 Numerical Flow Field Modeling 3

2.1 Overview on DNS, LES, RANS, and PDF . . . 3

2.1.1 Direct Numerical Simulation . . . 3

2.1.2 Large Eddy Simulation . . . 4

2.1.3 Reynolds Averaged Navier Stokes . . . 5

2.1.4 Probability Density Function . . . 7

2.2 Turbulent Closure Models . . . 8

2.3 Modeling of Combustion Chemistry . . . 14

3 Flame Response of Technically Premixed Flames 16 3.1 System Identification . . . 16

3.2 Entropy Transfer Function Modeling . . . 18

3.2.1 State Space Model . . . 19

3.2.2 Fast-Reaction-Zone Model . . . 20

3.2.3 Implementation in Fluent . . . 20

3.2.4 Diffusion Coefficient . . . 21

4 Kornilov Flame 24 4.1 Burner Geometry, Numerical Set-Up, and Boundary Conditions . . . 24

4.2 Flow Field without Combustion . . . 25

4.3 Diabatic Case . . . 25

4.4 Adiabatic Case . . . 27

5 BRS Burner 29 5.1 Burner Geometry, Numerical Setup, and Boundary Conditions . . . 29

5.2 Cold Configuration . . . 30

5.3 Hot Configuration . . . 32

5.4 Mesh Independance . . . 37

6 Entropy Transfer Function 39 6.1 Kornilov Case . . . 39

6.2 BRS Burner . . . 41

6.2.1 Standard Configuration . . . 41

6.2.2 Diabatic Boundary Condition . . . 43

6.2.3 Influence of Diffusion . . . 46

6.2.4 Validation of Method . . . 49

6.2.5 Mesh Independence . . . 50

7 Analytical Model 52 8 Conclusion and Outlook 57 Appendices 59 A 2S-CMS Reaction Mechanism in Chemkin Format 60 A.1 Mechanism Data File . . . 60

A.2 Transport Data File . . . 60

A.3 Thermodynamic Data File . . . 61

B BRS Scalar Transport 62 B.1 Diffusion Coefficient . . . 62

B.2 Spatial Distribution for Turbulent Diffusion Only . . . 63

B.3 Spatial Distribution for Laminar Diffusion Only . . . 64

B.4 Display of False Interpolation . . . 65

Bibliography 66

Roman Symbols

Q˙ Heat Release

〈sc〉s Average Flame Consumption Speed along the Surface Pr Prandtl Number

Sc Schmidt Number

a_i Constant for the NASA polynomials C_P Specific Heat Capacity

C_iε Constants for the k-ε-Model D_ω Cross-Diffusion Term Ea Activation Energy

F Blending Function for the SST-k−ωModel

f Volume Force

G_b Generation of Turbulent Kinetic Energy due to Buoyancy

G_k Generation of Turbulent Kinetic Energy due to Mean Velocity Gradient J_i,i Molecular Diffusion Flux Vector

Mi Molar Mass of Component i M_t Turbulent Mach Number

p_n Probability of Each Mode in the Eulerian Scheme r_k Reaction Rate for Specie k

Si j Strain Rate Tensor

u Velocity

Y_M Contribution of Fluctuating Dilation

A Pre-exponential Factor for the Arrhenius Equation

a Speed of Sound

b Factor in Arrhenius Equation g Gravitational Constant

h Enthalpy

k Turbulent Kinetic Energy

n Number of Units with Backmixing

O Operator To Connect An Input Signal With The Output Signal P Favre joint PDF of composition

p Pressure

Q Specific Heat Release

S Entropy

T Temperature

t Time

V Volume

x Spatial Variable Y Mass Fraction

Greek Symbols

β Thermal Expansion δ Delta Function

² Lennard-Jones Parameter η Dynamic Viscosity

γk Angular Velocity of Reference Frame for CalculatingΨi j

κ Kronecker delta function

ν Static Gain Ω Collision Integral

ω Specific Dissipation Rate

Ψi j Mean Rate of Rotation in a Rotating Reference Frame φ Equivalence Ratio

ψ Composition Space Vector

ρ Density

Σ Flame Surface Density σ Lennard-Jones Parameter θ Flow Variable

θm Conditional Mean Composition of Speciemin then^thmode ε Dissipation Rate

ϑ Arbitrary Scalar ξ Mixture Fraction

Acronyms

2S-CMS 2-Step Mechanism As Developed By Bibrzycki [16]

ARX Autoregressive Model with Exogenous Input BRS Beschaufelter Ringspalt

CFD Computational Fluid Dynamics CV Control Volume

DNS Direct Numerical Simulation EBU Eddy Break Up

ETF Entropy Transfer Function FRZ Fast Reaction Zone

IR Impulse Response LES Large Eddy Simulation

LOS Line of Sight

LTI Linear-Time-Invariant

MIMO Multiple-Input-Multiple-Output PDF Probability Density Function RANS Reynolds Averaged Navier Stokes SR Step Response

SST Shear Stress Transport UDS User Defined Scalar

The topic of lean premixed combustion has received increasing attention during the past years due to its lowered NO_x-emissions. However, these systems are very prone to instabilities [25, 39] which in turn generate noise. Fluctuations in the heat release cause direct combus- tion noise as pressure fluctuations are generated and indirect noise as temperature fluctua- tions are generated. The latter are known as temperature hot spots and in literature, they are also termed entropy waves [31]. Entropy waves themselves are silent [2, 32], however, when they are accelerated in the turbine, inhomogeneities in temperature and therefore inhomo- geneities in density result in a sound generating volume contraction.

Indirect combustion noise is subdivided into vorticity and entropy noise with the latter being thought of as its main component. Entropy noise is generated when temperature hot spots are accelerated in the turbine. Morgans et al. have separated five important stages when talk- ing about entropy waves [31]: 1) the generation of entropy waves, 2) the advection of entropy waves through the combustor, 3) the acceleration of entropy waves through a nozzle, 4) the passage of entropy noise through the turbine, and 5) the reflection of entropy noise back into the combustor. Entropy waves are generated due to e.g. fluctuations in the equivalence ra- tio (φ) [26, 36] directly at the flame front [45]. After generation, the entropy waves are advected through the combustor. Many studies have neglected any alteration of the entropy waves be- tween the flame and the combustor exit [7, 10, 11]. However, lately it has been found that, as a governing effect, the dispersion of the entropy waves in the combustion chamber leads to a significant modification [32, 36] and therefore their transport to the turbine inlet must be considered. The first to develop a model that considers the effects on a turbulent flow field in the combustion chamber on entropy waves was Sattelmayer [36]. He underlined the im- portance of realistically implementing the dispersion of equivalence ratio fluctuations before the flame and the entropy waves after the flame. Since then, many groups have analyzed the transport through the combustor [32, 38, 42, 45]. It has been found that it is necessary to con- sider the temperature fluctuations independent from the global heat release rate. Steinbacher et al. have only found recently that the entropy response shows a very different behavior than the heat release response using two sub-models [38]. A controversy in the literature remains on the significance of the turbulence. Morgans et al. have stated that the dispersion of the mean flow fields is the governing effect in the overall dispersion of the entropy waves. In con- trast, Xia et al. have found a significant difference between the time-averaged (mean) fields versus the time-varying (instantaneous) fields [45]. After the turbulent transport in the com- bustion chamber, the entropy waves undergo acceleration in the turbine nozzle causing the generation of entropy noise [10, 31]. Then, the entropy noise passes through the turbine cross-

ing multiple turbine blade rows until they appear at the exit where they can get reflected and while traveling upstream can pertubate the earlier stages and even the flame [31] creating a feedback loop.

There were many attempts to measure entropy waves and entropy noise experimentally, but many difficulties such as the harsh conditions in a burner or turbine exacerbate the experi- ments. Measuring the temperature fluctuations has been difficult since high frequency tem- perature oscillation data is hard to obtain with common measurement equipment [31]. When trying to measure the noise, it is complicated to distinguish between direct and indirect com- bustion noise [31]. Efforts have been made to correlate different frequencies to direct and indirect combustion noise [42], however simulation is still considered as a very necessary and useful approach to separate the effects.

This work will address the propagation of entropy waves in a turbulent swirl-stabilized com- bustor using computational fluid dynamics (CFD) combined with an approach known from system identification. Applying techniques from system identification in the domain of com- bustion modeling has been proven as an useful tool to describe the response of the flame to pertubation [34, 38]. By modeling the flame response as an multiple-input-multiple-output (MIMO) system, it is possible to correctly link the input variables such as the velocity (u) and the equivalence ratio (φ) to the output values of the heat release ( ˙Q) and the temperature (T).

This thesis is divided into three major parts with two chapters each. The first two chapter cover a short summary of the theoretical background on different modeling approaches of turbulent flow fields, the chemical reaction models, and system identification theory in com- bination with flame transfer functions. In the next two chapters, the simulation of the flow field of two burners is discussed where the multiple slit burner by Kornilov et al. served as a validation case while the "Beschaufelter Ringspalt" (BRS) burner was the case under more de- tailed investigation. The third important set of chapters deals with the analysis of the entropy transfer function (ETF) including the development of a model to predict the flame response.

This chapter covers the theoretical background behind the physical and chemical phenomenon which are necessary to explain the applied models regarding the flow field simulation. First, turbulent flow field modeling is explained using different simulation approaches such as Di- rect Numerical Simulation (DNS), Reynolds-Averaged-Navier-Stokes (RANS) Simulation, Large- Eddy-Simulation (LES), and Probability-Density-Function (PDF) Simulation. Turbulent clo- sure models such as the k-ε, the realizable k-ε, the k-ω, and the Shear-Stress Transport Model (SST) are presented.

2.1 Overview on DNS, LES, RANS, and PDF

When simulating turbulent combustion, many complex phenomena and processes have to be estimated. Therefore, quantifying different contributions and combining many physical effects into sets of equations has been a major accomplishment in the domain of CFD. Even though it would be desirable to describe all processes in as much detailed as possible, there is always a compromise between the complexity of the numerical approach and the compu- tational effort that must be undertaken. In general, four different methods are distinguished and explained in more detail in the following sections.

2.1.1 Direct Numerical Simulation

The first method is Direct Numerical Simulation (DNS). With DNS the continuity equation (Eq. 2.1 ), the Navier-Stokes equations, the energy equation, the species transport equations, and chemical reaction equations (Eq. (2.2) ) are solved directy, as the name implies, without applying any further turbulence model [37]. As the main differences between the models can be found in the continuity and the Navier Stokes equations, they will be the focus of this re- view.¹

Continuity Equation:

∂ρ

∂t +ρ∂u_i

∂x_i =0 (2.1)

with the densityρ, the timet, the velocity and the spatial coordinate in i directionui andxi

respectively.

1The equations in section 2.1 are displayed in the same representation as by Davidson [5] unless cited other- wise. Some variables were changed for simplicity and unity throughout this work without further explanation.

Navier-Stokes Equations for Newtonian Fluids:

ρ∂u_i

∂t = −∂p

∂xi

+ ∂

∂xj

µ

2ηS_{i j}−2 3η∂u_k

∂xk

κi j

¶

+ρf_i (2.2)

with the pressurep, the velocity in i directionui, the dynamic viscosityη, the strain rate ten- sorS_{i j} the volume force f_i =(0, 0,g) with the gravitational constant g and Kronecker delta functionκ, which is defined as

κi j=

(1 i=j

0 i6=j . (2.3)

The assumption of incompressible flow is made often since the Eq. (2.1) und (2.2) simplify significantly to:

Continuity Equation:

∂u_i

∂x_i =0 (2.4)

Navier Stokes Equation:

ρ∂ui

∂t = −∂p

∂x_i +η ∂²u_i

∂x_j∂x_j +ρf_i (2.5)

2.1.2 Large Eddy Simulation

The second method of modeling turbulent flow is the Large Eddy Simulation (LES) approach where only the large scale contributions are considered. The method applies a spacial filter on the instantaneous equations such that large eddies are kept, whereas the effects at smaller scales need to be modeled [8, 37]. There are several different techniques for filtering, but pro- grams operating with the finite volume method (such as the Fluent software) often apply a box filter in the form of Eq. (2.6):

θ(x,ˆ t)= 1

∆V Z

CVθ( ˜x,t)dV (2.6)

whereθis any flow variable,xis the spatial vector,t is the time,V is the volume, andCV is the defined volume on which the filter is applied. The circumflex accent represents LES filter- ing. The filtering length∆f in this case is defined asp³

∆V and is therefore equaling the length scale of each control volume [37]. Figure 2.1 shows graphically how the filtering method con- verts a fully turbulent flow field into a flow field containing only the large eddies.

As only filtered Continuity and Navier-Stokes equations are solved, the computational costs are reduced. A further decrease in computational effort is the third method termed as Reynolds- Averaged-Navier-Stokes (RANS) Simulation.

Figure 2.1:Display of LES Filtering adapted from Schwarze [37]: a) Flow field with large and small scale eddies; b) Flow field after LES filtering showing only the large scale eddies

2.1.3 Reynolds Averaged Navier Stokes

With the Reynolds-Averaged-Navier-Stokes approach, the equations are simplified as the fluc- tuating components are subdivided into the mean and the fluctuating component. There are several different types of averaging (time, spatial, ensemble), but mostly time averaging is used [37, 43]. In formulas, time averaging can be described as [37]:

θ(x)= lim

∆t→∞

1

∆t Z t+∆t

t θ(x, ˜t)d˜t (2.7) whereθis any flow variable,xis the spatial vector, andtis the time. The overline represents the time averaged value asθremains only a function of the location.

The continuity and Navier-Stokes equation as presented in Eq. (2.1) and Eq. (2.2) simplify significantly to:

Continuity Equation:

∂u_i

∂x_i =0 (2.8)

Navier Stokes Equation:

ρ∂u_iu_j

∂x_j = −∂p

∂x_i + ∂

∂x_j µ

η∂u_i

∂x_j −ρu⁰_iu⁰_j

¶

(2.9) Figure 2.2 visualizes the effect of Eq. (2.7). It shows how the RANS appraoch removes all large and small scale eddies leaving only the mean flow field: All field variables are averaged over an infinite amount of time removing all unsteady turbulent features. Therefore, closure models are required to model all time-dependent terms. [8]. Different closure approaches are pre- sented later in this work in Section 2.2.

Figure 2.2:Display of Reynolds Averaging: a) Turbulent flow field; b) Flow field after Reynolds averaging showing the averaged streamlines

Turbulence-Chemistry Interaction

Since for RANS and LES modeling the chemistry and the turbulence terms are unclosed, the interaction between the chemistry and the turbulence needs to be modeled. In literature, many different methods can be found including a one equation chemistry model for sim- ple combustion chemistry, models defining whether the chemistry is dominant (Arrhenius) or the turbulence is dominant (Eddy Break Up Model), statistical methods (Bray Moss Libby model), and flame surface density models [41]. This section covers in more detail the promi- nent Eddy Break Up model and the flame surface density model.

The Eddy Break Up model (EBU) assumes that the chemical reaction rate is only controlled by the turbulence. In formulas, the mean reaction rate is described as [41]:

r=C_{E BU}·ρ q

T_r⁰² tE BU

(2.10) with the model constantC_{E BU}, the densityρ, the turbulence time as estimated from the ki- netic energykand the dissipation rateε

t_{E BU} =k

ε, (2.11)

and the reduced temperatureT_r

T_r=C_p(T−T₁)

QY_F¹ (2.12)

with the specific heat capacityCp, the specific heat release of the fuelQ, and the initial mass fraction of fuel in the unburnt streamY_F¹.

The turbulence chemistry interaction is implemented into the set of equation as:

∂ρTr

∂t + ∂

∂x_i(ρTrui)= ∂

∂x_i µ

ρD∂Tr

∂x_i −ρu⁰_iTr0¶

+r (2.13)

Another possibility is to account for the chemistry turbulence interaction is the flame surface density model [41]. This approach describes the mean reaction rate by the flame consumption rate in correlation with the flame surface area. The basic assumption behind the model is if the flame density, meaning the surface area per volume, is large, the reaction rate is high.

Therefore, the mean reaction rate can be described as:

r=ρ0〈s_c〉sΣ (2.14)

with the density of the unburnt streamρ0, the average flame consumption speed along the surface〈sc〉s, and the flame surface densityΣ. The implementation in the equation set follows as with the EBU model Eq. (2.13).

2.1.4 Probability Density Function

The fourth approach uses a statistical method of modeling transport parameters by applying the Probability Density Function. It offers many advantages for chemically reacting flow fields as the chemical terms in the transport equation have not been subject to averaging or filtering and are thus closed. Therefore, it is especially suited for cases where the chemical reactions and molecular transport are tightly coupled and it is often used for finite-rate chemistry ef- fects in turbulent flames. PDF transport has been applied on premixed, non-premixed, and partially premixed flames [13].

One of the first and most influential persons in the development of the PDF transport ap- proach was Pope in the late 1980s with the most popular work from 1985 [35]. Pope intro- duced the relationship between particle models and PDF methods. Since then, many groups found different versions of applying the probability density function. In the following, the single-point, joint composition PDF approach as implemented in Fluent is further discussed.

The PDF has N+2 dimensions: one for each of the N species, the pressure and the tempera- ture. The transport equation for the PDF is presented in Eq. (2.15) [8, 15]:

∂(ρP)

∂t +∂(ρu_iP)

∂x_i +∂(ρr_kP)

∂ψk

= − ∂

∂x_i

hρ〈u_i⁰|ψ〉i + ∂

∂ψi

· ρ

¿1 ρ

∂J_i,k

∂x_i

¯

¯

¯

¯ψ À

P

¸

(2.15) where ρ is the density, P is the Favre joint PDF of composition,r_k is the reaction rate for specie k,ψis the composition space vector,u_i⁰ is the fluid velocity fluctuation, andJ_i,kis the molecular diffusion flux vector. The notation〈A|B〉signifies the conditional probability of A given that event B occurs.

The three terms on the left-hand side of equation describing the change of the PDF due to the unsteady rate, the mean velocity field, and to chemical reactions respectively are closed, whereas the two terms on the right-hand side describing the turbulent scalar flux require a closure model. In Fluent the gradient-diffusion assumption is made:

− ∂

∂x_i

hρ〈u⁰_i|ψ〉Pi

= ∂

∂x_i µρηt

Sc_t

∂P

∂x_i

¶

(2.16)

with Sct as the turbulent Schmidt number. A turbulence model as described later in Section 2.2 is needed to model the turbulent viscosityηt.

Since only a single-point PDF is solved, the interaction with the neighboring particles needs to be accounted for separately. This is done with either the Lagrangian particle or the Eulerian mesh method. The Lagrangian method requires more computational effort as it stochastically tracks Lagrangian particles through the domain. In constrast, the Eulerian method assumes a shape for the PDF in the form of a product of delta functions. This approach is shown in Eq.

(2.17):

P(ψ;x,t)=

Ne

X

n=1

p_n(x,t)

Ns

Y

m=1

δ£

ψm− 〈θm〉n(x,t)¤

(2.17) wherep_n is the probability of each mode,<θm >n is the conditional mean composition of speciemin then^thmode,ψmis the composition space variable of speciem, andδis the delta function.

2.2 Turbulent Closure Models

Viscous models link the turbulent viscosity to turbulent parameters and they offer a solution to the closure problems for RANS and PDF calculations. Different methods with a varying amount of equations ranging from one to 7 equations are used. Two-equation models have the advantage of computational simplicity, however, they have problems with predicting the onset and amount of separation in adverse pressure gradient flows [28]. This part covers sev- eral two-equation models such as the Shear-Stress-Transport (SST)-Model. Since they are in- corporated in the SST-Model, thek-ε- and thek-ω-Model are discussed in detail. The promi- nent realizablek-ε-Model is explained as an example on how the standard ideas have been adapted and optimized for certain applications. Finally, the implementation and the advan- tages of the SST-Model approach are discussed.²

Standard k-ε-Model

The standard k-ε-Model is one of the most widely used turbulent closure models. It includes two equations (Eq. (2.18) and 2.19) for the turbulent kinetic energykand the dissipation rate ε. As it is semi-empirical, it does not use the exact transport equations forε, but, for sim- plicity, assumes a proportial behavior of several terms [15, 37] (Eq. (2.19)). It offers a robust and reasonably accurate solution with low computational costs. As it has been applied many times, its strengths and weaknesses are well-known. It simulates reasonably simply flow fields such as the bulk of fully developed turbulent flow and is therefore often used when only the global flow field of a technical process has to be considered. But it struggles modeling bent

2The equations for all viscous models are taken from the Fluent documentation. Some equations were simpli- fied by directly inserting constants or simple terms without further explanation. For more detail, consult source [15]

and twisted streamlines and overestimates the viscosity in near wall regions due to the pres- sure rise [24, 37]. As a result, several improved models like the realizable k-ε-Model, which is described later in this chapter, were proposed.

The formulas for the standard k-ε-Model are introduced in the following. Equation (2.18) and (2.19) are the two major governing equations:

∂(ρk)

∂t +∂(ρkui)

∂x_i = ∂

∂x_j

h³η+ηt

´∂k

∂x_j i

+G_k+G_b−ρε−Y_M+B_k (2.18)

∂(ρε)

∂t +∂(ρεui)

∂x_i = ∂

∂x_j

h³η+ ηt

1.3

´ ∂ε

∂x_j i

+C_1εε

k(G_k+C_3εG_b)−C_2ερε²

k +B_ε (2.19) whereρis the density,ui is the velocity in i direction, xi is the spatial coordinate, ηt is the turbulent dynamic viscosity,G_kis the generation of turbulent kinetic energy due to mean ve- locity gradients (Eq. (2.20)),G_bis the generation of turbulent kinetic energy due to buoyancy (Eq. (2.21)),YM is the contribution of fluctuating dilation in compressible turbulence to the overall dissipation rate (Eq. (2.23)),C_1ε,C_2ε, andC_3ε are constants, andB_k andB_ε are user- defined source terms.

The generation ofkandεcan be described as:

G_k= −ρu_iu_j∂u_j

∂x_i (2.20)

G_b=βg_i η Prt

T xi

(2.21) where Pr_t is the turbulent Prandtl number andg_iis the component of the gravitational vector in the i^thdirection, andβis the thermal expansion as defined in Eq. (2.22).

β= −1 ρ

³∂ρ

∂T

´

p (2.22)

Y_M=2ρεM_t² (2.23)

with the turbulent Mach number

M_t = s

k

a² (2.24)

whereais the speed of sound.

The turbulent viscosity is calculated as:

ηt=ρC_ηk²

ε (2.25)

withC_ηbeing a constant.

Realizablek-ε-Model

The realizablek-ε-Model deviates from the standardk-ε-Model in two important points: It contains a different approach for calculating the turbulent viscosity and the equation for the dissipation rateεis derived from an exact equation for the transport of the mean-square vor- ticity fluctuation [15].

Eq. (2.25) for the turbulent viscosity remains valid, however, with the realizable k-ε-Model the termC_ηdoes not stay constant.

C_η= µ

4.04+p

6 cosζk· q

S_{i j}S_{i j}+Ψ˜i jΨ˜i j

ε

¶₋₁

(2.26) with ˜Ψi j being

Ψ˜i j=Ψi j−2εi j kγk (2.27)

Ψi j=Ψi j−εi j kγk (2.28)

whereΨi j is the mean rate of rotation tensor viewed a rotating reference frame with the an- gular velocityγk,ζdefined as:

ζ=1 3cos⁻¹

µp

6 Si jS_{j k}S_ki

¡pS_{i j}S_{i j}¢3

¶

, (2.29)

andSi j written out as:

S_{i j}=1 2

µ∂u_j

∂x_i +∂ui

∂x_j

¶

. (2.30)

Equation 2.19 for the dissipation rateεis exchanged by:

∂(ρε)

∂t +∂(ρεuj)

∂xi

= ∂

∂xj

h³η+ ηt

1.2

´ ∂ε

∂xj

i

+ρC₁S_ε−ρC₂ ε² k+q_η

ρε+C_1eε

kC_3εG_b+B_ε (2.31) whereρ is the density,u_i is the velocity in i direction, x_i is the spatial coordinate, ηis the dynamic viscosity with the indextrepresenting the turbulent parameter,C_1ε,C_2ε, andC_3εare constants, andB_εis a user-defined source term.

C₁=max

·

0.43, S^k_ε S^k_ε+5

¸

(2.32)

Standardk-ω-Model

The advantage of thek-ω-Model is often found in the improved performance in the near wall region where it is more accurate and robust than the k-ε-Model. The two equations model the transport of the turbulent kinetic energykand the specific dissipation rateω. The specific dissipation rate can be thought of as the ratio ofεto k.

∂(ρk)

∂t +∂(ρku_i)

∂xi = ∂

∂xj

·µ η+ηt

2

¶∂k

∂xj

¸

+G_k−Y_k+S_k (2.33)

∂(ρω)

∂t +∂ρωui

∂x_i = ∂

∂x_j

·µ η+ηt

2

¶∂ω

∂x_j

¸

+G_ω−Y_ω+S_ω (2.34) The turbulent viscosity is modeled as:

ηt=ρk

ω (2.35)

The production of turbulent kinetic energy remains modeled as in Eq. (2.20) while the pro- duction of the specific dissipation rate is implemented as decribed in Eq. (2.36)

G_ω=αω

kG_k (2.36)

α=α_∞ α^∗

µ¹

9+_2.95ηω^ρk 1+_2.95ηω^ρk

¶

(2.37) withα_∞=0.52 andα^∗equaling

α^∗=

βi

3 +_6ηω^ρk 1+_6ηω^ρk

(2.38) with the thermal expansion coefficientβi=0.072.

The dissipation of the turbulent kinetic energy k is described as:

Y_k=ρβ^∗f_β^∗kω (2.39)

f_β∗=











1 ¡₁

ω ∂k

∂xj

∂x∂ωj

¢≤0

1+680¡₁

ω ∂k

∂x j ∂ω

∂x j

¢2

1+400¡₁

ω ∂k

∂x j ∂ω

∂x j

¢2

¡₁

ω ∂k

∂xj

∂x∂ωj

¢>0 (2.40)

and

β^∗=0.09

µ4/15+¡_ρk 8ηω

¢4

1+¡ _ρk 8ηω

¢4

¶

£1+1.5F(2M_t)¤

(2.41)

with

F(2Mt)=

(0 2M_t ≤0.25

2M_t²−0.25² 2M_t >0.25 (2.42) whereM_t is the turbulent Mach number as described in Eq. (2.24).

The dissipation ofωis given by:

Y_ω=ρβf_βω² (2.43)

with

β=βi

µ 1−

0.09

µ4/15+¡_ρk

8ηω

¢4

1+¡_ρk

8ηω

¢4

¶

0.072 1.5F(2M_t)

¶

(2.44) withβi=0.072 and

f_β= 1+70

¯

¯

¯

1 4

³∂ui

∂x j−^∂_∂^{u j}_xi

´³_∂

u j

∂xk−^∂_∂^uk_{x j}

´

S_ki (0.09ω)³

¯

¯

¯

1+80

¯

¯

¯

1 4

³∂ui

∂x j−^∂_∂^{u j}_xi

´³_∂

u j

∂xk−^∂_∂^uk_{x j}

´

Ski

(0.09ω)³

¯

¯

¯

. (2.45)

SST-k-ω-Model

The SST-k-ω-Model was first introduced by Menter et al. in 1994 [15, 28]. It uses the Wilcox k- ω-Model as a baseline but in the freestream region switches to the k-ε-Model. Mathematically, the implementation of the two models is realized by using two blending functions F₁ and F₂, which are "1" when in the near-wall region activating the k-ω-Model and "0" in the free- stream region activation the k-ε-Model. However, the model differs from the standard k-ω/ε model: the SST incorporates a damped cross-diffusion derivative term in theωequation, the turbulent shear stress transport is accounted for in the definition of the turbulent viscosity, and the modeling constants are different [15].

The two blending functions are defined as:

F₁=tanh³ minh

max¡ pk

0.09ωy, 500η

ρy²ω, 4ρk σ_ω,2D_ω⁺y²

¢i´

(2.46) with the positive portion of the cross-diffusion termD_ω⁺

D⁺_ω=maxh 2ρ 1

σω,2

1 ω

∂k

∂x_j

∂ω

∂x_j, 10⁻¹⁰i

(2.47) and

F₂=tanh µ³

max h

2 pk

0.09ωy, 500η ρy²ω

i´2¶

. (2.48)

The transport equations are adapted and incorporate blending functionF₁:

∂(ρk)

∂t +∂(ρku_i)

∂xi

= ∂

∂xj

·µ

η+ ηt

¡_F1

2 +^(1−F_1.168¹⁾¢₋₁

¶∂k

∂xj

¸

+G˜_k−Y_k+S_k (2.49)

∂(ρω)

∂t +∂ρωu_i

∂xi = ∂

∂xj

·µ

η+ ηt

¡ _F1

1.176+(1−F₁)¢₋₁

¶∂ω

∂xj

¸

+G_ω−Y_ω+S_ω (2.50) with

G˜_k=min(G_k, 10ρβ^∗kω) (2.51)

withG_kas in Eq. (2.20).

The generation ofωis modeled as:

G_ω=ρα ηt

G˜_k (2.52)

with alpha as in 2.37, except withoutα_∞being a constant:

α∞=F₁α∞,1+(a−F₁)α∞,2 (2.53) and

α_∞,1=0.075

0.09 − 0.41² 2p

0.09 (2.54)

α_∞,2=0.0828

0.09 − 0.41² 1.168p

0.09 (2.55)

As already mentioned, the SST-k−ωModel incorporates the cross-diffusion. Its positive por- tion has already been explained earlier as part of one of the blending functions. For complete- ness, the overall cross-diffusion is defined as:

D_ω=2(1−F₁)ρσ_ω,2

1 ω

∂k

∂xj

∂ω

∂xj

(2.56) The turbulent viscosity is modeled as:

ηt =ρk ω

³ max£

(α^∗)⁻¹, SF₂ 0.31ω

¤´−1

(2.57) withα^∗as given in Eq.(2.38)

The dissipation ofkandωare calculated similar to the standardk−ωmodel:

Y_k=ρβ^∗kω (2.58)

Y_ω=ρβω². (2.59)

Additionally,βi is defined differently:

βi=0.075F₁+0.0828(1−F₁) (2.60)

2.3 Modeling of Combustion Chemistry

Modeling combustion chemistry is a complex phenomenon. First, it has to be decided how many species and which reaction pathways need to be considered. Then, all the thermody- namic and kinetic informaton, such as the heat of reaction and the equilibria data have to be implemented. Additionally, the coupling between the turbulence and the combustion has to be considered [41].

One of the most well accepted mechanisms for methane or natural gas as fuel is the GRI- Mech 3.0. This current release includes 53 species undergoing 325 elementary chemical re- actions. In addition, it contains information about the associated rate coefficient expressions and the thermodynamic parameters [12]. The mechanism has been experimentally validated on multiple occasions.

Due to the complexity of the GRI-Mech 3.0 mechanism, many simpler versions have been examined. As an example Bibrzycki et al. investigated the effect of a simpler reaction mech- anism on the laminar flame speed and found a good agreement for the two step chemistry (2S-CMS) with the more complex mechanisms for lean combustion [16]. The reaction mech- anism can be expressed as follows:

1^stReaction :C H₄+1.5O₂→CO₂+2H₂O (2.61)

2^ndReaction :CO+0.5O₂ *)CO₂ (2.62)

Methane (C H4) reacts irreversibly with oxygen (O2) from the air to form carbon monoxide (CO) and water (H₂O). The formedCOcan then further oxidize to carbon dioxide (CO₂) in a reversible reaction. The reaction order for the first step is 0.9 with respect to methane and 1.1 with respect to oxygen. The reaction order for the forward reaction of the second step is 1.0 and 0.5 forCOandO₂respectively. For the implementation of the Arrhenius equation as depicted in Eq. (2.63), the activation energyE_aof 35000 _mol^cal and 12000 _mol^cal for both reactions respectively was determined. Factorb was 1 in both cases and the pre-exponential factors were 2.0E+15 mol-cm-s-K and 2.0E+09 mol-cm-s-K respectively [30]. The implementation of the reaction file in Chemkin format is shown in A.1 in the appendix.

k=AT^bexp µ

− E_a RT

¶

(2.63)

In addition to the reaction mechanism, thermodynamic and transport information could be added to the Fluent software. Thermodynamic values were given as constants for NASA poly- nomials as a function of temperature. The NASA polynomials for the specific heat capacity (CP), the enthalpy (H) and the entropy (S) are defined with seven constants (a₁−a₇) as fol-

lows:

Cp

R =a₁+a₂T+a₃T²+a₄T³+a₅T⁴ (2.64) H

R =a₁+a₂T

2 +a₃T²

3 +a₄T³

4 +a₅T⁴ 5 +a₆

T (2.65)

S

R =a₁lnT+a₂T+a₃T²

2 +a₄T³

3 +a₅T⁴

4 +a₇ (2.66)

Premixed Flames

This chapter covers the basics of the theory behind system identification and how it is applied to find the flame transfer function and the entropy transfer function. The method of modeling the advection of fluctuations with a state space model and accounting for the combustion chemistry with the static gain is described in more detail. At the end, the determination of the diffusion coefficient is discussed.

3.1 System Identification

The principles developed by the field of system identification can be applied to combustion problems as transfer functions are developed. The fundamental concept relies on reducing the complex system to a black box with inlet and outlet variables. An example for such a black box system is given in Figure 3.1. Mean normalized fluctuations in velocity and equivalence ratio are input variables, whereas the output variables are mean normalized fluctuations in the heat release ˙Q and in temperatureT. The equivalence ratio is a measure of the fuel to air amount and is described in more detail later in this work. Both inlet pertubations can cause pertubations in the heat release connected by the flame transfer functionsF_u andF_φ. Fluc- tuations in the temperature originate from equivalence ratio pertubations and the transport and conversion is described by the entropy transfer functionE_φ.

In the following, the basic principles of system identification are summarized. To model the system, it is important to abstract all effects in the correct input and output variables. Fig- ure 3.1 gives one example of an multiple-input-multiple-output system. Even though there are also simpler methods with a single input and/or a single outlet or combination of sin- gle/multiple inlets/outlets [14, 44], modeling the flame as a MIMO system proved to be essen- tial when analyzing the entropy transfer function as the response of the two output variables is significantly different [38].

A crucial requirement for a simple and useful application of system identification techniques is a linear and time-invariant (LTI) system , which can be assumed for small pertubations [14]. The linearity includes the amplification and the superposition principle [1]. Given a out- put signaly(t) is connected via an operatorO with a input signalx(t), both being a function of the timet,

y(t)=O{x(t)} (3.1)

Figure 3.1:Scheme of the flame response as an MIMO system taken from Steinbacher et al.

[38]; the FTF due to velocity fluctuations and the FTF and ETF for equivalence ratio fluctuations are included; the FTFs result in changes in the heat release while the ETF gives temperature fluctuations as an output

the linearity (Eq. (3.2)) and the superposition (Eq. (3.3)) can be described as:

O{aw(t)}=aO{w(t)} (3.2)

O{x₁(t)+x₂(t)}=O{x₁(t)}+O{x₂(t)}. (3.3) Time-invariance means that a shift in the input signal byτ results in a shift in the output signal:

y(t−τ)=O{x(t−τ)} (3.4)

The advantage of a LTI system is that its complete response behavior can be described by its answer to the impulse response, which can be mathematically proven in the following as described in by Girod et al. [1]:

We assume again the responsey(t) to any given inputx(t). The impulse function, also known as the delta distribution, is defined as:

δ(t)=

(0 t6=0

∞ t=0 (3.5)

Since the delta function is 0 for all other values other than its argument (t), it can be used to describe any input functionx(t):

x(t)= Z _∞

−∞x(τ)δ(t−τ)dτ. (3.6)

The response tox(t) is:

y(t)=O

½Z _∞

−∞x(τ)δ(t−τ)dτ

¾

. (3.7)

Considering the linearity of the system leads to:

y(t)= Z _∞

−∞

x(τ)O{δ(t−τ)}dτ. (3.8)

This expression can be simplified by using the impulse responseo(t)=O{δ(t)}

y(t)= Z _∞

−∞

x(τ)o(t−τ)dτ (3.9)

and finally the definition of the convolution:

y(t)=x(t)∗o(t) (3.10)

with∗being the convolution operator. Eq. (3.10) shows directly how the response to any input signal is connected to the input signal itself and the response to the delta function.

Often also the step-function instead of the delta function is used as they are related by in- tegration/differentiation and therefore result in the same information. The advantage of the step function is that, depending on the system, it is often easier to practically implement.

To convert the input signal from the time to the frequency domain, either the Laplace or the Fourier transformation needs to be applied [1]. The response of a system is often summarized in a Bode diagram where the magnitude and the phase are plotted over the frequency.

3.2 Entropy Transfer Function Modeling

As described previously, system identification can be used model the response of pertuba- tions in the composition of the stream at the inlet. A common variable being used to describe the composition of the unburnt stream is the equivalence ratioφ. It is defined as the mass ratio of fuel to air compared to the stoichiometric ratio of fuel to air [15]:

φ=

mass of fuel mass of air

µ

mass of fuel mass of air

¶

st.

. (3.11)

With this definition, combustion systems with an equivalence ratio below 1 are considered as lean combustion.

Steinbacher et al. have developed a concept to simplify the modeling process of the transport of equivalence ratio fluctuations upstream of the flame to temperature fluctuations down- stream of the flame [38]. They could break down the complex processes to (1) the pure advec- tion of equivalence ratio pertubations from the combustor inlet over the flame to the outlet and (2) the combustion process. Due to the linearity of the problem, (1) and (2) can be applied in succession and the answers can be superpositioned. The advection of the equivalence ratio pertubation were modeled with a state space model and the combustion was implemented using a Fast-Reaction-Zone (FRZ) approach. Figure 3.2 visualizes in a one dimensional (1D)

Figure 3.2:Visualization of the changed order of transport/transfer; instead of modeling the transport through the unburt area, the conversion of reactants, and finally the transport through the burnt mixture, the transport through the whole domain is modeled and only at the end the FRZ is taken into account

approach the modification of the problem. In order the equivalence ratio pertubation are ad- vected through the unburnt mixture (index "u") to the flame front. After the length L_f the FRZ is reached and the reactants are chemically converted. Then the temperature fluctua- tions are convected downstream through the burnt mixture (index "b"). Due to the linearity of the problem, it is, however, a valid operation to change the order of the individual transfer functions. As Figure 3.2 demonstrates in green, the FRZ can be moved to the end of the prob- lem. Therefore, first the advection from the inlet to the outlet is modeled and afterwards the transfer function to convert from equivalence ratio fluctuations to temperature fluctuations is modeled.

3.2.1 State Space Model

The state space model simulates the advection of the equivalence ratio fluctuations through the combustor. It uses the steady state flow field data and the transport equation given in equation 3.12.

ρ∂³

ξ⁰ ξ

´

∂t +ρui

∂³

ξ⁰ ξ

´

∂xi

= ∂

∂xi

µ ρD∂³

ξ⁰ ξ

´

∂xi

¶

(3.12) Instead of using the equivalence ratio directly, the mixture fractionξ is being used as it al- lows the decoupling of advection and chemistry [38]. In contrast to the equivalence ratio, the mixture fraction is a conserved quantity over the flame front. The mixture fraction is defined as:

ξ= 1

1+^l_φ^st (3.13)

withlst =17.13 as the stoichiometric ratio from air to fuel (on mass basis) for methane-air flames.

For small fluctuations, the fluctuation of the mixture fraction is linearly connected to the equivalence ratio fluctuation:

φ⁰ φ = 1

1−ξ ξ⁰

ξ. (3.14)

This approach has been proven as a good and simple linear approximation instead of the elaborated transient non-linear calculation [38].

3.2.2 Fast-Reaction-Zone Model

The conversion at the flame front, which is assumed to be infinitely thin, from equivalence ratio fluctuations to temperature fluctuations due to chemical reactions is modeled through Eq. (3.15):

T_b⁰

T_b =νF R Z(φ_u)φ⁰_u

φ_u (3.15)

The static gainνF R Zis the value of the transfer function under steady conditions [17]. During the combustion process, the kinetics of the chemical reactions need to be taken into account.

The value for the static gain were taken from Steinbacher et al. [38], where they used the Can- terra software package [4] with the full GRI-Mech 3.0 mechanism [12]. They applied Eq. (3.16) with a spacing of∆φ=0.01 and plotted the curve as a function of the equivalence ratio in the red curve in Figure 3.3. They explain the constant static gain for an equvalence ratio between 0.6 and 0.8 as the achievement of equilibrium. Once the equivalence ratio approaches stoi- chiometric conditions, the temperature increases, which changes the kinetics of the chemical processes so a complete combustion cannot be reached. The fluctuations of the equivalence ratio and the temperature are in anti-phase when the static gain becomes negative since the combustion is incomplete and the specific heat capacity of the burnt mixture increases due to unburnt gases.

νF R Z=φ_u

T_b·T_b(φ_u+∆φ)−T_b(φ_u−∆φ)

2∆φ (3.16)

3.2.3 Implementation in Fluent

Via the User Defined Scalar (UDS)-Function in Fluent, it is possible to implement the analysis via a step function with the state space model. The transport of any arbitrary scalarϑcan be calculated with the following transport equation [15]:

∂ρϑ

∂t + ∂

∂x_i µ

ρuiϑ−D∂ϑ

∂x_i

¶

=S_ϑ (3.17)

Figure 3.3:Static gain as a function of the mean equivalence ratio of the unburnt mixtureφ_u by using the Canterra Software Package with the Gri-Mech 3.0 reaction mecha- nism (red) as adapted from Steinbacher et al. [38]; the other curves describe dif- ferent models or assumptions: for further information consult source [38]

with the densityρ, the velocity ini directionu_i, the spatial coordinatex_i, the diffusion coef- ficentD, and the source termS_ϑ.

The diffusion coefficient is calculation as in Section 3.2.4 and included in Fluent as a User Defined Function. The source term is set to zero.

The application of the step function was realized by setting the UDS inlet boundary condition to a constant value of 1. Unless stated otherwise all other boundary conditions were set to a flux of 0. All equations except the UDS transport equation are deactivated to only use the steady results of the flow field and limit the computational effort. The resulting scalar at the outlet is mass flux averaged:

SR(t)=

Rρuiϑ(t)dA

RρuidA (3.18)

with SR representing the step response as a function of the timet.

3.2.4 Diffusion Coefficient

The diffusion coefficient can contain laminar and turbulent contributions:

D=D_{l am}+Dt (3.19)

Laminar: Chapman-Enskog

The laminar diffusion coefficient for a gas mixture was modeled according to the Chapman- Enskog approach as described by Maier [27]. The method is based on the following equation [3, 27, 29]:

D_{l am}=1, 8583·10⁻³· T³²q

1 M1+_M¹₂ pσ²₁₂Ω12

(3.20) where T is the temperature, M_i is the molar mass of component i, p is the pressure, σ12 is the Lennard-Jones parameter from components 1 and 2, andΩis the collision integral. The indexes 1 or 2 label the parameter for component 1 or component 2, whereas the index 12 labels mixture parameter.

The Lennard-Jones parametersσ12and²12of a two component gas can be calculated from its containing gases

σ12=1

2(σ1+σ2) (3.21)

²12=p

²1+²2 (3.22)

Table 3.1:Lennard-Jones Parameters [3, 27]

σ[ ˚A] _k^²

B [K]

CH₄ 3.758 148.6 Air 3.711 78.6

The parameters are given in table 3.1. With this information and the Neufeld parameter from table 3.2, the collision integralΩ12can be calculated:

Ω12= A

T_r^B + C

exp(D·T_r)+ E

exp(F·T_r)+ G

exp(H·T_r) (3.23)

using the reduced temperature

T_r=k_BT

²12

. (3.24)