Tutorial 08 Custom FTF

The typical use case for thermoacoustic stability assessment will require case-speicific data. e.g, measured impedance functions or flame transfer functions. This tutorial explains how to specify custom functions in order to include them in your model.

Model set-up

The model is the same Rijke tube configuration as in Tutorial 01. However, this time the FTF is customly defined. Let's start with the definition of the flame transfer function. The custom function must return the FTF and all its derivatives w.r.t. the complex freqeuncy. Therefore, the interface must be a function with two inputs: a complex number and an integer. For an n-τ-model it reads.

function FTF(ω::ComplexF64, k::Int=0)::ComplexF64
    return n*(-1.0im*τ)^k*exp(-1.0im*ω*τ)
end

FTF (generic function with 2 methods)

Of course, it is sometimes hard to find a closed-form expression for an FTF and all its derivatives. You need to specify at least these derivative orders that will be used by the methods applied to analyse the model. This is at least the first order derivative for the standard mslp method. You might return derivative orders up to the order you need it using an if clause for each order. Anyway, the function may get complicated. The implicit method explained further below is then a viable alternative.

For convenience we can also specify a display signature, that is used to nicely display our function when the discretized model is displayed in the REPL. Let's say in our case we want the function to be named "ntau(z)" where "z" is a placeholder for whatever the name of the argument will be. We achieve this by adding a method to our FTF function that takes a symbol as input and returns a string.

function FTF(z::Symbol)::String
    return "ntau($z)"
end

FTF (generic function with 3 methods)

Now let's set-up the Rijke tube model. The data is the same as in Tutorial 01

using WavesAndEigenvalues.Helmholtz
mesh=Mesh("Rijke_mm.msh",scale=0.001) #load mesh
dscrp=Dict() #initialize model descriptor
dscrp["Interior"]=(:interior, ()) #define resonant cavity
dscrp["Outlet"]=(:admittance, (:Y,1E15)) #specify outlet BC
γ=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
x_ref=[0.0; 0.0; -0.00101]; #reference point
n_ref=[0.0; 0.0; 1.00]; #directional unit vector of reference velocity
R=287.05; # J/(kg*K) specific gas constant (air)
speedofsound(x,y,z) = z<0. ? sqrt(γ*R*Tu) : sqrt(γ*R*Tb);
c=generate_field(mesh,speedofsound);

btw this is the moment where we specify the values of n and τ

n=1; #interaction index
τ=0.001; #time delay

but instead of passing n and τ and its values to the flame description we just pass the FTF

dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref,FTF)); #flame dynamics

#... discretize ...
L=discretize(mesh,dscrp,c)

1006×1006-dimensional operator family: 

ω^2*M+K+ntau(ω)*Q+ω*Y*C

Parameters
----------
λ	Inf + 0.0im
ω	0.0 + 0.0im
Y	1.0e15 + 0.0im

and done! Note how the function is correctly displayed as "natu(ω)" if we wouldn't have named it, the FTF will be displayed as "FTF(ω)" in the signature of the operator.

Checking the model

Solving the problem shows that we get the same result as with the built-in n-τ-model in Tutorial 01

sol,nn,flag=mslp(L,340,maxiter=20,tol=1E-9,scale=2*pi)
# changing the parameters

Launching MSLP solver...
scale: 6.283185307179586
Iter   dz:     z:
----------------------------------
0		Inf	340.0
1		121.73944539642724	218.8954344765888 + 12.416794063505094im
2		59.03920596186603	171.84143313726085 + 48.075576972715im
3		11.172355692677408	170.75031345770205 + 59.19452429557703im
4		0.3944514828712027	171.14381739654743 + 59.2218486885704im
5		0.0004947274139142701	171.14332284343598 + 59.22183555715669im
6		7.767030819263656e-10	171.14332284266058 + 59.221835557201786im
Solution has converged!
...finished MSLP!
#####################
 Results
#####################
Number of steps: 6
Last step parameter variation:4.880169392400843e-9
Auxiliary eigenvalue λ residual (rhs):1.1691665198048127e-5
Eigenvalue:171.14332284266058 + 59.221835557201786im

Because the system parameters n and τ are defined in this outer scope, changing them will change the FTF and thus the model without rediscretization. For instance, you can completely deactivate the FTF by setting n to 0.

n=0

0

Indeed the model now converges to a purely acoustic mode

sol,nn,flag=mslp(L,340*2*pi,maxiter=20,tol=1E-11)

(####Solution####
eigval:
ω = 1709.429929657814 + 6.990634427877916e-13im

Parameters:
λ = 3.879158853133288e-9 + 5.058060765172243e-22im
Y = 1.0e15 + 0.0im
, 5, 0)

However, be aware that there is currently no mechanism keeping track of these parameters. It is the programmer's (that's you) responsibility to know the specification of the FTF when sol was computed.

Also, note the parameters are defined in this scope. You cannot change it by iterating it in a for-loop or from within a function, due to julia's scoping rules.

# Something else than n-τ...

