Sensitivity analysis of flame residence time

[Patrick Singal]

I am trying to run a sensitivity analysis of the residence time for an NH3/air flame. As far as I can tell, Cantera doesn't have any internal tools or examples for accomplishing this. While I attempted to create a sensitivity script adapted from an existing one for flame speed, I ultimately arrived at a dead end. Some pointers regarding the strategy needed to resolve this problem would be greatly appreciated.

Thanks!

Patrick 

Here is a minimal version of my flame residence time script:

import numpy as np

import cantera as ct

import matplotlib.pyplot as plt

species_of_interest = 'NO'

def getFlame(gas,P):

    flame = ct.FreeFlame(gas,width=0.05)

    flame.set_refine_criteria(ratio=3, slope=0.05, curve=0.05)

    flame.soret_enabled = True

    flame.transport_model = 'multicomponent'

    flame.solve(loglevel=0, auto=True)

    state = {

        'T [K]': flame.T,

        'vel [m/s]': flame.velocity,

        'dist [m]': flame.grid #list of distances [m]

    }

    for species in gas.species_names:

        state[f"X_{species}"] = flame.X[gas.species_names.index(species)]

    S_b = list(state['vel [m/s]'])[-1] # [m/s]

    Rjoule = 8.31446261815324 # [J/K/mol]

    full_time = [(distance)/S_b for distance in state['dist [m]']] # [s]

    conc = np.multiply(flame.X[gas.species_names.index(species_of_interest)], np.divide(P*ct.one_atm,np.multiply(Rjoule,state['T [K]'])))

    distance_index=np.argmax(conc)

    tau_f = np.array(state['dist [m]'][distance_index])/S_b

    state['tau [ms]']=[(ft - tau_f)*1000 for ft in full_time] # [ms]

    return state

gas = ct.Solution('gri30.yaml')

fuel = 'NH3'

phi=1.1

T = 450

P = 1

oxidizer = {'O2':1,'N2':3.76}

gas.TP = T, P*ct.one_atm

gas.set_equivalence_ratio(phi,fuel,oxidizer,basis='mole')

state = getFlame(gas,P)

tau = list(state['tau [ms]'])

X_NO = list(state['X_NO'])

plt.figure()

plt.loglog(tau,X_NO)

And here is what I wrote in an attempt at making a sensitivity script before I arrived at a dead end:

def getSensitivity(gas):

    f = ct.FreeFlame(gas,width=width)

    f.set_refine_criteria(ratio=3, slope=0.06, curve=0.1)

    f.solve(loglevel=0, auto=True)

    X0_NO = f.X[gas.species_names.index("NO")]

    print(X0_NO)

    # Su0=f.velocity[0] # m/s

    sensitivities = pd.DataFrame(index=gas.reaction_equations(), columns=["base_case"])

    dk = 0.1

    for r in range(gas.n_reactions):

        gas.set_multiplier(1.0)

        gas.set_multiplier(1+dk,r)

        f.solve(loglevel=0, refine_grid=False, auto=False)

        X_NO = f.X[gas.species_names.index("NO")]

        # Su = f.velocity[0]

        sensitivities.iloc[r,0]=(X_NO-X0_NO)/(X0_NO*dk)

        # sensitivities.iloc[r,0]=(Su-Su0)/(Su0*dk)

    gas.set_multiplier(1.0)

    sensitivities_subset = sensitivities[sensitivities["base_case"].abs() > threshold]

    return sensitivities_subset

[Daniel Thomas]

Hi Patrick,

Did you try doing this with the solve_adjunct method?  It sounds like you've already seen the definition of get_flame_speed_reaction_sensitivities() in onedim.py.  I also found this thread helpful.  Below is a starter function to get temperature or species sensitivity:  Can you define a gX function from f.velocity and the species that gets you the residence time definition you're after?

Hope this helps,

Daniel

import cantera as ct
import numpy as np
import time

def calculate_sensitivity_adjoint(f, S_var='T', i_grid=1, display=True, dk=1e-5):
        ''' Calculates the normalized sensitivity of flame temperature or species concentration w.r.t reaction rates.
            Uses the solve_adjoint method, following get_flame_speed_reaction_sensitivities():
            S_var:    Set to 'T' or species    
            i_grid:   Set to grid index position for analysis  '''
       
        def gX(sim):
            return sim.X[f.gas.species_index(S_var)][i_grid]

        def gT(sim):
            return sim.T[i_grid]
           
        def perturb(sim, i, dp):
            sim.gas.set_multiplier(1+dp, i)
               
        start_time = time.perf_counter()
        Nvars = sum(D.n_components * D.n_points for D in f.domains)  # domains are fuel_inlet, flame, oxidizer_inlet
        i_X   = f.inlet.n_components + f.flame.n_components*i_grid + f.flame.component_index(S_var)
        # i_X is the index in the global solution vector, where:
        #   grid point 1          grid point 2           grid point 3        ...
        #  __(inlet)__       ___(flame node 1)___   ___(flame node 2)___
        # |All components\      |All components\       |All components\
        # |u,V,T,etc.\
        dgdx  = np.zeros(Nvars)
        dgdx[i_X] = 1
        if S_var == 'T':
            T0 = gT(f)
            sensitivities = f.solve_adjoint(perturb, f.gas.n_reactions, dgdx, gT, dp=dk) / T0
        elif S_var in f.gas.species_names:
            X0 = gX(f)
            sensitivities = f.solve_adjoint(perturb, f.gas.n_reactions, dgdx, gX, dp=dk) / X0
        if display:
            print('  Calculated normalized sensitivies of {} at y = {:0.2f} mm (i={}) for {} reactions (Adjoint method in {:0.2f} s)'.format(
                  S_var, f.grid[i_grid]*1e3, i_grid, f.gas.n_reactions, time.perf_counter()-start_time))
        return sensitivities

gas = ct.Solution('gri30.yaml')

oxidizer = {'O2':1,'N2':3.76}

gas.TP = 450, ct.one_atm
gas.set_equivalence_ratio(1.1,'NH3',oxidizer,basis='mole')
f = ct.FreeFlame(gas,width=0.05)

f.set_refine_criteria(ratio=3, slope=0.06, curve=0.1)
f.solve(loglevel=0, auto=True)

S = calculate_sensitivity_adjoint(f,'NO',10,True)