Tutorial 01 Rijke Tube

A Rijke tube is the most simple thermo-acoustic configuration. It comprises a tube with an unsteady heat source somewhere inside the tube. This example will walk you through the basic steps of setting up a Rijke tube in a Helmholtz-solver stability analysis.

Note

The tutorial reproduces the validation case found in [1]. To run it yourself you will need the "Rijke_mm.msh" file. Download it here.

The Helmholtz equation

The Helmholtz equation is the Fourier transform of the acoustic wave equation. Given that there is a heat source, it reads:

\[\nabla\cdot(c^2\nabla \hat p) + \omega^2\hat p = - i \omega (\gamma -1)\hat q\]

Here $c$ is speed-of-sound-field, $\hat p$ the (complex) pressure fluctuation amplitude $\omega$ the (complex) frequency of the problem, $i$ the imaginary unit, $\gamma$ the ratio of specific heats, and $\hat q$ the amplitude of the heat release fluctuation.

(Note that the Fourier transform here follows a $f'(t) \rightarrow \hat f(\omega)\exp(+i\omega t)$ convention.)

Together with some boundary conditions the Helmholtz equation models thermo-acoustic resonance in a cavity. The minimum required inputs to specify a problem are

  1. the shape of the domain
  2. the speed-of-sound field
  3. the boundary conditions

In case of an active flame ($\hat q \neq 0$). We will also need

  1. some additional gas properties and
  2. a flame response model linking the heat release fluctuations $\hat q$ to the pressure fluctuations $\hat p$

Once you completed this tutorial you will know the basic steps of how to conduct a thermo-acoustic stability analysis

First you will need to load the Helmholtz solver. The following line brings all necessary tools into scope:

using WavesAndEigenvalues.Helmholtz

Mesh

Now we can load the mesh file. It is the essential piece of information defining the domain shape . This tutorial's mesh has been specified in mm which is why we scale it by 0.001 to load the coordinates as m:

mesh=Mesh("Rijke_mm.msh",scale=0.001)
mesh: Rijke_mm.msh
#################
points:     1006
lines:      0
triangles:  1562
tetrahedra: 3380
#################
domains: Cold, Flame, Flame_in, Flame_out, Hot, Inlet, Interior, Outlet, Walls

The displayed info tells us that the mesh features 1006 points as vertices 1562 (addressable) triangles on its surface and 3380 tetrahedra forming the tube. Note, that there are currently no lines stored. This is because line information is only needed when using second-order finite elements. To save memory this information is only assembled and stored when explicitly requested. We will come back to this aspect in a later tutorial.

You can also see that the mesh features several named domains such as "Interior","Flame", "Inlet" and "Outlet". These domains are used to specify certain conditions for our model.

A model descriptor is always given as a dictionary featureng domain names as keys and tuples as values. The first tuple entry is a Julia symbol specifying the operator to be build on the associated domain, while the second entry is again a tuple holding information that is specific to the chosen operator.

Model set-up

For instance, we certainly want to build the wave operator on the entire mesh. So first, we initialize an empty dictionary

dscrp=Dict()
Dict{Any,Any}()

And then specify that the wave operator should be build everywhere. For the given mesh the domain-specifier to address the entire resonant cavity is "Interior" and in general the symbol to create the wave operator is :interior. It requires no options so the options tuple remains empty.

dscrp["Interior"]=(:interior, ())
(:interior, ())

Boundary Conditions

The tube should have a pressure node at its outlet, in order to model an open-end boundary condition. Boundary conditions are specified in terms of admittance values. A pressure note corresponds to an infinite admittance. To represent this on a computer we just give it a crazily high value like 1E15. We will also need to specify a variable name for our admittance value. This will allow to quickly change the value after discretization of the model by addressing it by this very name. This feature is one of the core concepts of the solver. As will be demonstrated later.

The complete description of our boundary condition reads

dscrp["Outlet"]=(:admittance, (:Y,1E15))
(:admittance, (:Y, 1.0e15))

You may wonder why we do not specify conditions for the other boundaries of the model. This is because the discretization is based on Bubnov-Galerkin finite elements. All unspecified, boundaries will therefore naturally be discretized as solid walls.

Flame

The Rijke tube's main feature is a domain of heat release. In the current mesh there is a designated domain "Flame" addressing a thin volume at the center of the tube. We can use this key to specify a flame with simple n-tau-dynamics. Therefore, we first need to define some basic gas properties.

γ=1.4 #ratio of specific heats
ρ=1.225 #density at the reference location upstream to the flame in kg/m^3
Tu=300.0    #K unburnt gas temperature
Tb=1200.0    #K burnt gas temperature
P0=101325.0 # ambient pressure in Pa
A=pi*0.025^2 # cross sectional area of the tube
Q02U0=P0*(Tb/Tu-1)*A*γ/(γ-1) #the ratio of mean heat release to mean velocity Q02U0
2088.98730878897

We will also need to provide the position where the reference velocity has been taken and the direction of that velocity

x_ref=[0.0; 0.0; -0.00101]
n_ref=[0.0; 0.0; 1.00] #must be a unit vector
3-element Array{Float64,1}:
 0.0
 0.0
 1.0

And of course, we need some values for n and tau

n=0.0 #interaction index
τ=0.001 #time delay
0.001

With these values the specification of the flame reads

dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref,:n,:τ,n,τ))
(:flame, (1.4, 1.225, 2088.98730878897, [0.0, 0.0, -0.00101], [0.0, 0.0, 1.0], :n, :τ, 0.0, 0.001))