Of course, this is an academic example. In practice you will specify something very custom as FTF. Maybe a superposition of σ-τ-models or even a statespace model fitted to experimental data.

A fairly generic way of handling experimental FTFs was intoduced in [1]. This approach is not considering a specific FTF but parametrizes the problem in terms of the flame response. Using perturbation theory, it is then possible to get the frequency as a function of the flame response. This allows for closing the system with a specific FTF in a subsequent step and, hence, finding the eigenfrequency by solving a scalar equation! Here is a demonstration on how that works:

dscrp["Flame"]=(:flame,(γ,ρ,Q02U0,x_ref,n_ref)) #no flame dynamics specified at all
H=discretize(mesh,dscrp,c)
η0=0
H.params[:FTF]=η0
sol_base,nn,flag=mslp(H,340*2*pi,maxiter=20,tol=1E-11)

(####Solution####
eigval:
ω = 1709.429929657814 + 6.990634427877916e-13im

Parameters:
λ = 3.879158853133288e-9 + 5.058060765172243e-22im
Y = 1.0e15 + 0.0im
FTF = 0.0 + 0.0im
, 5, 0)

The model H has no flame transfer function, but the flame response appears as a parameter :FTF. We set this parameter to 0 and solved the model, i.e. we computed a passive solution. This solution will serve as a baseline. Note, that it is not necessary to use a passive solution as baseline, a nonzero baseline flame response would also work. Anyway, the baseline solution is used to build a 16th-order Padé-approximation. For convenience we import some tools to handle the algebra.

using WavesAndEigenvalues.Helmholtz.NLEVP.Pade: Polynomial, Rational_Polynomial, prodrange, pade
perturb_fast!(sol_base,H,:FTF,16)
func=pade(Polynomial(sol_base.eigval_pert[Symbol("FTF/Taylor")]),8,8)
ω(z,k=0)=func(z-η0,k)

ω (generic function with 2 methods)

Newton-Raphson for finding flame response

Now, we can deploy a simple Newton-Raphson iteration for finding the eigenfrequency for a chosen FTF as a postprocessing step. This is possible because we have a good approximation of $\omega$ as a function of the flame response $\eta$ and a FTF would link $\omega$ to a flame response. Hence, the eigenfrequency of the problem is implicitly given by the scalar (sic!) equation:

\[\eta = FTF(\omega(\Delta\eta))-\eta_0\]

Which we can first solve for $\eta$ using Newton-Raphson

\[Δη_{k+1}=Δη_{k}-\frac{FTF(ω(Δη_k))-Δη_k}{FTF'(ω(Δη_k))ω(Δη_k)-1}\]

and then for $\omega$ by plugging $\eta$ into $\omega=\omega(\eta)$.

n=1
τ=0.001
η=sol_base.params[:FTF]
for ii=1:10
    global omeg,η
    omeg=ω(η)
    η-=(FTF(omeg)-η)/(FTF(omeg,1)*ω(η,1)-1)
    println("iteration #",ii," estimate η:",η)
end
println(ω(η)/2pi," η :",η," Res:",FTF(omeg)-η)
sol_exact,nn,flag=mslp(L,ω(η)/2pi,maxiter=20,tol=1E-11,output=false) #solve again
ω_exact=sol_exact.params[:ω]
println(" exact=$(ω_exact/2/pi)  vs  approx=$(ω(η)/2pi))")

iteration #1 estimate η:-6.889689350013574 - 1.545586716417518im
iteration #2 estimate η:2.46764965679154 + 3.108439027260453im
iteration #3 estimate η:0.38931011361579326 - 0.6938411268487634im
iteration #4 estimate η:0.8329795102538471 - 1.294171967331538im
iteration #5 estimate η:0.6886569977998448 - 1.2734675411614125im
iteration #6 estimate η:0.6897696051382094 - 1.2763169125012945im
iteration #7 estimate η:0.6897680950356609 - 1.276316713108412im
iteration #8 estimate η:0.6897680950358203 - 1.2763167131087532im
iteration #9 estimate η:0.6897680950358202 - 1.2763167131087534im
iteration #10 estimate η:0.6897680950358197 - 1.2763167131087536im
171.14332064830836 + 59.22183491426868im η :0.6897680950358197 - 1.2763167131087536im Res:0.0 + 2.220446049250313e-16im
 exact=171.14332284265123 + 59.22183555720103im  vs  approx=171.14332064830836 + 59.22183491426868im)

Works like a charm!

Summary

You can analytically define custom FTFs. In order, for high-order perturbation theory to work, you will need to specify high order derivatives of your custom function, too. Sometimes this is cumbersome. As an alternative, you can first solve the problem parametrized in the flame response only and then closing the problem in a subsequent step using a custom FTF of which only the first derivative is known and a Newton-Raphson-Iteration.

References

[1] C. F. Silva, L. Prieto, M. Ancharek, P. Marigliani, and G. A. Mensah, Adjoint-Based Calculation of Parametric Thermoacoustic Maps of an Industrial Combustion Chamber, JEGTP, GTP-20-1520, doi:10.1115/1.4049295

This page was generated using Literate.jl.