Tutorial 07 Bloch periodicity
The Helmholtz solver was designed with the application of thermoacoustic stability assesment in mind. A very common case is that an annular combustion chamber needs assessment. These chamber types feature a discrete rotational symmetry, i.e., they are built from one constituent entity – the unit cell – that implicitly defines the annular geometry by copying it around the circumference multiple times. This high regularity of the geometry can be used to efficiently model the problem, e.g. only storing the unit cell instead of the entire geometry in a mesh-file in order to save memory and even performing the entire analysis thereby significantly accelerating the computations [1] This tutorial teaches you how using these features.
To do this tutorial yourself you will need the "NTNU_12.msh" file. Download it here.
Model Set-up
As a geometry we choose the laboratory-scale annular combustor used at NTNU Norway [2]. The NTNU combustor features unit cells that are reflection-symmetric. The meshfile NTNU_12.msh, thus, stores only one half of the unit cell. Let's load it...
using WavesAndEigenvalues.Helmholtz
mesh=Mesh("NTNU_12.msh",scale=1.0);Note that there are two special domains named "Bloch" and "Symmetry". These are surfaces associated with the special symmetries of the combustor geometry. "Bloch" is the surface where one unit cell ends and the next one should begin, while "Symmetry" is the symmetry plane of the unit-cell. As the mesh is stored as a half-cell both surfaces are boundaries. The code will automatically use them to create the full mesh or just a unit-cell. First, let's create a full mesh We therefore specify the expansion scheme, this is a list of tuples containing all domain names we like to expand together with a qualifyer that is either :full,:unit, or :half, in order to tell the the code whether the copied versions should be merged into one addressable domain of the same name spanning the entire annulus (:full), pairs of domains from two corresponding half cells should be merged into individually adressable entities (:unit), or the copied hald cells remain individually addressable (:half). In the current case we only need individual names for each flames and can sefely merge all other entities. The scheme therefore reads:
doms=[("Interior",:full),("Inlet",:full), ("Outlet_high",:full),
("Outlet_low",:full), ("Flame",:unit),];and we can generate the mesh by
full_mesh=extend_mesh(mesh,doms);The code automatically recognizes the degree of symmetry – 12 in the current case– and expands the mesh. Note how the full mesh now has 12 consecutively numbered flames, while all other entities still have a single name even though they were copied.
Discretizing the full mesh
The generated mesh behaves just like any other mesh. We can use it to discretize the model.
speedofsound(x,y,z)=z<0.415 ? 347.0 : 850.0; #speed of sound field in m/s
C=generate_field(full_mesh,speedofsound);
dscrp=Dict(); #model description
dscrp["Interior"]=(:interior,());
dscrp["Outlet_high"]=(:admittance, (:Y_in,0));# src
dscrp["Outlet_low"]=(:admittance, (:Y_out,0));# src
L=discretize(full_mesh, dscrp, C,order=:lin)22987×22987-dimensional operator family: ω*Y_out*C+ω^2*M+K+ω*Y_in*C Parameters ---------- λ Inf + 0.0im ω 0.0 + 0.0im Y_in 0.0 + 0.0im Y_out 0.0 + 0.0im
We can use all our standard tools to solve the model. For instance, Indlekofer et al. report on a plenum-dominant 1124Hz-mode that is first order with respect to the azimuthal direction and also first order with respect to the axial direction. Let's see whether we can find it by putting 1000 Hz as an initial guess.
sol,nn,flag=mslp(L,1000,tol=1E-9,scale=2pi,output=true);Launching MSLP solver... scale: 6.283185307179586 Iter dz: z: ---------------------------------- 0 Inf 999.9999999999999 1 123.12819008931804 1123.128190089318 2 0.4821754578897787 1123.6103655472077 3 0.00010345806850245335 1123.6102620891393 4 1.9541327962911847e-11 1123.6102620891197 Solution has converged! ...finished MSLP! ##################### Results ##################### Number of steps: 4 Last step parameter variation:1.227817847393453e-10 Auxiliary eigenvalue λ residual (rhs):1.7353565060630459e-6 Eigenvalue:1123.6102620891197
There it is! But to be honest the computation is kinda slow. Obviously, this is because the full mesh is quite large. Let's see how the unit cell computation would do....
Unit cell computation
Ok, first, we create the unit cell mesh. We therefore set the keyowrd parameter unit to true in the call to extend_mesh.
unit_mesh=extend_mesh(mesh,doms,unit=true)mesh: unit from NTNU_12.msh ################# points: 2329 lines: 12444 triangles: 1484 tetrahedra: 8446 ################# domains: Flame#0, Inlet, Interior, Outlet_high, Outlet_low
Voila, there we have it! Discretization is nearly as simple as with a normal mesh. However, to invoke special Bloch-periodic boundary conditions the optional parameter b must be set to define a symbol that is used as the Bloch wavenumber. The rest is as before.
c=generate_field(unit_mesh,speedofsound);
l=discretize(unit_mesh, dscrp, c, b=:b)1933×1933-dimensional operator family: ω*Y_out*C+ω*Y_out*exp(ib2π/12)*C+ω*Y_out*exp(-ib2π/12)*C+ω*Y_out*δ(b)*C+ω*Y_out*δ(b)*exp(ib2π/12)*C+ω*Y_out*δ(b)*exp(-ib2π/12)*C+ω^2*M+ω^2*exp(ib2π/12)*M+ω^2*exp(-ib2π/12)*M+ω^2*δ(b)*M+ω^2*δ(b)*exp(ib2π/12)*M+ω^2*δ(b)*exp(-ib2π/12)*M+K+*exp(ib2π/12)*K+*exp(-ib2π/12)*K+*δ(b)*K+*δ(b)*exp(ib2π/12)*K+*δ(b)*exp(-ib2π/12)*K+ω*Y_in*C+ω*Y_in*exp(ib2π/12)*C+ω*Y_in*exp(-ib2π/12)*C+ω*Y_in*δ(b)*C+ω*Y_in*δ(b)*exp(ib2π/12)*C+ω*Y_in*δ(b)*exp(-ib2π/12)*C+(1-δ(b))*D Parameters ---------- b 0.0 + 0.0im λ Inf + 0.0im ω 0.0 + 0.0im Y_in 0.0 + 0.0im Y_out 0.0 + 0.0im
Note how the signature of the discretized operator is now much longer. These extra terms facilitate the Bloch periodicity. What is important is the parameter :b that is used to set the Bloch wavenumber. For relatively low frequencies, the Bloch wave number corresponds to the azimuthal mode order. Effectively, it acts like a filter to the eigenmodes. Setting it to 1 will give us mainly first order modes. (In general it gives us modes with azimuthal order k*12+1 where k is a natural number including 0, but 13th order modes are very high in frequency and therefore don't bother us). In comparison to the full model, this is an extra advantage of Bloch wave theory, as we cannot filter for certain mode shapes using the full mesh. This feature also allows us to reduce the estimate for eigenvalues inside the integration contour when using Beyn's algorithm, as there can only be eigenmodes corresponding to the current Bloch wavenumber. Let's try finding 1124-Hz mode again. It's of first azimuthal order so we should put the Bloch wavenumber to 1.
l.params[:b] = 1;
sol,nn,flag = mslp(l,1000,tol=1E-9, scale=2pi, output=true);Launching MSLP solver... scale: 6.283185307179586 Iter dz: z: ---------------------------------- 0 Inf 999.9999999999999 1 127.20107543230797 872.798924567692 + 7.406416813969517e-14im 2 9.271457110544489 863.5274674571475 + 1.3191702956989754e-14im 3 0.04978223109540639 863.477685226052 + 4.341711409064263e-14im 4 1.4353235664278573e-6 863.4776837907285 - 9.6106619227306e-14im 5 1.0136616056359396e-12 863.4776837907295 - 1.2484850598797733e-13im Solution has converged! ...finished MSLP! ##################### Results ##################### Number of steps: 5 Last step parameter variation:6.369023706983804e-12 Auxiliary eigenvalue λ residual (rhs):6.580759984541837e-8 Eigenvalue:863.4776837907295 - 1.2484850598797733e-13im
The computation is much faster than with the full mesh, but the eigenfrequency is exactly the same!
Writing output to paraview
Not only the eigemnvalues much but also the mode shapes. Bloch wave theory clearly dictates how to extend the mode shape from the unit cell to recover the full-annulus solution.
The `vtkwrite function does this for you. You, therefore, have to provide it with the current Bloch wavenumber. If you provide it with the unitcell mesh, it will just write a single sector to the file.
data=Dict();
v=bloch_expand(full_mesh,sol,:b);
data["abs"]=abs.(v)./maximum(abs.(v));
data["pahase"]=angle.(v)./pi;
vtk_write("Bloch_tutorial",full_mesh,data);Summary
Bloch-wave-based analysis, significantly reduces the size of your problem without sacrificing any precission. It thereby saves you both memory and computational time.
References
[1] Efficient Computation of Thermoacoustic Modes in Industrial Annular Combustion Chambers Based on Bloch-Wave Theory, G.A. Mensah, G. Campa, J.P. Moeck 2016, J. Eng. Gas Turb. and Power,doi:10.1115/1.4032335
[2] The effect of dynamic operating conditions on the thermoacoustic response of hydrogen rich flames in an annular combustor, T. Indlekofer, A. Faure-Beaulieu, N. Noiray and J. Dawson, 2021, Comb. and Flame, doi:10.1016/j.combustflame.2020.10.013
This page was generated using Literate.jl.