Note that the order of the values in the options tuple is important. Also note that we assign the symbols :nand :τ for later analysis.

Speed of Sound

The description of the Rijke tube model is nearly complete. We just need to specify the speed of sound field. For this example, the field is fairly simple and can be specified analytically, using the generate_field function and a custom three-parameter function.

R=287.05 # J/(kg*K) specific gas constant (air)
function speedofsound(x,y,z)
    if z<0.
        return sqrt(γ*R*Tu)#m/s
    else
        return sqrt(γ*R*Tb)#m/s
    end
end
c=generate_field(mesh,speedofsound)
3380-element Array{Float64,1}:
 694.4379021914054
 694.4379021914054
 694.4379021914054
 694.4379021914054
 694.4379021914054
 694.4379021914054
 694.4379021914054
 694.4379021914054
 694.4379021914054
 694.4379021914054
   ⋮
 347.2189510957027
 347.2189510957027
 347.2189510957027
 347.2189510957027
 347.2189510957027
 347.2189510957027
 347.2189510957027
 347.2189510957027
 347.2189510957027

Note that in more elaborate models, you may read the field c from a file containing simulation or experimental data, rather than specifying it analytically.

Model Discretization

Now we have all ingredients together and we can discretize the model.

L=discretize(mesh,dscrp,c)
1006×1006-dimensional operator family: 

ω^2*M+K+n*exp(-iωτ)*Q+ω*Y*C

Parameters
----------
n	0.0 + 0.0im
λ	Inf + 0.0im
ω	0.0 + 0.0im
τ	0.001 + 0.0im
Y	1.0e15 + 0.0im

The return value L here is a family of linear operators. The shown info is an algebraic representation of the family plus a list of associated parameters. You might have noticed that the list contains all parameters that have been specified during the model description (n,τ,Y). However, there are also two parameters that were added by default: ω and λ. ω is the complex frequency of the model and λ an auxiliary value that is important for the eigenfrequency computation.

Solving the Model

Global Solver

We can solve the model for some eigenvalues using one of the eigenvalue solvers. The package provides you with two types of eigenvalue solvers. A global contour-integration-based solver. That finds you all eigenvalues inside of a specified contour Γ in the complex plane.

Γ=[150.0+5.0im, 150.0-5.0im, 1000.0-5.0im, 1000.0+5.0im].*2*pi #corner points for the contour (in this case a rectangle)
Ω, P = beyn(L,Γ,l=5,N=256, output=true)
(Complex{Float64}[1711.4187608614363 + 8.62316956628517im, 4370.048951426738 + 25.289648986111153im], Complex{Float64}[0.0003632488102475537 - 0.01834031619600166im 6.761239033017678e-5 + 0.04355322195158859im; 0.00036325300744217205 - 0.018407015391644795im 6.785964180127594e-5 + 0.04355337777869356im; … ; 0.0003403596292166756 - 0.025501769394082243im 9.418377340769333e-5 + 0.040765369621667644im; -0.0003194854554595479 - 0.05397508026360308im 0.00020041757690195927 - 0.038672353843390336im])

The found eigenvalues are stored in the array Ω. The corresponding eigenvectors are the columns of P. The huge advantage of the global eigenvalue solver is that it finds you multiple eigenvalues. Nonetheless, its accuracy may be low and sometimes it provides you with outright wrong solutions. However, for the current case the method works as we can verify that in our search window there are two eigenmodes oscilating at 272 and 695 Hz respectively:

for ω in Ω
    println("Frequency: $(ω/2/pi)")
end
Frequency: 272.38075549130394 + 1.3724200615938802im
Frequency: 695.5148921731194 + 4.024972645198529im

Local Solver

To get high accuracy eigenvalues , there are also local iterative eigenvalue solvers. Based on an initial guess, they only find you one eigenvalue at a time but with machine-precise accuracy. One of these solvers is based on the method of successive linear problems (mslp). You call it as follows:

sol,nn,flag=mslp(L,250*2*pi,output=true);
Launching MSLP solver...
Iter   dz:     z:
----------------------------------
0		Inf	1570.7963267948965
1		144.7512883689692	1715.5476151638657 + 7.516104496360754e-13im
2		6.1067775966380395	1709.4408375672276 + 6.992327766313931e-13im
3		0.01090787460429965	1709.4299296926233 + 6.990634438094383e-13im
4		3.4808635973604396e-8	1709.4299296578147 + 6.990634427843083e-13im
5		1.3642420526593924e-12	1709.4299296578133 + 6.990634427873366e-13im
6		1.5916157281026244e-12	1709.4299296578117 + 6.990634427876007e-13im
7		2.483965218359922e-26	1709.4299296578117 + 6.990634427875759e-13im
8		2.0194839173657902e-27	1709.4299296578117 + 6.990634427875779e-13im
9		9.087677628146056e-28	1709.4299296578117 + 6.990634427875788e-13im
10		2.2737367544323206e-13	1709.4299296578115 + 6.990634427876056e-13im
Warning: Maximum number of iterations has been reached!
...finished MSLP!
#####################
 Results
#####################
Number of steps: 10
Last step parameter variation:2.2737367544323206e-13
Auxiliary eigenvalue λ residual (rhs):8.556944110533481e-10
Eigenvalue:1709.4299296578115 + 6.990634427876056e-13im

The return values of this solver are a solution object sol, the number of iterations nn performed by the solver and an error flag flag providing information on the quality of the solution. If flag==0 the solution has converged. If flag<0 something went wrong. And if flag>0... well, probably the solution is fine but you should check it.

In the current case the flag is 1 indicating that the maximum number of iterations has been reached. This is because we haven't specified a stopping criterion and the iteration just ran for the maximum number of iterations. Nevertheless, the 272 Hz mode is correct. Indeed, as you can see from the printed output it already converged to machine-precision after 5 iterations.

Changing a specified parameter.

Remember that we have specified the parameters Y,n, and τ?. We can change these easily without rediscretizing our model. For instance currently n==0.0 so the solution is purely acoustic. The effect of the flame just enters the problem via the speed of sound field (this is known as passive flame). By setting the interaction index to 1.0 we activate the flame.

L.params[:n]=1
1

Now we can just solve the model again to get the active solutions

sol_actv,nn,flag=mslp(L,245*2*pi-82im*2*pi,output=true,order=3);
Launching MSLP solver...
Iter   dz:     z:
----------------------------------
0		Inf	1539.3804002589986 - 515.221195188726im
1		1019.4775205432949	1053.6470564007975 + 381.1032625593791im
2		23.472745200482013	1075.3252077129691 + 372.1017336373935im
3		3.361460350470734e-5	1075.3252115068922 + 372.10176703720964im
4		4.190347093584501e-11	1075.3252115068517 + 372.1017670371988im
5		7.41866506207056e-12	1075.3252115068453 + 372.1017670372026im
6		1.2736097495586632e-11	1075.3252115068326 + 372.1017670372029im
7		5.514690581383049e-12	1075.325211506838 + 372.1017670372021im
8		4.351399913731402e-12	1075.325211506834 + 372.1017670372006im
9		2.3437142008433856e-13	1075.3252115068337 + 372.10176703720055im
10		1.8410599683081325e-12	1075.325211506832 + 372.10176703720026im
Warning: Maximum number of iterations has been reached!
...finished MSLP!
#####################
 Results
#####################
Number of steps: 10
Last step parameter variation:1.8410599683081325e-12
Auxiliary eigenvalue λ residual (rhs):4.1541902042462e-9
Eigenvalue:1075.325211506832 + 372.10176703720026im

The eigenfrequency is now complex:

sol_actv.params[:ω]/2/pi
171.14332284265015 + 59.221835557199306im

with a growth rate -imag(sol_actv.params[:ω]/2/pi)≈ 59.22

Instead of recomputing the eigenvalue by one of the solvers. We can also approximate the eigenvalue as a function of one of the model parameters utilizing high-order perturbation theory. For instance this gives you the 30th order diagonal Padé estimate expanded from the passive solution.

perturb_fast!(sol,L,:n,16) #compute the coefficients
freq(n)=sol(:n,n,8,8)/2/pi #create some function for convenience
freq(1) #evaluate the function at any value for `n` you are interested in.
172.47570467361967 + 59.82427381825554im

The method is slow when only computing one eigenvalue. But its computational costs are almost constant when computing multiple eigenvalues. Therefore, for larger parametric studies this is clearly the method of choice.

VTK Output for Paraview

To visualize your solutions you can store them as "*.vtu" files to open them with paraview. Just, create a dictionary, where your modes are stored as fields.

data=Dict()
data["speed_of_sound"]=c
data["abs"]=abs.(sol_actv.v)/maximum(abs.(sol_actv.v)) #normalize so that max=1
data["phase"]=angle.(sol_actv.v)

vtk_write("tutorial_01", mesh, data) # Write the solution to paraview

The vtk_write function automatically sorts your data by its interpolation scheme. Data that is constant on a tetrahedron will be written to a file *_const.vtu, data that is linearly interpolated on a tetrahedron will go to *_lin.vtu, and data that is quadratically interpolated to *_quad.vtu. The current example uses constant speed of sound on a tetrahedron and linear finite elements for the discretization of p. Therefore, two files are generated, namely tutorial_01_const.vtu containing the speed-of-sound field and tutorial_01_lin.vtu containing the mode shape. Open them with paraview and have a look!.

Summary

You learnt how to set-up a simple Rijke-tube study. This already introduced all basic steps that are needed to work with the Helmholtz solver. However, there are a lot of details you can fine tune and even features that weren't mentioned in this tutorial. The next tutorials will introduce these aspects in more detail.

References

[1] F. Nicoud, L. Benoit, and C. Sensiau, Acoustic Modes in Combustors with Complex Impedances and Multidimensional Active Flames, AIAA, 2007, doi:10.2514/1.24933

This page was generated using Literate.jl.